Soft-gluon and Coulomb corrections to hadronic top-quark pair production beyond NNLO

We construct a resummation at partial next-to-next-to-next-to-leading logarithmic accuracy for hadronic top-quark pair production near partonic threshold, including simultaneously soft-gluon and Coulomb corrections, and use this result to obtain approximate next-to-next-to-next-to-leading order predictions for the total top-quark pair-production cross section at the LHC. We generalize a required one-loop potential in non-relativistic QCD to the colour-octet case and estimate the remaining unknown two-loop potentials and three-loop anomalous dimensions. We obtain a moderate correction of $1.5\%$ relative to the next-to-next-to-leading order prediction and observe a reduction of the perturbative uncertainty below $\pm 5\%$.


Introduction
The total top-quark cross section at hadron colliders is a key observable of the standard model that serves both as a probe of the heaviest known elementary particle as well as a benchmark process of perturbative Quantum Chromodynamics (QCD). Comparison of measurements at the Tevatron and LHC with precise theoretical predictions allows to measure the top mass in a well-defined scheme, where recent measurements have reached an accuracy of below 2 GeV [1,2]. Further applications include determinations of the strong coupling constant α s [3,4] and global PDF fits, see e.g. [5][6][7]. The experimental results from the LHC experiments have now reached an impressive accuracy, with an uncertainty of better than ±4% (see [8] for a recent review), challenging the best theoretical predictions based on next-to-next-to-leading order (NNLO) QCD [9][10][11][12] supplemented with next-tonext-to-leading logarithmic (NNLL) soft-gluon resummation [13][14][15][16][17][18].
Since a full computation of top-quark pair production at N 3 LO accuracy in QCD is currently out of reach, attempts to reduce the perturbative uncertainties further rely on resummation methods. In this article we explore prospects for a combined resummation of soft-gluon and Coulomb gluon effects at N 3 LL accuracy in the framework of [13,19] and provide expressions for the total cross section at N 3 LO in the partonic threshold limit β = 1 − 4m 2 t /ŝ → 0. Coulomb corrections arise from the exchange of gluons between the slowly moving top quarks and lead to corrections of the form (α s /β) k , which can be resummed to all orders using Green-function methods in non-relativistic QCD (NRQCD). Counting both the logarithmic soft-gluon corrections α s ln β and the Coulomb corrections α s /β arising at each order in perturbation theory as quantities of order one, a combined resummation of these effects rearranges the perturbative series of the partonic cross section into the schematic form [13] σ ∝σ (0) We also indicated a modified counting N n LL', where the fixed-order corrections in the second line are included at one order higher than in the "unprimed" counting. The combined resummation of soft and Coulomb corrections at NNLL accuracy [15] was implemented in the program topixs [18], whose current version includes the matching to the complete NNLO corrections [9][10][11][12]. NNLO+NNLL soft-gluon resummation with the Mellintransform method [16] with a fixed-order treatment of Coulomb corrections was implemented in the program top++ [17], which further includes the O(α 2 s ) constant terms in the resummation that are part of the NNLL' corrections in (1.1). Other NNLL resummations based on pair-invariant mass or single-particle-inclusive observables [20][21][22] also do not include the Coulomb corrections. In contrast, a combination of resummed Coulomb corrections with fixed-order soft corrections was performed in [23][24][25].
Top quarks are not dominantly produced at threshold at the LHC, so thresholdenhanced corrections do not necessarily require resummation and fixed-order perturbation theory is expected to be adequate for the total cross section. 1 To the extent that these terms nevertheless constitute a significant part of the full higher-order corrections, it is therefore justified to expand a resummed prediction to a fixed order. In this way, N n LL resummation predicts all N n LO terms that become singular for β → 0, while a constant correction relative to the Born cross sectionσ (0) is further included at N n LL'. Experience at NNLO [26][27][28][29] indicates that the singular terms provide a feasible approximation to the full result if corrections due to possibly sizeable non-singular terms are estimated in a sufficiently conservative way. For instance, the prediction σ +4.7 −4.7 pb was obtained in [15] using the partonic cross sections calculated in [27], where the first error refers to the scale uncertainty and the second one estimates terms beyond the threshold approximation. This result is consistent with the full NNLO calculation σ NNLO tt (7 TeV) = 167.0 +6.7 −10.7 pb [12] and improves the accuracy compared to the NLO prediction σ NLO tt (7 TeV) = 158.0 +19.5 −21.2 pb. This motivates the construction of approximate N 3 LO corrections, which should be based on N 3 LL accuracy to obtain all threshold-enhanced terms.
At present, a complete N 3 LL resummation for top-quark production is not feasible, since some three-and four-loop coefficients in the resummation function g 3 are not known for the colour-octet case. In addition, starting from NNLL logarithmic corrections arise also from Coulomb corrections, which are governed by renormalization-group equations in NRQCD [30,31], with anomalous dimensions only fully known for colour-singlet states. 2 Furthermore, the power corrections ∼ β 2 in the curly brackets in (1.1) must be controlled in order to achieve N 3 LL accuracy according to the counting α s ∼ β. This makes a complete resummation at this accuracy conceptually and technically challenging.
In this paper, we construct a partial N 3 LL approximation by including all currently known ingredients of N 3 LL soft-gluon resummation [32][33][34] and the Coulomb corrections from a recent calculation of e − e + → tt at N 3 LO in potential NRQCD [35][36][37]. We further include several sources of power-suppressed corrections, i.e. contributions of P-wave production channels and the combination of Coulomb corrections and a so-called nextto-eikonal logarithm. In this way we obtain all threshold enhanced N 3 LO terms for the colour-singlet state and are in a position to estimate the uncertainty due to missing ingredients for colour-octet states. We also include an O(α 3 s ) contribution to the Coulomb corrections that does not follow from the straightforward expansion of the well-known Sommerfeld factor for stable particles [38].
We compare our results to other recent works on N 3 LO effects in top-quark pair production. In [15] the expansion of the NNLL cross section was used, which is not sufficient to predict terms of the form α 3 s {ln 2,1 β, 1/β} exactly. In the momentum-space resummation approach [32,39,40] used in [15], this incomplete knowledge manifests itself in a residual dependence on unphysical hard and soft scales which results in a sizeable uncer-tainty. A prediction based on NNLL resummation in one-particle inclusive kinematics was made in [41]. No attempt is made to estimate the systematic uncertainties of this approximation. In [42] threshold resummation for the invariant-mass distribution is combined with information about the large-x limit. However, only the gluon-initiated subprocess is considered.
This paper is organized as follows: In Section 2 we outline the framework of our calculation and discuss consequences of collinear factorization and renormalization-group invariance for N 3 LO corrections. In Section 3 we discuss in detail the inputs of the partial N 3 LL resummation while the approximate N 3 LO results for partonic cross sections and predictions for the LHC are presented in Section 4. Explicit results for the renormalization group evolution of the two-loop soft and hard functions are given in Appendix A. The generalization of some contributions to the potential corrections to the required spin and colour states are discussed in Appendix B. Explicit results for the factorization-scale dependent contributions to the partonic cross sections are collected in Appendix C.

Factorization and resummation framework
In this section we collect results for the factorization and renormalization scale dependence of N 3 LO corrections, discuss the colour and spin states of top-quark production near threshold, and outline the resummation formalism used for the combined soft and Coulomb gluon resummation at NNLL [13,15,18,19].

