Resummation for (boosted) top-quark pair production at NNLO+NNLL' in QCD

We construct predictions for top quark pair differential distributions at hadron colliders that combine state-of-the-art NNLO QCD calculations with double resummation at NNLL' accuracy of threshold logarithms arising from soft gluon emissions and of small mass logarithms. This is the first time a resummed calculation at full NNLO+NNLL' accuracy in QCD for a process with non-trivial color structure has been completed at the differential level. Of main interest to us is the stability of the M_{tt} and top-quark p_T distributions in the boosted regime where fixed order calculations may become strongly dependent on the choice of dynamic scales. With the help of numeric and analytic arguments we confirm that the choice for the factorization and renormalization scales advocated recently by some of the authors is indeed optimal. We further derive a set of optimized kinematics-dependent scales for the matching functions which appear in the resummed calculations. Our NNLO+NNLL' prediction for the top-pair invariant mass is significantly less sensitive to the choice of factorization scale than the fixed order prediction, even at NNLO. Notably, the resummed and fixed order calculations are in nearly perfect agreement with each other in the full M_{tt} range when the optimal dynamic scale is used. For the top-quark p_T distribution the resummation performed here has less of an impact and instead we find that upgrading the matching with fixed-order from NLO+NNLL' to NNLO+NNLL' to be an important effect, a point to be kept in mind when using NLO-based Monte Carlo event generators to calculate this distribution.


Introduction
The top quark is the heaviest fundamental particle discovered so far. Because of its large mass, the top quark is the Standard Model particle which couples most strongly to the Higgs boson; as such, it plays a pivotal role in the study of the electroweak symmetry sector of the Standard Model. For this reason, and thanks to the large number of top-quark pairs produced at the Large Hadron Collider (LHC), accurate experimental measurements of several top-quark related observables are either already available or will become available in the next few years.
An exciting feature of top physics at the LHC is that the large collider energy enables the study of boosted top-quark production. In this context "boosted" refers to the kinematic regime where the energies of the produced top quarks are much larger than the top-quark mass. Boosted top quarks may appear in the study of many kinematic distributions: top quarks with very high pair invariant mass, top quarks with very large transverse momentum, or very forward top quarks (large rapidity). The LHC has already probed top quarks with transverse momenta around 1 TeV [1], and will extend the energy range up to a few TeV. While these boosted energy top quarks are more rarely produced than the low energy ones, and the high energy regions do not contribute significantly to the total cross section, they are phenomenologically important due to their potential to probe directly physics beyond the electroweak scale.
NNLO QCD calculations for differential cross sections in top-quark pair production were obtained in [23,24]. These calculations added to NNLO results for more inclusive quantities such as the total cross section [25][26][27][28] and the forward-backward asymmetry at the Tevatron [29]. In [23], distributions such as the tt (top-pair) invariant mass and the top-quark transverse momentum distribution were evaluated at the LHC with center-ofmass energies of √ s = 8 TeV and 13 TeV using renormalization and factorization scales µ varied around the top mass m t , and the typical scale uncertainties were estimated to lie below 10%. In [24], on the other hand, dynamical scale choices were investigated in order to determine which choice of scale is most appropriate for fixed-order studies of multi-TeV differential cross sections, based on the convergence of the fixed-order perturbative series.
A notable result of that study is that the high-energy tails of distributions can be quite sensitive to the parametric choice of factorization scale, even at NNLO and especially in the case of the top-pair invariant mass distribution. The fact that the application of fixed-order perturbation theory in the boosted regime is rather delicate is not entirely surprising. Indeed, in the case of boosted top quarks one encounters two potential difficulties. The first is that in the boosted regime the result of a fixed-order calculation contains mass logarithms of the form ln(E t /m t ) arising from quasicollinear gluon emissions. The second is that, due to the shape of the parton distribution functions (PDFs), the effective partonic center-of-mass energies in most events with boosted top quarks are not much larger than the invariant masses of the top quark pairs. This can lead to enhanced corrections from so-called soft logarithms or threshold logarithms, and indeed, in general the two types of logarithms multiply each other at a given order in perturbation theory, due to emissions which are simultaneously soft and collinear. At the LHC, these logarithmic corrections may become important numerically, to the point that higher order corrections are not generically much smaller than lower order ones, and while it may be possible to deal with this to some extent through a judicious scale-setting procedure in a fixed-order calculation, it is desirable to address the issue head-on by resumming these two types of logarithms to all orders in perturbation theory and adding them onto the fixed order result through a matching procedure.
A framework for the simultaneous resummation of threshold and small-mass logarithms in differential cross sections in top-quark pair production at the LHC was set up in [30], and applied to the phenomenological study of tt invariant mass and top-quark transverse momentum distributions in [31]. The results are valid in the soft gluon emission limit, but contain an extra layer of resummation for small-mass logarithms, so that overlapping soft and small-mass logarithmic corrections are properly taken into account. The factorization formalism underlying the resummation was derived using soft-collinear effective theory (SCET), and the resummation was carried out in Mellin space. While the focus of that phenomenological study was the boosted regime, a matching procedure with standard soft-gluon resummation results (i.e. without the simultaneous resummation of small-mass logarithms) at next-to-next-to-leading logarithmic (NNLL) accuracy [32] and to NLO fixedorder calculations was performed in order to expand the validity of the results to the full phase space.
In this paper we combine state-of-the-art results from analytic resummation and fixedorder perturbation theory in QCD in order to produce for the first time phenomenological predictions which match NNLO results with resummation of soft and small-mass logarithms at the level of differential tt invariant mass and top-quark transverse momentum distributions.
The paper is organized as follows. In section 2 we immediately address the technical aspect which is new in this work, namely the matching of the resummed calculation in the soft and boosted limit to the NNLO calculations. In section 3 we explain the resummation procedures used in two distinct kinematic limits. In particular, in section 3.1 we review the kinematics for top-pair production and the form of the differential cross section in Mellin space. In section 3.2 we discuss the soft limit, where threshold logarithms are large but mass logarithms are of generic size, while in section 3.3 we consider the boosted-soft limit, where both soft and mass logarithms are considered large. With all the analytic tools ready, in section 4 we study kinematic features of the top quark pair in the boosted region, and use these insights to determine appropriate choices for the scales appearing in the matching functions in the resummed formulas. In section 5 we present numerical results for the top-pair invariant mass and top-quark p T distributions, paying close attention to a comparison between pure NNLO results and those supplemented by resummation, and then present conclusions in section 6. We perform more comparisons among results at different orders in fixed-order and resummed perturbation theory in appendix A, and relegate some lengthy formulas for resummation exponents to appendix C.

