Top Quark Pair Production beyond NNLO

We construct an approximate expression for the total cross section for the production of a heavy quark-antiquark pair in hadronic collisions at next-to-next-to-next-to-leading order (N$^3$LO) in $\alpha_s$. We use a technique which exploits the analyticity of the Mellin space cross section, and the information on its singularity structure coming from large N (soft gluon, Sudakov) and small N (high energy, BFKL) all order resummations, previously introduced and used in the case of Higgs production. We validate our method by comparing to available exact results up to NNLO. We find that N$^3$LO corrections increase the predicted top pair cross section at the LHC by about 4% over the NNLO.


Introduction
The lack of discovery of any new physics signal so far has made of the Large Hadron Collider (LHC) even more a precision machine than it ever was. Because the LHC is a hadron collider, the largest uncertainties are related to nucleon structure, i.e. parton distributions, and higher order QCD corrections. This has consequently led to enormous progress in the computation of higher order corrections to QCD processes in the last several years, recently leading to the publication of first N 3 LO results for a hadron collider process [1].
In Ref. [2] some of us have proposed a general methodology for the determination of approximate expressions for higher-order corrections to QCD processes. The basic idea is that QCD corrections to many important hard processes are known to all orders in the strong coupling in two opposite limits: the high energy limit, in which the available centerof-mass energy is much greater than the invariant mass of the final state, and the soft limit, in which the invariant mass is close to threshold. By performing a Mellin transform, this knowledge can be turned into information on the singularities of the hard (partonic) Mellinspace cross section to any given finite order, viewed as an analytic function of the Mellin complex variable: the soft limit determines the singularity at infinity, and the high-energy limit determines the rightmost pole on the real axis. It is then possible to reconstruct an approximate form of the function by exploiting this knowledge: indeed, knowledge of all the singularities would determine the function completely.
The case of Higgs production, to which the methodology of Ref. [2] was first applied [2][3][4], is particularly simple in many respects: the leading-order process has fixed kinematics (only one scalar colorless particle in the final state), and the total cross section is known [5] to be dominated by its threshold limit. The validation of the methodology which comes from checking that it provides a good approximation to known results in the case of Higgs is thus not necessarily the most compelling.
In this paper we turn our attention to heavy quark production: a process which is rather more complicated in terms of kinematics and color structure, with the goal of performing a more stringent test of our methodology. Next-to-leading (NLO) QCD corrections to this process were computed long ago [6,7], though fully analytic expressions only became available much more recently [8], and were a first necessary step towards the full determination of the NNLO result which was very recently achieved [9][10][11][12]. These results revealed that the physical-space partonic cross section has a particularly complex singularity structure, which raises the question of how this relates to the Mellin-space singularities which are used for the approximation of Ref. [2], and also, that existing attempts at an approximate NNLO determination [13] were rather off the mark.
Here, we will present an approximate determination of the N 3 LO QCD corrections to the total cross section for heavy quark production. In Section 2 we will review the singularity structure of the heavy quark production partonic cross section both in physical and in Mellin space. Specifically, we will discuss the relationship between the total cross section and the invariant mass distribution, to which the methods of Ref. [2] apply, and show that the peculiar singularity structure of the physical-space cross section leads nevertheless to the standard Mellin-space singularities. The rest of our treatment closely follows that of Higgs production: in Sects. 3 and 4 respectively we present the derivation of the soft and high-energy singularities from the respective resummations. They are combined in Section 5 where we present our final result: first, we compare the approximation obtained through our procedure to the known exact result up to NNLO, and then we present our final approximate N 3 LO both at the parton and the hadron level.

Analytic structure of the partonic cross section
The analytic structure of the cross section for the production of a heavy quark pair potentially differs from that for Higgs production discussed in Ref. [2] for two different reasons, thereby potentially hampering, or requiring some adaptation of the methods used in that reference. First, the known [8] fact that the partonic cross section for heavy quark production already at NLO has unphysical singularities in the complex plane of its kinematic variables may suggest that the singularity structure of its Mellin transform also does not follow the general pattern discussed in Ref. [2]. Second, as we shall see in more detail below, the analytic structure discussed in Ref. [2] is generic for partonic differential cross sections, while the case we are presently considering is that of a total cross section. We will see that the first issue is of no concern, while the second requires a careful discussion of the relation of the total heavy quark production cross section to the corresponding invariant mass distribution.
We will address the two issues in turn, after a quick review of standard notation and general results.