Setup of the perturbative calculation
The total hadronic cross section for the production of a tt + X final state in collisions of hadrons N 1,2 with centre-of-mass energy s is obtained from the convolution of the partonic cross section with the parton luminosity, where the latter is defined in terms of the parton distributions functions (PDFs) The perturbative expansion of the partonic cross section in the strong coupling constant is conveniently expressed in terms of corrections relative to the Born cross section, pp are known exactly for all partonic channels [9][10][11][12]. The numerical predictions have been parameterized by fitting functions and implemented in the most recent versions of the programs top++ [17], HATHOR [43], and topixs [18]. The goal of this paper is to predict the partonic cross sections in the leading production channels pp = gg, qq up to O(α 3 s ) in the threshold limitŝ → 4m 2 t . The dependence of the cross section (2.3) on the renormalization scale µ r can be reconstructed from a calculation of the correctionsσ (n) pp (ŝ, m t µ f ) for µ f = µ r by re-expressing the result as an expansion in α s (µ r ). In this way, the renormalization-scale dependence of the N 3 LO cross section is obtained aŝ In a strict threshold expansion, the limitŝ → 4m 2 t of the partonic cross sections on the right-hand side is used and the terms in the last line are dropped, since they contribute to the terms O(β 0 ) relative to the Born cross section.
The factorization-scale dependence of the cross section can be obtained from results at lower perturbative order by exploiting the known factorization-scale dependence of the PDFs, which implies the evolution equation of the partonic cross section where P p/p (x) are the Altarelli-Parisi splitting functions for the splitting of a parton p into a partonp. In the threshold limit it is consistent to use the x → 1 limit of the splitting functions for a parton p in the colour representation r, 3 which is given in terms of the light-like cusp-anomalous dimension Γ r cusp and a subleading anomalous dimension γ φ,r The anomalous dimensions are all known at least up to three-loop level [44] and summarized in the conventions used here in [32,45].
In the higher-order corrections to the cross section, the dependence on the factorization scale can be made explicit by the decomposition Here some care has to be taken since subleading terms in the x → 1 limit can be enhanced by the Coulomb corrections. In can be checked, however, that the leading correction to the N 3 LO cross section from the O(1 − x) term in (2.6) is of order α 3 s β 2 , and therefore beyond N 3 LL.
with the so-called scaling functions f (n,m) pp and the variable ρ = 4m 2 t s . For S-wave production channels withσ (0) ∝ α 2 s (µ f )β it is useful to define modified scaling functions 4 g (n,m) pp In the threshold limit ρ → 1 we obtain the results for the N 3 LO scaling functions where the threshold limit of the lower-order scaling functions is consistently used on the right-hand side. The convolutions are defined as and the conventions for the coefficients β n and P

Top-quark production channels near partonic threshold
The resummation of soft-gluon and Coulomb corrections requires to decompose the partonic production cross sections into contributions with definite colour and spin states of the top-quark pair, σ (0),R, 2S+1 L J pp . For the colour representations, the familiar decomposition of the tensor product 3 ⊗3 = 1 + 8 of the SU (3) representations is used, where for gluon initial states the symmetric and antisymmetric colour channels with respect to gluon exchange, 8 s and 8 a , are distinguished. For the spin states, the spectroscopic notation for orbital angular momentum L = S, P, . . . , spin S and total angular momentum J is used.
In the quark-antiquark induced production channel, the top-antitop pair is dominantly produced in a colour-octet 3 S 1 state, while the dominant production channel in gluon fusion is given by a colour-singlet or symmetric octet 1 S 0 state. The threshold limit of the leading order S-wave cross sections is given by These definitions are related to [27] by g (n,m) pp = 2 n+m s (n,m) pp while our convention for the splitting function P (n) is 2 n+2 times the one used there.
where we have left an overall factor 4m 2 t /ŝ unexpanded. The spin label in these expressions will be dropped if no confusion can arise. For the counting α s ∼ β used in (1.1) also subleading terms of O(β 2 ) in the threshold expansion of the Born cross section must be taken into account at N 3 LL accuracy. For the S-wave production processes (2.13), these terms are given by a correction factor 1 − 1 3 β 2 . At the same order, there are contributions to the total cross section from P-wave production channels 3 P 0 , and 3 P 2 (see e.g. [46]) with the contributions to the total production cross section Since the second and third Coulomb corrections differ for S-wave and P-wave production channels (see e.g. [47]) it is necessary to distinguish the different angular momentum states contributing to the subleading Born contributions.
In addition to the colour-singlet and symmetric octet channels, there is also a kinematically suppressed contribution from the antisymmetric colour-octet channel As mentioned in [34] and discussed in the context of the violation of the Landau-Yang theorem in QCD [48,49], this suppression is not the signal of P-wave production but the result of an accidental cancellation in the Born S-wave matrix element.

Combined soft and Coulomb resummation at NNLL
Threshold resummation of soft logarithms [50,51] was first established for top-quark pair production at NLL [52,53] and more recently at NNLL for the total cross section [13][14][15][16] and pair-invariant mass or single-particle-inclusive observables [20][21][22]. The combination of soft-gluon and Coulomb-gluon resummation up to NNLL in the combined counting (1.1) was shown in [13,19] for pairs of heavy coloured particles produced in an S-wave state. This method was applied to top-quark production in [15,18]. The basis for the joint soft and Coulomb resummation is the factorization of the total partonic production cross section into a potential function J, a hard function H, and a soft function W [19]: which was derived using the leading-power Lagrangians of soft-collinear effective theory (SCET) and potential non-relativistic QCD (PNRQCD). Here E = √ŝ − 2m t is the energy relative to the production threshold. The sum is over the colour representations R of the final state top-pair system, i.e. the colour-singlet and octet states, and the total spin S = 0, 1. The index i denotes the colour basis for the hard scattering pp → tt [13]. The leading-order hard function for the S-wave scattering process pp → (tt) 2S+1 S J R for the colour state R and spin state S is related to the partonic Born cross sections in the corresponding channel at threshold (2.13) according tô In (2.16) we have made the conventional prefactor 4m 2 t /ŝ explicit that is implicitly included in the hard function in [19]. The soft function W R i in (2.16) is defined by a time-ordered expectation value of Wilson lines and corresponds to the 2 → 1 process of the production of a heavy particle in the colour representation R [13]. For resummation of threshold logarithms the momentum-space method [32,39,40] is used, where the soft and hard functions are evolved from a soft scale µ s ∼ m t β 2 and a hard scale µ h ∼ 2m t to the factorization scale µ f ∼ m t using renormalization-group equations (RGEs) summarized in Section 3 below. The potential function is given in terms of the imaginary part of the Coulomb Green function in non-relativistic QCD, which resums ladder diagrams with Coulomb-gluon exchange to all orders in α s . The resulting resummed cross section is of the form [19] σ res Here the Laplace transformation of the soft function, was introduced, where s = 1/(e γ E µe ρ/2 ). After carrying out the differentiations with respect to η in (2.18), this variable is identified with a resummation function which contains single logarithms, α s ln(µ s /µ f ), while the resummation function U i sums the Sudakov double logarithms α s ln 2 µ h µ f and α s ln 2 µs µ f . The precise definitions of these functions for the case of heavy-particle pair production are given in [19] and the expansions required for N 3 LL accuracy can be found in [32]. For the NNLL resummation carried out in [15], the NLO hard [23,54] and soft [13] functions have been used, together with the three-loop cusp-anomalous dimension [44] and the remaining anomalous dimensions in the evolution equations at two-loop accuracy. In the Coulomb sector, using the NLO potential function quoted in [15] resums all corrections of the form (α s /β) k and α s × (α s /β) k . This was supplemented by a leading resummation of logarithms by using a running Coulomb scale, and the inclusion of the leading so-called non-Coulomb correction [15,27,55], which give rise to a tower of terms of the form α 2 s ln β × (α s /β) k .