Matching fixed order and resummed calculations
We study the top-quark pair production process N 1 (P 1 ) + N 2 (P 2 ) → t(p 3 ) +t(p 4 ) + X(p X ) , (2.1) where N 1 and N 2 are the colliding hadrons (proton-proton for the LHC and future colliders, proton-antiproton for the Tevatron), and X is an inclusive hadronic final-state. In this work, the top quarks are treated as on-shell particles. QCD factorization allows us to write an arbitrary differential cross section in the schematic form The symbol ⊗ stands for a convolution over longitudinal momentum fractions of the initialstate partons i, j ∈ {q,q, g} in the (anti-)proton, the dσ ij are the differential partonic cross sections for the process i + j → tt +X, theX is a partonic final state, and φ i/N denotes the parton distribution function (PDF) of parton i in hadron N . While the PDFs are non-perturbative objects to be extracted from experiment, the differential partonic cross sections can be calculated in perturbative QCD. The aim of this paper is to provide state-of-the-art QCD calculations for two specific differential cross sections: namely the invariant mass of the top-quark pair, and the p T of the top-quark. The baseline for the calculation is NNLO in fixed-order perturbation theory, to which we add the two types of resummation mentioned in the introduction and to be described in detail in section 3. The first is performed in the soft limit of the differential partonic cross sections, where the top-quark pair carries almost all of the energy of the partonic collision. It can be obtained to NNLL m order using the results of [32], where we have labeled the logarithmic accuracy of the resummation with a subscript m to indicate that the result is obtained for arbitrary values of m t . As the energy of the top-quark pair becomes very large, this standard soft-gluon resummation itself develops logarithms which become large in the limit m t → 0. We call this the "boosted-soft limit", and perform a joint resummation of overlapping soft and small-mass logarithms using the formalism developed in [30]. In this case it is possible to increase the accuracy of the resummation to NNLL b order, where now the subscript b indicates that the results are valid in the boosted-soft limit, and thus neglect corrections which vanish in the limit m t → 0. The perturbative ingredients for these two types of resummation (anomalous dimensions and matching functions) and the order at which they need to be calculated to achieve a given resummation accuracy is summarized in table 1 in section 3.4.
We have just described three different calculational formalisms, each of which is tailored to a different kinematic situation. The NNLO calculation is optimal in regions of phase space where the top quarks are not highly boosted, and hard-gluon emissions are important. The NNLL m result is applicable when soft-gluon radiation dominates, but small-mass logarithms are unimportant. When soft-gluon radiation dominates, and the top-quarks are highly boosted, one would like to make use of the NNLL b results. To make optimal use of our results we would like to have a unified description over the whole phase space. For this purpose it is necessary to combine the different formulas in such a way that no contribution is counted more than once.
To understand such a matching procedure, we first consider matching NNLL m resummation with fixed-order results. In that case, the matching formula with (N)NLO reads where dσ NNLLm denotes the differential cross section evaluated to NNLL m accuracy. The first term in the above equation contains the all-orders resummation result in the soft limit, and the difference of terms in parenthesis contains subleading terms in that limit, which are taken into account by the fixed-order calculation. Explicit expressions for the fixedorder (N)NLO expansion of the resummation formulas needed in the matching are given in section 3.6. We next match the resummation formulas in the soft and boosted-soft limit with each other and with (N)NLO. To do so, we need to remove the overlap between the NNLL b and NNLL m results to all orders in α s . This is done by exploiting the fact that the boosted-soft resummation formula is the small-mass limit of the soft-gluon resummation formula at any fixed order in α s , so we must subtract out the leading term in the limit m t → 0 in order not to double count. The combined result, denoted NNLL b+m , is thus given by where the terms in the parenthesis account for contributions which are suppressed by m t /M in the boosted-soft limit and thus not included in the NNLL b result. Matching with fixed order then proceeds in analogy to eq. (2.3) resulting in Again, the terms in the parenthesis account for subleading terms in the soft limit, which are taken into account through a fixed-order calculation, but are not accessible to either of the resummation formulas. Calculating the subtraction term requires one to expand each term in eq. (2.4) to (N)NLO. The procedure for obtaining the different components of the above equation is described in more detail in section 3.6. Differential distributions obtained from the explicit evaluation of eq. (2.5) are a main result of this work, and can be found in section 5. Before going into numerical studies, we give details of the resummation procedure, including recipes for obtaining the different pieces used in the matching procedure, in the next section. These details can safely be skipped by a reader interested in purely phenomenological results. In section 4 we conduct a thorough analysis of the choice of scales for the matching functions which appear in the resummed results.

.1 Kinematics and differential cross sections
In this section we review the kinematics involved in describing the limits in which resummation is carried out. At Born level, and to leading order in the soft limit considered below, two partonic channels contribute to the partonic cross section: the quark-antiquark annihilation channel and the gluon fusion channel The momenta of the incoming partons are related to the hadron momenta according to p i = x i P i (i = 1, 2). The relevant Mandelstam invariants are defined as In fixed order perturbation theory, starting from NLO accuracy, a new 2 → 3 production channel, initiated by a quark and a gluon, opens up. At NNLO one needs to account for the contribution of 2 → 4 processes (at tree level), as well as the contribution of 2 → 3 processes (up to one loop) and the 2 → 2 production channels in eqs. (3.1) and (3.2) (up to two-loops). All of these channels are of course included in the NNLO results which we employ in this work. The soft emission region, which is of interest for the resummed calculation, is defined by the limit M 2 tt →ŝ (sometimes also referred to as the partonic threshold region). In this limit, the final state particles in addition to the top pair are soft. In order to describe the top-pair invariant mass distribution near the partonic threshold, it is convenient to introduce the following variables: The quantity β t is the 3-velocity of the (anti-)top quark in the tt rest frame, while β is often invoked to describe the partonic threshold for the total cross section [33][34][35][36]. In the soft limit z → 1, one has β → β t . Moreover, in that limit the scattering angle θ is related to the Mandelstam variables according to from which one can easily verify the usual relation M 2 tt + t 1 + u 1 = 0. We will perform resummation on the double differential cross section in the top-pair invariant mass and the scattering angle θ. Applying the generic QCD factorization formula eq. (2.2) allows one to write this differential cross section as where M tt and cos θ are in the ranges The hard-scattering kernels C ij are proportional to the partonic cross sections and can be calculated in perturbation theory at the factorization scale µ f , while L ij are nonperturbative parton luminosity functions, defined as In this work, we study the resummation in two limits where soft and collinear gluon emissions dynamically generate scales much lower than the scattering energy: In each of these limits the perturbative expansions of the hard-scattering kernels contain large logarithms of scale ratios, which can be resummed to all orders in α s . For the discussion of resummation that follows, it is convenient to study the cross section in Laplace or Mellin space. The Mellin transform and its inverse are defined bỹ (3.11) where in the inverse transform the real part of the contour c is chosen such that it lies to the right of all singularities in the functionf (N ). Convolutions such as the differential cross section in eq. (3.6) become simple products in Mellin space. Indeed, by performing the Mellin transform of eq. (3.6) with respect to τ , we find the differential cross section in Mellin space, which reads The soft limit z → 1 corresponds to N → ∞, as can be seen by taking the Mellin transform of the plus distributions appearing in the partonic cross section, which are generated by soft gluon emissions: and we have introduced the variableN = N e γ E in order to simplify the expressions, with γ E denoting the Euler constant. We note that the partonic cross section contains terms of the form α n s P k (z) where 0 ≤ k ≤ 2n − 1 at N n LO in its perturbative expansion. In Mellin space this becomes α n s L k where L = lnN and 0 ≤ k ≤ 2n, which will be important when we consider the tower of logarithms which can be resummed in section 3.4. In Mellin moment space, the soft and boosted-soft limits in eqs. (3.9) and (3.10) become Besides the invariant mass distribution, the distribution in the transverse momentum p T of the top quark is also interesting. In [37], a different formulation, dubbed "1PI", was employed to deal with such "single-particle-inclusive" observables in the boosted regime. In contrast, the formulation used in this paper was called "pair-invariant-mass" (PIM) kinematics. Near the partonic threshold, the PIM and 1PI formulations differ only by power-suppressed contributions. We can exploit this fact to express the differential cross section with respect to p T in terms of the double differential cross section in eq. (3.6), thereby avoiding the introduction of the 1PI formulation. In the threshold limit, the transverse momentum p T and the rapidityŷ of the top-quark in the partonic center-ofmass frame are given by and we can write the differential cross section with respect to p T andŷ in Mellin space as In the above equation, it is understood that M tt and cos θ should be expressed in terms of the integration variables according to where we have defined the transverse mass m T . The transverse momentum distribution can be obtained by integrating overŷ in the range