Notation and general results
The total cross section for the production of a heavy quark-antiquark pair can be written in factorized form as where ρ h = 4m 2 /s, s is the hadronic center of mass energy, L ij (ρ, µ 2 F ) are parton luminosities, m is the mass of the heavy quark, µ F , µ R are factorization and renormalization scales, and the scaling variable is ρ = 4m 2 /ŝ, withŝ the partonic center of mass energy. To simplify notations, in the following we will not indicate explicitly the dependence of parton luminosities on the factorization scale µ F and that of the partonic cross section on the strong coupling and on µ F , µ R .
The partonic cross sectionsσ ij admit a perturbative expansion in powers of the QCD coupling: where, thanks to the m −2 prefactor, the coefficientsσ (k) ij (ρ) are dimensionless. Eq. (2.1) is in the form of a convolution product, so its Mellin transform factorizes in terms of the Mellin space luminosity and partonic cross section function, defined respectively as If the Mellin transform integral has a finite convergence abscissa, the N -space partonic cross section is an analytic function of the complex variable N , defined by the integral representation Eq. (2.5) to the right of the convergence abscissa, and by analytic continuation elsewhere. Therefore, it is fully determined by the knowledge of its singularities.
In this paper, we will concentrate on the gg partonic channel, which is the most relevant (in the MS scheme) at LHC energies, while we leave the study of other parton subprocesses for future work. For this reason in the following we will drop the parton indices ij, with the understanding that i = j = g in all parton luminosities and partonic cross sections.
The singularity structure of a generic differential partonic cross section is relatively simple. The singularity at infinity is determined by soft gluon radiation, which leads to a growth of the cross section with increasingly high powers of ln N at higher perturbative orders. For finite N , singularities away from the real axis are not allowed, as they would have to come in pairs and would thus lead to oscillatory behaviour of the cross section at high energy. From Regge theory one expects that the rightmost singularity on the real axis is a multiple pole located at N = 1 [14], with further multiple poles along the real axis at N = 0, −1, −2, . . . . This expectation is confirmed by explicit fixed-order calculations.
While knowledge of the residues of all poles is required in order to fully determine the functionσ(m 2 , N ), its behaviour in the region 1 < Re N < ∞ is mostly controlled by the rightmost pole at N = 1, with poles further to the left having an increasingly small impact. Using the saddle point approximation one can show that the hadronic cross section is mostly determined by the behaviour of the partonic cross sectionσ(m 2 , N ) in the vicinity of a single (saddle) value of N on the real axis, to the right of the rightmost singularity [2,15] (see also Sect. 5.1). In Ref. [2] it was suggested that knowledge of the rightmost singularity on the real axis and of the singularity at infinity is sufficient to determine the partonic cross section in the region where the saddle-point is located with reasonable accuracy.
As mentioned in the introduction, however, in the specific case of heavy quark production, we encounter two difficulties. The first, is related to the fact that the physical momentum-space coefficient function has unphysical singularities. The second is due to the fact that the partonic cross sectionsσ ij (m 2 , N ) Eq. (2.6) vanish as N → ∞, rather than growing logarithmically.
We will address both issues in turn. The first issue turns out to be of no concern: the analytic structure in Mellin space is unaffected by the unphysical momentum-space singularities. The second issue instead is due to the fact that the reason why the partonic cross sections discussed in Ref. [2] grow as N → ∞ is that they are distributions, rather than ordinary functions: indeed, elementary properties of Laplace transforms imply that the Mellin transform of an ordinary function, if it exists, must vanish as N → ∞. Now, the integral of a distribution is an ordinary function, so it is clear that the behaviour of Ref. [2] can only hold for partonic cross sections which are sufficiently differential. As we shall see below, it is the invariant-mass distribution which behaves at large N in the way discussed in Ref. [2]. However, by exploiting the relation between total cross section and invariantmass distribution, it is possible to relate their respective large-N behaviors, and define a coefficient function whose singularity structure follows the general pattern discussed above.

Unphysical singularities in momentum space
The perturbative coefficientsσ (k) ij (ρ) display a class of spurious singularities in the complex ρ plane outside the physical region 0 < ρ ≤ 1. However, after Mellin transformation, even in the presence of such spurious singularities in ρ space, the result has the expected analytic structure, as we now show.
As an example, we present the case of NLO corrections of top pair production in the gg channel. In Ref. [8], the complete analytic form ofσ (1) gg (ρ) is computed, and it contains, in addition to the usual (physical) singularities in ρ = 1 and ρ = 0, four further singularities (branch points) located at ρ = −4, ρ = 4, ρ = −1, ρ = − 1 4 (see in particular the functions F i (ρ), i = 1, . . . , 4 of section 4 and the discussion of section 5 of Ref. [8]). For simplicity we focus our attention only on the first spurious singularity; similar considerations hold for all the others.
The function has a logarithmic branch cut starting at ρ = −4, because the integrand has a simple pole in However, the Mellin transform is (2.11) the ρ integral is convergent for Re N > 1, and the result can be analytically continued to the whole complex N plane, except the isolated points N = 1, 0, −1, −2, . . . , where the result has simple poles. This can be seen for example by repeatedly integrating by parts using z N −1 as a differential factor. This is the analytic structure one expects for the Mellin transform of a physical cross section. We conclude that the presences of unphysical singularities in the ρ space does not affect the general structure the partonic cross section in N space. Of course, it could be that the presence of new structures in the partonic cross section at higher perturbative orders affects the numerical size of the residues of singularities, for the same reasons why it makes standard scale-variation estimates of higher order terms unreliable.