Towards N LL resummation
In this section we collect the input required for N 3 LL resummation according to the systematics defined in (1.1) and identify missing ingredients, whose effect will be estimated in the phenomenological results.
For soft-gluon resummation alone, the required ingredients to increase the logarithmic accuracy are well-defined, but not fully available for N 3 LL resummation. For N n LL accuracy according to (1.1), the cusp-anomalous dimension has to be known at (n + 1)-loop order, the remaining anomalous dimensions in the evolution equations of the soft and hard functions to n-loop order. The fixed-order soft and hard functions have to be known to N n−1 LO (N n LO) accuracy for the N n LL (N n LL') approximation. For top-quark production, NNLO soft and hard functions can be obtained from results of [32][33][34]56] as discussed in Sections 3.1 and 3.2. The four-loop cusp anomalous dimension was recently computed for the quark case [57] but is still unknown for gluons. However, this affects only the full N 3 LL resummation but not the expansion to N 3 LO considered in this paper. Further, the three-loop soft anomalous dimension is known for the massless [58,59] but not the massive case.
Concerning the corrections in the Coulomb sector, the potential function for the coloursinglet, spin-triplet case is known at the required accuracy from a calculation of e − e + → tt [35][36][37]. This requires the insertion of O(α 2 s , α s β, β 2 ) suppressed potentials, which are not fully known for colour-octet states. In Section 3.3 we obtain the expansion of the potential function to O(α 3 s ). We partially generalize the result to general spin and the colour-octet case and estimate the remaining unknown colour-octet contributions.
However, the presence of Coulomb corrections near the total production thresholdŝ → 4m 2 t complicates the extension to higher accuracy compared to resummations for other kinematical threshold definitions such as the pair-invariant mass or single-particle-inclusive observables [20][21][22]. The enhancement of Coulomb corrections by negative powers of β requires to control the corrections involving positive powers of β in the curly brackets in (1.1) to achieve the desired accuracy using the counting α s ∼ β. In the effectivetheory framework of [13,19] used to derive the factorization (2.16), such corrections arise from power-suppressed interactions and production operators in SCET and PNRQCD and require an extension of the factorization formula (2.16) with generalized hard, soft and potential functions. In particular, the simplification of the colour structure to that of a 2 → 1 scattering process has only been shown at leading power. It is known that threeparticle colour correlations appear in infrared singular parts of the pp → tt scattering amplitudes [60,61], but do not contribute to the NNLO and NNLL cross section [19,27,33]. In the framework of [13,19] these corrections do not enter in the definition of the soft and hard functions and their evolution equations, but rather through the generalized soft and hard functions in the extended factorization formula. The complete treatment of the O(α 2 s , α s β, β 2 ) corrections required for a full soft and Coulomb resummation at N 3 LL is beyond the scope of this paper, but we will include several corrections of this form: • Corrections of the form α 3 s ln 2,1 β to the N 3 LO cross section arise from (ultra)soft corrections 5 due to the power-suppressed chromoelectric x· E s vertex in the PNRQCD Lagrangian. In the language of [13,19] these corrections arise from generalized soft and potential functions. We estimate these corrections using the known result for the colour-singlet case [62].
• Corrections due to subleading production operators are related to the P-wave production channels (2.14), which give rise to α 3 s × ( 1 β , ln 2 β, ln β) corrections relative to the leading S-wave channels since the β 2 suppression of the Born cross section can be compensated by the second Coulomb singularity. In Section 3.4 we apply the NLL resummation formula to P-wave channels, which is sufficient to obtain the threshold-enhanced terms at O(α 3 s ).
• The interference of the second Coulomb correction with subleading soft-collinear interactions suppressed by β 2 results in contributions ∼ α 3 s log β to the total cross section. While there is current interest in subleading soft-collinear effects using QCD factorization [63][64][65] and SCET [66][67][68][69] approaches, these effects have not been systematically studied for top-quark production. In Section 3.5 leading-logarithmic subleading corrections to the initial state are heuristically included using universal results from Drell-Yan and Higgs production.
In this way, subleading soft-collinear and soft-potential corrections are included. From the arguments used to exclude contributions from power suppressed soft-potential or softcollinear interactions at NNLO [19] it is expected that sub-leading corrections at N 3 LO separately affect the initial and final state and a cross-talk between soft-collinear and softpotential corrections only appears at higher orders. 6

Hard functions
The hard function satisfies the evolution equation Up to the three-loop level [44], the light-like cusp anomalous dimension satisfies the socalled Casimir scaling, i.e. Γ r cusp = C r γ cusp , where C r is the quadratic Casimir for the SU (3) representation r of the parton p. At the four-loop level, the recent result for the quark cusp-anomalous dimension [57] shows a violation of this property (see also [70,71]). The anomalous dimension γ V i can be written in terms of single-particle anomalous dimensions: where the anomalous dimensions for light partons, γ p , are known up to three-loop level. The results in the notation used here are collected in [72]. The anomalous dimension for massive particles γ R H,s is known up to two-loop level [13,73]. The anomalous dimension γ R,S J (µ) arises first at two-loop order because of additional IR divergences in the hard function, which are related to UV divergences of the potential function arising from insertions of non-Coulomb potentials. The RGE for the corresponding matching coefficient of the electromagnetic quark current was derived at NLL in the PNRQCD framework [30] and the fixed-order three-loop anomalous dimension was computed in [74], see also [75]. For the colour-octet case, the O(α 2 s ) result for γ R,S J can be obtained from the non-Coulomb corrections to the potential function obtained in [27,55], and is given in (3.34). The three-loop result is not known yet.
The hard functions are required up to the two-loop level for NNLL and N 3 LL resummation. The one-and two-loop solution of the evolution equation (3.1) is given explicitly in Appendix A.1. The initial conditions for the evolution can be obtained from comparing the scaling functions obtained from the expansion of the resummation formula to the threshold limit of the result of a diagrammatic calculation of the total NLO and NNLO cross section in the corresponding colour and spin state. Subtracting constant contributions to the potential and soft functions, this relates the hard function at some initial scale to the constant (i.e. ρ-independent) term in the scaling function f (n,0) pp,(R) , which will be denoted by C

R(n) pp
(the spin-dependence of the scaling functions will be left implicit for notational simplicity). The resulting relations of the NLO and NNLO constants C to the one-and two-loop hard functions are given in A.7 and A.8. In order to extract the hard functions, we choose the initial scale µ h = m t , which simplifies the relation to the C R(n) pp compared to the more natural scale 2m t since the factorization-scale dependence is conventionally expressed in terms of ln( µ f mt ) in (2.7). The one-loop hard coefficients are obtained from the comparison to the threshold expansion of the NLO cross section [23,54] where the notation for the perturbative expansion of the hard function is given in (A.4).
The results for the two-loop constants in the threshold limit have been obtained in [34], In the expression for h [15] the non-logarithmic terms with colour factor C A should have the opposite sign. The numerical implementation is not affected by this typo.
We denote these constants by a bar to denote that the correction factors in the threshold expansion in [34] are defined relative to the leading cross sections (2.13) atŝ = 4m 2 t , whereas we keep the factor 4m 2 t /ŝ = 1 − β 2 unexpanded. Because of the enhancement due to the second Coulomb correction, this leads to a difference in the expression for the two-loop constants in the two conventions, see Eq. (A.9). In addition, since no spin decomposition is performed for the NNLO results [34], contributions from sub-leading Pwave production channels as given in (3.48) below have to be subtracted. 8 A further contribution to the NNLO constant arises from the subleading antisymmetric colour-octet gluon channel (2.15) since the accidental suppression of the Born matrix element does not persist at one loop [34,48,49]. Therefore, the one-loop squared contribution in the 8 a channel leads to an O(α 2 s /β 2 ) correction relative to the leading order cross section (2.15) [34] which contributes at the same order as the constants C

R(2) gg
in the dominant coloursymmetric production channels. Since the top-quark pair is produced in an S-wave state in the 8 a channel and the soft-gluon and Coulomb corrections do not distinguish between the 8 s and 8 a representations, the contribution to the resummed cross section can be correctly taken into account by adding the constant (3.9) to the one of the symmetric octet (3.8).
We therefore definē where the normalization of C 8a(2) gg in eq. (4.8) of [34] was taken into account and N c = 3 was used in the second step. The constantsC R(2) pp are related to the two-loop coefficients of the hard function by (A.10). For the number of light flavours n f = 5 and using N c = 3 and T F = 1/2 we obtain

Soft function
The RGE of the Laplace-transformed soft function (2.19) reads 14) The anomalous dimension γ W is related to that of the hard function according to with the anomalous dimensions γ φ,p entering the evolution of the parton distribution functions in the x → 1 limit (2.6). In the second equality of (3.15), soft anomalous dimensions of massless partons p in the colour representation r, γ r s ≡ γ p + γ φ,p , have been introduced. At least up to three-loop level for the massless partons and the two-loop level for massive particles, the soft anomalous dimensions satisfy Casimir scaling The known results for the soft anomalous-dimension coefficients are summarized in (A. 16) and (A.17). The unknown three-loop soft anomalous dimension for massive particles, γ H,s , will be kept explicitly in the results to estimate the resulting uncertainty.
The one-and two-loop solution of the evolution equation (3.14) is given in Appendix A.2. As initial conditions, we require the functionss R i (0), which are available for arbitrary colour representation [13] at one-loop level, while the two-loop coefficients have been calculated for the colour singlet [32,56] and octet [33] for identical initial state representations (i.e. C r = C r ),