Resummation in the soft limit
Resummation of top-quark hadroproduction cross sections in the partonic threshold limit z → 1 was considered in [38]. More recently, top-pair production in the soft limit was reanalyzed by means of SCET methods 1 . With this technique, it was possible to study the resummation of soft gluon emission corrections both in PIM and 1PI kinematics [32,40] in momentum space. Here we follow closely the discussion and notation of the NNLL calculation of [32], although we perform the resummation in Mellin rather than momentum space. This approach was adopted in [31] for top-pair production. Subsequently, it was also applied to the evaluation of the soft emission corrections to ttW ± [41], ttH [42], and ttZ [43] to NNLL accuracy. The starting point is the factorization of the partonic cross section in the soft limit. In Mellin space, this takes the form The hard functions H m ij and the soft functions s m ij , referred to generically as matching functions, are matrices in color space -explicit results for the two-by-two matrices in the qq channel and three-by-three matrices in the gg channel up to NLO in α s can be extracted from [32], where the derivation of a momentum-space factorization formula analogous to eq. (3.22) is given in detail. Comparing to the notation employed in [32], here we have introduced a superscript m on these functions, indicating that they contain full dependence on the top quark mass m t . This is to distinguish them from the corresponding functions with m t = 0 used later in the boosted-soft limit.
Given the factorized form of the partonic cross section in Mellin space, one can derive and solve renormalization group (RG) equations for the component functions. This allows one to evaluate the hard and soft functions at an arbitrary hard scale µ h and soft scale µ s , where large logarithms are absent. One then uses RG evolution to obtain the hard scattering kernels at the factorization scale µ f . These RG-improved hard-scattering kernels can be written as where we have suppressed the dependence of all functions in eq. (3.23) on the variables M tt , m t and cos θ. The evolution matrices U m ij contain all large logarithms in an exponentiated form, and thereby resum them to all orders in α s . The explicit form of the evolution matrices was derived in [32], and in Mellin space reads (suppressing for the moment the subscript ij and the argumentN ) The explicit exponential in the above equation is color-diagonal, and contains the evolution functions where γ stands for an anomalous dimension such as Γ cusp or γ φ , and β(α s ) = dα s (µ)/d ln µ is the QCD β-function. The matrix u m in the second line of eq. where γ h,m is the color non-diagonal part of hard anomalous dimension with full mass dependence, and P denotes path-ordering. The definition and the explicit expressions for the various anomalous dimensions can be readily found in the appendix of [32]. Although we have dropped indices indicating the partonic channel in eqs. (3.24) and (3.26), one should keep in mind that the anomalous dimension γ φ , (related to the PDFs), the cusp anomalous dimensions Γ cusp , and the non-color diagonal anomalous dimension γ h,m are different in the quark-annihilation and gluon fusion channels. The above integrals appearing in the evolution matrices can be evaluated and truncated to a given logarithmic order. The results can be expressed in terms of the strong coupling constant evaluated at the various scales, α s (µ f ), α s (µ h ) and α s (µ s ), as was done in [32]. This is convenient if the soft scale µ s is chosen directly in momentum space, and is thus a real number. However, in the current paper, we employ a more conventional choice of the soft scale in Mellin space, µ s ∼ Λ/N for some mass scale Λ, which is now a complex number. It is therefore more convenient to re-express α s (µ s ) in terms of α s (µ h ), keeping in mind that ln(µ h /µ s ) is a large logarithm. Actually, it is conventional to express all α s (µ i ) in terms of α s (µ h ), where µ i is any scale other than µ h appearing in the formula. This can be done by using the perturbative evolution of the strong coupling, which up to 3-loop order is given by (see, e.g. [44]) Logarithmic orders are then determined according to powers of α s (µ h ) with ln(µ h /µ i ) counted as order 1/α s (µ h ). For this reason we introduce two O(1) parameters: λ s and λ f defined by The part of the exponent proportional to the identity matrix in eq. (3.24) can now be expanded as a series in α s (µ h ), and the resulting evolution matrix can be written in the form Explicit expressions for the RG exponents g m i are given in appendix C.1, while a method for evaluating u m is detailed in [32] using techniques from [45].
The all-order resummed hard scattering kernel eq. (3.23) is formally independent of the matching scales µ s and µ h . However, the truncation of the resummed formula to a given logarithmic order introduces residual dependence on these scales. In order to perform the resummation, one must choose these scales in such a way that the fixed-order expansion of the hard and soft functions are free of large logarithms. The explicit form of the oneloop hard and soft functions reveals that the leading logarithmic terms are α s ln 2 (M 2 tt /µ 2 h ) and α s ln 2 (M 2 tt /N 2 /µ 2 s ), respectively, which motivates the naïve choices µ h ∼ M tt and µ s ∼ M tt /N . However, we will study the analytic form of the hard and soft functions in greater detail in section 4 in order to make a more informed choice of appropriate scales. It should be noted however, that in picking the soft scale µ s directly in Mellin space (as we do in this work) the resummed hard-scattering kernel contains a branch cut at large N due to the Landau pole in the running of α s . This leads to ambiguities in the choice of contour for the inverse Mellin transform in eq. (3.11). We come back to this issue when discussing the numerical implementation of our results in section 3.5. As mentioned earlier, one could also choose the soft scale directly in momentum space, which is then independent of N , as was done in [32]. With this choice the inverse Mellin (or Laplace) transform for the soft function can be carried out analytically and is free from the Landau pole problem. On the other hand, this comes at the price of resumming a different tower of logarithms compared to the pure partonic threshold ones, as discussed in [40] and explored in more detail in, e.g. [46,47].

Resummation in the boosted-soft limit
Resummation of top-quark hadroproduction cross sections in the boosted-soft limit eq. (3.10) was considered in [30]. We first collect the main results from that paper in the absence of perturbative corrections involving closed top-quark loops. In that case, the functions H m ij and s m ij in eq. (3.22) can be subfactorized in the m t → 0 limit as In the above formulas, the hard functions H ij and soft functions s ij without the superscript m are independent of the top-quark mass m t . They were calculated to NNLO in [48] and [49], respectively, and can also be applied to di-jet production. All m t -dependence is factorized into the two functions C D and s D , which are related to the perturbative heavy-quark fragmentation function [50] and were extracted at NNLO in [30]. After this refactorization, the result for the partonic cross section in the boosted-soft limit reads As for the soft limit, the resummed hard-scattering kernel can be obtained by deriving and solving RG equations for the component functions in the above factorization formula. We write the result in the form In the l.h.s. of eq. (3.34) we dropped the dependence on the arguments M tt , cos θ and m t . Similarly, U ij , H ij and s ij also depend on M tt and cos θ.
As mentioned with the massive hard and soft functions, we postpone discussion of the most appropriate scale choices for the matching functions H ij and s ij until section 4. However, since the functions C D and s D depend only on the scales m t and m t /N , we can immediately make the assignment µ dh = m t and µ ds = m t /N . The evolution of the hard and soft functions is encoded in the functions U ij , which are matrices in color space, while the evolution of C D and s D is in the color-diagonal function U D . Suppressing for the moment the channel labels ij, the explicit expressions for the evolution matrices are given by The functions S A , a γ , and u are defined analogously to those in eqs. (3.25) and (3.26). Here A is given by 2Γ q cusp in the qq channel and Γ q cusp + Γ g cusp in the gg channel. Definitions of the various anomalous dimensions and their explicit perturbative expansions can be found in [30]. In order to write the perturbative expansion of the evolution functions in the same form as eq. (3.30) we introduce two additional O(1) parameters, λ dh and λ ds in complete analogy to λ s and λ f , as defined earlier in eq. (3.29). Following the same procedure as in the soft limit, the evolution matrices can then be written as and We list the (lengthy) expressions for the RG-exponents g i and g D i in appendix C.2. In the presence of heavy-quark loops the factorization of the partonic cross section in the boosted-soft limit is more involved. It is necessary to introduce additional coefficients related to matching six-flavor PDFs, heavy-quark fragmentation functions, and α s onto five-flavor ones. To ease notation we cluster these contributions from heavy quarks together in coefficientsc ij t . Such corrections are proportional to powers of n h , the number of heavy flavors, and introduce additional m t dependence into the formula via logarithms of the form ln n (m t /M tt ). Our factorized hard scattering kernel is then written as It is not entirely clear whether logarithms of the form ln n (m t /M tt ) appearing through heavy-quark loops 2 can be systematically resummed. Therefore, we add them onto the resummation formula eq. (3.34) using fixed-order perturbation theory. At NNLO, these mass logarithms come from two sources: those from the interference of two-loop amplitudes with tree-level ones, and those from one-loop amplitudes squared. The one-loop squared contributions can be extracted from the results in [51,52], while the two-loop terms were calculated in [53,54]. Numerically, we have found that the contributions of these n h terms to the differential cross sections are almost negligible.