Total cross section and invariant-mass distribution
We will first show that the invariant-mass distribution for heavy-quark production has the large N singularity structure discussed in Ref. [2] and dominated by Sudakov radiation, then we will prove that the coefficient function C(N ), defined by factoring out the leading order contribution in the total cross section: exhibits the same Sudakov enhancement at large N . The invariant-mass distribution Σ(m 2 , ξ, z) is defined as which is a function of two dimensionless ratios, due to the presence of an extra energy scale M 2 , the invariant mass of the quark-antiquark pair. We choose to express it as a function of z = M 2 /ŝ, and of ξ = 4m 2 /M 2 = ρ/z, and we insert the factor ofŝ in Eq. (2.14) for later convenience. The total cross section is related to the invariant-mass distribution bŷ Using this last relation, we will show that C(N ) implicitly defined in Eq. (2.13) has the same Sudakov singularities in the large N limit as the invariant mass distribution Σ(m 2 , ξ, N ) in the limit ξ → 1. This limit is called the absolute threshold limit, and it corresponds to the double limit ξ → 1 and z → 1 of Σ m 2 , ξ, z in z space. Indeed, we see from Eq. (2.15) that, when ρ → 1, only the region z → 1 (and hence ξ → 1) contributes to the integral.
The behaviour of Σ(m 2 , ξ, z) for z close to one is governed by Sudakov resummation. In analogy with Drell-Yan or resonant Higgs production, the O(α n s ) perturbative coefficient for Σ(m 2 , ξ, z) is a linear combination of the distributions with 0 ≤ k ≤ 2n − 1, plus contributions which are less singular in the threshold limit z → 1. Correspondingly, the perturbative coefficients of The resummation of these large logarithmic contributions has been performed in recent years, both in momentum space with Soft Collinear Effective Theory (SCET) techniques [16][17][18] and in Mellin space [19][20][21]. In the latter case one finds [20] up to corrections which vanish in the large-N limit. Here, bold symbols are used for matrices in color space. In the case of heavy quark pair production in gluon-gluon fusion, these are 3 × 3 matrices because there are three independent colour configurations [21]. The function G(N ) is given by where A(α s ) [22,23] and D(α s ) [24] have perturbative expansions in powers of α s , which are known up to O α 3 s (see for instance [25][26][27]). The soft anomalous dimension Γ S (ξ, α s ) and the soft (matrix) function S(α s ) originate from soft emissions in the presence of a heavy quark pair. Finally, the matrix function H(ξ, α s ) represents the hard contribution to the cross section. NNLL accuracy is achieved by expanding the cusp anomalous dimension A up to three loops, the characteristic function D and soft anomalous dimension Γ S up to two loops, and the soft and hard functions S, H to one loop. Moreover to predict the term proportional to the delta function δ(1 − z) at O(α 2 s ), we just need the sum of the two loop contributions to the hard and soft functions, H and S, which can be determined by matching with NNLO calculation.
However, what we are interested in is the absolute threshold limit ξ → 1 of Eq. (2.21). We now show that in this limit Eq. (2.21) greatly simplifies. For ξ ∼ 1 we can replace M 2 by 4m 2 in the argument of α s , the difference being suppressed by powers of 1 − ξ. Furthermore, it was noted in Ref. [28] that in this limit the matrix Γ S (ξ, α s ) is diagonal in the singlet-octet basis defined e.g. in Ref. [21]. To order α 2 s it takes the form [28] 24) and the matrix Π 8 is the projector over the octet subspace: We now turn to the soft matrix S(α s ). By a suitable definition of H, it can be chosen to be the unity matrix I at leading order. Its NLO expansion, sufficient to achieve NNLL accuracy, is given by Thus, in the absolute threshold limit, Eq. (2.21) reduces to The calculation of the color trace in Eq. (2.27) is greatly simplified by noting that, for any complex number a, As a consequence, the resummed invariant mass distribution splits into the sum of an octet and a singlet component: The matrix H(ξ, α s ) cannot be expanded in powers of ξ − 1, because it is proportional to the phase-space factor √ 1 − ξ. However, one may factorize the leading order coefficient as in Ref. [28], to obtain The coefficients H (i) (ξ), i ≥ 1 are now analytic around ξ = 1. Thus where h (i) = H (i) (1) are constant matrices. A further simplification arises from the particular structure of the matrix H (0) (ξ, α s ). Indeed, one finds whereσ LO I (m 2 , ξ) are the leading-order total cross sections in each color configuration. Since we have It follows that are independent of N . Eq. (2.30) may be cast in a more familiar form, by reabsorbing the logarithmic terms in the factor 1+C A α s (4m 2 /N 2 )/π in a redefinition of the function Γ S (α s ). This also generates N -independent terms, which can be absorbed in a redefinition of the constantsḡ (i) 8 [28].
Putting everything together, we obtain Finally, recalling the relation Eq. (2.18) between total cross section and invariant-mass distribution, we get is the Mellin transform of the Born level total cross section. By comparing Eq. (2.42) with the definition of the coefficient function Eq. (2.13), we see that has the same singularity structure as the Mellin-transformed invariant-mass distribution in the N → ∞ limit.

Analytic structure of the coefficient function
We now focus on the coefficient function C(N ), implicitly defined in Eq. (2.13). We have shown in the previous section that in the N → ∞ limit C(N ) has the same singularity structure as the Mellin-transformed invariant mass distribution, whose behaviour in turn is determined by soft-gluon emissions. On the other hand, we note that C (N ) has the same leading singularity in N = 1 asσ(m 2 , N ), becauseσ LO (m 2 , N ) is subleading in the high-energy regime. Therefore, we can construct an approximation to C(N ) according to the procedure of Ref. [2], as where C soft contains the terms predicted by Sudakov (soft) resummation, and C h.e. the terms predicted by BFKL (high-energy) resummation. Explicit expressions for both components are given in the following Sections.
3 Large-N contributions

