Non-resonant and electroweak NNLO correction to the $e^+ e^-$ top anti-top threshold

We determine the NNLO electroweak correction to the $e^+ e^-\to b\bar{b}W^+W^- X$ production cross section near the top-pair production threshold. The calculation includes non-resonant production of the final state as well as electroweak effects in resonant top anti-top pair production with non-relativistic resummation, and elevates the theoretical prediction to NNNLO QCD plus NNLO electroweak accuracy. We then study the impact of the new contributions on the top-pair threshold scan at a future lepton collider.


Introduction
The precision study of the top pair production threshold is among the main motivations for the construction of a high-energy e + e − collider [1]. About 100 fb −1 of integrated luminosity spread over ten center-of-mass energies distributed around √ s ≈ 345 GeV can provide a measurement of the top-quark MS mass with an experimental uncertainty of about 50 MeV [2][3][4]. This must be compared to the ultimate precision possible at the LHC, which is constrained to O(1 GeV) due to the limited understanding of the relation between the MS mass and the mass parameter in the calculation and simulation of the final state from which the top mass is directly reconstructed. There has been some progress in the quantification of this relation when the mass is reconstructed from two-jettiness in e + e − collisions in the boosted top regime [5], but the extension of this approach to hadron collider processes requires the consideration of additional effects [6]. In addition, the top width, the strong coupling constant and the top Yukawa coupling can be extracted from the threshold scan to varying degree of accuracy.
The threshold region is defined as the kinematic regime where the top quarks have a small three-velocity v = ( √ s/m t − 2) 1/2 of the order of the strong coupling constant α s . Thus, the top quarks are non-relativistic and are subject to the colour Coulomb interaction, that would facilitate the formation of toponium bound states if the top quarks were stable. The sizeable top decay width caused by the electroweak interaction also prevents hadronization. Therefore, the top threshold dynamics is governed by the colour Coulomb interaction, which must be treated non-perturbatively, while the strong coupling α s ≪ 1 is still small. This interplay between the strong Coulomb attraction and the large top decay width has first been realized in [7,8]. A significant effort has since been invested into providing high-precision predictions for top pair production near threshold. The major focus has naturally been the strong interaction effects, which have now been computed to next-to-next-to-next-to-leading order (NNNLO) accuracy [9] in an expansion where α s ≪ 1 and v ≪ 1, but α s /v = O(1). The effective field theory formalism and ingredients that underlie this calculation are summarized in [10], to which we refer for more details on the QCD aspects of the calculation. 1 The NNNLO QCD result has finally settled the issue of the poor convergence of the perturbative expansion up to NNLO [23]. The NNNLO corrections are well behaved and the remaining scale uncertainty of the QCD result is at the level of ±3%. Similarly, it has been observed that the RG-improved prediction at the (almost) next-to-next-toleading logarithmic order [24] stabilizes the scale uncertainty at the level of ±5%.
In the present work we are concerned with electroweak effects and non-resonant production of the observable final state bbW + W − + X of the decayed top anti-top pair. An analysis of various electroweak effects [25] has demonstrated that they are as large as 10%. Thus, the full NNLO non-resonant and electroweak contributions must be included to salvage the precision of the prediction. Even more importantly, as will be discussed below, they are required to obtain a well-defined result, since the pure QCD cross section by itself contains divergences proportional to the top-quark decay width [15], which are cancelled only once the non-resonant production is included [26,27].
The main result of this work is the NNLO calculation of all electroweak and nonresonant effects. We also provide an implementation of initial-state radiation in a scheme consistent with Coulomb resummation and the inclusion of O(α) electromagnetic corrections, following a similar treatment as for the W + W − threshold [28,29]. To define the precise meaning of "NNLO" for electroweak effects, we note that they introduce the electromagnetic (α em ), SU(2) electroweak (α EW ) and top-quark Yukawa (λ t ) coupling. For the purpose of power counting we do not distinguish between α em and α EW and count that is, an electroweak coupling counts as two powers of the strong coupling, which is consistent with counting Γ t ∼ m t α EW ∼ m t v 2 , which is always adopted in the pure QCD calculation. The pure QCD calculation up to NNNLO then accounts for all terms in the total cross section σ of the form where the global factor α 2 EW v accounts for the phase-space suppression of the cross section near the threshold and the electroweak production in e + e − collisions. The electromagnetic, electroweak, Yukawa and non-resonant terms are of the parametric form where the first line refers to resonant and the second to non-resonant production. We note the absence of phase-space suppression and Coulomb resummation for the non-resonant part. The non-resonant contribution is known at NLO [26], but only partial results are available at NNLO [27,30,31]. On the resonant side, the (α em /v) k terms arise from the QED Coulomb potential. 2 These as well as all Yukawa coupling effects have already been included up to NNNLO in [25]. This result together with the NLO non-resonant and the NNNLO QCD calculation has been made available in the QQbar threshold code [32]. The NNLO non-resonant and the remaining NNLO electroweak contributions are computed in this work, thus elevating the precision at the top-pair threshold to complete NNNLO QCD+Yukawa and NNLO EW+non-resonant. The ellipses in (1.3) denote third-order electroweak and non-resonant terms that remain unknown. The outline of the paper is as follows. In Section 2 we describe how the calculation is split into resonant and non-resonant contributions, such that no double-counting occurs and the divergences are cancelled consistently. We also discuss the implementation of an invariant mass cut. For the practical calculation we split the total cross section into three separately finite parts, which are computed, each within its own computational scheme, in Sections 3, 4 and 5, respectively. Section 6 describes a consistency check we performed for our results and the comparison with some previous results. In Section 7 we analyze the importance of the various contributions for the threshold scan including initial-state radiation. We conclude in Section 8. Several appendices collect technical results, in particular the implementation of the new results into the QQbar threshold code.
2 Setup of the computation 2.1 Resonant and non-resonant separation in unstable particle EFT Precision calculations of top pair production near threshold are most conveniently done in potential non-relativistic effective field theory (PNREFT) [33,34], which describes the dynamics of slowly moving particles with three-momentum m t v coupled to ultrasoft radiation/massless particles with energy m t v 2 after hard and soft effects have been integrated out. The computation contains uncancelled divergences proportional to the top-quark width, which start at NNLO in dimensional regularization. The top-pair production cross section is thus an ill-defined quantity. Instead one must consider the final state of the decay products bbW + W − + X. The narrow-width approximation is not applicable since the top width is not small compared to the top kinetic energy E = √ s − 2m t ∼ m t v 2 . 3 The above final state can also be produced non-resonantly, i.e. without an intermediate non-relativistic top pair. The resonant and non-resonant production mechanisms cannot be distinguished physically and must be summed. Only the sum is well-defined and finite-width divergences must cancel [15]. This cancellation has already been demonstrated up to NNLO [27], and will be reproduced in the computation of the full NNLO correction in this paper.
To account for the non-resonant production mechanism, one must embed the effective theory framework for the QCD result [10] into Unstable Particle Effective Theory [35,36]. The complete NNLO cross section can be written as the sum of a resonant and a nonresonant contribution σ NNLO (s) = σ NNLO res (s) + σ NNLO non-res (s). (2.1) The resonant contribution has the form σ NNLO res (s) ∼ Im k,l It is understood that the imaginary part refers only to discontinuities of the forward amplitude that correspond to a bbW + W − X final state. 4 The production operators O (l) annihilate the incoming e + e − states and produce a nearly on-shell top and anti-top quark with small relative velocity. The matrix element is evaluated within PNREFT, appropriately generalized from QCD to account for electroweak effects and top decay.
In addition one must consider the interactions of the energetic initial-state electrons.
The C (l) are the hard matching coefficients of the production operators. They also receive electroweak corrections and furthermore acquire an imaginary part from diagrams involving cuts corresponding totbW + and tbW − final states. The imaginary part arises, for example, from the interference of the process e + e − → W W * , where the off-shell W decays totb with the process e + e − → tt, where the on-shell t decays to W b. In unstable particle theory this contribution appears in the resonant term, since the separation into resonant and non-resonant is done strictly on the basis of the virtuality of the top propagators, which in this example is small for both t andt. It originates from cuts over hard propagators that correspond to the physical final state bbW + W − X. Hard cuts over the tt final state are not possible kinematically near threshold. Thus, the leading corrections are fromtbW + and tbW − cuts and are of the order α 3 EW , which constitutes a NLO contribution to the cross section σ LO ∼ α 2 EW v. The nonresonant term arises from expanding the full-theory diagrams in E. Since both E and α EW count as two orders in the expansion, the NNLO contribution is given by the QCD O(α s ) corrections to the process e + e − →tbW + + tbW − , computed directly at the threshold √ s = 2m t , while actual bbW + W − cuts as well as electroweak and E/m t corrections are of the order α 4 EW and only contribute at NNNLO. The construction implies that the poles of internal top propagators in the non-resonant contribution are not regulated by a finite-width prescription, since any width terms would have to be expanded out. This leads to singularities at phase-space boundaries (p b + p W + ) 2 → m 2 t , which must be regulated dimensionally. The 1/ǫ poles cancel exactly the finite-width divergences that appear in the resonant contribution. The computation of the QCD correction to the process e + e − →tbW + + tbW − with this specific prescription, required for consistency with the resonant PNREFT calculation in dimensional regularization, is the major result of the present work.

Organization of the computation
We now discuss the structure of the phase-space endpoint divergences in more detail. The clarification of their diagrammatic origin allows us to divide the sum of resonant and non-resonant NNLO contributions into several separately divergence-free parts, and this separation determines the organization of the actual calculation. The cross sections of the processes e + e − →tW + b and e + e − → tW −b are equal by CP symmetry, hence we shall only consider the final statetW + b below and multiply the result by two in the end. In unitary gauge the NLO non-resonant contribution is given by the diagrams shown in Figure 1 [26]. At NNLO real and virtual gluon corrections must be considered. While this appears to be a standard NLO QCD correction computation to a 2 → 3 process, existing automation tools can nevertheless not be employed due to the endpoint divergences, which are present in addition to the usual UV and IR singularities.
To illustrate this issue, we consider the phase-space integral of a virtual diagram such as h ix below, where the integrand f ix is a Lorentz scalar, i.e. it only depends on scalar products of its arguments. This allows us to define The Heaviside function accounts for the optional cut on the invariant mass of the top quark as will be discussed in Section 2.3. Since the bottom quark mass can be safely neglected for this calculation, for the total cross section y = m 2 W /m 2 t . The real corrections can be brought into the same form as (2.4) with the variable t * ≡ (p W + + p b + p g ) 2 /m 2 t instead of t. The endpoint divergences originate from the region t → 1, where the integrand becomes singular due to negative powers of (1 − t) = (m 2 t − (p W + + p b ) 2 )/m 2 t , which stem from top-quark propagators becoming resonant. In [27] the leading terms in an expansion around t = 1 of the integrands g ix (t) were obtained using the expansion by regions approach [37,38]. The remaining t-integration for the expanded result is trivial, The divergent integrals with a ≥ 1 are regulated dimensionally by the bǫ in the exponent, which is inherited from the d − 1 dimensional phase-space integral. At NNLO endpointdivergent integrals with a = 1, 3/2, 2 are present, but only those with a = 1 manifest as 1/ǫ poles. This is related to the well-known property of dimensional regularization, that it renders some power-like divergent integrals finite for ǫ → 0. It is obvious from (2.6) that the integrands g ix (t) must not be expanded in ǫ, because it would spoil the dimensional regularization of the endpoint divergences. This implies that the loop integrals in the virtual corrections cannot be expanded in ǫ, since even the tree-level phase-space integration is divergent. Expressions for scalar one-loop integrals in general d dimensions with up to four external legs were obtained recently [39], but Figure 2: Gluon corrections to the tree-level diagram h 1 . This set of endpoint divergent diagrams is UV and IR finite and will be denoted as the squared contribution. Symmetric diagrams and diagrams with tW −b (g) cuts are not displayed.
a simpler strategy is to take the results for the endpoint divergent terms from [27] as subtractions to the complete integrand. The integrals (2.4) are decomposed as follows: (2.7) where the required coefficientsĝ (a,b) ix of the series expansion in (1 − t) are available up to order O(ǫ 0 ) from [27]. This renders the t-integration on the right-hand side finite and allows us to expand the subtracted expression in the square bracket in ǫ. Thus, the Figure 3: Additional endpoint singular diagrams for the NNLO non-resonant part. This set of endpoint divergent diagrams also contains UV divergences and will be denoted as the interference contribution. Symmetric diagrams and diagrams with tW −b cuts are not displayed.
integral can be performed numerically. Additionally, we require the O(ǫ) contributions toĝ (1,b) ix , because the coefficients with a = 1 are multiplied with a 1/ǫ pole in (2.7). In total the NNLO non-resonant correction requires the evaluation of the order of 100 diagrams obtained by attaching one gluon to the diagrams in Figure 1 in all possible ways. Fortunately only about 15% of those contain endpoint divergences. They have been identified in [27] and are shown in Figures 2 and 3. They are computed manually by applying the subtractions (2.7). The remaining large number of finite diagrams is computed in an automated fashion using suitably edited MadGraph code. This latter contribution will be referred to as the automated part σ aut .
The endpoint divergent diagrams are divided into two parts. The first is given by the QCD corrections to the diagram h 1 , shown in Figure 2, and is denoted as the squared contribution σ sq . It is UV and IR finite, because it includes the complete virtual, real and counterterm contributions to h 1 . The remaining endpoint divergent diagrams, shown in Figure 3, are referred to as the interference contribution σ int . In addition to the endpoint divergences, the interference part contains UV divergences, which are cancelled by endpoint-finite counterterm contributions contained in the automated part. We disentangle the two types of divergences by performing the subtraction (2.7) and obtain is endpoint-divergent but UV-finite, and is endpoint-finite but UV-divergent. In total, this allows us to split the non-resonant . . .