Resummation accuracy
Having obtained resummed hard scattering kernels in the soft eq. (3.23) and boosted-soft eq. (3.34) limits, we now examine what level of resummation can be achieved given the current status of perturbative calculations. At this point, it should be pointed out that there exist two naming schemes for the logarithmic accuracies of resummed results, as discussed in [46] and summarized in table 1 of [55]. While they are purely conventions and one is free to choose either, it is important to have internal consistency with the earlier works [30-32, 37, 40]. We therefore adopt the so-called "Notation " outlined in table 1 of [55] to denote the accuracies of our resummed results. In table 1, we list the perturbative orders at which the matching functions and anomalous dimensions need to be evaluated in order to achieve resummation at a given logarithmic accuracy.
As highlighted in section 3.1 in the discussion following eq. (3.14), in Mellin space the perturbative expansion of the resummed cross section gives corrections of the form α n s L k where L = lnN . The power k of logarithms included in the expansion of the resummed result at a given logarithmic accuracy is also indicated in the last column of table 1. As can be seen there, the difference between the NNLL and NNLL accuracies amounts to a single logarithm at each order in perturbation theory.
The cusp anomalous dimension is fully known to three-loop order [56], results for the other anomalous dimensions to NLO can be found in [57][58][59][60][61][62][63], and the massive hard H m ij and soft functions s m ij have been extracted to NLO [32]. We can therefore perform resummation in the soft limit up to the NNLL accuracy 3 . On the other hand, the matching functions in the massless limit H ij , s ij , C D and s D are all known to NNLO [30,48,49], enabling Table 1. Our naming scheme for the logarithmic accuracies. We list the perturbative orders at which the cusp anomalous dimension, the QCD β-function, all other anomalous dimensions and matching functions need to be evaluated in order to obtain resummation at a given logarithmic order.
resummation to NNLL accuracy in the boosted-soft limit. In terms of the evolution factors in eqs. (3.30), (3.37), and (3.38), we need to keep the first three g-functions and compute the evolution matrices u m and u to second order to obtain NNLL or NNLL accuracy. Keeping only the first two g-functions and the LO u m and u matrices results in NLL resummation as can be seen by the lower perturbative order of the anomalous dimensions in the first line of table 1.

Mellin inversion
Given the results for Mellin-space resummed hard-scattering kernels in eqs. (3.23) and (3.34), we must perform the inverse Mellin transform of eq. (3.11) in order to get the differential cross section in momentum space. For the invariant mass distribution we have 40) and the transverse momentum distribution follows similarly. As indicated in eq. (3.40), there are two ways of carrying out the inverse transform. The first line utilizes the parton luminosity functions in Mellin space (N -space), while the second line uses the PDFs in momentum space (x-space). The x-space PDFs are easier to obtain, but the z-integration is numerically unstable due to the singular behaviour of L(y) around y ∼ 0, and one needs to use some tricks to improve the convergence [66]. The use of the N -space PDFs, on the other hand, avoids these instabilities and leads to a fast numerical implementation. However, public libraries such as LHAPDF [67] only give x-space PDFs and one needs to construct the N -space ones by performing the Mellin transform and analytically continuing to complex N . To this end we employ the methods from [55,68], namely, we approximate the x-space luminosity functions in terms of Chebyshev polynomials, from which the Mellin transform can be carried out analytically. At fixed order in α s , the inverse Mellin transform in eq. (3.40) is well-defined (in the sense of the delta function and plus distributions) and independent of the integration contour as long as it lies to the right of all (physical) singularities of the hard-scattering kernels c ij . However, after resummation the c ij develop dependence on the running coupling α s (µ) at the soft scale µ s with the canonical choice µ s ∼ M tt /N (and likewise with the soft-collinear scale µ ds ∼ m t /N in the boosted-soft case). This introduces an unphysical Landau pole singularity at large N in c ij , and the inverse transform is ambiguous against the choice of contour. In this paper we adopt the so-called Minimal Prescription (MP) [69], in which the contour is chosen to be to the left of the Landau pole but to the right of all other singularities. It is because of this prescription that the first integration in z runs from τ to ∞, not τ to 1. Our hard scattering kernelc ij no longer vanishes for z > 1, but does asymptote to zero quickly enough so that parton model assumptions are not violated. This point is addressed in detail in appendix B of [69]. Other prescriptions such as the Borel Prescription exist in the literature [70], and we have checked that in our case there is no large numerical difference between the two.

Ingredients for matching with fixed-order
We outlined the procedure for matching fixed-order calculations with soft (and boostedsoft) gluon resummation formulas in section 2. In the following, we make more precise some of the definitions introduced in the matching equations.
In eq. (2.3), we need the NLO and the NNLO expansions of the NNLL m resummed result. For the NLO expansion, we can make use of the following result: The reason that the NLO expansion of the NNLL formula is so simple is that µ h and µ s dependence in the NLO matching functions cancels against factors that come from expanding the RG evolution factors in eq. (3.23). The end result can thus be obtained directly by setting these scales to µ f at the beginning as indicated in eq. (3.41), which turns off the RG evolution and leaves the NLO matching functions evaluated at those scales. A similar result holds for the NNLO expansion of the NNLL b result, namely The NNLO expansion of the NNLL m resummed result, which we write as is not as simple, because the NNLO matching functions are absent. As a result, the µ h and µ s dependence does not completely cancel at NNLO. We can express dσ NNLLm,(2) by inserting the following c (2) into eq. (3.12) or (3.18) in place of c ij : where we have suppressed all arguments in the expansion coefficients of the hard and soft functions with the exception of the scale at which they should be evaluated. The expansion coefficients are defined through The form of c (2) in eq. (3.44) requires a bit of explanation. First of all, eq. (3.41) tells us that c (2) should vanish if we set µ h = µ s = µ f , since the NNLL m calculation includes only the NLO matching coefficients. Secondly, note the fact that the NNLL m resummed result dσ NNLLm would be independent of µ s and µ h up to NNLO if the NNLO contributions from the hard and soft functions were known and included (which would upgrade the resummation accuracy to NNLL m ). Therefore, in a fixed order expansion of the resummed result to NNLO, one should recover the first line of eq. (3.44) by adding the NNLO contributions from the hard function (evaluated at µ h ) and the soft function (evaluated at µ s ). The form of eq. (3.44) follows directly from these two facts. Evidently, we must also explain how to evaluate eq. (3.44) without knowing the massive two-loop hard and soft functions H (2) m and S (2) m . The logarithmic terms of these two functions can be determined from their RG equations and in fact this is all one needs. The non-logarithmic "constant" terms which could also appear in H (2) m and S (2) m (and which are not determined by the RG equation) cancel between the first and second lines in eq. (3.44) and therefore do not contribute to c (2) . This concludes the matching of the NNLL m results with fixed-order calculations.
We now turn to eqs. (2.4) and (2.5), which describe the matching between (N)NLO, NNLL m , and NNLL b results. For this we need the NNLO expansion of dσ NNLL b given in eq. (3.42), as well as the NNLO expansion of the NNLL resummation formulas given in eq. (3.43). A further ingredient is the m t → 0 limit of the NNLL m formula appearing in eq. (2.4). To evaluate that we exploit the fact that the boosted-soft resummation formula is the small-mass limit of the soft resummation formula at any fixed order in α s . This leads to the result This follows because setting µ dh = µ h and µ ds = µ s in the boosted-soft result removes RG evolution between the functions H and c D in eq. (3.31), and s and s D in eq. (3.32), thus leaving behind the leading contributions from threshold resummation in the limit m t /M tt → 0. It is worth mentioning some additional subtleties concerning the virtual top-quark loops at NNLO (we mentioned in the footnote on page 14 that we do not consider real topquark pair emissions). As discussed in the last paragraph of section 3.3, it is not entirely clear how to resum these contributions to all orders in α s , and we choose to add them in fixed order. For the NNLO+NNLL result this is automatically taken into account by the matching formula while for the NLO+NNLL result we have to add them manually. A complication arises from the fact that the soft resummed result, dσ NNLLm , generates some (but not all) of the α 2 s n h terms through RG running. We need to subtract these terms out before adding back the full NNLO heavy quark contributions introduced at the end of section 3.3 in order to avoid double counting. Again, we only consider these contributions for completeness, practically their numerical impact is negligible.

