Threshold expansion at order alpha_s^4 for the t-tbar invariant mass distribution at hadron colliders

We calculate the leading O(alpha_s^4) contributions to the invariant mass distribution of top-quark pairs produced at the Tevatron and LHC, in the limit where the invariant mass of the t-tbar pair approaches the partonic center-of-mass energy. Our results determine at NNLO in alpha_s the coefficients of all singular plus distributions and scale-dependent logarithms in the differential partonic cross sections for q-qbar, gg ->t-tbar + X. A numerical analysis showing the effects of the NNLO corrections on the central values and scale dependence of the invariant mass distribution is performed. The NNLO corrections are found to significantly enhance the cross section and reduce the perturbative uncertainties compared to the NLO calculation.


Introduction
The top quark is the heaviest known particle in the Standard Model (SM) of fundamental interactions. Because of its large mass, it is expected to couple strongly with the fields responsible for electroweak symmetry breaking, and the detailed study of top-quark properties is likely to play a key role in elucidating the origin of particle masses. In addition to observables such as the top-quark mass and total inclusive top-quark pair production cross section, differential cross sections are also of interest. For instance, the tt invariant mass distribution can be used to measure m t [1], and in searches for physics beyond the SM. The presence of bumps in the smoothly decreasing tt invariant mass distribution would be a clear signal of an s-channel heavy resonance [2], which is predicted in many beyond-the-SM scenarios [3].
To date, thousands of top-quark events have been observed by two different experiments at the Fermilab Tevatron, and the top-quark mass has been extracted at the percent level. Searches for narrow width particles decaying into top-quark pairs have been pursued [4], and results for the top-pair invariant mass distribution were recently obtained from data collected by the CDF collaboration [5]. The experiments at the LHC are expected to observe millions of top-quark events per year already in the initial low-luminosity phase, bringing the study of top-quark properties into the realm of precision physics. In particular, the total inclusive top-quark pair production cross section is expected to be measured with a relative error of 5% to 10% at the LHC [6]. To make optimal use of the data requires equally precise theoretical predictions for the measured observables.
Theoretical calculations in QCD rely on the factorization formula for the differential cross section, which is of the form where the symbol ⊗ stands for a convolution. The hard-scattering kernels C ij are related to the partonic cross section and can be calculated as a series in α s , whereas the f i are parton distribution functions (PDFs) for the partons i = q,q, g in the incoming hadrons and must be extracted from data. Current theoretical predictions are based on next-to-leading order (NLO) calculations of the total cross section [7][8][9][10] and differential distributions [11,12] in fixed-order perturbation theory, supplemented with soft gluon resummation to next-to-leading-logarithmic (NLL) order [13][14][15][16][17]. The NLO computations suffer from theory uncertainties larger than 10%, both for Tevatron and LHC center-of-mass energies. These uncertainties are due to our imperfect knowledge of the parton distribution functions, and also to the truncation of the perturbative series in the strong coupling constant, which introduces a dependence on the unphysical renormalization and factorization scales into physical predictions. This theoretical uncertainty is typically reduced by including more terms in the perturbative series, and for this reason the calculation of the differential partonic cross section to next-to-next-to-leading order (NNLO) has been an area of active research [18][19][20][21][22]. However, due to the complexity of the diagrammatic calculations, complete NNLO results are not yet available. In the absence of full results, there has been much recent activity in NNLO calculations and next-to-next-to-leadinglogarithmic (NNLL) resummation for the total cross section near production threshold [23][24][25][26][27], but the analogous results for differential cross sections remain unknown. The goal of this Letter is to present a subset of the NNLO corrections to the hard-scattering kernels for the tt invariant mass distribution at hadron colliders. In particular, we focus on the region where the invariant mass M of the top-quark pair approaches the partonic centerof-mass energy √ŝ , and compute NNLO corrections to the hard-scattering kernels for the qq and gg channels at leading order in 1 − z, where z = M 2 /ŝ → 1 in the threshold region. The leading-order term in this threshold expansion for the differential cross section is equivalent to the virtual-soft approximation, and is written in terms of singular plus distributions and delta functions in the variable 1−z. Our results determine the coefficients of all plus distributions of the form [ln n (1 − z)/(1 − z)] + , as well as all µ-dependent pieces multiplying the δ(1 − z) term, whereas a remaining piece of the delta-function coefficient is left undetermined. The basis for our calculations is a factorization formula for the hard-scattering kernels in the partonic threshold region, and is described in Section 2. In this region, the hard-scattering kernels factorize into products of hard functions, related to virtual corrections, and soft functions, related to real emission in the soft limit [13]. These functions satisfy certain renormalization-group equations, which determine their dependence on the scale µ. By calculating the hard and soft functions at one-loop order, and using results for the two-loop anomalous dimensions recently derived in [28][29][30][31] (see also [32][33][34][35][36]), the µ-dependent logarithms in the hard and soft functions can be determined exactly to NNLO using the renormalization group. The logarithms in the soft function are of the form ln(ŝ(1 − z) 2 /µ 2 ), and uniquely determine the coefficients of the [ln n (1−z)/(1−z)] + distributions. This is similar in spirit to the calculations of [16] for the soft corrections to the NNLO differential cross section using threshold resummation techniques in Mellin space, but goes beyond those results by completely determining the coefficient of the [1/(1 − z)] + distribution, which is sensitive to process-dependent two-loop anomalous dimensions. In Section 3 we perform a short numerical analysis illustrating the impact of the NNLO corrections on the central values and scale dependence of the tt invariant mass distribution.