Potential corrections
The potential function for S-wave production in the colour channel R α and spin state S is given in terms of the imaginary part of the zero-distance Green function of the Schrödinger equation, where we use the abbreviation v = E/m t . The LO Green function is obtained as solution of the Schrödinger function with the Coulomb potential and sums up all correction of the form (α s /β) k . For higher-order corrections we make use of the method described in [35] where insertions of higher-order potentials are treated perturbatively. The N 3 LO Green function for the colour-singlet, spin-triplet case computed in [35][36][37] includes all corrections up to terms of order ( to the S-wave production processes (2.13). The expansion of the NNLO Green function reads ImG where the kinematical correction ∼ v 2 arises from the relativistic kinetic energy correction term ∂ 4 /8m 3 t in the NRQCD Lagrangian. Furthermore, the relation of the variable v used in the PNRQCD Green function to the non-relativistic velocity β reads at this accuracy Taking these corrections into account leads to the value familiar from the treatment of top-quark pair production in e − e + collisions. The coefficient d J receives corrections at O(α s ) and becomes scale-dependent. However, these corrections only contribute to the constant term at O(α 3 s ) and are therefore beyond N 3 LL accuracy. Starting from NNLO the potential function is explicitly scale dependent, where the anomalous dimension is the same one that appears in the evolution equation of the hard function (3.1). In a fixed-order expansion to O(α 3 s ) it has the form The three-loop anomalous-dimension coefficients γ R,S(2,i) J are known for the colour-singlet case [74]. 9 The results for the colour-octet case are currently not known. However, in the fixed-order expansion of the cross section to O(α 3 s ) these terms cancel against corresponding contributions from the evolution of the hard function (3.1).

Potentials
For the computation of the Coulomb Green function we use the formalism of PNRQCD in the conventions of [35]. After performing a colour and spin decomposition (see Appendix B), the potential for a quark-antiquark pair can be written in the form The relation to the notation of [74] is γ . Note that the coefficient β 0 in the convention of [74] is 1/4 that in our convention. A modification of γ (3) v due to the d-dimensional treatment of spin was pointed out in [77].
with the colour factor where C R is the quadratic Casimir of the representation R. For the computation of the O(α 3 s ) corrections to the potential function, we require the Coulomb potential V R C and the 1/m potential up to the two-loop level and the 1/m 2 potentials V R,S 1/m 2 and V R p as well as the so-called annihilation correction ν R,S ann at one-loop. The Coulomb potential at the relevant accuracy is known both for the colour-singlet [78] and octet states [79] and is given in (B.6).
The 1/m 2 potential is spin dependent and reads at leading with the explicit values The one-loop result for the colour-singlet case is given in [35] and generalized to the octet case in Appendix B, see also [80]. The annihilation contribution arises from local four-fermion operators in NRQCD [81] and was omitted in the calculations of [15,27]. The leading correction arises for the colouroctet, spin-triplet case, ν 8,1 ann (0) = 1, (3.29) and results in a contribution to the NNLO threshold expansion of the top-quark pair production cross section in the quark-antiquark channel first obtained by an explicit twoloop calculation in [34]. The annihilation corrections to the one-loop matching coefficients of the four-fermion operators [81] also arise in the spin-singlet channel and therefore also enter in the gluon-fusion initial state. The corresponding one-loop coefficients and the resulting corrections to the potential function are given in Appendix B.2. Therefore, all ingredients of the potential are known at the required accuracy, apart from the 1/m potential, which is known at the one-loop level [27,55] for the singlet and octet cases, but only for the singlet case at the two-loop level [82]. In addition, at O(α 3 s ) also ultrasoft corrections due to the power-suppressed x· E s vertex in the PNRQCD Lagrangian become relevant, which are also only known for the colour-singlet case [62]. We will estimate the effect of the 1/m potential and the ultrasoft corrections in the colour-octet case by a naive replacement of colour factors.

Expansion of the potential function
Given the above results for the potentials, we can obtain the expansion of the potential function to O(α 3 s ) using the expressions from the N 3 LO calculation of e − e + → tt [35][36][37]. The annihilation corrections are computed as discussed in Appendix B.2. A brief description of the methods used for the expansion of the Green function in α s is given in the following. The required corrections are expressed in terms of nested harmonic sums and sums over (poly-)gamma functions. They depend on α s through the parameter where Γ t is the top-quark decay width. In the former case, the summation limits depend on λ, but we can always transform the sums in such a way that this dependence is shifted to the summand. In both cases, we then expand the summand in the limit λ → 0 before performing the sum. Afterwards, the coefficients of the expansion are given by so-called multiple zeta values, which we take from the program summer [83]. A simple example for the series expansion of a single harmonic sum is where ζ(n) = S n (∞) is Riemann's zeta function with integer argument. Since polygamma functions are related to harmonic sums, after expanding in α s we obtain similar expressions in that case as well. The infinite series can then be truncated at the required order in α s . Finally, we can take the limit Γ t → 0 + (with E > 0), which yields the imaginary part of the Green function from the discontinuity of the square root in λ. Bound-state effects arising for E < 0 as taken into account in the resummed calculation [15] are not included in this naive expansion. In the O(α 3 s ) corrections, these bound-state effects give rise to a contribution localized at E = 0 [38], which is added separately below.
The resulting expansion of the potential function in the strong coupling constant is written as where the scale dependence enters only through the variable The constant a 1 arises from the NLO Coulomb potential and is given in (B.7). The O(α 2 s ) coefficient of the anomalous dimension γ R,S J entering the RGE (3.22) can be read off the expansion of the potential function and is given by For ν R,S ann (0) = 0 and the colour-singlet case D 1 = −C F this result agrees with the NLL anomalous dimension of the electromagnetic quark current in eq. (14) of [30] for the fixedorder values of the potential coefficients D 1,s = α s (µ), as well as with [74]. The constant term in (3.32) is given by All terms in (3.32) apart from the constant were obtained in [27], with the exception of the annihilation contribution which was added in [34,55]. The computation of the O(α 3 s ) corrections to the potential function requires the use of the N 3 LO Green function, which includes all corrections of the form α 3 s × (α s /β) k . The result can be split into several contributions: The expansion of the LO-Green function includes a term ∼ α 3 s /E whose imaginary part is nonvanishing for unstable particles (see e.g. [84]) but does not yield a correction of the form α 3 s /v 3 if the decay width is neglected. Treating the imaginary part of the Green function carefully in the distributional sense, it was shown that there is instead a delta-function contribution at O(α 3 s ) [38]: As discussed in [38] this contribution is taken into account in the resummed calculation [15,18], which includes the resummed Coulomb corrections above threshold and the boundstate corrections below. The naive expansion of the NNLO Green function in α s is fully known for all required colour and spin states and gives contributions enhanced by one inverse power of β. The scale-dependence again enters only through the variable L E : with the constant where the coefficient a R 2 of the NNLO Coulomb potential is given in (B.8). All contributions of (3.38) apart from the constant c (0,1) J,3 are included in the implementation of NNLL resummation in [15] and the corresponding approximate N 3 LO prediction. For the coloursinglet case and the spin states S = 1 and S = 0, this result reproduces the imaginary part of the threshold expansion of the vector and pseudo-scalar current correlators given in the Appendix of [85] if the corresponding matching coefficients [86][87][88] are taken into account. 10 The expansion of the N 3 LO correction to the Green function to O(α 3 s ) gives rise to purely logarithmic corrections α 3 s ln 2,1 β and scale-dependent terms α 3 s ln 2,1 ( mt µ ) governed by the RGE (3.22). The result can be written in the form The result of [35][36][37] and the generalization of the 1/m 2 potential to general spin in (B.14) allows to compute the coefficients c (i,0) J,3 and the anomalous-dimension coefficients γ R,S(2,i) J exactly for the colour-singlet state. For the colour-octet case, we use the result for the 1/m 2 potential (B.14) and estimate the unknown 1/m potential and the ultrasoft corrections by a naive replacement of colour factors. The ultrasoft corrections to the Green function in [62] are presented in a form where logarithmic terms ln α s , ln µ mt are given explicitly while a function δ us (Ê) of the variableÊ = E mtα 2 s is only available numerically. Using the fact that ln α s -terms cannot appear in a fixed-order calculation and thus must be canceled by logÊ-terms in the remainder function δ us (Ê), it can be checked that all logarithmic O(α 3 s )-contributions can be reconstructed with the replacement ln α s → ln v. We have checked that the expansion of the Green function computed in [35][36][37] reproduces the three-loop anomalous-dimension coefficients γ R,S(2,i) J for the colour-singlet case [74]. 11 The result for the remaining coefficients is given by

41)
10 In the pseudo-scalar case the so-called singlet contributions to the two-loop matching coefficient [88] have to be set to zero since they are not included in [85]. 11 In this comparison the scale-dependence of the coefficient d J in (3.19) must be taken into account, see e.g. Eq. (3.45) in [35] for the case of the vector current.
The terms given explicitly in these expressions hold both in the colour-singlet and octet cases. The contributions from the 1/m potential and ultrasoft corrections in the coloursinglet case are given by where the terms containing only the single colour factor C F arise from the ultrasoft corrections. 12