Figure 4:
The middle panel shows the diagrams accounting for the bare absorptive contribution to the matching coefficients C (k) . The lower panel sketches the respective contribution to the cross section, where the ladder exchanges of gluons cause the Coulomb singularities (α s /v) k and are, therefore, of the same order. Only the diagram with a single gluon exchange contains a 1/ǫ pole, which implies a scheme dependence of the finite part. We obtain the same diagrams (up to symmetry) as in Figure 3 by restoring the full theory graphs in place of the effective operators, i.e. by replacing C contribution into the following parts Only the first two terms contain endpoint divergences. The third term, enclosed in square brackets, is finite. The endpoint divergences cancel with the resonant contribution. Specifically, the endpoint divergence σ (EP div) int of the interference contribution is compensated by σ C (k) Abs,bare from the bare absorptive parts C (k) Abs,bare of the hard matching coefficients C (k) appearing in (2.2). The C (k) Abs,bare are given by the diagrams in the upper and middle panel of Figure 4, which have a direct correspondence to the diagrams h ia in Figure 3. Following this observation we split the resonant contribution into two parts, Abs,bare + σ res, rest , (2.12) where the remainder σ res, rest contains various contributions described in detail in Section 3. Here, we only point out that σ res, rest cancels the endpoint divergence of the squared contribution. Thus, we can now split the cross section into three separately finite parts . (2.13) The finiteness allows us to evaluate each of the parts (I), (II) and (III) in a different computational scheme. They will be computed in Sections 3, 4 and 5, respectively. An overview over the divergences that appear in the individual parts (I), (II) and (III) is given in Table 1. The computations are performed in the top-quark mass pole scheme.
The results are then converted to an IR renormalon-free mass scheme for the numerical evaluations performed in Section 7. The bottom-quark mass is neglected in all contributions except σ aut , where the default value m b = 4.7 GeV of MadGraph is used unless indicated otherwise.

Implementation of a "top invariant mass cut"
The main result of this work is the non-resonant NNLO correction to the full cross section σ(e + e − → bbW + W − X), but we also present results with loose cuts on the invariant mass of the top and anti-top quark where p t (pt) denotes the momentum of the (anti-) top quark. The implementation of cuts in the effective field theory framework has been discussed in [29] and depends on the scaling of the cut parameter ∆M t with respect to the power counting parameters of the EFT. The cut is termed "loose" when ∆M t ≫ Γ t . Thus, a loose cut never affects the resonant contribution to the cross section, where the off-shellness of the tops is parametrically of order m t Γ t . Since the physical final state is bbW + W − X without reference to whether it was produced through an intermediate top or anti-top, it is necessary to define what is meant by the (anti-) top momentum. An invariant mass cut of the form (2.14) was already implemented in the NLO non-resonant calculation [26], but since at this order the partonic final state is always bbW + W − , one simply defines p t = p W + + p b , pt = p W − + pb. At NNLO the final state may contain an additional gluon and a definition of the observable is required at the hadronic and the partonic level. Any sensible definition of an observable called "top momentum" should be equal up to an amount of order Γ t to the momentum of the top quark in the hypothetical limit that the top quark were a stable particle. It should also lend itself to simple, infrared-finite theoretical computations. On the other hand, the assignment of top momentum to the final state of a non-resonantly produced bbW + W − X event is rather arbitrary and a matter of definition subject to UV finite IR finite EP finite correspond to the virtual plus counterterm and real contributions to σ sq , respectively. With ⋆ we indicate contributions that are endpoint divergent by power counting, but finite in dimensional regularization. The contributions that make up σ res, rest are defined in Section 3. the previous two requirements. Because non-resonant production is a sub-leading effect, different definitions differ only by small amounts.
In the following we describe an algorithm that satisfies these requirements. The algorithm is most likely not optimal and serves only as a proof of concept. In the first step we cluster any hadronic or partonic event into the objects W + , W − , b-jet,b-jet and other jets. For the purpose of this discussion energetic leptons and photons are also among the "other jets". We require that the event contains exactly one W + , one W − and at least one b-and oneb-jet. 5 Any jet algorithm can be used to define these objects.
In a second step we group the above pre-defined objects into exactly two clusters. The momenta of these two clusters define the (anti-) top momentum. This fulfills the abovementioned requirement, since close to threshold in an event with a resonant top and anti-top, the momentum of any other particle can be at most of order Γ t .
To implement the second step, assume that the event contains N other jets and let S = {p Ji , i = 1 . . . N} be the set of jet momenta. A partitioning of S consists of two disjoint sets A,Ā such that A ∪Ā = S. The momentum of the "top cluster" and the "anti-top" cluster for a given partitioning are defined as If there is more than one b-jet (b-jet) in the event, p b (pb) refers to the most energetic b-jet (b-jet), and the remaining ones are considered to be part of the set S. We then define p t and pt by the value of p tA and ptĀ, respectively, of the partitioning A, which minimizes the product An event passes the cut (2.14), if the so-defined top momenta satisfy the inequality (2.14) . We now apply this definition to the partonic NNLO calculation. The partonic final states are tt only at LO, tt and tW −b ,tW + b at NLO, and at NNLO the previous and tW −b g,tW + bg. Here t (t) means a final state that can originate from on-shell (anti-) top decay at the given order, with invariant mass within m 2 t by an amount of order m t Γ t . Consider first the tt final state. One might think that the t (t) decay products are automatically clustered into the correct parent p t (pt) and hence the event always passes the loose cut as desired. However, this is not always true, since an energetic gluon from top decay may be radiated collinear to theb from anti-top decay, in which case the jet algorithm merges it into theb-jet rather than the b-jet. In this case the loosecut condition (2.14) might not be satisfied, and the event is missed even though both tops were produced resonantly. 6 Since our calculation does not include the kinematics of resonant top decay, we cannot correct it for this misassignment. However, the probability for such a misassignment is at most of order where R is the half-opening angle of theb-jet. The first factor represents the suppression of energetic, large-angle gluon radiation and the second the jet area on the unit sphere.
practice, this will be inefficient although a Monte Carlo study for the top forward-backward asymmetry at the ILC concluded that the discrimination between bottom and anti-bottom jets can be achieved with a purity of 80% at about 60% efficiency [40]. We also do not discuss the non-trivial combinatorial problem of reconstructing the W bosons from their hadronic decay in the presence of additional jets, since in the logic of our discussion the W bosons are considered as stable particles. 6 If the gluon is not energetic but ultrasoft with momentum of order Γ t , the misassignment is irrelevant and the loose-cut condition remains satisfied.
The numerical estimate is obtained for R ≈ 0.3. We can therefore safely neglect this error.
Next we consider the non-resonant final state t W −b . (With obvious modifications the following discussion applies to the CP-conjugate final state.) To NNLO accuracy, on-shell top decay must be taken in NLO, and may contain an additional gluon. Whenever there is no gluon or the gluon is merged with the b-orb-jet, the set of partitionings is empty and the definition of p t and pt is the sum of the appropriate W and b-jet momenta. The loose cut is passed except when the gluon is misassigned to theb-jet as above, but in this case the probability for this to occur is even further suppressed due to the suppression of non-resonant production in the first place. If, on the other hand, the jet algorithm returns an additional (gluon) jet, there are two partitionings, one where the gluon jet momentum is (correctly) added to the top decay, i.e. to p W + + p b , and the other, where it is not. The first, correct, possibility will almost always minimize χ in (2.16) and then satisfy the loose-cut condition, whenever the invariant mass of the non-resonant W −b pair is larger than ym 2 t , where y ≡ (m t − ∆M t ) 2 /m 2 t . Hence, imposing the cut yields a single Heaviside function θ((p W − + pb) 2 − ym 2 t ) in the phase-space integral, as in (2.4). The other partitioning where the additional gluon jet is incorrectly combined with the non-resonant W −b to form the anti-top momentum can minimize χ only if the invariant mass of the W − ,b-jet and gluon jet accidentally adds up to m 2 t within an amount m t Γ t . This possibility is suppressed by the NNLO probability for the process to happen in the first place times the small phase-space fraction where the kinematic requirements for misassignment are satisfied, and hence can be neglected at NNLO.
Finally we discuss the t W −b g final state, which appears at NNLO in the non-resonant part. At NNLO, it is sufficient to consider the resonant top quark decay in the tree approximation. Hence, the discussion of the W + W − bbg final state from above can be repeated, except that now the partitioning that minimizes χ with overwhelming probability is the correct combination of the gluon jet momentum with the non-resonant W −b pair. Hence, up to a negligible error, the loose cut (2.14) is implemented in the realemission phase-space of the NNLO non-resonant contribution as the Heaviside function θ((p W − + pb + p g ) 2 − ym 2 t ). We have not implemented other cuts, but it is in principle straightforward to do so as long as they are loose. A general cut is a function c(p i ) of the external momenta, which evaluates to one if the event passes the cut and to zero otherwise. We define the complementary cut asc(p i ) = 1 − c(p i ). Assuming that c(p i ) is loose, the complementary-cut cross section σ(c) is purely non-resonant and free of endpoint divergences. It can therefore be computed with automated NLO parton level event generators such as MadGraph [41]. The non-resonant contribution with the original cut is then given by subtracting σ(c) from the total W + W − bbX cross section, where the cancellation of divergences between the resonant and non-resonant parts has already been taken care of. This approach will also be exploited in Section 6.1 to perform a powerful check of our computation.
A generalization to arbitrary cuts would affect the resonant contributions and is beyond the scope of this work. Recently, first results of an implementation of the fully differential cross section with NLL accuracy near threshold matched to the NLO fixed order result have been presented [42], but the generalization of this method to NNLO accuracy as discussed here is not straightforward.

Part (I)
The scheme for part (I) as defined in (2.13) is fixed by the existing QCD results for σ NNLO res . The resonant QCD cross section factorizes into a leptonic tensor L and the correlation function of two top-quark currents, Π(q 2 ). The former is evaluated in four dimensions; the latter completely in d dimensions in the naive dimensional regularization scheme (NDR). The squared contribution contained in part (I) factorizes into the same leptonic tensor L and a hadronic tensor H, and the same conventions must be applied. We compute this part in Section 3.2. The electroweak NNLO corrections to the resonant part must also abide by this scheme (except for σ C (k) Abs,bare contained in part (II)), and we consider them first.