NNLO corrections at threshold
In this section we discuss our method for determining NNLO corrections to the differential pp (pp) → tt + X cross section in the threshold region. To describe the threshold region we introduce the variables where M is the invariant mass of the tt pair, √ s is the hadronic center-of-mass energy, and √ŝ is the partonic center-of-mass energy. We define the threshold region as the limit z → 1, with ρ = 4m 2 t /M 2 a generic O(1) variable, which is however not too close to unity. Note that this is different from the threshold limit ρ → 1 often used for the total cross section.
The calculation of the differential cross section in the threshold region is greatly simplified, since there is little phase space available for real gluon emission. Hard emissions are suppressed by powers of 1−z, and the partonic scattering process is dominated by hard virtual corrections and the real emission of soft partons, which can be calculated in the eikonal approximation and exponentiate into Wilson lines. The phase space for these processes is effectively twobody, and the fully differential cross section can be written in terms of the kinematic variables entering the Born-level result.
At the Born level the differential cross section receives contributions from the quarkantiquark annihilation and gluon fusion channels The partonic cross section is a function of the kinematic invariantŝ and momentum conservation impliesŝ + t 1 + u 1 = 0. The hadronic cross section is obtained from the partonic one by convoluting with PDFs. The fully differential cross section involves three variables, which can be chosen, for instance, as the invariant mass and rapidity of the tt pair, and the scattering angle θ between p 1 and p 3 in the partonic center-of-mass frame. The kinematic invariants in (4) can be written in terms of these variables according tô Less differential results are derived by integrating the triply differential rate over the appropriate phase space. In this Letter we will be interested in the invariant mass distribution of the tt pair. We shall write this distribution at threshold as where we have introduced the parton luminosity functions For simplicity, we set the factorization and renormalization scales equal and refer to them as µ. The luminosities for qq are understood to be summed over all species of light quarks, and we have made explicit that mixed channels such as qg are power suppressed by 1 − z in the threshold region [13]. The hard-scattering kernels C ij factorize into hard functions, related to virtual corrections, and soft functions, related to soft real emissions [13]. The result is We have used a label i = qq, gg to distinguish the different partonic channels, but will suppress it in the formulas below, referring instead to a single coefficient function. Note that the hard and soft functions are matrices in a color decomposition of the QCD amplitudes. We have usedŝ = M 2 /z ≈ M 2 everywhere except in the first argument of the soft function, since we expect that the exact results contain logarithms of the form ln(ŝ(1 − z) 2 /µ 2 ), as in the case of Drell-Yan production [37]. In [38], we will give details on the derivation of the factorization formula (8) in soft-collinear effective theory. The perturbative expansion of the hard-scattering kernel can be written as The purpose of this Letter is to present new results for the NNLO coefficient C (2) . As mentioned in the introduction, our derivation relies on the renormalization-group equations for the hard and soft functions. The soft function contains singular plus distributions, and its renormalization-group equation is non-local. It is convenient to avoid this complication by working instead with the Laplace transform of the soft function [39,40], defined as Whereas the soft function contains distributions, the Laplace-transformed function is a polynomial in its first argument and satisfies a local evolution equation. To rewrite the hard-scattering kernel in terms of this function, we introduce the notation where ∂ η is a differential operator with respect to an auxiliary variable η. The factorization formula for the hard-scattering kernel then takes the form [38] To evaluate the above formula in terms of distributions in the variable z, one must use analytic continuation to regulate the divergence at z → 1, take derivatives with respect to the auxiliary variable η, and then set η → 0. This procedure yields results in the form of plus distributions and delta functions given in (21) below. From the above discussion, we conclude that it is sufficient to focus on the calculation of c at NNLO. Defining its perturbative expansion in analogy with (9), and using the expansions for the hard and soft functions (8), the result at NNLO reads In the second line we have defined expansion coefficients multiplying explicit powers of ∂ η . We see that obtaining a complete result for the threshold expansion at NNLO would amount to calculating the hard and soft matrices to this order. Although the NLO coefficient c (1) can be extracted from known results for the virtual corrections and real emissions in the soft limit, the hard and soft matrices themselves are so far unknown, except for in the absolute threshold limit ρ → 1 [25,26]. A main result of our work is the calculation of these functions at NLO, for both the qq and gg channels; details of the calculation will appear in [38]. As for the NNLO corrections, we will show in what follows that it is possible to extract all of the coefficients except for c (2) 0 by using the renormalization-group equations for the hard and soft functions, along with the one-loop results and the two-loop anomalous dimensions from [29][30][31]. The renormalization group also determines the µ-dependent logarithms in this remaining term, but to obtain the scale-independent piece would require to evaluate at NNLO the virtual corrections and soft real emissions.
We now sketch how to calculate the coefficient in (14) using the evolution equations for the hard and soft functions. The hard function satisfies [13] where Γ H is given by the matrices Γ qq or Γ gg from [31]. The anomalous dimensions to two-loop order can be decomposed as where the object γ h is defined through a comparison with the explicit results of [31], and Γ cusp is equal to C F γ cusp for qq and C A γ cusp for gg, with γ cusp the universal cusp anomalous dimension. The evolution equation for the soft function follows from that for the hard functions and PDFs, along with RG-invariance of the cross section. The result is [38] where the soft anomalous dimension is given by The object γ φ is defined through the large-x limit of the Altarelli-Parisi splitting functions, which reads Since the anomalous dimensions for the hard function and PDFs are known to two-loop order, we can calculate that of the soft function to this order.
Given the evolution equations, the anomalous dimensions, and the one-loop hard and soft matrices, it is a simple matter to deduce a general expansion of the form and to calculate all of the coefficients up to h in (14), as well as the µdependent piece of c (2) 0 . The hard-scattering kernel C is obtained by evaluating (12). In doing so, powers of ∂ η in c are converted into delta functions and plus distributions in the variable z. The final result can be written as Explicit results for the coefficients in both the gg and qq channels as functions of the variables (M, m t , cos θ, µ) are given in the computer program attached with the electronic version of this article. We stress that the coefficients D 0 , . . . , D 3 are completely determined from our calculations. The coefficient C 0 , on the other hand, is correct only in its µ-dependence. It also receives scale-independent contributions related to the translation from c to C, and from the product of NLO hard and soft functions, which we discard for consistency. To determine the remaining scale-independent piece would require a complete NNLO calculation of the hard and soft matrices. Note that, because this piece is unknown, the expansion of C 0 in terms of ln(M 2 /µ 2 ) considered so far is not unique: expansions in terms of ln(m 2 t /µ 2 ) (or other choices) still satisfy the same evolution equations. In our phenomenological study in the next section we shall consider these two different expansions as one indication of systematic errors related to the missing, scale-independent piece of C 0 . The function R(z) is finite in the limit z → 1 and is not determined from the leading-order threshold expansion. However, it receives a contribution related to the z −η factor in (12). We choose to keep these contributions when evaluating the NLO corrections, but drop them from the NNLO ones. Our results for D 1 , D 2 , D 3 , and the µ-dependent terms agree with [16], while the exact result for D 0 contains process-dependent two-loop effects and is new.
The coefficients D i and C 0 do not depend on the variable z. It is therefore possible to convolute the plus distributions and delta function with the luminosities as in (6). For the phenomenologically interesting values of 0.03 < τ < 0.3 discussed in the next section, we can then make a few general comments concerning the hierarchy of the different plus-distribution and delta-function terms. The contribution of the leading [ln 3 (1−z)/(1−z)] + distribution is the largest in all cases. The [ln 2 (1 − z)/(1 − z)] + distribution is the second largest, its contribution being about 1/3 that of the leading one. The contribution of the delta-function is about 10-20% of the leading plus distribution, and is of comparable size to the [ln(1 − z)/(1 − z)] + and [1/(1 − z)] + distributions.