P-wave contributions
The gluon-fusion initial state to top-quark pair production receives contribution from the 3 P 0 and 3 P 2 production channels (2.14), which contribute at the same order in the threshold expansion as the kinematic correction to the leading S-wave production channels taken into account in (3.19). The relative P-wave contribution to the total cross section is parameterized by the quantity where σ (0),R pp is the threshold limit of the leading S-wave production cross section (2.13). For top-quark pair production, P-wave production contributes only to the gluon initial state in the symmetric colour representations with The leading-order P-wave potential function for stable particles above threshold is given by [89,90] with the perturbative expansion The α 3 s correction agrees with explicit results of the threshold expansion of the current correlation function [85,91].
The combination of the Born suppression ∼ β 2 of the P-wave cross section and the two-loop Coulomb enhancement leads to a constant α 2 s contribution relative to the leading S-wave production channel and threshold-enhanced N 3 LO contributions of order α 3 s /β and α 3 s ln β. In order to obtain these corrections, it is sufficient to consider soft-gluon resummation at NLL accuracy. In [47] the factorization (2.16) has been shown at least at NLL for partonic channels with leading P-wave production. Here we heuristically apply this result also for the case of subleading P-wave contributions to an S-wave dominated process. 13 We therefore insert the P-wave potential function in the factorization formula (2.16) and use the resummed soft and hard functions at NLL accuracy, setting the soft-and hard scales appearing in the momentum-space resummation formalism to µ s = m t ln 2 β and µ h = m t . Leading logarithms of potential origin are taken into account by using the scale µ C = m t ln β in the potential function.
The constant contribution at NNLO obtained in this way is given by ∆σ R (2) pp ,P = σ with the scaling functions Note that the form of these corrections relative to the "constant" NNLO term (3.48) differs from the usual NLO threshold corrections for an S-wave production process (see e.g. [19]) by the coefficient of the Coulomb correction and the presence of the β 0 -term that arises from the LL-running of the potential function.

Next-to-eikonal correction
Contributions of the order α 3 s ln β can arise from the interplay of power-suppressed NLO corrections α s ln β × β 2 and the NNLO Coulomb correction α 2 s /β 2 . A full analysis of these corrections in the EFT framework is beyond the scope of this paper. However, for the N 3 LO threshold approximation only the LL next-to-eikonal corrections are relevant that arise only from initial-state radiation and therefore can be obtained from results for the Drell-Yan process and Higgs production (see e.g. [63,65,92]).
These corrections can be incorporated in our framework by including an additional term in the one-loop soft function obtained by the replacement in the position-space soft function W R i (2m t (1 − z)) given in eq. (C.5) in [13]. This leads to the LL next-to-eikonal contribution where the factor (2m t ) −1 arises from the normalization convention of the LO soft function, W i,n.e. (2m t (1 − w)) , (3.54) where we have performed a transformation of variables ω = 2m t (1 − w) and approximated The only threshold-enhanced contribution at N 3 LO arises from the second Coulomb correction in the NNLO potential function (3.32). We therefore obtain the logarithmic correction to the N 3 LO scaling function as (3.55)

Approximate N 3 LO results
We are now in a position to obtain approximate N 3 LO predictions based on the partial N 3 LL resummation constructed in the previous section. In Section 4.1 we present explicit results for the scaling functions and discuss our estimate of the currently unknown colouroctet coefficients. In Section 4.2 we present results for the total hadronic top-quark pair production cross section and discuss the scale dependence and the estimate of the remaining theoretical uncertainty.