Resonant electroweak effects
Electroweak corrections to the resonant cross section are computed in the non-relativistic EFT framework extended from QCD to the full Standard Model. We consider them to NNLO in the counting scheme (1.1).
For ease of reference and to set up notation, we briefly recapitulate the well-known expressions for the LO cross section [32] where σ 0 = 4πα 2 /(3s) is the high-energy limit of the photon-mediated muon pair production cross section at leading order, E = √ s − 2m t and Γ is the on-shell top-quark width as defined below. At LO the top pair is produced via s-channel exchange of a photon or Z boson. The couplings of the fermions to the Z boson are given by where T f 3 is the third component of the weak isospin of fermion f , e f is the fermion electric charge measured in units of the positron charge, and s w and c w are the sine and cosine of the Weinberg angle, respectively. The S-and P-wave production operators are given by 7 with leading-order matching coefficients (3.8) Here ψ (χ) is the non-relativistic top (anti-top) field and e c i denotes the effective field (as defined in soft-collinear effective theory (SCET)) of an energetic electron moving in the light-like direction n µ i . In the present context, the directions n 1 and n 2 are set by the electron and positron beams, respectively. The collinear electromagnetic Wilson lines have been introduced to make the operators invariant under collinear gauge transformations in SCET, as well as the light-like vectorsn i with n i ·n i = 2. The factor of 4πα/s has been absorbed into the operators to render the coefficients dimensionless and of order one. The P-wave production operators and their Wilson coefficients will be required below. Note that because the cross section is constructed as an expansion in E, the energy-dependence of the s-channel photon and Z boson propagators could be expanded around s = 4m 2 t , in which case the short-distance matching coefficients would be truly energy-independent. However, we apply a convention where we keep the full s-dependence in the s-channel propagators, which therefore appears in (3.3) to (3.8).
The renormalization scheme for the electroweak parameters adopted here is the (m W , m Z , α(m Z )) scheme. The Weinberg angle is then given by c 2 . The electromagnetic coupling α em from now on is denoted by α, where α refers to the scale dependent electromagnetic coupling α(µ α ) defined through the photon vacuum polarization, which interpolates between the fine-structure constant α 0 = α(0) and the input parameter α(m Z ).
In (3.1) G 0 (E+iΓ) denotes the non-relativistic zero-distance Coulomb Green function in dimensional regularization [43,44], which describes the propagation of the top-anti-top pair at LO in the non-relativistic EFT. It is expressed through the variable where ψ is the logarithmic derivative of the gamma function. Furthermore, the logarithm L λ = ln(λµ/(m t α s C F )) = − 1 2 ln(−4mE/µ 2 ) appears. 8 After separating σ C (k) Abs,bare from the NNLO resonant contributions as explained around (2.12), the remaining parts are The pure QCD S-wave contribution σ QCD has been obtained in the formalism employed here in [34]. Top-pair production in a P-wave state σ P-wave was computed in [22]. Higgs contributions σ H that only involve the top Yukawa coupling have been computed already up to NNNLO [25,45]; similarly the effect σ δV QED of the LO QED Coulomb potential δV QED = −4παe 2 t /q 2 [25]. At NNLO, top decay introduces additional contributions to the bilinear part of the PNREFT Lagrangian, which contribute σ Γ to the resonant cross section (Section 3.1.1). While there are no electroweak contributions to the nonrelativistic potential at NNLO (Section 3.1.2), there are electroweak corrections to the hard matching coefficients C (k) . The contribution σ C (k) EW from the real part of the hard matching coefficients is given in Section 3.1.4. Contrary to the QCD case the electroweak hard matching coefficients contain an imaginary part from cuts over all possible final states. ThetW + b (tW −b ) cuts contribute to the e + e − → bbW + W − cross section [46]. The imaginary part is split into a bare contribution σ C (k) Abs,bare (Section 4.1) and a contribution from field renormalization σ C (k) Abs,Z t (Section 3.1.3), because the two parts are treated in different schemes. Partial results for the mixed-QCD-electroweak corrections to the hard matching coefficients C (k) are available [45,47], but they only contribute at NNNLO and will not be considered here. Finally, we consider effects from initial-state radiation (ISR), σ conv IS (Section 3.1.5).