Choosing the kinematics-dependent matching scales
The (massive) hard and soft functions depend on the Mandelstam variables t 1 , u 1 andŝ, and the appropriate choice of the matching scales µ h and µ s appearing in the resummation formalism depends on the kinematic regime probed by the differential cross section under study. The focus of this section will be to identify a well-motivated set of these kinematicsdependent matching scales for the top-pair invariant mass and p T distributions, which can be used all the way from the low-energy regime where all kinematic invariants are on the order of the top-quark mass, up to the boosted regime where m t is small compared to M tt or p T . We begin our analysis by highlighting and comparing some important kinematic features of the low-energy and high-energy regions of the M tt and p T distributions in section 4.1. We then devote section 4.2 to scale choices for the top-pair invariant mass distribution and section 4.3 to those for top-quark p T distribution.
Throughout the rest of the paper we use the following numerical inputs. We fix the LHC collider energy to √ s = 13 TeV, take m t = 173.3 GeV, and use the NNPDF3.0 PDF sets with α s (M Z ) = 0.118 [71] in conjunction with LHAPDF6 [67]. We use NNLO PDFs for all predictions unless otherwise indicated. In the numerical evaluation of the resummed formulas we have made use of the CUBA integration library [73].

Some kinematic considerations that underpin the choice of scales
In this section we point out some important kinematic features of the M tt and top-quark p T distributions which are instrumental in determining appropriate values of the matching scales µ h and µ s . In particular, we study the high-energy and low-energy regimes for both distributions, and explain the differences between them that impact their description in fixed-order and resummed perturbation theory. The main idea is that soft gluon resummation works in the limit where the kinematic features of a given observable resemble those of the LO process. This can be studied by introducing kinematic variables sensitive to higher-order hard emissions. A kinematic variable we find particularly useful is and R T ≤ 1, with R T = 1 corresponding to θ = π/2 (central scattering). Moreover, one finds that the Jacobian factors arising from rewriting the double differential cross sections in terms of R T and either M tt or p T have integrable singularities proportional to 1/ 1 − R 2 T . For instance, Beyond LO, however, H T is only constrained to be smaller than √ŝ . Hard emissions generate non-vanishing cross section in the region R T > 1, in addition to contributing to the R T < 1 region, and the singularity at R T = 1 is resolved into a Jacobian peak. The quantity R T thus offers a useful kinematic discriminant: the more sensitive an observable to the region R T ≤ 1, the greater the potential for an improved prediction through resummation. Conversely, observables characterized by R T > 1 are inaccessible to soft kinematics and dominated by hard emissions.
In figure 1, we show the distributions of R T in four kinematic regions, normalized to the integrated cross section in each region. The first is M tt ∈ [380, 420] GeV, representative of "low M tt ", M tt 2m t ; the second is M tt ∈ [2500, 3000] GeV, representative of "high M tt ", M tt 2m t ; the third and fourth regions are p T ∈ [50, 100] GeV and p T ∈ [1200, 1400] GeV, representative of "low p T " and "high p T " respectively. We will also refer to these as "low-energy" and "high-energy" bins of M tt or p T . An interesting observation is that the differences between the LO and NLO distributions are quite distinctive in the four cases. It is the task of this subsection to explore the explanations and implications of these facts.
We first discuss the low-M tt bin of the top-pair invariant mass distribution, whose R T distribution is shown in the top-left plot of figure 1. At both LO (the green curve) and NLO (the black curve), the distribution is peaked at R T ∼ 1, a fact following from eq. (4.3). However, at NLO, hard emissions generate a non-negligible fraction of the distribution in the R T > 1 region. In terms of a perturbative description beyond NLO, this indicates that while soft-gluon resummation can potentially improve the description in this bin, matching with NNLO in order to describe the R T > 1 region is needed for a precise prediction.
The R T distribution in the high-M tt bin is displayed in the bottom-left plot of figure 1. A remarkable feature of this bin is that the R T distribution is peaked at values close to R T ∼ 0.2, both at LO and NLO. This feature can be traced to a property of the partonic cross section for gg → tt scattering, namely that the Born diagrams possess tand uchannel singularities in the limit m t → 0. These singularities are cut off by the finite value of the top-quark mass. However, we show below with analytic arguments that their presence implies that in the boosted regime where m t M tt the integral over cos θ used to calculate the differential cross section is dominated by the regions of small R T . One implication of this feature is that the high-invariant mass region is especially amenable to soft-gluon resummation, as we shall see in section 5. Another is that, in spite of what one might expect, the most relevant perturbative scale in the high-energy tail of the top-pair invariant mass distribution is H T rather than M tt itself. As a result, the partonic cross section, as well as the hard and soft functions appearing in the factorized form of the cross section in the soft limit, are better evaluated at scales µ ∼ H T in this region. This was indeed a main result of the fixed-order study of the top-pair invariant mass distribution to NNLO [24], which identified µ f = H T /4 as the choice at which fixed-order perturbation theory converges best.
We now turn to the R T distributions in the two bins of top-quark p T . The low and high p T regions are shown in the top-right plot and the bottom-right plot of figure 1, respectively. The Jacobian peaks at R T ∼ 1 are still present, similar to the M tt case. The distributions exhibit long tails toward lower values of R T . These come from the tand u-channel singularities, the same effect at work in the high M tt bin discussed in the last paragraph. On the other hand, the behavior of the NLO distribution in this case is rather different from the M tt case. The measured value of the top quark p T can only constrain the combined p T of the anti-top-quark and the extra parton, but puts no constraints on their separate transverse momenta. It is therefore kinematically allowed to have the tt pair recoiled against a separate hard parton, such that one ends up with a large H T and a small M tt . This explains the long tail at the NLO towards R T > 1, particularly in the high p T bin. In fact, at µ f = m T /2, one finds that 47% of the distribution lies at R T > 1. From these facts, we expect that it is important to incorporate the effects of hard emissions for the p T distribution, especially in the high-p T region.