Partonic cross sections
Our default implementation of the partial N 3 LL resummed cross section is obtained from the resummation formula (2.18) by inserting the N 3 LO potential function, the resummation functions U i and η at N 3 LL accuracy, and the two-loop hard and soft functions obtained in Sections 3.1 and 3.2. In addition, the P-wave contribution (3.49) and the contribution from the interference of next-to-eikonal logarithms and the Coulomb corrections (3.55) are added. The resummation formula is expanded to O(α 3 s ), setting the soft and hard scales introduced in the resummation formalism to µ s = k s m t v 2 and µ h = k h m t , where the unknown constants k X are of order one. For N 3 LL resummation, all thresholdenhanced N 3 LO terms are independent of these unphysical parameters.
We will present the results for the approximate N 3 LO corrections for the different partonic production channels in the form where the threshold limit of the partonic production cross sections for a given colour channel (2.13) was factored out. The renormalization-scale dependence can be obtained using (2.4). We do not perform a decomposition into states of definite spin, i.e. the scaling functions for the gluon initial states include the P-wave contribution (3.49). The result of the ω-integral over the potential function in the resummation formula (2.18) is expressed in terms of the energy variable E = √ŝ − 2m t ≡ m t v 2 . To write the approximate N 3 LO corrections as functions of the customary variable β = √ 1 − ρ, the terms ∼ β 2 in the relation of the two variables (3.20) must be included. Taking the overall factor v in the potential function (3.31) into account, the relevant corrections to the thresholdenhanced terms at O(α 3 s ) can be obtained by the replacements where the O(1) terms contribute to the constants C pp in the scaling functions. Note that in contrast to the NNLL implementation [15], at this order one should not replace the threshold limit σ R(0) pp in (4.1) by the full Born cross section, since this would treat these kinematic corrections in an inconsistent way. In order to estimate kinematic ambiguities of the threshold approximation, we will also consider the results obtained with the unexpanded scaling functions in terms of the variable v.
The results for threshold expansion of the scale-independent scaling functions for the different production channels are for the explicit values N c = 3, n f = 5,  The coefficients of the terms ln 6,...3 β, 1 β ln 4,...1 β and 1 β 2 ln 2,...0 β reproduce the results of the expansion of the NNLL cross section in [15] (up to shifts 41839.4 → 41558.7 in the coefficient of the ln 3 β term and 6528.61 → 6506.96 in the ln β/β term in the quarkantiquark channel due to the annihilation contribution, which was omitted in [15]). The threshold approximation of the scaling functions f (3,i) for i = 0 related to the factorizationscale dependence are given in Appendix C. These results were obtained both from the expansion of the resummation formula and by exploiting the known factorization scale dependence of the PDFs as discussed in Section 2.1, finding agreement with the results obtained using the two methods.
The results for the colour-octet channels in (4.3) depend on the unknown three-loop coefficient of the massive soft anomalous dimension γ (2) H,s and the coefficients δc J,3 parameterizing the unknown corrections to the potential function from the 1/m potential and ultrasoft corrections. In our default predictions we set the three-loop massive soft anomalous dimension to zero and estimate the uncertainty based on the known one-and two-loop results: As default value for the potential coefficients in the colour-octet case, we take the estimate obtained by the naive replacement C F → (C F − C A /2) in the known colour-singlet result (3.43). This amounts to the numerical values We estimate the uncertainty due to this estimate as ±2δc (n,0) J,3 | octet approx. . It is seen from (4.3) that the numerical effect of the unknown potential corrections is expected to be larger than that of the three-loop soft-anomalous dimension.
The constants C pp (R) in the threshold expansion of the scaling functions (4.3) can only be obtained from a three-loop computation and will be set to zero in our default approximation. The expansion of the N 3 LL resummation formula yields expressions for these constants that depend on the unphysical scale parameters k s , k h , and k C introduced in the resummation formalism and are given in Appendix C. We will use these results to estimate the uncertainty of our result due to effects beyond the threshold approximation. Choosing as a default k X = 1 and varying the scales independently in the interval 1 2 ≤ k X ≤ 2 with the constraint 1 2 ≤ ks k h ≤ 2 we obtain the estimates We further define the colour-averaged constant in the gluon channel, In Fig. 1 we compare the result (4.3) for the approximate N 3 LO corrections with an earlier approximation (denoted as N 3 LO B ) based on NNLL resummation [15], where the terms ∼ ln 2,1 β and ∼ 1 β beyond NNLL accuracy were dropped. In the figures we plot the integrand of the convolution of the approximate N 3 LO corrections (4.1) with the parton luminosity in the formula for the hadronic cross section (2.1) as a function of β, including the Jacobian ∂τ /∂β: The gray band indicates the uncertainty of our default approximation (4.3) as estimated from the dependence of the constants (4.6) and (4.9) on the unphysical scales k s , k h , and by using the kinematic variable v instead of β. In the figure the MMHT2014NNLO PDFs [6,93] are used to compute the parton luminosity.
From the comparison of the curves for our best prediction N 3 LO app to the earliear approximation N 3 LO B it can be seen that the effect of the corrections beyond NNLL is small in the limit β → 0 but becomes significant for β 0.3. As already observed in [15], the corrections to the integrand dσ/dβ in the gluon-fusion channel become very large for β → 0 and could lead to corrections of the order of several hundred picobarn to the total cross section at the LHC with 13 TeV. In [15] it was speculated that this might indicate a poor convergence of the fixed-order expansion, despite an overall small numerical effect of resummation on the total cross section. In the total cross section, these large positive corrections are cancelled to a large extent by negative corrections for β > 0.2, where, however, the threshold approximation is less reliable. The integrand approaches zero rapidly for β → 1, which indicates that we are at least not including spuriously large corrections from the region where the threshold approximation breaks down. The gray band in the figures indicates the growing uncertainty of the threshold approximation for β 0.4. The relative uncertainty is larger in the quark-antiquark channel, however, the impact of both channels on the uncertainty of total cross section is similar as seen in (4.11) below.

Phenomenological results
In order to obtain predictions for the total hadronic top-quark pair production cross section, we use as default the MMHT2014NNLO PDFs [6,93]  where "kin." denotes the uncertainty due to the expansion in terms of v instead of β. It is clearly seen that the uncertainties of the threshold expansion as estimated by the variation of the constants and the kinematic ambiguity dominate over the uncertainties due to the unknown soft anomalous dimension and potential coefficients. In addition to the corrections included in the scaling functions (4.3), there is the distributional contribution arising from the third Coulomb correction (3.37). Inserting this correction to the potential function into the factorization formula (2.16) gives the correction 15 The resulting correction at the LHC with √ s = 13 TeV is given by ∆σ C 3 tt (13TeV) = 0.60 pb, (4.14) 15 The present notation is related to the one in [38] by H R pp ,i = in agreement with [38]. This contribution to the N 3 LO cross section is therefore of a similar size as the uncertainties due to the unknown coefficients in N 3 LL resummation. Our default prediction is obtained by adding the correction ∆σ N 3 LOapp tt from the threshold expansion of the scaling functions and the correction (4.13) to the NNLO cross section [9][10][11][12], σ In order to estimate the scale uncertainty, we vary both renormalization and factorization scales independently in the interval 1 2 ≤ µ f , µ r ≤ 2 with the constraint 1 2 ≤ In addition to the uncertainty of the perturbative calculation estimated by the scale variation, the prediction of the hadronic cross section also relies on the value of the strong coupling constant and the PDFs, with their respective uncertainties. For the MMHT2014NNLO PDF, the combined PDF+α s uncertainty of the NNLO top-quark pairproduction cross section at the LHC is about ±2.5% [93], and therefore below the perturbative uncertainty of the approximate N 3 LO result (4.16). Following [94] we note that the evolution of the NNLO PDFs is performed including the three-loop splitting function, so the use of NNLO PDFs in an N 3 LO calculation is formally consistent, although the use of N 3 LO predictions in a PDF fit might lead to a non-negligible effect on the PDFs, in particular in the case of top-pair production [95]. In the following we discuss only the estimate of the perturbative uncertainties of our calculation.
It is interesting to compare the approximate N 3 LO result to the prediction of NNLL resummation matched to NNLO. Using the programs topixs 2.0.1 [18] and top++ 2.0 [17], which implement the combined soft-gluon and Coulomb resummation in momentum space [15] and soft-gluon resummation in Mellin space [16], respectively, we obtain The resummed Coulomb corrections beyond NNLO and bound-state corrections included in topixs only have a small effect, +2.8 pb for the bound-state corrections and less than one picobarn for higher-order Coulomb corrections. Furthermore, top++ includes the twoloop constants C (2) pp , which are part of the NNLL' approximation. Setting these constants to zero, one obtains the top++ result of 812.20 pb, which is closer to the result of topixs. The remaining difference indicates the size of sub-leading corrections, which are treated differently in the different resummation methods. The top++ result is also closer to the approximate N 3 LO prediction (4.16), which indicates the numerical relevance of the interplay of the two-loop constants and the soft threshold corrections. The total scale uncertainty of the two NNLL results is similar to each other and to that of the approximate N 3 LO correction. The topixs result also includes an estimate of the ambiguities of the resummation procedure, which includes the variation of the soft scale used in the momentum-space approach and the comparison of the expansion in v and β. Since our approximate N 3 LO cross section is based on a partial N 3 LL resummation, the dependence on the soft-scale is reduced compared to the NNLL resummation and enters only in the estimate of the constant term. No attempt to estimate the ambiguities of the resummation procedure is made in top++. The fact that the approximate N 3 LO corrections are close to the resummed results indicates that the perturbative expansion is well-behaved, despite the large corrections to the partonic cross sections for small β observed in Figure 1.
In Figure 2 we compare the factorization-scale dependence of the approximate N 3 LO cross section (red dash-dotted) to the NNLO result (blue, long-dashed) and the two resummed predictions from top++ (green, dashed) and topixs (magenta, dotted). The renormalization scale is set equal to the factorization scale. In the topixs calculation, all other scales are fixed to their default value. It is seen that the inclusion of the approximate  N 3 LO corrections further flatten the scale dependence compared to the NNLO prediction and brings the curve closer to the resummed top++ result. In contrast, the topixs curve follows more closely the NNLO result. In Figure 3 we consider instead the factorization scale dependence for fixed renormalization scale and vice-versa. No results from topixs are shown since the renormalization scale is not separated from the factorization scale. Again it is seen that the inclusion of the approximate N 3 LO corrections stabilize the scale dependence and brings the result closer to the resummed prediction.
Approximate N 3 LO results were also published in [41,42]. In these papers, no analytical results for the partonic cross sections are given so we can only compare numerical results for the total cross section. In addition, different PDFs are used for the results. The result of [41] is based on NNLL resummation for one-particle inclusive kinematics. At the LHC with 13 TeV, a correction of 2.7% relative to the NNLO result is found with a scale uncertainty of +2.9% −2.0% . The inherent uncertainty of the constructed approximation is not estimated. Ref. [42] constructs a partial N 3 LL soft resummation in Mellin space using similar input as our approximation, but without the contributions of the N 3 LO potential function. Furthermore, they also include 1/N -suppressed contributions in Mellin space and information on the large-β behaviour. As a result, they find larger corrections of +4.2% at 13 TeV. The scale uncertainty of ±2.7% is comparable to that of our result (4.16). The scale-dependence of our result in Figure 3 is qualitatively similar to that shown in Figure  9 of [42]. An approximation corresponding to the standard Mellin-space resummation [16,17] dubbed "N -soft" yields smaller corrections of +2.3% and is therefore closer to our prediction. The ambiguities of the approximation are estimated as ±1.9%, which is larger than our estimate of +1.2% −0.8% in (4.16). Therefore, while our default implementation results in a smaller correction than the one found in [41,42], the results are consistent if the estimate of the uncertainties due to the threshold approximation are taken into account. In particular, using the expansion in terms of v instead of β, we obtain a larger N 3 LO effect of 2.3% which is in good agreement with the "N -soft" prediction of [42] and with [41].
It is interesting to compare our uncertainty estimate in more detail to the one of [42]. The "approx" uncertainty estimate in [42] includes contributions from two sources: the treatment of the subleading contributions in Mellin space and the parametric uncertainties due to the unknown three-loop constant term and the massive soft anomalous dimension for the colour-octet case. The estimate of the three-loop constant in Eqs. (5.3) and (5.5) of [42] amounts to an uncertainty of the cross section of ±2.6 pb at the LHC, which is smaller than our estimate in (4.12). In contrast, the estimate of the colour-octet softanomalous dimension in [42] amounts to an uncertainty of ±4.6 pb which is much more conservative than our estimate. 16 Nevertheless, combining the two sources of parametric uncertainty in quadrature results in a similar parametric uncertainty as in our result. The dominant contribution to the uncertainty estimate of [42] arises from the ambiguities due to the treatment of the subleading contributions in Mellin space, which contribute an uncertainty of ±1.8%, compared to our smaller estimate of the kinematic ambiguities using the expansion in v instead of β. In addition, Ref. [42] estimates the uncertainty from omitting the N 3 LO approximation for the quark-antiquark and quark-gluon initial states as 1%. Instead, our result (4.11) includes the quark-antiquark channel with a contribution of only 0.2%.