Finite-width corrections to the NNLO Green function, σ Γ
Additional terms appear in the PNREFT Lagrangian due to the instability of the top quark and its coupling to photons. The coupling of the top quarks to ultrasoft photons must be multipole expanded in the spatial component, just like the interactions with the ultrasoft gluons. Only the leading term is relevant at NNLO. However, its contribution vanishes, because the multipole-expanded field can only resolve the net electric charge of the top anti-top pair, which is zero. In analogy to QCD, the couplings in the Lagrangian (3.13) can be removed by a field 8 When the 1/ǫ pole is related to a finite-width divergence, we set µ = µ w and write G (w) 0 (E) and L (w) λ to distinguish the finite-width scale-dependence of the resonant contribution from the µ r scaledependence due to the strong coupling, cf. [14,22]. transformation involving an ultrasoft Wilson line (cf. [48]). The generalization of the bilinear part of the PNREFT Lagrangian is [35,36] where ψ(χ) is the non-relativistic top (anti-top) field and ∆ is a hard matching coefficient. It can be determined by matching the top propagator in the effective theory to the full theory. In the pole mass scheme we obtain where Γ is the pole width of the top quark defined through the gauge-invariant position of the pole of the top propagator in the complex p 2 plane. We note that with this convention (3.14) contains the term , which has the form of a mass shift. It can be absorbed into the definition of the pole scheme by adding −Γ 2 /4 to the right hand side of (3.16), which completes the square and defines a different convention used e.g. in [49]. Electroweak corrections to the top-pair production cross section near threshold have also been considered in [50]. The absence of the Γ 2 correction to the Green function in [50] implies that this different convention is also adopted there. Thus, one must be careful to account for this difference in the definition of the top pole mass when comparing their results to ours. The term (iΓ/2)(ψ † ψ − χ † χ) in (3.14) belongs to the LO Lagrangian and must be treated non-perturbatively. It leads to the replacement E → E + iΓ, which defines the QCD contribution, and makes the argument of the Green function in (3.1) complex. The two remaining terms in (3.14) that contain the width are of NNLO and can be treated perturbatively. Only two simple single insertions are required. We denote the correction to the Green function G 0 (E) from the terms ( , respectively. They are given by where ψ 1 (x) = ψ ′ (x) is the first derivative of the polygamma function. The corresponding NNLO contribution to the Green function is In the implementation of the cross section in the QQbar threshold code [32] the top width is treated as a parameter. This implies that higher-order corrections to the tree-level width Γ 0 are also treated non-perturbatively through the replacement E → E + iΓ. A subtlety arises at electroweak NNLO when this result is combined with the non-resonant contribution, which is computed in dimensional regularization. The pole part of the NNLO non-resonant contribution is proportional to Γ 0 with a finite part that follows from expanding diagrams up to O(ǫ 0 ). For consistency, the tree-level contribution to the width in (3.14) must be treated as a d-dimensional hard matching coefficient. Hence the O(ǫ) terms in the d-dimensional tree-level expression of the top width contribute finite terms to the resonant part from their multiplication with the finite-width 1/ǫ poles. These finite parts are not included when Γ is treated as a fourdimensional numerical parameter, and must be added separately. 9 The LO pole width, which is required in d dimensions, is given by where x W = m 2 W /m 2 t and µ w denotes the scale related to the finite-width divergences as discussed in [14,22]. The contribution from the O(ǫ) terms of (3.20), which multiply the finite-width divergence contained in (3.19) and the one in the pure QCD result, to be added to the cross section is where Γ 0 is the ǫ → 0 limit of (3.20). On the whole, we obtain In the numerical evaluation we resum the perturbative corrections to the would-be toponium bound-state poles. Due to the instability of the top quarks, we are dealing with a non-Hermitian Lagrangian, cf. (3.14). The implications have been discussed in [53]. The positions of the would-be toponium poles are the complex eigenvalues of the non-Hermitian Hamiltonian, where E n is the bound-state energy assuming stable top quarks and Γ n = 2Γ + δΓ n is the total inclusive width of the bound state. In accordance with the earlier discussion, the top-quark width Γ is treated as a parameter. The corrections δΓ n describe the effects of time dilatation on the top decays due to the residual movement of the top quarks and the annihilation of the would-be toponium state through strong (e.g. tt → ggg) or electroweak (e.g. tt → l + l − ) interactions. The eigenstates of a non-Hermitian operator do not form an orthogonal basis of the Hilbert space [53]. This implies, that the completeness relation must be modified. We consider the sets of right and left eigenstates 10 Assuming that the Hamiltonian transforms as under time reversal, and that the eigenvalues are non-degenerate, the eigenstates can be normalized such that they form a bi-orthogonal set [53] m | n = δ mn , (3.26) which implies that the completeness relation takes the form The property (3.25) implies that the state |ñ is exponentially growing at the same rate as |n is decaying, which facilitates the normalization (3.26). After applying (3.27) the Green function takes the following form near the poles which generalizes the expression for the QCD result [14]. The pole position and residue of (3.28) have the following perturbative expansion with the LO expressions n (x). This holds, because the non-Hermitian part of the LO Hamiltonian is proportional to the identity operator and thus only affects the eigenvalues E (0) n but not the eigenstates |n (0) .
The results for the non-relativistic Green function in perturbation theory do not take the form (3.28), but contain higher-order poles This allows us to read off the NNLO correction to the bound state parameters from the contribution (3.19) to the Green function As discussed above, the Γ 2 term in (3.14) has the form of a mass shift and therefore leads to an n-independent correction (3.35) to the position of the pole, while it does not affect the residue. The iΓ∂ 2 term in (3.14) accounts for time dilatation, which reduces the total width of the would-be toponium resonance by (3.36). Since it is non-Hermitian, it also makes the residues complex due to (3.37). The corrections to the bound states from QCD effects as well as the QED Coulomb and Higgs potentials can be found in [14] and [25], respectively. Using this input we can resum the higher-order poles in the expanded Green function by the replacement where the expanded term has the form (3.34). In the actual implementation [32] we apply the pole resummation procedure to G(E) alone in the electroweak contributions, but to the entire current correlation function of vector and axial vector currents in the pure QCD contributions. Figure 5: Cancellation of the electroweak gtt vertex correction and the QCD γtt vertex correction in the on-shell scheme. The vector current itself is conserved and therefore not renormalized.

Mixed QCD-electroweak NNLO corrections in PNREFT
In addition to the kinetic terms (3.14) the PNREFT Lagrangian contains potential interactions. We now show that there are no mixed QCD-electroweak corrections at NNLO. The construction of PNREFT proceeds by first integrating out fluctuations at the hard scale, which yields NREFT, and then integrating out fluctuations at the soft scale. In the first step, one must consider electroweak corrections to the hard matching coefficients of the QCD vertex, and vice versa. The relevant diagrams are shown in Figure 5. The loop momenta in Figure 5 are hard, while the momenta of the external particles can be either soft, potential or ultrasoft and must be expanded out of the loop integral. Therefore the vertices are effectively evaluated at zero external momentum, and the corresponding contribution is exactly cancelled by the on-shell external field renormalization factor.
The potentials are determined in the matching procedure between NREFT and PN-REFT. The diagrams that contribute to the 1/q 2 potential at order α s α are shown in Figure 6, where the momenta of the external top quarks are potential and the loop momentum is soft. The contributions of the first and second diagram are identical, and are exactly opposite to those of the third and fourth diagram, which implies that the sum of the diagrams in Figure 6 vanishes. We have not drawn the four diagrams that involve soft vertex corrections to the tree-level potentials, because these corrections are scaleless and vanish in dimensional regularization. Last, but not least there are no contributions from insertions of the one-loop corrections to the hard matching coefficients in the treelevel potential, because these coefficients vanish as argued above. We conclude that no mixed QCD-electroweak potentials appear at NNLO.
Furthermore, we demonstrate that the potential does not receive any pure electroweak NNLO corrections from the exchange of Z bosons. We count the mass of the Z bosons as hard and therefore have to integrate out the Z boson in the hard matching to NREFT. This implies, that all interactions that are mediated by the Z boson in the full theory + + + become local in PNREFT. The Z-boson exchange potential corresponds to the full-theory diagrams shown in Figure 7 and is proportional to α EW /m 2 t in momentum space. Thus, it is suppressed by (α EW /α s ) × (q 2 /m 2 t ) ∼ v 3 compared to the LO colour Coulomb potential and only contributes at NNNLO.
Finally, we comment on the so-called 'jet-jet' interactions that were considered in [54]. These are corrections involving gluon emission from the final-state bottom quarks and it was demonstrated in [54] that they vanish at NLO. In their calculation, the authors of [54] first consider the subgraph I µ , which corresponds to the cut to the right of the gluon of the third diagram in Figure 5. Their result for I µ scales as √ α s Γ 0 /|k|, where k is the gluon three-momentum, which is either potential or ultrasoft in their case. They then show by explicit computation that all NLO corrections that involve the subgraph I µ and/or its CP conjugate J µ vanish. In our approach, where loop integrals are strictly expanded according to the scaling of the momentum regions, the only non-vanishing contribution to the subgraph I µ comes from the region of hard loop momentum, where no inverse powers of |k| can appear because the external momenta are expanded out. Therefore, the absence of any 'jet-jet' interactions at NLO is a matter of simple power counting, which implies that corrections can first appear at the relative order α EW , i.e. at NNLO. We have already proven that there are also no 'jet-jet' interactions at NNLO by demonstrating that electroweak corrections to the QCD vertex in Figure 5 vanish.

Absorptive part from field renormalization, σ C (k)
Abs,Z t The hard matching coefficients C (k) become complex at NNLO due to bW + loop corrections. The imaginary part contributes to the finite-width divergence of the resonant cross section σ res, rest in (3.12) and, thus, it has to be determined in d dimensions in accordance with the scheme used to evaluate the other components of part (I). The bare absorptive contribution to C (k) on the other hand, is contained in part (II) and there-fore has to be computed in a different scheme. Since the two parts are treated using different conventions, we find it convenient to separate them in notation. The matching coefficients up to NNLO are expanded as The real part of the electroweak corrections, C EW , does not yield a finite-width divergence at NNLO. Thus, it is not necessary to split it as well. The absorptive part (3.40) of the hard matching coefficient is available in four dimensions [46]. We have repeated the calculation of the individual contributions in the schemes described above. In four dimensions we reproduce the result of [46], where the normalization factor is necessary because the definition of the hard matching coefficients in [46] differs from (3.39).
The bare part C (k) Abs,bare will be given in Section 4.1. From the field renormalization, we obtain The contribution to the NNLO cross section is given by where the finite terms from the multiplication of the 1/ǫ divergence in the real part of the Green function (3.10) with the O(ǫ) parts of (3.42) and (3.43) are included.

Electroweak contributions to the hard matching coefficient
The real part of the electroweak contributions to the NNLO matching coefficients C (k) has been computed in [55][56][57]. Pure QED corrections have been neglected there. Therefore, we split C The hard QED vertex correction to the γe + e − and Ze + e − vertices contains divergences that cancel among initial-state radiation (ISR) contributions (see Section 3.1.5). We therefore assign it to σ conv IS to render both, σ C (k) EW and σ conv IS , separately finite. There is no contribution from the box diagram involving two photons, since only its interference with the production of the top pair through the vector component of the s-channel γ or Z boson is of NNLO and the correlator of three vector currents vanishes [58]. The box diagram with a photon and Z boson is considered to be a non-QED correction to the photon-exchange contribution and is therefore already part of C (k) WZ . Thus, the only pure QED effects are the hard photon vertex correction to the γtt and Ztt currents and the photon self energy, which yield As in [59], the renormalized photon self-energy Π AA R (s, µ 2 α ) is defined in the scheme of [60], and will be discussed below. The non-QED contributions are where C ew V,A (ν = 1) is given in [57], α 0 is the fine-structure constant and Π AA R (4m 2 t , 0) coincides with the expression for Π AA R from [57]. The subtraction terms are present because Higgs effects which only involve the top Yukawa coupling have already been included separately as part of σ H in [25] and the photon self-energy is contained in the QED contribution (3.46). Corrections that involve Higgs couplings to gauge bosons or Goldstone bosons remain in (3.47). We note that the photon self-energy terms in (3.46) and (3.47) differ, because we use a renormalization scheme which is different from [55][56][57]. The matching coefficients given in [55][56][57] are expressed in terms of the fine-structure constant α 0 . This scheme suffers from a large spurious dependence on the light fermion masses, that cancels explicitly with the self-energy corrections to the matching coefficients when the fine-structure constant is expressed in terms of a less infrared-dependent definition of the electroweak coupling constant. Therefore, we write the cross section in terms of the running on-shell coupling α ≡ α(µ α ) from [60]. In this scheme the renormalized photon self-energy takes the form where the bare self energy Π AA is taken from [55]. The explicit factor 1/µ 2 α appears, because [55] defines the photon vacuum polarization Π AA (s) as a dimensionful quantity and does not imply a power-dependence of the cross section on the scale µ α . In the limits µ α → 0 and s → 4m 2 t the scheme of [60] converges to the scheme of [55][56][57], i.e. α → α 0 and Π AA (µ 2 α )/µ 2 α → Π ′,AA (0), and the self-energy terms in (3.46) and (3.47) coincide with each other and with the respective expression in [57].
The cross section receives the contribution from the electroweak corrections to the hard matching coefficients of the production operators.

IS
In this section we take into account effects from QED initial-state radiation. As such we count all corrections that involve an additional photon attached only to the external e ± states relative to the LO cross section. ISR was already considered in the 1980s [8] at the leading logarithmic order, i.e. summing corrections of the form (α ln(s/m 2 e )) k to all orders, whereas later works concentrated on the 'partonic' tt cross section. We extend the treatment of ISR to NNLO+LL accuracy below. The non-resonant part is only affected by QED radiation at NNNLO and will not be considered. With the exception of the effects from the hard momentum region, all contributions are universal and our treatment closely follows the one for W pair production near threshold in [28,29]. In fact, the equations in this section can often be obtained directly from those in [29] by substituting c (1,fin) p,LR → −4 + π 2 /12, and by adapting the different tree-level process. 11 When the electron mass is neglected the ISR contribution involves the hard, k ∼ m t , and ultrasoft, k ∼ m t v 2 , momentum regions, and in addition two hard-collinear momentum regions,n i · k ∼ m t , n i · k ∼ m t v 2 , k i⊥ ∼ m t v (i = 1, 2) familiar from SCET, where n i ,n i are pairs of light-like vectors with n i ·n i = 2 defined by the electron (i = 1) and positron (i = 2) momentum. Real collinear emission is kinematically forbidden in the resonant part, because it carries away a hard momentum fraction, which pushes the top pair off-shell. Virtual collinear corrections are scaleless. Hence, the hard-collinear regions vanish [28], and we are left with the hard and ultrasoft contributions. We evaluate these separately. A hard photon cannot be exchanged between the incoming and outgoing electrons, since this would also push the top pair off-shell. Thus the only correction from the hard region is the QED γ/Zee vertex correction which contributes to the hard matching coefficients C (v,a) . We find As it should be, this agrees with the QCD analogue of the hard matching coefficient of the vector current to SCET, first obtained in this context in [61] in DIS kinematics. We only kept the real part, because the imaginary part comes from cuts that do not correspond to the final state bbW + W − . The correction to the cross section from hard ISR is σ The contributions from the ultrasoft momentum region are shown in Figure 8. Virtual ultrasoft corrections are scaleless. The diagram with the photon attached to incoming and outgoing electron vanishes, because it is proportional to the square of the light-like direction n 1 of the electron beam. No ultrasoft corrections that couple to the collinear and non-relativistic sector occur at NNLO, because the leading ultrasoft photon coupling to the final state vanishes, as discussed in Section 3.1.1. Thus, the contribution to the cross section from the ultrasoft region is due to the right diagram in Figure 8 and reads (3.52) When the small electron mass is neglected, the photon radiation corrections are given by the sum of (3.51) and (3.52). We observe that the 1/ǫ 2 pole cancels, but a collinear divergence remains, because the cross section is not infrared safe for m e = 0. Subtracting this divergence defines a scheme-dependent 'partonic' cross section.
The divergence is regularized by the non-zero electron mass, which in turn yields large logarithms ln(s/m 2 e ). They can be resummed into an electron distribution function Γ LL ee (x), which describes the probability of finding an electron with momentum xp in the "parent electron" with momentum p. The cross section with resummed ISR from the electron and positron is given by Expressions for the structure function at leading-logarithmic (LL) accuracy can be found in [62][63][64][65], where LL implies that all terms of the form α n ln n (s/m 2 e ) are summed to all orders. The resummation of the next-to-leading logarithms (NLL) α n+1 ln n (s/m 2 e ) is crucial for the precision program at a future lepton collider, but the structure functions are presently unknown at this order.
At LO, the 'partonic' cross section σ conv (s) is given by (3.1). At higher orders it depends on the scheme used to regularize and subtract the collinear divergence. The scheme dependence cancels in the convolution with the structure functions. This implies that we have to adapt the results (3.51) and (3.52), which correspond to a minimal subtraction scheme, to the conventional scheme in which the structure functions Γ LL ee (x) are defined. This procedure has been described in detail in [28,29]. First, we need to convert the dimensional regulator of the collinear divergences to a finite electron mass regulator. Then, the O(α) terms that appear in the convolution of the structure functions with the LO cross section have to be subtracted from the fixed order NNLO partonic cross section to avoid double counting.
The first point is accomplished by noting that the presence of the additional scale m e ≪ m t v 2 introduces the additional momentum regions hard-collinear: with k 2 ∼ m 2 e and k 2 ∼ m 2 e v 4 , respectively. The soft-collinear region contributes in the diagrams shown in Figure 8. As before, the left diagram vanishes, and one finds The collinear 1/ǫ poles cancel in the sum of the hard and hard-collinear, and ultrasoft and soft-collinear contributions, separately. The collinear sensitivity is instead expressed through the large logarithms ln(4m 2 t /m 2 e ). The remaining singularities cancel in the sum over all regions. To make the cancellation explicit, one can expand the factor 1/k 1+2ǫ in the distribution sense: where a > 0 is arbitrary and we have introduced the modified plus-distribution We obtain which is finite, such that the four-dimensional expression (3.10) for the LO Green function can be used. The a-dependence cancels. We determine the subtraction terms by expanding the convolution of the LO cross section with the structure function in the coupling constant. We take the expression for the electron structure function from [65] with β ≡ β exp = β S = β H = 2(α/π)(ln(s/m 2 e ) − 1), given by The perturbative expansion of the structure function can be written as For the determination of the subtraction term only the limit x → 1 is important, The O(α) term in the convolution of the leading order partonic cross section with the structure functions is  x). We also set a = m t and neglected the difference between s and 4m 2 t in the argument of the logarithm. The initial-state QED correction to the partonic cross section in the conventional scheme for the electron structure function is given by (3.60) with (3.63) subtracted, resulting in The imaginary part of the Green function is neglected for E < −m t outside the nonrelativistic regime. The photon radiation contribution (3.65) to the cross section in this scheme is finite and free of large logarithms of the electron mass.