Threshold resummation and analyticity
In this subsection, we will extend the procedure outlined in Ref. [2] for the Higgs production cross section to the coefficient function C(N ) of Eq. (2.44). The cusp anomalous dimension A(α s ) has been computed up to three loops [23], and the colored characteristic anomalous dimensions D I (α s ) are known completely at two loops [29], allowing us to achieve NNLL accuracy. Hence, we are able to determine all terms of order α n s ln m N , with 2n − 3 ≤ m ≤ 2n, in the coefficients of the perturbative expansions of C(N ). The inclusion of α 2 s contribution in the constant termsḡ I enables the extension of our prediction to 2n − 4 ≤ m ≤ 2n. We are thus able to predict all large-N non-vanishing contributions to C(N ) up to O α 2 s , and all logarithmically enhanced contributions except the single log and the constantsḡ The single log at O(α 3 s ), formally a N 3 LL contribution, is produced by the third order of the colored characteristic anomalous dimensions D I . At this order only the singlet is fully known [27], but we do not know the contributions to D I = 0 in the following. The extra uncertainty induced by this assumption on our results will be discussed in Sect. 5.3 below.
As pointed out in Ref. [2,5], the quality of the soft approximation to the full cross section significantly depends on the choice of subleading terms which are included in the resummed result. Since the resummation procedure only fixes the coefficients of logarithmically divergent terms, there is a certain freedom in defining how the soft approximation is constructed, by including contributions that are suppressed in the limit N → ∞. The idea proposed in Ref. [2,3] is to include subleading terms that are known to be present in the exact calculations, in order to preserve the known analytic structure at small N .
In principle, the exponent G I (N ) Eq. (2.39) is ill-defined, because the integration range includes the Landau pole of the strong coupling. This problem is usually avoided by computing the integral in the large-N limit: the resummed N -space result is then well-defined as a function of ln N . The logarithm, however, has an unphysical cut at N = 0, where the cross section should only have poles.
However, as pointed out in Ref. [2], if we are only interested in the expansion of Eq. (2.39) in powers of α s (m 2 ) to finite order, the problem of the Landau pole does not arise. The Mellin transform in Eq. (2.39) may then be computed exactly, provided the integrand is understood to be expanded in powers of α s (m 2 ) to some finite order. The result is a function of N which has the correct logarithmic behaviour at N → ∞, but is free of unphysical cuts on the real negative axis from N = 0.
An explicit calculation yields where D k (N ) are the Mellin transforms of the distributions and the coefficients b n,k,I up to order n = 3 are given in the Appendix A, Eqs. (A.1) and (A.2), and depend respectively on A(α s ) and D I (α s ) only. Explicit forms of the functions D k (N ) are given e.g. in Ref. [2]; they are seen to only have poles on the real axis.
Following the argument of Ref. [2], we now observe that the logarithmic enhancement at threshold arises from the integration over the transverse momentum of the emitted gluons, which has the form where P gg (z) is the gluon-gluon splitting function for z < 1 and 3)). Thus, it is natural to include subleading terms in Eq. (2.39) by restoring the factor of 1/ √ z in the upper integration bound, and by replacing A (1) by the expansion of A g (z) about z = 1 up to some finite order (keeping the full expression of A g (z) is not advisable, because an unphysical singularity in z = 0 would appear; see Ref. [2]).
We therefore replace G I (N ) witĥ which differs from G I (N ) by non-logarithmically enhanced terms. The factor of z in the first terms arises after expansion of A g (z) in powers of 1 − z to first order; indeed, it turns out that An extra power of z simply amounts to a shift of N by one unit in the Mellin-space result. This clearly does not affect the logarithmic behaviour at N → ∞. This extra z factor multiplies also all higher order contributions to the cusp anomalous dimension A(α s ).
The constants c I have been introduced by requiring When expanding in powers of α s as in Eq. (3.1), we can effectively get rid of the constants by simply defining distributions whose Mellin transform differ by 1/N suppressed terms at large N , namelyD In this way we find n,k,I are the same as in Eq. (3.1), and the argument of the Mellin transform of the distributions associated with the cusp term are shifted by one. Explicit expressions for the Mellin transformsD k (N ) are given in Appendix A. Note that here, unlike in Ref. [2], we do not shift the D terms; however, the difference is subleading, and the impact is negligible.
It is important to observe that, unlike in the case of Higgs production, the inclusion of the second term in the expansion Eq. (3.5), does not lead to the full inclusion of all subdominant contributions of the form α n s N −1 ln 2n−1 N . This is because contributions of this order to C(N ) may arise both from soft emission, but also from interference of powers suppressed terms in the leading order cross section. We will use these terms (by turning on and off the shift in N ) as a way to estimate the uncertainty in our procedure.

Coulomb singularities
In Sect. 3.1 we have studied the large-N terms originated by soft gluon emission. In the case of heavy quark pair production, however, there is another class of contributions which, order by order in perturbation theory, do not vanish in the large-N limit, altogether unrelated to soft emission. It was pointed out many years ago [30] that pair interaction dynamics and bound state effects may be relevant in the threshold regime. This multiple exchange of "Coulomb gluons" between the two heavy quarks in the final state leads to corrections of order (α s /β) m with β = √ 1 − ρ, which are usually referred to as Coulomb singularities. For m large enough, such contributions may compete with, and even dominate over, Sudakov logarithms in the threshold limit. It turns out, however, that these contributions can be resummed to all orders (see for example Ref. [31] and references therein), and the result of resummation vanishes in the large-N limit. Nevertheless, these contributions must be taken into account in the absolute threshold limit ρ → 1.
Coulomb terms have been studied in Refs. [29][30][31][32][33][34]. In particular, it was shown in Ref. [32] that Coulomb singularities and soft singularities factorize in Mellin space in the N → ∞ limit and in the singlet-octet basis:  It can be shown that the effect of Coulomb terms on the physical cross section is very small. In Ref. [29] the impact of the contribution J (2) I (N ) on the hadronic cross section was estimated to be about 0.5% of total NNLO correction, and even less at higher orders. Similar conclusion can be drawn from inspection of Fig. 1, where we compare, in N space, the coefficient function Eq. (3.9) with and without the Coulomb factors J I (N, α s ), both at NLO and NNLO. The effect at N ∼ 2, which is the relevant region for Mellin inversion (see Section 5.1), is very small.
A complete resummed expression of Coulomb contributions has been obtained in the context of potential non-relativistic QCD (pNRQCD) [36][37][38]. This calculation can be used to test whether the pattern observed at NLO and NNLO continues at higher orders. In Fig. 2 we present the same comparison as in Fig. 1 at N 3 LO, with J