Scales for the M tt distribution
The main goal of this and the next subsection is to identify the optimal choices of the hard scale µ h and the soft scale µ s , based on the kinematic features of the hard and soft functions. Since the hard and soft functions may be evaluated analytically, this will also help to understand the kinematic features of the R T distribution given in the previous subsection.
The philosophy of RG-improved perturbation theory is to choose the matching scales such that the fixed-order expansion of the hard and soft functions is well behaved. The massless hard and soft functions depend on the kinematic invariants M 2 tt , t 1 and u 1 . As long as all of these scales are of the same size numerically, the choices µ h ∼ M tt and µ s ∼ M tt /N free these functions of potentially large logarithmic corrections and thus ensure good perturbative convergence. The situation becomes more subtle in the boosted regime, where m t M tt . In that case, when the top quark is produced in the very forward region, cos θ → 1, the kinematic invariants develop a hierarchy |t 1 | M 2 tt ∼ |u 1 |. An analogous hierarchy develops for very backward production, with |t 1 | ↔ |u 1 |. Both of these situations correspond to the region H T M tt , since for Born kinematics On the other hand, the region where m t ∼ M tt corresponds to H T ∼ M tt , irrespective of the value of cos θ. The top-pair invariant mass distribution is calculated from the double differential cross section in eq. (3.6) by integrating over the scattering angle in the region −1 < cos θ < 1. At large M tt , the results for the R T distribution shown in the bottom-left panel of figure 1 and discussed in the previous subsection make clear that the integral is dominated by the region where | cos θ| ∼ 1 and R T is significantly smaller than unity. This kinematic feature is explained by the fact that at Born level the gg partonic cross section has tand u-channel singularities, related to the hard function. For example, in the limit t 1 → 0, the LO hard function reads where x t ≡ −t 1 /M 2 tt and N c = 3 is the number of colors in QCD. The expression in the limit u 1 → 0 is obtained by replacing x t → 1−x t . Therefore, at fixed M tt m t , the cross section gets large contributions from the region H T M tt , due to tand uchannel enhancements in the gg channel. The qq channel is free of such tand uchannel enhancements.
The dynamical enhancement of the forward and backward scattering regions at large M tt has important implications for the choice of the matching scales µ h and µ s . We can study this issue analytically by expanding the higher-order corrections to the hard and soft functions in the limit R T → 0. The gg partonic cross section is symmetric under the exchange of t 1 and u 1 , so the R T → 0 limit can be easily obtained from the x t → 0 results, which we focus on for concreteness. While the soft function itself has no 1/x t singularities in this limit, it enters the factorization formula in a matrix product with the hard function, so we must deal with both functions at once in order to take the correct limit of the differential cross section. This leads us to study the higher-order perturbative corrections at the level of the objects and where we have suppressed the dependence of the matching functions on all parameters and scales other than µ h and µ s . Explicit results for the NLO corrections can be written as and those for the NNLO corrections as where we have set N c = 3 and the number of light quarks to N l = 5 in the above equations. 4 An important feature of eqs. (4.9) and (4.10) is that both the NLO and the NNLO corrections depend on the two physical scales −t 1 and M tt (through the ratio x t ). Therefore, any choice of µ h and µ s will lead to corrections of the form α n s ln m (x t )/x t in the x t → 0 limit. However, the structure of such corrections is rather different for the hard and soft functions. For the hard function, the choice µ 2 h = −t 1 frees the NLO corrections of such logarithmic corrections in x t → 0 limit, and reduces the logarithmic terms in the NNLO corrections to a single power of ln x t . On the other hand, the choice µ h = M tt generates a double logarithmic series whose corrections have the form α n s ln 2n (x t )/x t . Using that √ −t 1 = H T /2 in the x t → 0 limit along with the symmetry of the gg channel under t 1 ↔ u 1 , one thus expects the perturbative corrections from the hard function to the M tt distribution to be well behaved across phase space with the choice µ h = H T /2.
For the soft corrections S gg , double logarithmic corrections of the form α n s ln 2n (x t )/x t are generated even when µ s = √ −t 1 /N is chosen. To understand this result, we note that NLO corrections to the massless soft function can be written in the form [30,49] where s 12 = s 34 = M 2 tt , s 13 = s 24 = −t 1 , s 14 = s 23 = −u 1 , and the sum over (I, J) excludes the terms where I = J. The matrices w IJ differ for the qq and gg channels, and can be found in, for instance, [49]. The NLO corrections to the soft functions are thus characterized by the three different scales s IJ /N , the relative importance of each scale being determined by the properties of matrix elements w IJ , which are pure color factors.
Taking the x t → 0 limit of S (N)NLO gg then leads to the double logarithmic series mentioned above, irrespective of the choice of µ s .
The analytic results above give very useful insight into the nature of perturbative corrections arising from the hard and soft functions. In particular, they hint at the use of a H T -based scale for the hard function, which is also a reasonable choice for the soft function. However, they are derived in the formal limit R T → 0, which for the top-pair invariant mass distribution is relevant because of a dynamical enhancement from tand uchannel singularities in the high M tt region. It is therefore useful and necessary to supplement these analytic arguments with a numerical study. To do this, we define NLO and NNLO K factors for the hard and soft functions in the following way. For the hard functions, we evaluate the differential cross section with respect to M tt using H (N)NLO ij for the hard-scattering kernel in eq. (3.12), and define NLO and NNLO K factors by dividing these results by those found using H LO ij . These ratios are indicated by K H,(N)NLO ij (M tt , µ h ). The parton luminosity cancels in the ratio since the soft function is independent of the Mellin parameterN at LO, so we can write Although not immediately apparent from the results above, one finds that the real parts of H Ratios K S,(N)NLO (M tt , µ s ), which take into account corrections from the soft function, are defined similarly. The important difference in the case of the soft function is that it depends onN and one must take the product with the Mellin-transformed parton luminosities before performing the inverse Mellin transform in order to get the contribution to the differential cross section. For simplicity, we use these luminosities evaluated at the scale µ f = M tt . We have checked that the luminosity dependence nearly completely cancels in the ratio defining the soft K factors, so that the results depend very little on the exact choice of µ f and also the collider energy. and µ s = rM tt /N (blue lines). At M tt = 400 GeV there is little difference between the K factors with the two different parametric choices, and both the hard and the soft K factors are moderate as long as the proportionality factor r is not too small. At the higher value of M tt , the corrections with the two different parametric scale choices differ by quite a large amount. For the hard function, the corrections remain moderate for µ h ∼ H T /2, as anticipated from the analysis above. For the soft function, the NLO corrections cannot be made smaller than about 50%. This happens, for instance, at µ s ∼ H T /N . The K factor at this scale is flat with respect to changes of µ s around this value, at both the low and high values of M tt and the NNLO corrections are also moderate.
The above analysis leads us to identify µ h = H T /2 and µ s = H T /N as a well-motivated choice of matching scales across the full range of M tt in the gg channel. We have checked that the soft and hard K factors in the qq channel (which gives a considerably smaller contribution to the cross section at both large and small M tt ) are also well behaved for these choices. Therefore, we will use these choices by default in all further numerical analysis of the M tt distribution. Of course, at the level of differential cross sections the dependence on these scales cancels against that in the RG evolution factors, so that resummed results are independent of the exact choice at the order at which one is working. The unphysical, residual dependence can be reduced by calculating more orders in the logarithmic series, but with a proper choice of matching scales these higher-order terms are expected to be small corrections.
It is worth noting that [31], as well as all other works on soft-gluon resummation with PIM kinematics, used M tt -based matching scales instead of H T -based ones. For moderate M tt there is little difference, but at higher M tt the perturbative uncertainties estimated from µ h and µ s variations are larger with the M tt -based choice. Some numerical results with M tt -based choices are given at the end in appendix A.

Scales for the p T distribution
We now turn our attention to the top-quark p T distribution. We again motivate a suitable choice of µ h and µ s by studying K factors for the hard and soft functions analogous to those for the M tt distribution, but this time obtained by substituting the hard-scattering kernel in eq. (3.18) with the appropriate terms in eqs. (4.7) and (4.8) before inverting the Mellin transform to compute the cross section as a function of p T . We show the results in figure 3 for a high and a low value of p T , examining both M tt and m T based scales. Note that while K H,(N )N LO (M tt , µ h ) displayed explicitly in eq. (4.12) can be calculated without reference to the parton luminosities, this is not the case for K H,(N )N LO (p T , µ h ). We used µ f = m T in the luminosities in calculating both the hard and soft K factors, and have checked that varying µ f to other values produces only a small effect in the K factor ratios.
The soft and hard K factors exhibit only a mild hierarchy between the two types of scale choices in the low-p T region, which is even smaller in the high-p T region. This can be understood from figure 1, which shows that in both p T regions the cross section sits mainly at 2m T ∼ M tt due to the Jacobian peak at R T = 1. The K factors are moderate when µ h = m T and µ s = 2m T /N , and we shall use these as the default choices in our resummed  calculations. Recall that in the soft limit the p T of the top is equal to that of the anti-top, so that m T = H T /2 and these are the exact same choices as for the M tt distribution after a trivial renaming.