Phenomenological applications
In this section we study the phenomenological implications of our results. The main goal is to show the impact of the NNLO corrections on the central values and perturbative uncertainties of the tt invariant mass distributions at the LHC and Tevatron. However, given that our NNLO results are incomplete, it is important to clarify under what conditions they are expected to give a good approximation to the full results. While it is difficult to assign a systematic uncertainty to an unknown correction, we can examine in parallel how well the same approximation works at NLO, and use this as an indication of its expected validity at NNLO.
At NLO we can compare the following three calculations of the invariant mass distribution: 1. The full result, given by the Monte Carlo program MCFM [41].
2. The full threshold expansion, which is equivalent to the full result at leading order in 1 − z.
3. The approximate threshold expansion, which is equivalent to the full result at leading order in 1 − z for the plus distributions and µ-dependent terms, but not in the µindependent piece of the delta-function term. Here we also have the choice of writing the logarithms in the delta-function coefficients in terms of ln(m 2 t /µ 2 ) or ln(M 2 /µ 2 ).
4. The approximate threshold expansion, where the delta-function term is neglected entirely.
To the extent that the second, third, and fourth options agree well with the full result at NLO, it would seem plausible that our NNLO result, which is limited to the third and fourth options, might also agree well with the full NNLO result, which is unknown. We thus investigate these different approximations to the invariant mass distributions for a range of M and µ, and study under what circumstances they show a reasonable agreement. The naive expectation would be that they agree only at very high values of the invariant mass, since then τ = M 2 /s → 1 and the integrand in (6) is needed only in the z → 1 limit. Moreover, the integrals over the plus distributions in (21) are more singular at τ → 1 than the unknown portion of the correction proportional to the delta function. However, due to the rapid fall-off of the PDFs at large x, the distribution at very high τ is difficult to measure. The most interesting region for phenomenology would be from M ∼ 2m t to around 1 TeV at the Tevatron and up to several TeV at the LHC; this region of the invariant mass corresponds to τ < 0.3. In order for our results to be a good approximation in this region, it is necessary that the luminosity functions ff ij (τ /z, µ) fall off so fast for τ /z → 1 that only the largest values of z give significant contributions to the integrand in (6). In that case, a dynamical enhancement of the partonic threshold region may occur, even if τ is not close to unity [37,42,43]. The agreement of the three different approximations at moderate values of τ is thus a measure of whether this dynamical enhancement actually occurs. With these points in mind, we begin by comparing the exact result and the complete threshold expansion at NLO. The results are shown in Figure 1. Here and throughout the analysis we use the MSTW2008 NNLO PDFs [44], take α s (M Z ) = 0.117 with three-loop running in the MS scheme with five active flavors, and use m t = 173.1 GeV in the pole scheme. We use the same set of PDFs at each order of perturbation theory so as to better illustrate the size of the different perturbative contributions to the hard-scattering kernels C ij . We choose µ = 2m t as the central value and vary it between m t and 4m t . For the case of the Tevatron, we see that the threshold expansion is in good agreement with the full result even at values M ∼ 2m t , where τ ≈ 0.03, and gets progressively better at higher values of M. The situation is only slightly worse for the LHC with √ s = 10 TeV, even though the value of τ for a given M is about 25 times smaller. We also note that at LHC energies the threshold expansion matches the exact result better at the higher value of µ = 4m t than at µ = m t . In what follows, we shall also explore the scale choice µ ∼ M. Since such a choice is not practical in the Monte Carlo program MCFM, from now on we shall use the exact threshold expansion to represent the full NLO result, keeping in mind the good agreement seen in Figure 1. We next investigate in more detail the approximate threshold expansions (options 3 and 4 above) at NLO and NNLO. In Figures 2 and 3 we compare the full threshold expansion at NLO with the different approximate expansions at the same order, and also show the approximate expansions at NNLO (these are generated by adding the NNLO corrections to the full threshold expansion at NLO). In Figure 2 we show the differential cross section as a function of µ for two different values of M. We notice that the different approximations at NLO show better agreement with the full threshold expansion for higher values µ ∼ M, especially at LHC energies, where the NLO approximations at lower values of µ differ greatly from the full results. At NNLO, the different approximations in options 3 and 4 do not differ substantially from one another, when one considers their highest and lowest values in the range M/2 < µ < 2M, which is how perturbative uncertainties are typically estimated. In Figure 3 we show results for the same approximations at NLO and NNLO, but this time as a function of M. Given the better agreement at higher µ observed in the previous two figures, we have made the default choice µ = M. We observe that at NLO the approximate threshold expansions recover a significant portion of the exact NLO correction, both at the Tevatron and LHC. It is not unreasonable to expect that the same is true of our approximate NNLO results, although they are clearly no substitute for a complete NNLO computation.
In Figure 4 we show the invariant mass spectrum as a function of M, with bands repre- senting the uncertainty associated with scale variations in the range M/2 < µ < 2M. For the NNLO result, we have displayed only the approximate threshold expansion in which the logarithms in the delta-function coefficient are written in terms of ln(M 2 /µ 2 ); those in the other two NNLO approximations are similar, and very nearly contained within this band. At both the Tevatron and the LHC we note a reduction of scale uncertainty upon including the NNLO corrections, which is the expected benefit of including more terms in the perturbative expansion. However, we must also observe from Figure 2 that at NLO the perturbative uncertainty estimated by varying the scale in the range M < µ < 2M in the approximate threshold expansions is actually smaller than in the complete threshold expansion at both M = 400 GeV and M = 1000 GeV at the LHC. We thus cannot rule out that the very small uncertainty of the NNLO result at small values of M is also an underestimate.