Final prescription for the soft contribution
We conclude this section by giving our final prediction for the soft-emission contribution to the total cross section to N 3 LO. We define where we have made explicit the dependence on α s = α s (m 2 ) and it is understood that the exponentials should be expanded in powers of α s up to order α 3 s . We will take the average between the two approximations, as the central value of our prediction, and the difference as an estimate of the uncertainty of our procedure. Hence, our soft approximation to the coefficient function will be finally given by We also consider the result which corresponds to the standard resummation, as given to NLL in Ref. [33], and extended to NNLL in Ref. [12] (and in the associate public code top++), which is expressed as a function of positive powers of ln N , which we will call N -soft, following Ref. [2]. This is given by

Small-N contributions
In order to extract the leading small-N singularity of the partonic cross section, we will use the so-called high-energy or k T factorization technique, first described in Ref. [40] for the total cross section, and more recently extended to rapidity distributions [41]. We follow the resummation procedure developed by Altarelli, Ball and Forte (ABF) [42]. For a detailed derivation of resummation, which affects both coefficient functions and evolution of the parton densities, we refer the reader to the original literature (e.g. [40,[42][43][44]); here we only summarize the main results which are relevant for our discussion.
In the k T factorization formalism, small-N singularities are obtained by computing the leading-order partonic cross section for the relevant process, with off-shell incoming gluons. We therefore define an off-shell partonic cross section α 2 s m 2σoff-shell (ρ, ξ 1 , ξ 2 ), which is a function of the scaling variable of the process ρ, and of the transverse momenta of the two off-shell gluons: ξ i = k 2 Ti /m 2 . Resummed results can be obtained through the determination of the so-called impact factor The prefactor R(M 1 )R(M 2 ) accounts for factorization scheme dependence [43]; in the MS scheme The impact factor Eq. (4.1) for the production of a heavy quark pair was computed in Ref. [45]. Here, we are interested in its expansions in powers of M 1 , M 2 , in the vicinity of N = 1: The coefficients h i 1 ,i 2 relevant at N 3 LO can be easily obtained by performing the expansions of the previous formula and they are given in the Appendix B, Eq. (B.1). Following Refs. [42,46], the resummation is performed by identifying the Mellin variable M i with the resummed DGLAP anomalous dimension (together with sub-leading running-coupling effects): By expanding the anomalous dimension to fixed perturbative order we have all the ingredients to construct our N → 1 approximation, which is simply found by substituting the expansion (4.7) in Eq.  4). Now we define a resummed coefficient function by factoring out the leading order contribution, evaluated in the high-energy limit, i.e. N = 1 where, in the second equality, we use the fact that and we defineh We have omitted the dependence of the coefficients on µ 2 F for simplicity. Since we are going to combine the small-N approximation with the large-N approximation to obtain an estimate of the full coefficient function, we require that the two limiting behaviors do not interfere with each other. In particular, we require that the high-energy contribution vanishes when N → ∞. Manifestly, Eq.  to all orders in α s . As pointed in Ref. [2], this choice for the subtraction is a compromise between the contrasting goals of not changing the small-N singularities structure and of damping strongly enough as N increases. In momentum space the subtraction of Eq. (4.11) corresponds to damping the ρ → 1 behaviour of the partonic cross section through a multiplicative factor (1 − ρ) 2 . One additional modification is needed to define our prediction for the cross section in the high-energy regime. We note that the anomalous dimensions vanishes at N = 2 due to momentum conservation. This implies that the small-N approximation Eq. (4.8) vanishes in N = 2. This value of N marks the transition between the small-N approximation (not accurate if N 2) and the large-N approximation (not accurate when N 2).
However, Eq. (4.8), differs from 0 in N = 2 because the small-N limit of the anomalous dimension γ + , Eq. (4.7), contains only the leading and the next-to-leading singularities in N = 1, and not a full fixed order expression. Momentum conservation can be enforced [42] directly in Eq. (4.8), by adding to C ABF (N ) a function f mom ∝ 1/N . Following Ref. [2], we construct our small-N momentum-conserving approximation as where the constant k mom (α s ) is fixed by requiring that C h.e (2) = 0 to all orders in the perturbative expansion. The enforcement of momentum conservation in our small-N approximation allows us to estimate the contribution of subleading poles in N = 0, −1, −2, . . . , which are not controlled by the high-energy approximation. Indeed, the full partonic cross section does not vanish at N = 2: only the contribution to it driven by hard radiation from external legs does. Hence we can estimate subleading small-N effects by allowing the partonic cross section to deviate slightly from zero at N = 2 due to contributions from subleading poles. We thus fix k mom to be k mom (α s ) = C ABF-sub (2) ± ∆k mom (α s ), (4.14) and we take the variation ∆k mom (α s ) to reach as its maximum value 15% of the size of soft contribution C A-soft (N ), Eq. (3.13) at N = 2. This is a somewhat more conservative estimate of the uncertainty in comparison to the corresponding one adopted in Ref. [2], due to the fact that we expect the contribution from subleading poles to be potentially somewhat larger in this case, based on the behaviour of known orders, and possibly also due to the issues discussed in Sect. 2.2. Namely, we choose in Eq. (4.14), In conclusion, this means that the small N contribution, rather than being completely switched off at N = 2, is small at this point, and gets switched off somewhere in its vicinity.