Results and discussion
In this section we give our main results for the top-pair invariant mass and (anti) top-quark p T distributions, as well as the total cross section, with a focus on comparing NNLO results with NNLO+NNLL ones. Some further comparisons across different perturbative orders and between standard and threshold resummed PDFs are presented in appendices A and B respectively. Although we present only a limited set of results for the LHC operating at a center-of-mass energy of 13 TeV, distributions with alternate binning and at different collider energies can be produced on request from the authors.
Results for the absolute (normalized) M tt distribution are shown in left (right) panel of figure 4. The NNLO results use µ f = H T /4 by default (we shall always set the renormalization scale appearing in the NNLO calculation to µ r = µ f unless otherwise specified), which is the scale favored by the analysis of perturbative convergence of the fixed-order series performed in [24]. The NNLO+NNLL results are obtained from the matching relation so that this method amounts to adding the uncertainties from independent scale variations in quadrature. A remarkable feature of figure 4 is that the NNLO+NNLL and NNLO results are in close agreement when µ f = H T /4 is chosen. To add context to this result, we compare in figure 5 the ratio of the NNLO and NNLO+NNLL results with µ f = M tt /2 to the NNLO result with µ f = H T /4, using the same set of matching scales and method of estimating perturbative uncertainties as in figure 4. These two figures deliver a couple of important messages. Firstly, the NNLO+NNLL result is rather stable against switching the factorization scale between H T -based and M tt -based schemes. This implies that the even higher order corrections to the NNLO+NNLL result are not so important. On the other hand, the NNLO result changes drastically when switching the schemes. In particular, higher order contributions beyond NNLO encoded in the resummation produce a very large effect for the choice µ f = M tt /2, as already forseen in [31]. Given these observations, the close compatibility between the NNLO+NNLL result (with either scale choice) and the NNLO result with µ f = H T /4 is a highly non-trivial fact. This provides an important confirmation of the result of [24], which favors the choice µ f = H T /4 for the fixed-order calculation of the M tt distribution. The overall picture emerging from the above analysis is that the perturbative description of the top-quark pair invariant mass distribution is under good control.
Results for the absolute (normalized) average top/anti-top (p T,avt ) distribution at NNLO and NNLO+NNLL are shown in the left (right) panel of figure 6. The NNLO checked that the uncertainties estimated this way differ little from those obtained by varying µ f and µr with the 7-point method. The NNLO+NNLL calculation varies four resummation scales and also µ f = µr independently and adds the uncertainties in quadrature, so a direct numerical comparison with the 7-point method is not straightforward. However, we have experimented with a 7-point scan over µ f and µ h , and found that the uncertainty estimates change only marginally compared to the quadrature method. results (with which resummation is matched) have been calculated using the definition dσ dp T,avt = 1 2 dσ dp T,t + dσ dp T,t , where p T,t (p T,t ) denotes the transverse momentum of the top (anti-top) quark, and we have labeled the distributions in figure 6 accordingly. The p T distribution is calculated using the scale choice µ f = m T /2 (where m T refers to the transverse mass of either the top or anti-top quark depending on the distribution under consideration), which is favored by the study [24]. The resummed results use µ h = m T and µ s = 2m T /N by default, as justified in the previous section. The bands refer to perturbative uncertainties estimated through scale variations using the same procedure as for the M tt distribution above. We see that the NNLO+NNLL result is consistent with the NNLO one. On the other hand, we show in appendix A that upgrading matching with fixed-order from NLO+NNLL to NNLO+NNLL is an important effect for the p T distributions, especially in reducing the scale uncertainties in the high p T region. This is an important fact to keep in mind when using NLO-based Monte Carlo event generators to model p T distributions. Finally, in figure 7 we show results for the total cross section, obtained in several different ways. The NNLO and NNLO+NNLL results with µ f = H T /4 are obtained by integrating the top-pair invariant mass distribution in figure 4, while those with µ f = m T /2 are obtained by integrating the p T distribution in figure 6. In these results with dynamical scales, perturbative uncertainties are estimated through the same procedure of scale variations used for the distributions, and are displayed as error bars in figure 7. These are compared to the "standard" results for the total cross section, which are calculated using fixed scales with µ f = µ r = m t by default. We obtain them from the Top++ program [74], which implements both the NNLO results from [28] as well as a soft-gluon resummation in the absolute threshold production limit β t → 0 [75]. In these fixed scale results, pertur-NNLO NNLO + NNLL 700 750 800 850 σ(pb) bative uncertainties are estimated in both the NNLO and the NNLO+NNLL results by varying µ f and µ r up and down by a factor of two using the seven-point method. Evidently, while this resummation result for the total cross section is also labelled NNLO+NNLL in figure 7, one should keep in mind that it uses a different framework than the current work, including the treatment of resummation scales and how they are varied as just described. From Figure 7 we see that the integral of both the NNLO+NNLL M tt distribution with µ f = H T /4 and the NNLO+NNLL p T distribution with µ f = m T /2 yield nearly the same total cross section as the widely quoted result from the Top++ program. This shows that the results of this work not only provide the most precise QCD results for the top pair invariant mass distribution and the top quark transverse momentum distribution across phase space, but also give the correct normalization for these distributions.

Conclusions
In this paper we combined state-of-the-art results from soft-gluon resummation (NNLL ) and fixed-order perturbation theory (NNLO) in order to produce NNLO+NNLL predictions for the top-pair invariant mass and the average top/anti-top quark p T,avt distributions at hadron colliders. These results represent the most complete QCD calculations of these observables to date. They are also the first instance where an NNLO calculation has been supplemented with resummation in a process where the Born-level cross section contains four partons and thus has non-trivial matrix structure in color space.
The resummation formalism used here contains several elements which have not appeared in the literature so far. Some of these involve the details of weaving together three different kinds of calculations to obtain results optimized throughout phase space. In this procedure, it is crucial to avoid counting the same contribution more than once. In particular, in section 2 we presented a matching procedure which allows us to combine NNLO results in fixed order [24], NNLL b results in a joint resummation of overlapping soft and collinear logarithms [30], and NNLL m results in pure soft-gluon resummation [32], in order to achieve what we have called NNLO+NNLL accuracy. All ingredients required to implement this matching procedure, as well as the resummation itself, carried out in Mellin space, were given in section 3 and appendix C.
Our analysis in section 4 revealed some important kinematic features of the M tt and p T distributions. These can not only be understood analytically using results from soft-gluon resummation, but also affect its implementation, and are useful to keep in mind when interpreting the numerical results for the M tt and p T distributions given in section 5. In the case of the M tt distribution, the analysis of section 4 demonstrates that in the boosted regime where m t M tt , the most relevant hard scale is not M tt itself but rather H T . This fact is related to the dynamical enhancement of the forward and backward scattering regions due to the tand u-channel diagrams appearing in the Born level partonic cross section in the gg initiated production process. We used this feature to identify a wellmotivated set of matching scales for the kinematics-dependent hard and soft functions appearing in the resummation formalism, and also to argue that the high M tt region is particularly amenable to soft gluon resummation. On the other hand, we observed that the high-energy region of the p T distribution is rather sensitive to hard emissions, so that matching resummation with NNLO is an essential improvement. This is certainly the case for the analytic resummation performed here, but should also be kept in mind when predicting the high-energy tail of the p T distribution with NLO-based event generators.
In section 5 we presented numerical results for the absolute and normalized M tt and p T,avt distributions, as well as the total cross section, valid to NNLO+NNLL . For the p T,avt distribution the resummation effects are mild, especially at the scale µ f = m T /2 favored by the NNLO analysis of perturbative convergence in [24]. For the M tt distribution, an interesting outcome of our analysis is the stability of the NNLO+NNLL results under parametric changes in µ f , as shown in figures 4 and 5. Given the large shifts in the NNLO calculation under such µ f changes, it is an important result that the NNLO+NNLL results stabilize the differential cross section close to the NNLO prediction with µ f = H T /4, which is the setting favored by [24] and currently being used in all NNLO phenomenology. The consistency between the NNLO and NNLO+NNLL results gives us confidence that even higher-order corrections are under good control.