The squared contribution
In this section we discuss the calculation of the squared contribution σ sq in (2.13), which is given by the diagrams in Figure 2. The computer programs Package-X [66], Feyn-Calc [67,68] and LoopTools [69] have been employed for certain steps of the computation. The result for the scalar four-point integral in the diagram h 1b was taken from [70].
The individual diagram contributions to the hadronic tensor H are evaluated in d dimensions and written in the form (2.7). The numerical t (or t * ) integral contains all terms with positive integer or half-integer powers of (1 − y). With the exception of h 1b , the subtracted integrands were all obtained in analytical form. The integrand for h 1b contains an additional numerical angular integral. The expressions for the integrands are rather lengthy and will not be given explicitly. The numerical integrals are plagued by integrable singularities involving 1/ √ 1 − t and ln(1 − t) terms, that cause numerical instabilities in the evaluation of some diagrams. As a remedy, we computed additional terms in the expansion in (1 − t) analytically and used them as further subtractions.
The contributions corresponding to the second term on the right-hand side of (2.7) are given by the sum of the respective expressions in [27] and terms from the O(ǫ) contributions toĝ (1,b) ix . The latter encapsulate the dependence of the squared contribution on the computational scheme and are therefore specified below. In the notation of [27], we obtain 12 where H

(EP fin) 1x
is the contribution from the first term on the right-hand side of (2.7). The prefactor is defined as The other diagrams in Figure 2 do not contain 1/ǫ poles from the endpoint divergence and, therefore, no terms of this type are present, i.e. H 1x = H 1x | from [27] + H (EP fin) 1x . The contribution of an individual diagram h ix , g i to the non-resonant part is where n s is a symmetry factor, that is either two for diagrams which are symmetric with respect to the cut, or four for diagrams which are not symmetric with respect to the cut. The photon couplings are v γ f = −e f and a γ f = 0, where e f is the fermion charge measured in units of the positron charge, and couplings of the fermions to Z bosons are given by (3.2). The photon mass obviously vanishes, m γ = 0. In (3.67), O(ǫ) terms in the leptonic tensor have been discarded, as discussed at the beginning of this section. We have checked explicitly that IR and UV divergences cancel in the sum over the diagrams in the squared contribution.