Conclusions
We have constructed a combined resummation of threshold logarithms and Coulomb corrections to the total top-quark pair-production cross section at hadron colliders at partial N 3 LL accuracy and used the result to compute the threshold limit β = 1 − 4m 2 t s → 0 of the partonic N 3 LO cross sections qq → tt and gg → tt. Our result makes use of the state-of-the-art N 3 LO results for the non-relativistic Green function from a calculation of e − e + → tt. We have also included contributions from sub-leading P-wave production channels and power-suppressed so-called next-to-eikonal corrections that give rise to threshold-enhanced N 3 LO corrections due to Coulomb enhancement. We have generalized the so-called 1/m 2 one-loop potential in PNRQCD to the spin-singlet and colour-octet state. However, the three-loop massive soft-anomalous dimension, the two-loop 1/m potential in PNRQCD and the ultrasoft corrections to the non-relativistic Green function are currently not known for the production of top-quark pairs in a colour-octet state. We have estimated the uncertainty due to these missing ingredients on our result. The corrections from the N 3 LO Green function and the missing contributions for the colour-octet state affect the coefficients of the terms α 3 s ln 2,1 β, and can have a larger numerical effect than the unknown three-loop massive soft-anomalous dimension.
Using the results for the threshold expansion of the partonic N 3 LO cross sections, we obtain moderate corrections +1.5% relative to the NNLO predictions for the LHC with √ s = 13 TeV and a reduction of the factorization-and renormalization-scale uncertainty 16 The conventions of [42] are related to our conventions according to (4π) 3 (D H,s +. . . , where additional terms arise due to the two-loop soft function [13]. The uncertainty estimate of (D to the ±3%-level. We estimate the uncertainty of our prediction due to the threshold approximation to be about ±1%. The total perturbative uncertainty therefore approaches a similar level as the PDF+α s uncertainty of ±2.5%. Our default prediction is smaller than other approximate N 3 LO predictions [41,42] but consistent with these results if the estimate of the uncertainties due to the threshold approximation are taken into account. The difference between these predictions indicates the possible size of subleading terms in the threshold expansion.
We have explicitly provided all necessary ingredients to implement our results in a numerical program. Since the numerical result of the approximate N 3 LO corrections and their scale uncertainty are close to the resummed NNLL predictions, this provides a simple, computationally less expensive way to obtain the leading corrections beyond NNLO. We plan to include these results in a future version of topixs, as well the implementation of the fully resummed threshold corrections at partial N 3 LL accuracy.

Acknowledgments
We would like to thank Martin Beneke for useful discussions and for making a copy of Ref. [80] available. The work of CS is supported by the Heisenberg Programme of the DFG. We acknowledge support by the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe" during parts of this work.

Note added
In this version v3 of the arXiv submission of the paper, typos in the coefficients of the ln β-terms in the scaling functions f (3,0) gg(1) in (4.3), f (3,1) gg (8) in (C.12), and the 1/β-term in f (3,2) gg(1) in (C.7) have been corrected. These typos do not affect the implementation used to compute the numerical results of the cross sections. However, in addition the scale dependence in the results (4.16) and (4.18)-(4.20) has been corrected.
In the latest version 3.0 of the program topixs, which can be obtained from http://users.ph.tum.de/t31software/topixs/ the approximate N 3 LO corrections computed in this paper are available as an option by setting APPROX NNNLO=1 in the input file topixs.cfg.

A Renormalization group functions
In this appendix we give explicit results for the NNLO expansion of the hard and soft functions. The various anomalous dimensions appearing in the expressions are expanded in the strong coupling constant according to For the expansion of the β-function and the splitting functions we use the conventions which corresponds to β 0 = 11 3 C A − 4 3 n l T f .

A.1 Two-loop hard functions
The perturbative expansion of the hard function in the QCD coupling can be written as where we have extracted the leading-order hard function H R(0) pp (µ h ) that includes a factor α 2 s (µ h ). We also left the spin dependence implicit. The relation of the one-and two-loop coefficients at a scale µ to the initial conditions at µ h = m t follow from the evolution equation (3.1): Here the factor α 2 s (µ h ) in the leading-order hard function was taken into account. We have also defined the colour factor C rr = C r + C r .
Expanding the resummation formula for the NNLL and NNLL' approximations to O(α s ) and O(α 2 s ), respectively, one obtains the relation between the constant terms C pp (m t ) C rr 4(8 + 9 ln 2 2) − 48 ln 2 − − β 0 C rr −24(3 ln 2 2 − 6 ln 2 + 8) ln 2 + 11π 2 ln 2 − 112ζ 3 3 − 22π 2 3 + 128 Although the natural hard scale is of the order 2m t , these relations simplify for the choice µ h = m t made here. In particular, they become independent of the anomalous dimensions γ V i , which do not obey Casimir scaling. In the NNLO coefficient, the constant terms in the two-loop soft function (3.18) and the NNLO potential function (3.35) enter. Additionally, the coefficient d J from kinematic corrections to the matching of the potential function (3.19) and the relation (3.20) of the variables v and β have been taken into account. The constants C R(2) pp in (A.8) are related to the results in the conventions of [34], here denoted by a bar, according to Here the first correction term arises from multiplying the prefactor 4m 2 t /ŝ = 1 − β 2 with the second Coulomb correction. Further, since no spin decomposition is performed in the results for the constantsC R(2) pp in [34], the P-wave contribution (3.48) with the coefficient (3.45) was subtracted to obtain the hard function for the S-wave process.
For the number of light flavours n f = 5 and using N C = 3 and T F = 1/2, but keeping the colour representations arbitrary, the NNLO result becomes

A.2 Two-loop soft function
We consider the perturbative expansion of the soft function in Laplace space, The coefficients of the variable ρ are fixed by the anomalous dimensions in the evolution equation (3.14), so the only required input are constant coefficientss H,s +s Here it was used that the soft anomalous dimensions of the light partons γ s vanish at the one-loop level. The initial conditionss (2)R i (0) given in (3.18) were obtained from twoloop calculations of the soft function in the literature by noting that the soft function in Laplace space is obtained from the position-space result in [32] by the simple replacement L = 2 ln iz 0 µe γ E 2 → −ρ and from the Mellin-space result in [33] by the replacement The explicit results for the cusp anomalous dimension are given by T F n f , The one-and two-loop coefficients of the soft anomalous dimensions of the light partons are given in [32] while the three-loop coefficient can be obtained from the results for γ (2)p and γ (2)φ,p in [32,45,72] The heavy particle anomalous-dimension coefficients are only known up to two loops T F n f . In this appendix we provide the projection of the PNRQCD potential onto spin-singlet and triplet states and obtain the one-loop 1/m 2 potential for the colour-octet case. The conventions of the PNRQCD Lagrangian follow [35]. The general colour dependence of the quark-antiquark potential can be decomposed in the two equivalent forms 17 where the potential coefficients in the two conventions are related by with the colour factor (3.25). The potential for arbitrary spin in a given colour channel R is of the form [35] Here the spin-dependence is written in a tensor-product notation, a ⊗ b, where a (b) refers to the spin-matrix on the quark (antiquark). In addition to (B.3), also the so-called annihilation contribution from local NRQCD four-fermion operators must be included for hadronic top-quark production. This is discussed in App. B.2. The potential (B.3) can be further decomposed into spin-singlet and triplet contributions according to where the spin projection in d − 1 space dimensions can be performed as discussed in Section 4.5 of [35]. In this way, one obtains the colour-and spin-projected potential Additionally, the term ∂ 4 /8m 3 t in the PNRQCD Lagrangian is also treated as a potential (cf. Eq.(4.60) of [35]).
We require the Coulomb-potential up to two-loop accuracy, which reads in four dimensions with the one-loop coefficient The two-loop coefficient depends on the colour representation [78,79], The 1/m potential for the colour-singlet case is quoted in Section 4.4.2 of [35]. The colouroctet one-loop result can be found in [55], while the corresponding two-loop result is currently not known. The resulting uncertainty on our predictions is estimated as discussed in Section 4.1.
The spin-dependence in (B.5) enters only through the 1/m 2 potential At tree level, the potential coefficients in the 1/m 2 potential are The one-loop result for the colour-singlet 1/m 2 potential is given in Section 4.4.3 of [35]. We argue now that the colour-octet result can be obtained from the singlet case after an appropriate adjustment of the colour factors. Following [35], the one-loop 1/m 2 potential can be calculated from an NRQCD computation consisting of a "hard" contribution of tree diagrams with one-loop matching coefficients, and a "soft" contribution of one-loop NRQCD diagrams. The soft contributions receive contributions from the diagrams in Fig. 9 of [35]. It can be seen from the NRQCD Feynman rules that the colour factors are either the same as the box or crossed-box topologies in the full QCD diagrams, or can be written as linear combinations thereof. The projection of these colour factors on the singlet and octet states can be written in the form C box → D 2 R and C c-box → D R D R + C A 2 so that in both cases the result for the general representation is obtained from the singlet result by the replacement (−C F ) → D R . There are further contributions from the abelian vertex correction, self-energy insertions and charge-renormalization counterterms. These contributions contain explicit factors of C F but are scaleless, so they do not contribute to the finite part of the potential. The hard corrections receive contributions from one-loop matching coefficients of four-fermion operators, the NRQCD quark-gluon vertex and the gluon two-point functions. Since the former arise from box-diagrams in full QCD, the same replacement of colour factors as in the soft contribution can be applied, while the colour factors arising from the coefficients d i of the vertex and two-point functions are unmodified. 18 As result, we obtain the 1/m 2 potential for the spin-singlet/triplet 18 Formally this result is obtained by adjusting the colour projection of the contribution of the Wilson coefficients of the four-fermion operators in the expressions for V and colour-singlet/octet states with v S,(1) where the O( ) contributions do not contribute to the logarithmic corrections considered here. This result agrees with an explicit calculation of the singlet and octet cases [80]. The potential coefficient V R,(1) p is the same for both the singlet and the octet case.

B.2 Annihilation contribution
The contribution of local four-fermion operators to the NRQCD Lagrangian is given by δL 4f = −δV 4f ab;cd ψ † c χ b χ † a ψ d , (B. 16) where the labels collectively denote the spin and colour quantum numbers. This convention agrees with [55]. The four-fermion Lagrangian can be decomposed into contributions with definite spin and colour quantum numbers, analogously to (B.1) and (B.4) The contributions to the four-fermion operators arising from "scattering" diagrams (i.e. diagrams present if ψ and χ correspond to different quark flavours) are already included in the PNRQCD 1/m 2 potential discussed in appendix B.1. For the "annihilation" contributions (i.e. diagrams only present for identical quark flavours) we use the notation Here we have neglected imaginary parts, which contribute to toponium decay into light hadrons and should not be taken into account for the total cross section for pp → tt → bbW + W − . The corrections to the NNLO Green function due to an insertion of an annihilation correction is [55] δ ann G   R (0, 0; E) is the Green function with the insertion of one NLO Coulomb potential. The second term has the same form as the double insertion of the NLO Coulomb potential and a delta potential treated in [96] and included in the calculation of [35][36][37], so the annihilation correction can be obtained from these results. The logarithmic O(α 3 s ) corrections to the potential function are given by

C Scaling functions
In this appendix we collect the numerical expressions for the threshold expansion of the scaling functions f (3,i) pp (R) with i ≥ 1 parameterizing the factorization-scale dependence, and the constants C qq (8) in the scaling functions f (3,0) pp (R) . The results of the scaling functions with i = 2, 3 are already fully determined at NNLL accuracy. We include the expressions given already in Appendix C of [15] for completeness. The constant terms in the scaling functions depend on the parameters k s , k h , and k C in the unphysical soft, hard, and Coulomb scales µ s = k s m t v 2 , µ h = k h m t , mu C = k C m t v.

C.1 Quark-antiquark channel
The constant in the threshold expansion in the quark-antiquark channel is given by The constant terms in these functions are given for completeness but are not included in our default approximation. These results agree with those obtained at NNLL, with the exception of the coefficient of the log β-term and the constant in f (3,1) qq (8) , which are a new result.

C.2 Gluon fusion colour-singlet channel
The constant in the scaling function for the colour-singlet gluon-gluon initial state is given by The Coulomb-scale dependence arises here because the P-wave contributions are only included with the LO potential function. The NLL treatment of the P-wave corrections leads also to some spurious contributions due to incomplete cancellations between the P-wave contribution and contributions arising from the term involving h (0),R pp ,P in the two-loop constant (A.9), which are entering in the S-wave contribution that is treated at N 3 LL. Since we only use the constant as an uncertainty estimate, we do not attempt to fix this issue.
The scaling functions of the scale-dependent contributions read