Approximate cross section up to N 3 LO
We are now ready to present our results for top pair production cross section at N 3 LO, by using Eq. (2.45) to combine the large-N terms Eq. (3.15) and the small-N terms Eq. (4.13). We first recall how N -space parton-level results can be related to physical hadron-level results using the saddle point methods, and then we present the parton and hadron level results in turn.

Saddle point approximation of Mellin inversion integrals
The physical hadronic cross section can be related to the underlying partonic cross section by viewing the former as the inverse Mellin transform of its factorized expression Eq. where c is to the right of all singularities of the integrand. It can be shown on general grounds that f (N ) has a unique minimum at N = N 0 on the real N axis, which allows us to estimate the integral by the saddle-point technique.
The position N 0 of the saddle typically depends very weakly on the partonic cross section, and is mainly determined by the parton luminosity; the saddle-point approximation, with the inclusion of quadratic fluctuations around the saddle, turns out to be generally quite accurate [47]. It is then possible to infer properties of the hadronic cross section from the behaviour of the partonic cross section at the value of N which corresponds to the saddle point for given hadronic kinematics.
The position of the saddle point N 0 as a function of the collider energy √ s, for m = m t = 173.37 GeV, computed using NNPDF 3.0 [48] parton distributions, is shown in Fig. 3. The value of the saddle turns out to be around N = 2.5 at LHC energies. The value of N 0 is a very slowly decreasing function of the total energy.

Parton-level results
We can now combine the results of Sections 3 and 4 to obtain an estimate of the coefficients of the perturbative expansion for the total production cross section of heavy quark pairs. We first test our procedure against exact fixed-order calculations, which are available at NLO [6,8] and NNLO [12]. At NLO we specifically compare to the fit to the analytic result of Ref. [8] presented in Ref. [49] . In Fig. 4 we show our approximation to the NLO and NNLO in Mellin space (labeled approx), compared to the exact results of Ref. [12,49] soft terms (green band), and the uncertainty Eqs. (4.14), (4.15) on the high-energy terms (blue band).
The agreement is excellent at NLO in the whole range displayed in Fig. 4. At NNLO there is a slight discrepancy in the region between 1 < N < 1.6. This is due to the fact we are only including the LL contribution in this region, while it was noted in Ref. [8,12] that NLL contributions close to N = 1 for top pair production are sizable. This region of N , however, would only be relevant for very high energy colliders ( √ s 100 TeV).
At the N values which are relevant for LHC energies, the agreement between the approximate and exact results is excellent; for lower collider energy, as the saddle point moves towards larger values, the uncertainty on the approximate prediction is smaller.
We now turn to the N 3 LO contribution. We recall that the constantsḡ 1 , are not known. They will be set to zero for the time being; the impact of this missing information on our prediction for the physical cross section will be discussed in Section. 5.4 below. The function C (3) approx (N ) is plotted in Fig. 5, together with the soft and the high-energy contributions that go into it according to Eq. (2.45).

Hadron-level results
We come now to the hadronic cross sections. It is convenient to define the gluon-channel K-factors where α s = α s (m 2 ), σ gg (m 2 , ρ h ) is the contribution to the hadron-level cross section Eq. (2.1) from the gluon-gluon subprocess, and σ LO gg (m 2 , ρ h ) the corresponding leading-order approximation. All results will be obtained using the partonic cross sections of Sect. 5.2, with factorization and renormalization scale µ R = µ F = m, and NNPDF 3.0 NNLO parton distribution functions [48], with α s (M 2 Z ) = 0.118 and n f = 5. Scale uncertainties on our final results will be discussed in Sect. 5.4 below.

√s[TeV]
exact N-soft approx Figure 6. Comparison between the exact result and our approximation for the NLO contribution to the K-factor Eq. (5.2) from the gluon channel, plotted as a function of the collider energy √ s. The result obtained expanding the standard NNLL resummation (N -soft) is also shown.
The NLO and NNLO K-factors at a pp collider are shown in Figs. 6 and 7, respectively, as functions of the center-of-mass energy √ s, and compared to the exact results of Ref. [12,49]. We also show the values obtained by expanding out to the given order the standard resummed result of Refs. [12,33], i.e. using the N -soft approximation Eq. (3.16). The main result of this work, namely the N 3 LO contribution to the K-factor in the gluon-gluon channel as a function of the collider energy is shown in Fig. 8. Numerical results for the K-factors at LHC energies are collected in Tab. 1.