Part (II)
It would be a natural choice to use the same scheme for part (II), given by (see (2.13

Absorptive contribution to the matching coefficient
The bare absorptive part of the matching coefficients C (k) Abs,bare is given by the diagrams shown in the second row of Figure 4. We define the scheme as follows: The coefficients C (k) Abs,bare are calculated in four dimensions, but the loop integrations in the third row of Figure 4, i.e. the ones related to the non-relativistic Green function, are performed in d dimensions. The Dirac algebra is completely treated in four dimensions. We describe in Section 4.2 how the interference contribution must be treated to achieve consistency with this scheme.
Our results for C (k) Abs,bare in four dimensions are The contribution to the cross section is given by Abs,bare Re We recall that at LO the Dirac structure of the top anti-top pair becomes trivial in the non-relativistic regime and only yields a prefactor 3 − 2ǫ. By introducing the factor 3/(3−2ǫ) in front of the Green function in (4.4), we adapted the expression to the scheme described above, which involves four-dimensional Dirac algebra. The contribution (4.4) is not affected by loose cuts.

Endpoint divergence of the interference contribution
The endpoint-divergent part of the interference contribution has been computed in [27]. The expression (2.9) also contains an endpoint-finite term from the O(ǫ) terms inĝ (1,2) ia multiplying the 1/ǫ pole. This term carries the dependence on the computational scheme and must, therefore, be treated in the same scheme as the contribution (4.4). We evaluate it using the expansion by regions approach described in [27]. For each of the diagrams in Figure 3, we treat the loop contained in the corresponding diagram in Figure 4, i.e. the right loop in h 2a and the left loop in h 3a and h 4a , in four dimensions. The Dirac algebra is also done in four dimensions, but the remaining loop integrations are performed in d dimensions. In the notation of [27] we obtain with I L WW = 1 for diagram h 3a with a photon attached to the W W vertex, and −c w /s w for the W W Z vertex. The endpoint divergent contributions of h 2a and h 3a follow from equation (3.67) with n s = 4. The contribution of h 4a is given by with symmetry factor n s = 4.
To verify that the treatment of the scheme is consistent we computed the finite sum of the contributions from the diagrams h 2a and h 3a and the contributions from the corresponding diagrams in C (k) Abs,bare also in the scheme of part (I) and found perfect agreement with the results presented above. Applying the scheme of part (II) simplifies the computation, especially for h 4a , since it avoids the more complicated integration of the left loop in h 4a in d dimensions.

Part (III)
The part (III) σ (EP fin) int + σ aut contains the automated part σ aut , which is evaluated with MadGraph. The automated part is UV divergent and therefore scheme dependent. The divergence and the corresponding scheme dependence cancel with the endpoint-finite part σ (EP fin) int of the interference contribution. We first describe the implementation of the automated part in MadGraph in Section 5.1. This fixes the scheme in which σ (EP fin) int is computed in Section 5.2.

The automated part
We first recall some aspects of MadGraph, which are relevant to our definition of the computational scheme.
1. The subtraction of IR singularities is performed automatically using the FKS method [71,72]. The IR singularities in the real corrections are subtracted before the phase-space integration and the subtraction terms are then integrated over the phase space of the real emission and added to the virtual corrections, where they cancel the IR singularities that arise in the loop integrals. The phase-space integration is then always done in four dimensions.
2. In the virtual corrections, MadGraph uses rational R 2 terms [73] to absorb the (−2ǫ)-dimensional parts of the numerators. For a given diagram with amplitude C the decomposition takes the form where D i = (l + p i ) 2 − m 2 i , quantities with a bar are (4 − 2ǫ)-dimensional and quantities without are 4-dimensional. The non-R 2 term can be written as a sum over 4-dimensional coefficients multiplying d-dimensional tensor integrals. The (−2ǫ)-dimensional parts related to the implementation of the 't Hooft-Veltman scheme [74] in MadGraph are all contained in the R 2 terms.
3. The amplitudes for the non-R 2 terms, the R 2 terms, the UV counterterms and the FKS subtraction terms are written as separate lists, each of them containing the coefficient of the 1/ǫ 2 pole, the 1/ǫ pole and the finite part. Afterwards, only fourdimensional operations are performed, i.e. the multiplication with the conjugated four-dimensional born amplitude and the four-dimensional phase-space integration.
Given the way that σ aut is defined, we never have to modify the construction of an amplitude A i , but we have to remove certain contributions A i A * j in the squared amplitude |A| 2 = i,j A i A * j . All the contributions associated with the diagrams in Figure 2 have to be removed, i.e. also the R 2 parts, the UV counterterms and the FKS subtraction terms.
There is however an ambiguity in the subtraction of the contributions in Figure 3, which determines the scheme in which σ (EP fin) int must be computed. We choose to only subtract the non-R 2 terms of h ia with i = 2, 3, 4. Following the discussion of the items 1 and 2 above, this implies that σ (EP fin) int has to be computed by using dimensional regularization for the tensor integrals. All other steps in the computation of σ (EP fin) int are then performed in four dimensions.
In the following, we describe the steps we performed in MadGraph to obtain the automated part in the scheme defined above. It is obvious that this cannot be achieved by modifying the process generation, because the automated part does not correspond to a squared amplitude. We therefore first generate the full process e + e − →tW + b including QCD corrections. By not invoking the complex mass scheme, we make sure that the selfenergy insertions are treated perturbatively. Hence, the cross section diverges rapidly for center-of-mass energies approaching √ s = 2m t from below. We remove the contribution from the endpoint divergent born diagram h 1 , the diagrams shown in Figures 2 and the non-R 2 terms from Figure 3 by editing the code generated by MadGraph. Finally, we have to deactivate some checks inside the code, that are invalidated by the modifications. MadGraph checks if the 1/ǫ 2 and 1/ǫ poles vanish for a number of phase space points. Here, this is not the case because the automated part is UV divergent. We have addressed this issue with in two ways -by deactivating the check or by performing a minimal subtraction of the UV divergence -and found agreement of both approaches. The minimal subtraction was also used to verify the cancellation of the UV divergence with the endpoint-finite part of the interference contribution σ (EP fin) int . Furthermore, due to the subtractions, the tree-level cross section and the real corrections are no longer the squared absolute value of an amplitude and, thus, no longer positive for all phase-space points. The positivity of these expressions is not necessary to make the code run properly, but is only used as an internal check [75]. Therefore, we can safely switch it off. The code can now be evaluated directly at the threshold √ s = 2m t . The contribution σ aut is given by the difference of fixed-order runs at NLO and LO, multiplied by a factor two to account for the tbW − contributions. Further details on the implementation and modifications in MadGraph are provided in Appendix B. The evaluation of the automated part in the code QQbar Threshold relies on a precomputed grid as described in Appendix A.5. Since the contribution σ aut is rather small, we do not aim for more precision than about 10% in the automated part. The resulting error of the cross section is less than one per mille. To reach this target precision we set up MadGraph to generate an integration grid from four iterations with 15000 points per integration channel and perform the actual integration using six iterations with 100000 points for each point of the QQbar Threshold grid. More precise results are possible at the cost of a considerably increased computing time for the generation of the grid.

Endpoint-finite part of the interference contribution
We recall that the endpoint-finite part of the interference contribution has the form (2.10). As detailed in Section 5.1 it must be evaluated by taking only the tensor integrals of the virtual loop in d dimensions and then performing all other steps in the computation strictly in four dimensions. Within this scheme, we have determined the endpoint subtracted integrands for the diagrams h 2a , h 3a analytically and for h 4a as a one-dimensional angular integral. We refrain from giving the lengthy results for the integrands. As described in Section 3.2 we use additional terms in the expansion of the amplitudes in 1 − t as subtractions to deal with integrable divergences that appear in the limit t → 1 of the t-integration. The result for σ (EP fin) int is given by applying the same prefactors and symmetry factors as for the endpoint divergent part of the interference contribution σ (EP div) int in Section 4.2.
6 Checks and implementation 6

.1 Consistency checks
Having performed the computation of the non-resonant part in the presence of the invariant mass cut (2.14), denoted by c ∆Mt (p i ), allows us to perform a very powerful numerical consistency check. The non-resonant cross section σ non-res (c ∆Mt ) in the presence of the complementary cutc ∆Mt (p i ) = 1 − c ∆Mt (p i ) is finite. Therefore, we can evaluate it using unedited MadGraph code. On the other hand, it can be obtained from our result by taking the difference σ non-res − σ non-res (c ∆Mt ). The comparison for various values of the cut ∆M t numerically tests the whole non-resonant result in (2.11), with the exception of the contributions from the O(ǫ) parts of theĝ (1,b) ix terms in (2.7), which originate from the t ( * ) → 1 region and are independent of the value of the cut, i.e. are not present in σ non-res (c ∆Mt ).
The result of our check is shown in Figure 9. Here, we have rearranged the contributions as follows, where σ h 1 (c ∆Mt ) is the contribution to the non-resonant part from the diagram h 1 at NLO (Figure 1) in the presence of the complementary cut. The line in Figure 9 shows our semi-analytical result for σ check (c ∆Mt ) obtained by means of (6.1). The points show the same quantity determined by evaluating (6.2) using MadGraph. The contribution from diagram h 1 is included in σ check (c ∆Mt ), because our edited MadGraph code, described in Section 5.1, corresponds to the combination σ aut (c ∆Mt ) − σ h 1 (c ∆Mt ) that appears in (6.2). We performed the same check for the individual contributions from the diagrams h ia with i = 2, 3, 4. In particular, this provides very welcome reassurance that the scheme dependence within part (III) has been treated consistently. Within estimated numerical errors we find good agreement, if the bottom-quark mass m b is neglected, as is done in our calculation.

Implementation in QQbar Threshold
All of the aforementioned NNLO corrections have been implemented in the new version 2 of the public code QQbar threshold [32]. A summary of the code changes and some code examples for the new functions are provided in Appendix A. QQbar threshold can be downloaded from https://www.hepforge.org/downloads/qqbarthreshold/. An updated online manual is available under https://qqbarthreshold.hepforge.org/.

Comparison to other approaches
While a complete calculation of NNLO electroweak and non-resonant contributions to the top-pair threshold as reported here has never been done before, NNLO non-resonant corrections have been evaluated in certain approximations in [50] and [31]. We briefly comment on these approximations and their limitations here.

Comparison to [50]
The leading NNLO non-resonant contributions for the case of "not-too-loose" cuts satisfying Γ t ≪ ∆M t ≪ m t were determined in [50] within the so-called phase-space matching (PSM) approach. This result captures the first terms in the expansion in Λ/m t (Λ 2 ≡ 2m t ∆M t − ∆M 2 t ) of the full non-resonant result, namely the terms of order m 2 t /Λ 2 , m t /Λ and (m t /Λ) 0 log Λ. The latter correspond to the endpoint-divergent terms computed in [27], which give the approximate result labelled "aNNLO" in Figure 14 below. Because of the "not-too-loose" cut condition, the PSM approach does not allow the calculation of the bbW + W − X total cross section near the top anti-top threshold.
The agreement between the PSM result and the full non-resonant computation of the non-analytic terms in the expansion in the invariant-mass cut parameter Λ/m t can be understood as a consequence of the cancellation of singularities between adjacent regions of loop momentum [27]. The Λ/m t non-resonant terms are obtained in the PSM approach by computing the ultraviolet behaviour of the resonant amplitude where the cut on the invariant mass of the top and anti-top quark has been implemented. Therefore Λ effectively acts as a regulator of the ultraviolet singular behaviour of the resonant part of the amplitude, that is obtained assuming on-shell top quarks, when the latter is further taken into the off-shell limit, i.e. for |p t | ≫ Γ t . On the other hand, the endpointdivergent terms are obtained from the non-resonant part of the amplitude, that assumes off-shell tops with |p 2 t − m 2 t | ≫ Γ t , upon going to the (infrared) on-shell limit within a distance regulated by Λ. The fact that both expansions provide the same divergent terms in Λ/m t is thus a consequence of the cancellation of the dependence on the cut-off Λ that separates the resonant and non-resonant regions. For the limitations on the PSM result to provide higher-order terms in the Λ/m t expansion we refer the reader to [27].

Comparison to [31]
Another approach, introduced in [31], aims at the computation of the non-resonant contribution to the total cross section in an expansion in ρ 1/2 , where ρ = 1 − m W /m t is treated as a small parameter. Even though ρ is not small in reality, one may hope that with sufficiently many terms in the expansion, a good approximation might be obtained.
Indeed, the exact NLO non-resonant result from [26] was reproduced by combining a deep expansion with Padé approximants.
Our concern here is the computation of the first two terms in the ρ 1/2 expansion of the NNLO non-resonant contribution. In the present notation, the first term in the expansion given in [31] reads explicitly σ (1),nr [31] 3) It is immediately clear from this expression that the meaning of "non-resonant" is different from ours, in which case the non-resonant contribution is analytic in energy and has a 1/ǫ pole. It appears that [31] does not distinguish between what we call non-resonant and absorptive matching coefficient contribution to the resonant part and directly constructs the expansion of the diagram in ρ 1/2 , such that (6.3) gives the sum of all contributions at order O(α s /ρ).
It is instructive to construct the O(α s /ρ) terms from the results in [30] and in the present paper. We find that they arise only from Abs,Z t (6.4) in part (I) and specifically from diagrams h 1a and h 1b in σ sq . Each of the two terms contains a 1/ǫ pole, which cancels in the sum. This holds separately for the two diagrams h 1a and h 1b plus their corresponding resonant counterparts 13 that contribute to σ C (k) Abs,Z t . We note that the leading term (6.3) from [31] originates only from diagram h 1a . Our result for this diagram including its resonant counterpart indeed agrees with the above except for the constant term +1 (see (C.4), (C.5) in Appendix C). However, as was already mentioned in [27,30], contrary to the statement made in [31] there is a nonvanishing contribution from h 1b at the same order. We computed the O(α s /ρ) from this diagram explicitly, and find that the complete O(α s /ρ) contribution to the total cross section reads Note the absence of a logarithmic dependence on ρ in the sum of all contributions (see (C.4)-(C.7) for the individual results). This can be traced to the cancellation of 1/ǫ divergences and the scaling of the leading momentum regions that contribute to the 1/ρ enhanced term. Furthermore, the coefficient of ln |E + iΓ| differs by a factor of two, which is related to the contribution of the diagram h 1b as described in Appendix C. We therefore disagree with the NNLO non-resonant result given in [31] already from the leading term in the ρ 1/2 expansion.
The authors of [31] did not actually attempt the calculation of diagram h 1b but referred to [54] to claim that it must not contribute. However, as already discussed in Section 3.1.2, the purported vanishing of h 1b , called "jet-jet" contribution in [54], refers to a different order in the non-relativistic expansion, namely NLO, and is reflected in the present framework as the non-renormalization of the coupling of the top quark to a potential gluon and the Coulomb potential by electroweak effects. Moreover, when the ρ 1/2 expansion is constructed from momentum regions, the leading 1/ρ term arises from a momentum region that was missed in [31]. Further details on the comparison and diagram h 1b can be found in Appendix C.

Discussion of results
Recent experimental studies [2,3] concluded that the statistical uncertainties of the top threshold scan at a future e + e − collider can be very small in realistic running scenarios. Thus, when discussing the impact of the electroweak and non-resonant corrections in this section, we focus on the theoretical uncertainties. An experimental analysis based on the theory prediction available in QQbar threshold [32] is in progress [4] and will combine statistical and systematic experimental errors with theory uncertainties.
To avoid the IR renormalon ambiguities, we exclusively employ the PS shift (PSS) mass scheme defined in [14,32,59]. For the numerical evaluation we adopt the input values where m PS t is the top-quark PS mass [76] and the running electroweak coupling is taken from [60], see the discussion in Section 3.1.4.

Size of the electroweak effects
We define a reference QCD prediction by adding the small P-wave contribution [22] to the S-wave result of [9]. The result is shown by the grey hatched band labelled "QCD" in Figure 10, which is spanned by variation of the renormalization scale µ r between 50 GeV and 350 GeV. Figure 10 also shows the net effect of all the corrections discussed above, excluding ISR, which will be considered below. These non-QCD effects slightly increase the height of the peak and move it towards smaller center-of-mass energies. Above the peak the cross section is slightly decreased by about 3.0 − 3.6%. Overall, the effect of the non-QCD corrections is to make the resonance more pronounced. The largest effect is observed below the peak, where the absorptive parts of the matching coefficients and the non-resonant contribution dominate the non-QCD corrections. Here, the bands cease to overlap at around √ s = 341.8 GeV. The size of the uncertainty band is somewhat increased and now reaches up to ±5.2% directly below the peak, where the uncertainty estimate for the QCD result is ±3.8%. In the remaining regions it is about ±3%. The increase of the scale uncertainty is mainly due to the Higgs potential insertion as was already observed in [25]. The size of the individual contributions is shown in Figure 11. We have already discussed the Higgs, QED Coulomb and NLO non-resonant corrections in [25], but briefly recapitulate the results here to give a complete overview over the non-QCD correction  up to NNLO. In the top-panel we show the relative effect of the Higgs contribution σ H at NNLO and NNNLO. At NNLO there is an almost constant relative shift, because only the hard-matching coefficient c (2) vH is present. At NNNLO, there is also a contribution from the local Higgs potential, which modifies the position of the peak. Due to the attractive nature of the potential, the binding energy is increased and the peak is shifted to the left. At the same time the Higgs corrections increase the cross section by 3 − 8%, depending on the value of √ s, and make the peak more pronounced. The comparison of the dashed and solid curves demonstrates that the inclusion of the NNNLO corrections is important for correctly capturing the energy dependence of the Higgs effects, which is crucial for a reliable measurement of the top Yukawa coupling. The remaining electroweak contributions to the 'partonic' resonant cross section, Abs,bare+Z t , are shown in the middle panel. The dashed line corresponds to the correction from the QED Coulomb potential only. It is attractive and therefore leads to an increase of the cross section by 2 − 8% and a shift of the peak towards smaller center-of-mass energy. The solid line shows the full correction. The width contribution σ Γ decreases the cross section by 0 − 1.5% depending on the energy. Including the real part of the electroweak matching coefficient leads to an almost constant relative shift of about −3.3%. The absorptive part of the matching coefficient multiplies the real part of the non-relativistic Green function, which has a broad peak, roughly centered around the point where the imaginary part has its maximal slope, on top of a smooth background. Thus, the absolute contribution has only a mild energy dependence and is of the order of −3% near and above the peak. However, it becomes even more important below the peak, where the cross section is small and modified by up to −15%.
The lower panel illustrates the behaviour of the non-resonant contribution to the total cross section. Its absolute size is nearly energy-independent. Thus, the shape of the curves is given by the "inverse" of the resonant cross section. At NLO, the effect is of the order −(3 − 4)% near and above the peak and reaches up to −22% for low center-ofmass energies, where the resonant cross section becomes small. The NNLO corrections compensate about 40% of the NLO result. This is in contrast to the findings of [27,30], where an enhancement of the negative non-resonant correction from an approximate NNLO result was observed. The apparent discrepancy is entirely explained by the very different choice made in [27,30] for the finite-width scale (µ w = 30 GeV) compared to the present (µ w = 350 GeV). The dependence of the full result on µ w is very mild as discussed below and, thus, mainly the size of the individual contributions is affectedmost notably the non-resonant correction and the one from the absorptive part of the hard matching coefficients.
We recall that the bands in Figure 10 only include the variation of the renormalization scale between 50 GeV and 350 GeV, while the scale µ w = 350 GeV is kept fixed. The dependence on the scale µ w cancels exactly between all contributions of a given order. We show the µ w dependence of the resonant cross section and the full cross section in Figure 12. For the resonant-only cross section, it is mild near and above the peak, but is significantly larger than the renormalization scale dependence below the peak. The sensitivity to µ w is greatly reduced for the full cross section, where the variation between  20 and 700 GeV considered in the plots only yields a ±(0.2 − 0.3)% effect near and above the peak and only a mild ±1.8% below the peak. The remaining µ w dependence is of NNNLO, where the full QCD corrections, but only a few electroweak effects are included and therefore no full cancellation is achieved. The central value µ w = 350 GeV for the finite-width scale is chosen near the hard scale to make the corresponding logarithms in the non-resonant part small. The logarithms of µ w are introduced by the separation into different momentum regions and are therefore spurious in nature. Explicitly, some of the 'large' logarithms ln v contained in the full cross section are split as follows where the first logarithm to the right of the equality sign is part of the non-resonant contribution and the second one of the resonant. Choosing µ w ∼ m t captures the 'large' logarithms present at NNNLO in the resonant part and renders the logarithms contained in the non-resonant part small. While the NNNLO resonant contributions are already partially known, the NNNLO non-resonant corrections are beyond the present computational limits. Thus our scale choice minimizes the uncertainty from the missing NNNLO contributions. 14 Variation of the scale µ w can be used to estimate the size of the missing NNNLO nonresonant corrections. The corresponding bands for the resonant plus NLO non-resonant and full cross section are shown in Figure 13, where we have varied µ w between 20 and 700 GeV. We observe that the inner (red) band is entirely contained in the outer (grey) one and much narrower. Thus, the chosen range of the finite-width scale variation provides a reasonable estimate of the NNLO non-resonant correction. However, an estimate of the missing NNNLO non-resonant correction based on the width of the "full" (red) band in Figure 13 is potentially less reliable, because the leading NNNLO terms might not cause any µ w dependence. This would be similar to the situation at NLO, where the leading non-resonant effect arises, yet there is no µ w dependence of the resonant contribution at this order at all, since the divergence from factorizing resonant and non-resonant contributions is purely linear.
We discussed the possibility of imposing loose cuts, which affect only the non-resonant part of the cross section, in Section 2.3. The dependence on the cut defined in (2.14) is shown in Figure 14, where the dotted and solid lines denote the NLO and NNLO nonresonant contribution. Very loose cuts with ∆M t ≥ 30 GeV have only a mild influence on the cross section. Tighter cuts ∆M t = (30, 20, 10, 5) GeV reduce the cross section by (0.007, 0.014, 0.037, 0.084) pb. We observe that for ∆M t around 4 GeV the NNLO nonresonant contribution becomes as large as the NLO one. Here, the assumption that the cut is loose is no longer appropriate and our description breaks down.  in Figure 14 shows the approximate NNLO result [27], which includes only the endpointdivergent terms as ∆M t → 0, for comparison. It describes the dependence on the cut very well, since the endpoint-divergent terms are most sensitive to it, but it is shifted by -0.004 pb for the full cross section and up to -0.013 pb including invariant mass cuts. In the absence of any cuts the exact result corresponds to a 46% correction with respect to the approximate NNLO result. We note, however, that for the scale choice of [27] and in the range of loose, but not too loose cuts Γ t ≪ ∆M t ≪ m t the approximation is much better.
We finally discuss the effects of initial-state QED radiation, which have so far only been taken into account in the experimental studies. Figure 15 shows the partonic cross section σ conv and its convolution with the electron structure functions. The QED contribution σ conv IS to the partonic cross section from (3.65) is a small effect of the order −(0.6 − 1.3)%. The convolution, however, reduces the cross section by 28 − 44%. The black band is spanned by four different implementations of the convolution (3.53) of the full NNNLO QCD plus NNLO EW cross section with the structure functions. This involves an extrapolation of the cross section for energy values outside of the range of the grids available in QQbar Threshold [32]. We either use the shape of the LO cross section below √ s = 328 GeV, rescaled to match the full result at √ s = 328 GeV, or an alternative implementation that interpolates linearly between σ( √ s = 320 GeV) = 0 pb and our result at √ s = 328 GeV. Numerically, we find a small difference of 0.1% near and above the peak, which goes up to 0.8% at √ s = 340 GeV. 15 For both extrapolations we consider the convolution (3.53) with the structure functions as defined in (3.60) and a purely LL approximation where we set β = (2α/π) ln(s/m 2 e ) in the structure function and accordingly modify the non-logarithmic ISR contribution (3.65) for the different subtraction term (3.63). The difference is formally a NLL effect and provides a rough estimate of the overall size of NLL ISR corrections. It amounts to about 1.4% above the peak and reaches up to 2.1% in the region where the slope is large.
For comparison we furthermore show as the dashed line the expression where the ISR resummation is only applied to the LO cross section. Since the LL resummation modifies a N k LO correction by order one, the difference between ISR and ISR 0 is formally a NLO effect. This emphasizes that it is mandatory to perform the convolution with the full partonic result. We see, as it is of course expected, that ISR is a huge effect, reducing the cross section by 28 − 44%. It also leads to a significant modification of the shape. The peak is shifted by almost 200 MeV to the right and smeared out considerably. Its height is reduced by about 40%. This emphasizes the need for a full NLL treatment of ISR and a proper analysis of the convergence and remaining uncertainty, which is of universal importance for high-energy e + e − collider processes, but beyond the scope of this work. We further note that at the level of NNLO electroweak accuracy the partonic cross section depends on the scheme employed for the electron structure function, and a phenomenological convolution as often applied in experimental studies in an unspecified scheme is no longer adequate.

Sensitivity to Standard Model parameters
Since the non-QCD effects computed in this paper cause substantial corrections to the cross section we provide an update of the discussion in [9,25] of the sensitivity of the top threshold scan to Standard Model parameters. Figures 16, 17 and 18 estimate the sensitivity by comparing the effects of parameter variations to the scale uncertainty in terms of the relative variation to a reference cross section.
All electroweak and non-resonant effects discussed in this paper are included in the figures, in particular also the ISR corrections. There are small quantitative differences with respect to [9,25], such as a small reduction of the height of the peaks present in the top-mass variation curves near 344.5 GeV, but the essence of the results and the associated conclusions remain unchanged. It is especially noteworthy that the huge ISR 15 While the grid could technically be extended to smaller values of √ s, the PNREFT and unstableparticle EFT breaks down far below the threshold. Improving the accuracy in this region would require matching the EFT description to the fixed-order calculation of the full non-resonant process as discussed for a single-particle resonance in [35,36] Figure 10). The prediction is normalized to the full cross section.
correction discussed above does not degrade the sensitivity. Since the bulk of the ISR correction is produced by the convolution with the luminosity function, we expect that the additional convolution of the cross section with the collider-specific beam function will not dilute the sensitivity to the parameters, either. From Figures 16 and 17 we expect the threshold scan to be sensitive to variations of about ±40 MeV for the top-quark PS mass, ±60 MeV for the top-quark width, +20 −25 % for the top-quark Yukawa coupling and ±0.0015 for the strong coupling constant α s (m Z ), when only a single parameter is varied at a time. These numbers are obtained from comparing the width of the band for the parameter variation with the one from the The predictions are normalized to the full cross section and the uncertainty band is the same as in Figure 16.
theoretical uncertainty, and requiring that the former is larger than the latter for a sufficient range in energy. This leaves open the question of how well the corrections from the simultaneous variation of several parameters can be disentangled from their energy dependence, which particularly concerns the Yukawa and the strong coupling, where variations lead to similar effects as seen in Figure 17 for the energy dependence and in Figure 18 for the position and height of the peak in the cross section. This needs to be addressed within realistic simulations, which include experimental uncertainties as well.