A Perturbative stability across orders
In this section we perform some comparisons of results for the top-pair invariant mass and average top/anti-top-quark p T distribution across different perturbative orders. dσ/dp T,avt (pb/GeV) p T,t/t ∈ [50, 100] GeV to be drawn from the figure is that the NNLO+NNLL results at different µ f congregate near the NNLO one with µ f = H T /4. Figure 9 shows results for the p T,avt distribution at different perturbative orders in the two sample bins considered in section 4.1. The resummed results use the same default matching scales as in section 5 and are completely analogous to those in figure 8. Compared to the M tt distribution, where the parametric difference between the H T and M tt -based µ f is large, there is no such hierarchy of scales to consider for the p T distribution. Therefore, we have shown results for three different m T -based choices, ranging from µ f = m T /2 to µ f = 2m T by default. While the resummation is of some benefit in stabilizing the (N)NLO results in the low-p T bin, the picture is less clear in the high-p T bin. For instance, there is a dramatic reduction in µ f dependence in the NNLO+NNLL results compared to the NLO+NNLL ones. This is an indication that the high-p T region is more sensitive to hard radiation than the high-M tt region. We have given some qualitative explanations for why this should be the case when discussing the R T distribution in section 4.1. Numerically, we have found that the NLO results for the high-p T region of the distribution are quite sensitive to both the qg channel and the R T > 1 region, in a strongly µ f -dependent fashion. Softgluon resummation cannot stabilize such µ f dependence, which explains the importance of matching to NNLO in fixed order.
The uncertainties associated with each of the distributions presented here result from  Table 2. Cross section in "high" and "low" energy bins. The first uncertainties refer to those from µ f variation only while the second set are generated by the variation of the matching scales. a combination of the uncertainties generated from the variation of each scale in accordance with eqs. (5.1) and (5.2). It is interesting to decompose the source of these uncertainty bands in terms of the contributions which arise from varying the factorization scale µ f compared with the other matching scales. We present such a sample decomposition in table 2 for the two "low" and "high" energy bins used throughout section 4. We show the cross section in each bin with two sets of uncertainties, the first refers to those obtained through variations of µ f alone while the second refers to the combined uncertainty generated by varying each of the matching scales. In most instances, the dominant contribution to the uncertainty arises from the variation of µ f . For the M tt distribution at high M tt the uncertainty from the matching scales is larger than at low M tt while across the p T spectrum this source of uncertainty remains constant. Finally, it is interesting to compare the results of the present work with the earlier resummation results at NLO+NNLL accuracy in [31]. To this end it is important to note that the differences between the current work and the work in [31] are two-fold: 1) we have matched with the exact NNLO calculation here compared to an NLO matching in [31]; and 2) we have employed different settings of the factorization scale and the matching scales than [31]. Concerning the scale choices, we remind the reader that in this work we by default use µ h = H T /2, µ s = H T /N , µ f = H T /4 for the M tt distribution, and µ h = m T , µ s = 2m T /N , µ f = m T /2 for the p T distribution. We refer to this scale setting as the "new scales". On the other hand, [31] by default uses µ h = M tt , µ s = M tt /N , and µ f = M tt /2 (for the M tt distribution) or µ f = m T (for the p T distribution). We call this scale setting the "old scales". In order to quantify the effect of switching to the new scales and the effect of matching with NNLO, we compare the 3 kinds differential cross sections: 1) NLO+NNLL with old scales from [31]; 2) NLO+NNLL with new scales from the present work; and 3) NNLO+NNLL with new scales which are the best predictions of the present work. Such a comparison for the absolute as well as normalized M tt distributions is shown in figure 10. It can be seen from the plots on the left side that, by changing to the new scales, the apparent scale uncertainties are reduced. This is a hint that with the new scale choice the higher order corrections are indeed smaller. The effect of matching with the NNLO results is shown in the plots on the right side of figure 10. We see that the matching changes the differential cross section most significantly in the low M tt bins, which can be expected. A similar comparison in the case of p T distributions is shown in figure 11. One can see that the change of default scales does not reduce the scale uncertainties (actually leading to larger uncertainties in the high p T bins). Matching with the NNLO is important here, which stabilizes the scale variation. Overall, the NNLO+NNLL results are consistent with the NLO+NNLL results when the new scales are used.

B Comparison with resummed PDFs
In the main part of this work, all computations have been carried out using the (N)NLO NNPDF3.0 PDF set [71]. However, one can also use threshold resummed PDFs [72], into which the effects of soft gluon radiation are incorporated. Including these contributions in the PDFs can produce non-negligible differences in threshold resummed cross section predictions when compared to those obtained with regular PDFs. An example of this can be found in [76], which considered sparticle pair production. In particular, for a large range of sparticle masses considered it was found that while K-factors (defined as the ratio of the resummed result to the fixed order one) obtained with regular PDFs were almost always greater than unity, this was not always the case when resummed PDFs were employed. Although the resummed PDFs are obtained with a reduced data set and are therefore not yet suitable for precision calculations, we consider this an appropriate opportunity to explore their implications for top-quark pair production at the LHC. In order to compare results obtained with threshold resummed PDFs to those with standard PDFs we compute the M tt and p T,avt distributions using the same settings as in section 5. To enable a fair comparison the NNPDF collaboration also provides PDFs which do not include threshold resummation but which are compiled from the same data set as the resummed ones. Specifically we use the NNPDF30 nnlo disdytop PDFs as the benchmark for PDFs without threshold resummation and the NNPDF30 nnll disdytop PDFs for those with threshold resummation.
In figure 12 we compare the central values of the NNLO and NNLO+NNLL predictions for the p T,avt and M tt distributions, obtained using µ f = m T /2 and and µ f = H T /4, respectively. The NNLO predictions are calculated using the PDFs without resummation, while the NNLO+NNLL predictions are computed using the PDFs with (labelled NNLL PDF) and without (labelled FO PDF) threshold resummation. In this manner, the NNLO and NNLO+NNLL predictions with the fixed order PDFs act as an equivalent set of predictions to those in the main part of our paper, but now with a PDF set which can be compared to those incorporating threshold resummation. The left plot in figure 12 shows the M tt distribution, where the effect of using resummed PDFs is significant. Here the use of resummed PDFs produces a suppression in the cross section at high M tt compared to predictions produced with regular fixed order PDFs. The plot on the right, showing the p T,avt distribution, also displays noticeable changes. Here again the resummed PDFs produce a greater suppression of the cross section in the tail of the distribution than when regular PDFs are used. It will be interesting to perform further studies using resummed PDFs in the future as the fits improve and experimental measurements in the tails of these distributions become more accurate.

C The RG exponents
In this appendix we collect explicit expressions for the RG exponents appearing in eqs. (3.30), (3.37), and (3.38). In order to ease notation we introduce the following shorthand Note that since the g-functions are derived from the part of the evolution functions in eqs. (3.24) and (3.35) which are proportional to the identity matrix in color space, the factors of iπ which appear there cancel out in the resummed cross section. As such, we do not retain these iπ factors in our expressions for the g-functions.

C.1 Soft limit
First, we present the g m i functions appearing in the evolution factor eq. (3.30) for the threshold resummed result:

C.2 Boosted-soft limit
Here we present the g-functions which appear in the evolution factors eqs. (3.37) and (3.38) for the boosted-soft resummation formula. The functions g i are simply given by their massive counterparts (g m i as above) using the replacement, for each instance of γ φ in the g m i in appendix C.1. Each instance of Γ cusp (α s ) in g m i is replaced with A(α s ) (this is due to the presence of S A and a A rather than S Γ and a Γ in eq. (3.35)) with A(α s ) given by A(α s ) → 2Γ q cusp , qq-channel Γ q cusp + Γ g cusp , gg-channel .
We decompose each of the g D i , which are functions of three arguments into two twoargument functions g D i,dh and g D i,ds as follows g D i (λ dh , λ ds , λ f ) = g D i,dh (λ dh , λ f ) + g D i,ds (λ ds , λ f ) .
Using this decomposition, we present below the functions as used in this work.