√s[TeV]
N-soft approx We note that our approximation reproduces the exact result within the estimated uncertainty both at NLO and at NNLO, in the whole energy range displayed in the plots. In comparison to the case of Higgs production, the uncertainty is larger, both because threshold resummation is known to a lower accuracy, and also because, as already mentioned, we have less control over ln k N N terms. Our result is seen to differ by varying amounts at each order from the N -soft one, simply obtained by expanding out the resummed result. At NLO and NNLO the origin of this difference is twofold: first, our approximation also includes the high-energy contribution and second, the functional form of the contribution obtained by expanding out the N -soft result differs from that adopted in our approximation, as explained in Sect. 3.3. At N 3 LO two further differences are, first, that our result includes the single logarithmic term, even if its coefficient is only partially known, which  is absent in a NNLL resummation. However, the numerical impact of this contribution is small. Second, that the function of α s which multiplies the resummed exponent is given by g I (α s ) in the N -soft resummed result, and byḡ I (α s ) when constructing our approximation, the relation between the two being given by Eq. (3.18). We will specifically discuss the impact of this choice in the next Section (see in particular Tab. 2).

N 3 LO top pair production cross section at LHC
We collect now final results for top pair production at the LHC energies and its uncertainty. Our prediction for the total cross section with µ F = µ R = m is LHC7: σ N 3 LO approx ρ h , m 2 = 177.43 ± 2.99 + 0.10ḡ (3) − 0.10 D The constantḡ (3) is not known. The difference between the coefficients D parametrizes the missing information about the single log at N 3 LO, as discussed in Section 3.1. Eq. (5.3) is obtained using exact expressions for the NLO [8] and NNLO [9][10][11][12], combined with our approximation for the N 3 LO in the gluon channel. We now consider various sources of uncertainty. First, we consider the uncertainty related to missing coefficients. In Tab. 2, we list the values of the coefficientsḡ (i) , and of D (i) 1 for the two known orders, as well as the coefficients g (i) defined using the analogue of Eq. (5.4) but for the expansion coefficients of g I (α s ) Eq. (3.18) which is used in the standard resummed results of Refs. [12,33], and the difference between the two, which is known to a higher perturbative order. Based on these values, we note that the perturbative behaviour of the coefficientsḡ (i) appears to be rather more stable than that of the coefficients g (i) , and we reasonably expect the unknown coefficients,ḡ (3) and D 1 , to be both of O(1). We conservatively estimate the uncertainty related to these unknown coefficients asḡ (3) = 0 ± 5, D which we add in quadrature to the approximation uncertainty of Eq. (5.3). We also note in passing that the N -soft approximation assumes g (3) = 0, while our estimate corresponds to a value of order ten for this constant, thereby partly explaining the difference observed in Fig. 8 between our result and the N -soft approximation, as discussed in the previous Section. A second source of uncertainty is due to the fact that we only predict the contribution of the gluon channel. In Tab. 3, we show the contributions to the top-pair production cross section coming from different partonic channels at LO, NLO, NNLO, for the specific case of LHC at 8 TeV. We note that gluon fusion contribution is always dominant, and largely dominant at high perturbative orders. This remains true for other LHC energies (while at the Tevatron the qq component becomes dominant). Therefore, we expect that our estimate for the N 3 LO contribution based on the gg channel only is only mildly affected by the lack of knowledge of the other channels. However, in order to take into account the uncertainty due to the missing partonic channels, we consider the dependence of our result on the factorization scale. Indeed, all the scale dependent terms at N 3 LO are available (for all channels), but the inclusion of all of them would consistently make the residual scale dependence of relative order α 4 s . On the other hand, if we include µ F dependent terms only in the gg channel, the residual scale dependence is of order α 3 s , since it misses the compensation between channels. The difference between these two ways of varying the  scale can be thus taken as an estimate of the size of the contribution from the missing quark channels.
Finally, uncertainties related to missing higher-order terms can be estimated by varying the renormalization and factorization scales in the usual way. The dependence of the cross section at LO, NLO, NNLO and approximate N 3 LO on the factorization and renormalization scales, is shown in a wide range in Fig. 9; the NNLO+NNLL result of Refs [12,33] is also shown. The factorization scale variation at N 3 LO with µ F is shown both retaining contributions from all channels, and from the gg channel only. The factorization scale dependence is rather mild already at NNLO, and even milder at approximate N 3 LO when all channels are included. When only the gg channel is included, a stronger scale dependence is observed, particularly at the extremes of the range, where sizable contributions from the missing channels are generated. The renormalization scale dependence is somewhat stronger than the factorization scale dependence, but at N 3 LO it has flattened out almost completely, thereby indicating a good perturbative convergence. It is interesting to observe that the NNLO+NNLL result has a milder scale dependence than the NNLO, but still a stronger scale dependence than the approximate N 3 LO. This is what one may expect based on the observation that effectively, because the process is far from threshold, the resummation is providing some approximation mostly to N 3 LO, which is however less complete than the approximation which is constructed here.
Our where the "channels" uncertainty has been computed as (± half) the difference between the µ F scale variation evaluated with only the gg channel or with all the channels (NNNLO gg and NNNLO curves in Fig. 9), in the range m/2 ≤ µ F ≤ 2m with µ R = m. The "scales" uncertainty is instead obtained through a canonical seven-point variation, namely m/2 ≤ µ R , µ F ≤ 2m with 1/2 ≤ µ R /µ F ≤ 2, computed with all channels. We observe that the approximation uncertainties, though conservatively estimated, are rather smaller than the scale uncertainty, and in fact adding in quadrature scale and approximation uncertainties we end up with an overall theoretical uncertainty on our N 3 LO result of 3.5%, not much larger than the scale uncertainty itself. This uncertainties can be compared to the PDF uncertainty, which at the LHC √ s = 13 TeV (with NNPDF 3.0 PDFs) is of order 2%. Additional uncertainties come from the values of α s and m t : see Ref. [50] for a more detailed discussion.
We observe that the uncertainty due to scale variations at NNLO is about 5% at the collider energies we are considering, which is larger than the overall uncertainty on our N 3 LO estimate even accounting for approximation uncertainties. The inclusion of our approximate N 3 LO contribution appears thus to be advantageous, and it leads to a decrease in theoretical uncertainty which is of the size one would expect when going ftom NNLO to N 3 LO.
We may compare our approximate N 3 LO result to that which would be obtained using the N -soft result shown in Fig. 8. The latter leads to a N 3 LO contribution which corresponds to an increase of 3%, 2.8%, 2.3%, 2.4% in comparison to the NNLO at LHC √ s = 7, 8, 13, 14 TeV, respectively. This is rather lower than our approximate result, which corresponds to an increase of 4.3%, 4.5%, 4.2%, 4.3% respectively. The reasons for this have been discussed in the previous section. As discussed above, and as it is apparent from Fig. 9, the scale uncertainty on this result is somewhat smaller than that on the NNLO result, but larger than that on our result.
Finally, we note that an approximate N 3 LO result is also presentes in Ref. [51]. This result is obtained by considering logarithmic terms enhanced at partonic threshold in the differential cross section, inclusive in one of the two heavy quarks produced. The result is then integrated to obtain the inclusive cross section. The enhancements over the NNLO results were found to be 4%, 3.6%, 2.7%, 2.6% at √ s = 7, 8, 13, 14 TeV respectively. We leave a more detailed comparison to the results of Ref. [51], as well as to the resummed result obtained in SCET (e.g. Refs. [32,35] and Refs. [16][17][18]), for future work.