Conclusions
The recent advance in the QCD calculation of the top anti-top threshold [9] has motivated the consideration of non-QCD effects of potentially similar size to the third-order QCD correction. While Higgs/top-Yukawa coupling effects up to the third order were already obtained in [25], the present work completed the calculation of NNLO electroweak corrections and in particular the NNLO non-resonant contribution to the e + e − → bbW + W − X process near the top-pair production threshold. This elevates the theoretical prediction to NNNLO QCD plus NNLO electroweak accuracy, including for the first time initialstate radiation in a scheme consistent with one-loop QED corrections. The new effects are indeed non-negligible compared to the ±3% accuracy estimated for the pure QCD calculation and are therefore essential for accurate top and Standard Model parameter determinations from the threshold. They have been implemented in the new version 2 of the public code QQbar threshold [32]. Despite the level of sophistication already achieved, further improvement could be considered or might be necessary, such as the combination of the NNLL summation of logarithms of E/m t in the QCD part [24] with the NNNLO fixed-order calculation [9], the inclusion of already known NNNLO electroweak corrections (see [45,47]) or the oneloop correction to the Higgs potential (a N4LO effect) together with the terms required to make these additions factorization-scheme independent. To cancel the finite-width µ w scale dependence of the NNNLO QCD result completely, the non-resonant part is needed to the same accuracy, which appears prohibitive at present. Finally, a consistent implementation of QED initial-state radiation with next-to-leading logarithmic accuracy seems to be a general prerequisite for accurate predictions of scattering at a future highenergy e + e − collider.

]
As before, the complete nonresonant contribution can be disabled by setting the option resonant_only to true ( ResonantOnly −> True in Mathematica).

A.2 Initial-state radiation
Initial-state radiation requires the computationally expensive convolution with structure functions. Therefore, this correction is not included automatically.
After defining the luminosity function with the electron structure functions Γ ee from (3.53) the cross section after initial-state radiation is given by Here, σ conv is the partonic cross section including the non-logarithmic initial-state radiation correction σ conv IS (see (3.65)). The non-logarithmic correction can be included with the option setting In principle, the convolution integral in (A.2) covers the whole energy range from zero to the nominal center-of-mass energy. However, our prediction for the cross section is only valid in the vicinity of the threshold. Sufficiently below the threshold the actual cross section becomes negligible. This implies that we can introduce a lower cut-off x min in the integral (cf. Section 7.1). In the following we choose x min = (330 GeV) 2 /s. A further, purely numerical problem arises from the integrable divergence of the luminosity function for x → 1. In order to eliminate this divergence, we can change the integration variable to t = (1 − x) β and write the cross section as 3) with the modified luminosity functionL(t) and a cut-off t max = (1−x min ) β . The function β = −2α(µ α )/π[log(m 2 e /s)+1] is available in QQbar threshold as ISR_log(sqrt_s, alpha) in C++ and ISRLog in Mathematica.
Finally, version 2 of QQbar threshold provides an integrate function in the header integrate.hpp, which can be used to compute the convolution integral as shown below.
The following C++ code prints the cross section σ = 0.591736 pb after initial state radiation for a center-of-mass energy of Numerically, the structure function under the replacement β → 2β is very close to the luminosity function. The same holds for the modified versions of both functions. Indeed, substituting modified_luminosity_function(t, beta) with modified_structure_function(t, 2 * beta) in the example changes the result for the cross section to σ = 0.59169 pb, i.e. by less than 10 −4 . This observation can be used to somewhat accelerate the computation of the convolution at the cost of accuracy.

A.3 Width corrections
Among the resonant NNLO electroweak corrections listed in (3.12), only σ conv IS and σ Γ are not already available in version 1 of QQbar threshold. In version 2, the correction σ Γ proportional to the top-quark width is included by default in the prediction for the cross section. Its components (cf. (3.19), (3.22)) can be controlled individually with the new contributions options v_width_kinetic (Eq. (3.18)), v_width2 (Eq. (3.17)), and width_ep (Eq. (3.21)). The respective Mathematica contribution names are vwidthkinetic, vwidth2, and widthep.
To incorporate the width corrections to the quarkonium energy levels from (3.35) the function ttbar_energy_level(n, mu, {m, width}, order, opts) can now take both the mass and width as arguments. Similarly, the ttbar_residue function can now take both arguments. In this way, this function includes the width corrections to the wave functions from (3.37). The corresponding Mathematica functions TTbarEnergyLevel and TTbarResidue are similarly extended. Finally, the toponium width including the corrections in (3.36) can be computed with the new function ttbar_width(n, mu, {m, width}, order, opts) (TTbarWidth in Mathematica). 16

A.4 Note on backwards compatibility
Disabling the new corrections in version 2 via the contributions option will produce results that are similar, but not identical to version 1 of the code. There are two causes for the difference.
First, as detailed in Sections 3.1.3 and 4.1, the calculational scheme for the NNLO electroweak corrections to the resonant cross section has been changed. Consequently, predictions for the cross section at or beyond NNLO that include electroweak corrections will differ between the two versions, even if all new corrections are disabled. The numerical differences are typically less than 1%, but can amount to almost 10% for a small renormalization scale and far below the threshold, where the cross section is already very small and the (new and old) NNLO electroweak corrections are sizeable.
Second, in contrast to the original code, version 2 now captures the full dependence on the scale µ w . The numerical effect of this change is of the order of a few per mille for low energies and significantly less than one per mille in the peak region.