Conclusions
We calculated O(α 4 s ) contributions to the tt invariant mass distribution at hadron colliders. The calculation was based on an effective field-theory approach to the factorization of the perturbative hard-scattering kernels into matrix-valued hard and soft functions in the partonic threshold regionŝ M 2 . At the technical level, it involved calculating the hard and soft functions at NLO, and then determining at NNLO the logarithmic µ-dependence of these functions using the renormalization group and recent results for two-loop anomalous dimensions. Our computations yield exact results for the coefficients of all singular plus distributions and µdependent logarithms in the differential partonic cross section at NNLO, in the limit where the invariant mass of the tt pair approaches the partonic center-of-mass energy. They agree with those previously obtained in [16] using a Mellin-space approach to soft gluon resummation, but go beyond them by uniquely determining the coefficient of the [1/(1 − z)] + distribution, which contains process-dependent two-loop effects related to anomalous dimensions of the soft and hard functions. To obtain the missing scale-independent piece of the delta-function term would require the exact two-loop calculation of these functions.
In Section 3 we performed a numerical study of the tt invariant mass distribution at Tevatron and LHC energies. In both cases, we noted the expected decrease in scale uncertainty after including the NNLO corrections, and also noticeable changes in the central values, especially at higher values of the invariant mass. Given that our NNLO result is still incomplete, even in the threshold region, part of our analysis was focused on determining how well the analogous approximation works at NLO. We observed that at NLO the full result in QCD and its exact threshold expansion agree fairly well. The same was true of the approximate threshold expansion at NLO, where the scale-independent pieces of the delta-function term are neglected. This provides some evidence that the steep fall-off of the PDFs at large x induces a dynamical enhancement of the partonic threshold region, which would mean that the terms we have calculated at NNLO provide a useful approximation to the full result.
Our finding that the NNLO perturbative corrections, which we have calculated in this Letter, have quite an important impact in the region of large invariant masses, as seen in Figure 4, provides a strong motivation for improving this estimate with more sophisticated calculations. We expect that the scale dependence indicated by the width of the various bands in the figure can be reduced significantly by resumming the partonic threshold logarithms to all orders in perturbation theory. Based on the results presented here and in [30,31], such a resummation is now possible at NNLL order. It will be discussed in a forthcoming paper [38]. Ultimately, however, only a complete NNLO calculation of the spectrum will provide us with a complete picture of the true residual uncertainty in the theoretical predictions.