Conclusions
We have constructed an approximate expression for the N 3 LO contribution to the production cross section of a heavy-quark pair at hadron-hadron colliders. We have focused on the gluon-gluon initiated subprocess, which gives the largest contribution at the LHC.
We have obtained our result by extending the method developed by some of us for the case of Higgs production in gluon-gluon fusion [2], and based on reconstructing the Mellin-space partonic cross section from its known singularities. Heavy-quark production is a process with a rather more complicated kinematic and color structure compared to inclusive Higgs production. Furthermore, fixed-order coefficient function for this process are known to have a rather non-trivial singularity structure in physical space, at first sight unrelated to the threshold or high-energy limits. Therefore, application to this case of the the technique suggested in Ref. [2] and applied there to Higgs production in gluon fusion provides a rather stringent test of this methodology. In this study we have shown that the method of Ref. [2] provides excellent approximations to known results up to NNLO, thereby validating he methodology.
Having established the reliability of the methodology even in this more subtle case, we have used it to produce a N 3 LO approximate partonic cross section for heavy quark production. We have then focused on the tt cross section, which is of great interest at the LHC. We have found that the approximate N 3 LO correction amounts to an increase in comparison to the NNLO prediction of 4.3%, 4.5%, 4.2%, 4.3%, for pp collisions at √ s = 7, 8, 13, 14 TeV, respectively. Inclusion of this correction reduces the scale uncertainty to 3%, with a combined uncertainty on the approximation itself of comparable size or smaller. Our final overall conservatively estimated uncertainty is thus somewhat smaller than the uncertainty on the exactly known NNLO result, and inclusion of our approximate result appears to be advantageous. Distributions in the Higgs Boson Era".

A Coefficients in the large-N contribution
The coefficients b 8 is unknown. To complete our resummed formula for the soft emission, we need explicit expressions for the functionsD k (N ), Mellin transforms of the distributionŝ where the plus distribution is defined by The full calculation is presented for example in Ref. [2,52]. Here we give the final result: where we have defined The Coulomb functions J I (N, α s ) are computed by taking a Mellin transform of the resummed momentum-space results obtained in the context of pNRQCD [32,34]. For more details about the procedure of this Mellin transformation, see Ref. [39]. The explicit expression for the Coulomb functions, up to N 3 LO, is the following: 93 − 10n f − 432 β 0 ζ 3 π + 36πβ 0 (γ E − 2 ln 2 + 2 (ψ 0 (2N ) − ψ 0 (N ))) − C 2 F π 108 + C 2 F (C A + 2C F ) 93 − 10n f − 432 β 0 ζ 3 π + 36πβ 0 (γ E − 2 ln 2 + 2 (ψ 0 (2N ) − ψ 0 (N ))) with B (a, b) = Γ(a)Γ(b) Γ(a+b) the Beta function. We now turn to the functionsḡ I (α s ). They must be computed by matching with exact results, which are available up to NNLO. Moreover, as pointed out in Ref. [12], at NNLO we can only infer from the exact result the sumḡ (2) =ḡ (2) 1 +ḡ (2) 8 , but not the individual contributionsḡ (2) I . The latter have been estimated by the same procedure adopted in Ref. [12], givingḡ I (α s ) = 1 +ḡ It can be checked that the uncertainty due to this additional guess is completely negligible in comparison to our error bands. For simplicity we do not give the explicit factorization scale dependence ofḡ(α s ).

B Coefficients in the small-N contribution
The only ingredients which are needed for the high-energy contribution are the coefficients h i 1 ,i 2 of the expansion of the impact factor, up to N 3 LO: The coefficientsh i 1 ,i 2 are simply obtained dividing each h i 1 ,i 2 by h 0,0 . The factorization scale dependence can be easily restored following the procedure explained for example in Ref. [2]. The leading and next-to-leading singular contributions to the anomalous dimension γ + are