A.5 Calculation of the non-resonant correction
The dynamic numeric evaluation of the NNLO non-resonant corrections is computationally prohibitively expensive. Hence, QQbar threshold internally uses interpolation of a precomputed grid. For reference purposes, a copy NNLO nonresonant grid.tsv of this internal grid is provided in the directory given by the function grid_directory in C++ and the variable GridDirectory in Mathematica. The coordinates of the grids are given by x W = m 2 W /m 2 t , accounting for variations of the top-quark mass, and y w = (1 − y)/(1 − x W ), which covers changes in the invariant mass cut discussed in Section 2.3. The remaining two grid entries Σ automated (x W , y w ) and Σ manual (x W , y w ) parametrize the automated and the manual part of the non-resonant cross section for µ w = m t . To obtain their contribution to the cross section, these entries have to be multiplied by a factor of α s (µ r )σ 0 Γ t . The complete NNLO correction to the non-resonant cross section for arbitrary µ w is then given by where the coefficient of the logarithm reads Abs ) . (A.5) As mentioned before, the dependence on µ w has to cancel exactly against the dependence in the resonant cross section. Like in the resonant part, we therefore do not expand out the energy dependence of the s-channel propagators in Σ log . The last term in (A.4) is required because we have expressed the non-resonant cross section in terms of the allorder width Γ t . The NLO non-resonant part is proportional to Γ t and therefore implicitly width of the top quark itself as opposed to the width of a toponium bound state.
contains the NNLO correction (δΓ 1 /Γ 0 )σ NLO non-res , where δΓ 1 is the NLO QCD correction to the top-quark width. The same contribution appears in the NNLO calculation of the non-resonant part and we must include the last term in (A.4) to subtract this double counting.
Note that highly unphysical top-quark masses, i.e. x W < 0.15 or x W > 0.3 and extremely tight invariant mass cuts y w < 0.01 are not supported. Furthermore, the default values are assumed for the remaining Standard Model parameters, such as the values of m W and m Z .
Our code for producing the NNLO non-resonant grid depends on a number of software packages, including MadGraph5 aMC@NLO [41], FastJet [77], and Cuba [78]. In order to facilitate reproducing our results without having to install all dependencies, we provide an image that can be executed using the Docker virtualisation software. After installing Docker and downloading the image nnlo nonres grid entry.tar.gz from https://www.hepforge.org/downloads/qqbarthreshold/, it can be imported with docker load −i n n l o _ n o n r e s _ g r i d _ e n t r y . tar . gz and run with docker run −it amaier / n n l o _ n o n r e s _ g r i d _ e n t r y <xw > <yw > with xw and yw replaced by the respective values for the parameters introduced above. The last line of the output corresponds to a grid entry in the same format as in the reference grid.
For convenience we provide a Mathematica interface to the calculation of the NNLO non-resonant grid entries. After importing the Docker image as described above and loading the QQbarGridCalc Mathematica package distributed with QQbar threshold, a grid entry can be computed with Q Q b a r C a l c N N L O N o n r e s o n a n t G r i d E n t r y [ xw , yw , V e r b o s e −> True ] which returns a list {Σ automated , Σ manual }. Setting the Verbose option to False will suppress intermediate output.
Especially for large values y w ∼ 1 the calculation of the automated contribution can fail, in which case a slight change in the input parameters may help. In practice, this is not a severe problem as the automated contribution becomes essentially constant in this region. The precision of the automated calculation is not very high, and the values obtained can easily deviate from the ones in the reference grid by around 10%. Since the NNLO non-resonant contribution itself is not very large and typically dominated by the manual and logarithmic contributions, this translates to an error of at most one per mille in the final cross section. One way to reduce this error further would be to calculate the grid entries several times and average over the results.
In principle it is possible to compute an entirely new grid with the Docker container. In practice it is computationally much more efficient to calculate Σ automated in the absence of an invariant mass cut, i.e. for y w = 1, and derive the entries for all other values of y w exploiting complementary cuts as discussed in Section 2.3. The entry for some y w with 0 < y w < 1 is then given by Σ automated = Σ automated yw=1 − Σ automated , where the phase space integral in Σ automated is restricted to the complementary region 0 ≤ t ≤ 1 − (1 − x W )y w . In this way, the numerically problematic endpoint region t → 1 is only computed once for each value of x W .

B Implementation of the subtractions in MadGraph
We briefly describe the implementation of the subtractions of the diagram h 1 in Figure 1 and all the diagrams in Figures 2 and 3 in MadGraph5 aMC@NLO 2.5.0.beta2. The MG5 aMC code for the computation of the process e + e − → tbW − at NLO in QCD is created in the directory TBWsubtractions by entering the commands MG5_aMC > g e n e r a t e e + e− > t b~w− [ QCD ] MG5_aMC > output T B W s u b t r a c t i o n s in the MadGraph5 prompt. First, we subtract the diagram h 1 from the code in the directory ∼/SubProcesses/P0 epem wmtbx/. 17 To this end, we identify the corresponding diagram numbers in MadGraph as 3 and 4 using born.ps. The subtraction is achieved by removing the terms proportional to AMP(I) * DCONJG(AMP(J)) with I,J= 3, 4 in the squared matrix element. This affects the function BORN in born.f and BORN_HEL in born hel.f. We note that one should avoid first adding and then subtracting the terms to avoid numerical instabilities since these contributions are divergent at threshold.
To subtract the real corrections g i in Figure 2 we remove the respective terms in the squared real amplitude given by the function MATRIX_1 in matrix 1.f where the corresponding set of I,J values can be determined from matrix 1.ps. To maintain separate IR finiteness of the real and virtual contributions we also have to edit the FKS subtraction terms in the files b sf 001.f, b sf 002.f, and b sf 003.f accordingly. This is done by removing the terms containing the product of the tree-level amplitudes 3 and 4 in the functions B_SF_00i. In the folder ∼/SubProcesses/P0 epem wmtbx/V0 epem wmtbx for the virtual corrections we first apply the usual subtractions for the squared tree-level amplitude to the function MATRIX in born matrix.f. The interference of a given one-loop diagram with the tree-level diagrams is evaluated by the function CREATE_LOOP_COEFS in polynomial.f. We create a copy called CREATE_LOOP_COEFS_h1bcd and which is modified by removing the interference with the tree-level diagrams 3 and 4. This allows us to remove the diagrams h 1b , h 1c , h 1d and h ia with i = 1, . . . , 4 by modifying the calls to CREATE_LOOP_COEFS in coef construction 1.f for the loop diagrams corresponding to the left-hand sides of the cuts. We either add the suffix h1bcd or comment out the calls. Again the relevant diagram numbers in MadGraph can be identified from the graphical representation loop matrix.ps. The same changes are applied to the multiple precision version of the virtual corrections given in mp compute loop coefs.f and mp coef construction 1.f.
Last but not least one needs to modify the counterterms given in loop matrix.f. The identification of the Madgraph IDs of the counterterm diagrams is more complicated since they are not drawn but must be inferred from the code. The counterterm amplitudes AMPL(K,I), where the first index K= 1, 2, 3 denotes the finite part, the 1/ǫ pole and the 1/ǫ 2 pole, are defined in helas calls ampb 1.f and helas calls uvct 1.f. The first file contains R 2 -terms and mass renormalization counterterms which are attributed to the loop diagrams in the same order in which they appear in loop matrix.ps. The second file contains the multiplicative wave function renormalization counterterms for the tree-level diagram. To subtract the diagrams h 1e , h 1f , h 1g and the R 2 contributions of the remaining diagrams in the squared contribution we remove the terms proportional to AMPL(K,I) * DCONJG(AMP(J)) with I= 11, 12, 16-23, 28-31 and J= 3, 4 in the function SLOOPMATRIX in loop matrix.f.
In addition we can implement the minimal subtraction of the UV divergences in h ia with i = 2, . . . , 4 by modifying the interference of the divergent part of the wave function renormalization of the tree-level diagrams 3 and 4 with the other tree-level diagrams. Explicitly, we multiply AMPL(2,I) with I= 29, 31 in the subroutine HELAS_CALLS_UVCT_1 in helas calls uvct 1.f with a factor 2/3. As discussed in Section 5.1 the R 2 -terms and finite parts of the wave function renormalization contributions for the interference contribution are not modified. Alternatively, one can simply deactivate the check for UV finiteness which yields the same results.
Following the discussion in Section 5.1 we have to deactivate an internal MadGraph consistency check for the positivity of the squared real amplitude to make the modified code run without producing error messages. This is done by removing the code block if ( wgt . lt .0. d0 ) then ... endif in ∼/SubProcesses/fks singular.f. Older MadGraph versions also require modifications in the file ∼/SubProcesses/P0 epem wmtbx/BinothLHA.f to allow for negative values of the squared Born amplitude.
Our modified version of the MadGraph code is shipped with the grid generation routines in the new version of QQbar Threshold. We have checked our procedure by applying similar modifications to the process e + e − →tbW + generated with the older MadGraph version 2.4.3. Furthermore, we have verified that an analogous set of changes correctly removes the Z-boson exchange contribution to the process e + e − → tt at NLO by comparing the results to the ones obtained by excluding the Z-boson exchange already in the process generation.
C Further details on the comparison to [31] We extend in this appendix the discussion about the discrepancy with the result for the cross section at leading order in the ρ 1/2 expansion of [31] and its connection to diagram h 1b . First we explain why the cancellation of finite-width and endpoint divergences requires a non-vanishing contribution from diagram h 1b at leading order. A similar argument was already put forward in [27,30]. Then we show that the leading-order term in h 1b comes from a loop-momentum region which is not among those considered to construct the unstable top EFT formulated in [31]. The NNLO resonant corrections related to h 1a and h 1b are given by the diagrams displayed in Fig. 19, where the loop-momenta in the top anti-top loops are expanded according to the potential (p) scaling, while the loop momentum in the bW loop is hard (h). It should be understood that in the diagram with the self-energy insertion in the propagator one has to consider only the NNLO piece (the cut self-energy in the on-shell limit gives the top-quark width, which is a LO contribution since p 2 − m 2 t ∼ Γ in the potential region; such terms are already accounted for by the replacement E → E + iΓ in the non-relativistic propagator). The diagram with the cut self-energy contributes to σ Γ (3.22) and σ C (k) Abs,Z t (3.44), once the symmetric diagrams are considered. In particular, the latter arises from field renormalization due the absorptive part of the electroweak one-loop self-energy. We have isolated that contribution in the second line of Figure 19, and split it such that one half of it can be attributed to field renormalization of the top quark leaving the production vertex, and the other half to the renormalization of the top-quark field entering the ttg vertex. The latter contribution is exactly cancelled by the electroweak correction to the Coulomb potential (third diagram in the second line of Figure 19), because upon expanding out the external (potential) momenta from the self-energy and vertex loops, these diagrams are equivalent to the renormalized vertex in the on-shell scheme for zero transfered momentum (see Figure 5). Therefore, the resonant counterpart of h 1b is equal to minus one half of the diagram with the field Figure 20: Forward-scattering diagram whose imaginary part is related to cut diagrams h 1b .
renormalization of the top quark leaving the production vertex, which is proportional to the coefficient C (k) Abs,Zt written in (3.44). It can be easily checked that C (k) Abs,Zt behaves as Γ 0 /(1 − x W ) ≃ Γ 0 /ρ, and that the contribution to the cross section σ C (k) Abs,Z t contains a α s /ǫ divergence from the real part of the Green function. Therefore both resonant diagrams in the first line of Figure 19 contain α s /ǫ × Γ 0 /ρ divergences, that are cancelled with endpoint divergences from the corresponding non-resonant diagrams h 1a and h 1b as was shown by explicit computation in [30]. In [31] only the non-resonant diagram analogue to h 1a is considered, while it is argued that the analogue to h 1b must vanish quoting results from [54]. As already explained in Section 6.3, the results from the latter refer to the vanishing of resonant contributions related to the top-quark instability at NLO, while the contribution under discussion here is of NNLO.
The dominant terms in the ρ expansion can also be obtained upon application of the method of regions [37,38]. Let us consider the forward scattering diagram in Figure 20, whose imaginary part corresponds to the sum of cut diagrams h 1b and g 5 (plus left-right symmetric ones), since no other cuts are kinematically possible. The authors of [31] discuss the contributions from the regions that are obtained by replacing v → ρ 1/2 in the hard, soft, potential and ultrasoft region. However, for the case of the diagram in Figure 20, which is related to h 1b discussed above, it can be shown that the leading order contribution comes from an additional region with parametrically smaller virtuality or order ρ 2 m 2 t , that was not considered in [31]. With the momentum assignment of Figure 20, it corresponds to k 0 i ∼ m t ρ 2 , k i ∼ m t ρ and l ∼ m t ρ. We obtain for the relevant scalar integrals in this region: = (−1) η π(2ρ) 1−2η−6ǫ 2m 2+2η+6ǫ t e 6iπǫ Γ 1 2 + ǫ Γ 1 2 − ǫ Γ (η + 2ǫ) 2 Γ (1 − η − 2ǫ) Γ (2 − 2ǫ) Γ (2η + 4ǫ)