P-wave contribution to third-order top-quark pair production near threshold

The next-to-leading order (NLO) P-wave Coulomb Green function contributes at third-order to top-pair production in e+ e- collisions near threshold. In this paper we compute the NLO P-wave Green function in dimensional regularization, as required for a consistent combination with non-resonant production of the W+ W- b\bar{b} final state, and present a phenomenological analysis of the P-wave contribution. We further briefly discuss squark production near threshold and top-pair production in gamma gamma collisions, where no S-wave contribution is present, and the P-wave thus constitutes the dominant production process.


Introduction
A future high-energy electron-positron collider will allow a very precise measurement of the top-antitop production cross section near threshold. From the threshold scan several standard model parameters, like top-quark mass, width, and Yukawa coupling, can be extracted with high precision. It was found in several studies [1][2][3] that the top mass can be determined with an uncertainty well below 100 MeV. Contrary to direct reconstructions at hadron colliders, there is no ambiguity in relating the result to a precisely defined mass parameter. To achieve this level of accuracy requires precise theoretical predictions for the threshold cross section. The challenge is that conventional perturbation expansions in the strong coupling α s fail for threshold production, since it involves multiple scales. In terms of the mass m t and velocity v of the top-quark these are the hard scale m t , the soft scale m t v, and the ultrasoft scale m t v 2 . At the ultrasoft scale, the colour-Coulomb force is non-perturbatively strong. In terms of Feynman diagrams, this implies that when the velocity is of order of the strong coupling, Coulomb singularities of the form (α s /v) n have to be summed to all orders. This can be achieved by successively integrating out the hard and soft scale leading to the non-relativistic QCD (NRQCD) [4][5][6] and potential NRQCD (PNRQCD) [7][8][9][10][11] effective field theories, respectively. Within this framework, described in detail for top-quark pair production near threshold in [12], the dominant contribution from the S-wave correlation function has been computed at next-to-next-to-next-to-leading order (NNNLO) [13][14][15][16].
The axial-vector coupling of the top quark to the Z boson gives rise to a P-wave contribution to the top-pair production cross section. In this work we compute the corresponding P-wave PNRQCD correlation function. Being suppressed by v 2 relative to the S-wave it contributes only starting from NNLO. The complete NNNLO calculation of the threshold correction therefore requires a NLO calculation of the P-wave correlation function, which we perform here in dimensional regularization. Some results for the P-wave Green function were already obtained in [17][18][19][20], but none of these computations were performed in dimensional regularization. Dimensional regularization is, however, required for the following reason: The imaginary part of the P-wave Green function, which is relevant for the cross section, is divergent already at leading order in the nonrelativistic expansion, if the finite width of the top quark is included. This divergence and the resulting scheme dependence cancel only when non-resonant corrections to the process e + e − → W + W − bb are added. The separation of resonant and non-resonant contributions can be performed consistently in unstable-particle effective field theory [21,22], in which the non-resonant terms appear as a hard region. The corresponding diagrams are computed as usual in dimensional regularization as has already been done in [23,24]. Consistency then requires that the non-relativistic, resonant part is also computed in dimensional regularization. We emphasize that to determine the top-quark mass precisely it is necessary to compute the process e + e − → W + W − bb including non-resonant terms, since the top mass is ultimately determined from the rise of the cross section near threshold, where non-resonant effects are important [23].
The outline of the paper is as follows. In Sec. 2 we give a brief overview of the framework of the calculation. A detailed discussion can be found in [12,16], however, we repeat some of the formulas given there in order to make the present paper selfcontained. The P-wave Green function up to NLO is computed in Sec. 3 with some technical details relegated to the appendix. In Sec. 4 a numerical analysis of the P-wave contribution to top-pair production at threshold is presented. The absolute size of the P-wave contribution is rather small due to the small axial couplings of the top quark. We therefore also discuss the P-wave dominated (s)top threshold production processes γγ → tt with different photon helicities and e + e − →tt in Sec. 5 and Sec. 6, respectively, where P-wave production is the dominant production mechanism. We conclude in Sec. 7.

Effective theory setup
The production of a top pair in e + e − annihilation is mediated by photons and Z bosons. While the coupling of photons to fermions is purely vector-like the Z boson couples to vector currents and axial-vector currents with the respective strengths where θ w is the Weinberg angle, e f the electric charge and T f 3 the third component of the weak isospin of the fermion f . For the vector current j (v) µ =tγ µ t and the axial-vector current j (a) µ =tγ µ γ 5 t we define the two-point functions We denote by R = σ Xtt /σ 0 the inclusivett production cross section σ Xtt = σ(e + e − → ttX) normalized to the high-energy limit of the µ + µ − production cross section σ 0 = 4πα 2 em /(3q 2 ). It can be related to the imaginary part of the two-point functions by the optical theorem ) .

(2.3)
Near the production threshold s ≡ q 2 ≈ 4m 2 t the usual perturbation theory in α s breaks down, because the velocity v of the top quark is of the same order as the strong coupling and contributions that scale as (α s /v) k have to be summed to all orders. Instead, in non-relativistic perturbation theory, one expands in α s and v around the nonperturbative solution that sums these terms. Explicitly, the re-organized expansion takes the form The computation of the vector current spectral function Π (v) (q 2 ) to NNNLO is summarized in [12,16]. This work focuses on the axial-vector current spectral function Π (a) (q 2 ). In a first step the hard modes are integrated out. Matching of the axial-vector current to NRQCD yields an expansion where bold face letters and latin indices refer to d −1 dimensional vectors, and d = 4 −2ǫ is the space-time dimension. The hard matching coefficient c a is given by [25] 1 We observe that (2.5) contains a covariant derivative. Thus the axial-vector current spectral function is suppressed by the well-known P-wave factor of v 2 compared to the S-wave and contributes to R only starting from NNLO. The two-point function takes the form A further matching to potential NRQCD integrates out the soft modes and the potential gluons and light quarks. To the required order the PNRQCD Lagrangian reads where ψ (χ) denotes the potential quark (antiquark) field. The coupling to the ultrasoft gluon field A 0 (t, 0) can be removed by a field redefinition, which, in general, modifies the external current [26], but cancels for the colour-singlet currents relevant to production through photons and Z bosons. Henceforth, ultrasoft gluons can be ignored, since they contribute only at higher order. Since the top pair is produced in a colour-singlet state, only the colour-singlet projection of the general potential is required. To the considered order it consists purely of the Coulomb potential where the d dimensional coefficients of the LO and NLO Coulomb potential are They are required up to order ǫ for reasons that will become clear later. The LO Coulomb potential contributes to the same order as the leading kinetic terms in the PNRQCD Lagrangian and thus has to be treated non-perturbatively, while the NLO correction V (1) C is a perturbation. Thus the LO Lagrangian describes the propagation of quark-antiquark pairs, where ladder diagrams with exchange of an arbitrary number of potential gluons between the quark-antiquark pair have been resummed. The quarkantiquark pair propagatorG 0 (p, p ′ ; E) satisfies the d-dimensional Lippmann-Schwinger equation where E = √ s − 2m t denotes the energy relative to the production threshold, m t the top-quark pole mass, and the scaleμ = µ [e γ E /(4π)] 1/2 has been chosen such that the minimal subtraction of 1/ǫ poles corresponds to the MS rather than the MS scheme.
Here and in the following the tilde is used to indicate that the Green function is given in momentum space. Its Fourier transform is the solution to the Schrödinger equation in four dimensions. An expression for general d is not available, but for d = 4 several representations are known [27][28][29][30]. We find it convenient to use the integral representation from [27], (2.17) with r ′ < r and p = √ −m t E, and P l (z) the Legendre polynomials. We also introduced the variable The Feynman rules required for higher-order computations have been derived in [12]. As also discussed there the soft matching of the two-point function is trivial and hence where now is the P-wave Green function at the origin. In passing to the second equation, we used that due to the spin-independence of the Coulomb potential (2.10) the (d−1)-dimensional spin algebra in (2.20) can be evaluated once and for all. In PNRQCD perturbation theory up to NLO, the P-wave correlation function reads where we have inserted the NLO Coulomb potential. As yet we have neglected the sizable width of the top quark Γ t = 1.33 GeV. For the leading-order S-wave contribution the width effect is accounted for by the replacement E → E +iΓ t in the spectral function [31,32]. We follow this prescription to define the pure QCD calculation of the pair-production cross section and therefore assume henceforth that E can take complex values. Due to the p · p ′ factor already the leading order result for the P-wave contains ultraviolet divergences of the form E/ǫ. Following the above description yields poles proportional to Γ t /ǫ in the imaginary part, and a scale-dependence related to these poles. To make this explicit, we distinguish between the scale µ r at which we evaluate the running coupling α s = α s (µ r ) and the scale µ w that arises from the finite-width divergences. While the residual µ r scale-dependence must always be of higher-order, the divergences and the dependence on the scale µ w have to cancel against non-resonant electroweak corrections involving a W b-loop correction to the off-shell top propagator [15]. So far only the 1/ǫ pole corresponding to the LO P-wave correlation function is known [24,33,34]. The finite term is required to cancel the corresponding scheme dependence. Since the finite term as well as the NNNLO non-resonant terms related to the NLO P-wave correlation function are presently unknown, we keep the dependence on µ w (and the associated poles) explicit in our analytical result.
3 Computation of the P-wave Green function

Leading order
Simple power counting shows that ladder diagrams with up to four loops are ultraviolet divergent, as compared to two loops for the S-wave, due to the additional factor p · p ′ in (2.21). These diagrams, the sum of which is denoted by G P (≤3ex) 0 , therefore have to be computed in d dimensions. We have used FIRE [35,36] to perform the reduction to a small set of master integrals. Results for these master integrals are available in the literature [37] and are in agreement with special cases of the more general calculation for the NLO contribution presented below. We obtain where we have given the contribution from the one-loop, zero-gluon exchange diagram explicitly. It exhibits the characteristic 1/λ 3 ∼ E 3/2 ∼ v 3 threshold behaviour of P-wave production. The higher-loop integrals I (n−1)0 P [1], corresponding to diagrams with n gluon exchanges, are given in App. A. (The notation is explained in more detail in the context of the NLO calculation after (3.16).) The remaining part G P (≥4ex) 0 is finite and can be calculated in d = 4 dimensions. We perform this calculation in position space. Eq.
where . . . denotes the angular mean of the respective expression. This projects out the P-wave, l = 1, component of the Green function G 0 (r, r ′ ; E). The expression (3.2) can be computed using the representation (2.17) with appropriate subtractions for the parts with less than three gluon exchanges. The sum G gives the expression for the correlation function in dimensional regularization: with a ∈ {r, w}, γ E the Euler-Mascheroni constant, and ψ(z) the logarithmic derivative of the Gamma function. In (3.3) subtracting the 1/ǫ poles gives the result in the MS scheme. Note that due to the overall 1/λ 3 factor only the term proportional to λ in square brackets results in a finite-width divergence. We have checked that this divergence in the imaginary part of the LO Green function (3.3) agrees with the result given in [24,38]. Neglecting the width of the top, the imaginary part is finite, and reads where v ≡ E/m t , and E n = −(m t C 2 F α 2 s )/(4n 2 ) with n ≥ 2 are the l = 1 bound state energies. Eq. (3.5) agrees with [17].

Next-to-leading order
Analogous to the S-wave computation [16], we define the single-insertion function where q ij = p i −p j . In terms of this the NLO correction to the Green function contained in (2.21) is given by As for the LO Green function divergences only occur in diagrams with up to four loops. We therefore split the NLO correction into a divergent (a) and a finite part (b) as indicated in Fig. 1. Since the top quark width Γ t cannot be neglected, the imaginary part of (a) contains divergences of the type Γ t /ǫ arising from poles of the form E/ǫ. Thus, contrary to the S-wave case, where no such divergences are present in the computation of corrections from the Coulomb potential, the NLO Coulomb potential (2.12) cannot be expanded in ǫ prior to integration in the computation of the loop integrals of part (a). However, the momentum integrals of part (b) are finite and the potential can be expanded before integration here. To deal with part (b), we start with the computation of the complete I P [1 + u] in position space and later perform the necessary subtractions to obtain I (b) P [1 + u]: We proceed by first solving the integral over r: (3.9) Taking the derivatives and the limit x, y → 0, we can perform the integrations over t 1 and t 2 and obtain (3.10) After the necessary subtraction the remaining part (b) is finite and can be computed in d = 4 dimensions. We obtain Here the last three terms are the first three terms in the expansion in α s , which corresponds to up to two gluon exchanges to the left or right of the NLO Coulomb potential insertion. This is precisely part (a), hence the above, subtracted expression is part (b) as desired. A strategy to solve this kind of integral is presented in [16]. We obtain 4 F 3 (1, 1, 4, 4; 5, 5, 3 − λ; 1). (3.14) Here ψ n (z) is the nth derivative of the psi function. We provide some useful formulas for the evaluation of the generalized hypergeometric function 4 F 3 in App. B. In terms of j P (u) and j ′ P (u) from (3.13) and (3.14) we can write part (b) of (3.7) as where we have expanded the NLO Coulomb potential (2.12) and written the logarithm of q 2 as a derivative at zero u.
The divergent part (a) is given by where the I nm P [1 + u] denotes the contribution to the single insertion function from the diagram with n potential gluon exchanges to the left and m to the right of the potential insertion. Part (a) then takes the form of (3.7) with I P replaced by (3.16), and requires the calculation of some four-loop diagrams in dimensional regularization. The results for the I nm P [1 + u] needed are given in App. A. The complete NLO correction to the Green function in dimensional regularization given by the sum of parts (a) and (b) reads: The dependence on the two scales µ r and µ w has been obtained with the procedure described in App. A. Alternatively, the dependence on µ r can be obtained using oneloop running of α s in the LO Green function. The remaining logarithms of µ must then be assigned to µ w . We note that the dependence on the scale µ w is polynomial in E and cancels in the imaginary part for Γ t = 0, which provides a useful consistency check.

Pole resummation
At negative energies the P-wave Green function contains poles corresponding to quarkantiquark bound states with angular momentum l = 1. In (3.17) they appear as poles of polygamma and hypergeometric functions for positive integer λ ≥ 2. Near these bound states the exact Green function has the form where E P n is the energy of the nth P-wave bound state and ψ ′ n (0) the derivative of the corresponding wave function at the origin. Both take on a perturbative expansion in the strong coupling constant 19) The NLO Green function expanded in α s therefore takes the form The singular terms near E P (0) n can be resummed into a single pole to all orders by subtracting (3.20) from the Green function and adding (3.18) with the energies and derivatives of the wave function in NLO approximation [10]. The leading-order expressions can be read off from the imaginary part (3.5) of the LO result: To obtain the NLO corrections we expand (3.17) for λ near positive integer n. We find (n + 2) − 2nψ 1 (n + 2) + 2a 1 + n 2 − 1 n(n − λ) 2 2β 0 (L r n +ψ(n + 2)) + a 1 + regular . With this it is straightforward to obtain: e P 1 =2a 1 + 4β 0 L r n +ψ(n + 2) , (3.23) where L r n = ln (nµ r /(m t C F α s )) . We have checked that this agrees with the results of [18,19].

P-wave top-pair production cross section
In this section we discuss the phenomenological aspects of the P-wave contribution to the top-pair production cross section near threshold. All expressions so far have employed the pole mass definition of the top quark. Since the pole mass suffers from an infrared renormalon ambiguity [39][40][41], which is not present in the top-pair cross section itself, we show results using the PS mass definition [42], which eliminates this spurious infrared sensitivity. This has been implemented in the PS Shift (PSS) and PS Insertion (PSI) schemes [16]. Denoting by δm t the difference between the pole mass and the PS mass, the former is defined by where the value of δm t is order-dependent. That is, in LO, we only use the leading-order expression δm t ∝ µ f α s , whereas in NLO, δm t includes the µ f α 2 s term that contains the a 1 correction to the Coulomb potential [42]. 2 The PSI scheme is obtained by reexpanding the right-hand side in α s . We find, however, that the difference between the two schemes is very small and would not be visible in the figures below. We therefore only show the results in the PSS scheme. For the numerics we adopt the parameter values α s (M Z ) = 0.1184, m PS t ≡ m PS t (µ f = 20 GeV) = 171 GeV, which corresponds to a LO (NLO) pole mass of m t = 172.025 (172.433) GeV, the top-quark width Γ t = 1.33 GeV and the Weinberg angle sin 2 θ w = 0.23. The strong coupling is evolved to the scale µ r in the four-loop approximation.
Since G P (E) is divergent there exists an ambiguity in (2.19) whether G P (E) or Π (a) (q 2 ) should be minimally subtracted, since the factor that relates both depends on d. This ambiguity is only resolved when the finite term of the non-resonant contribution is computed. Since the non-resonant calculation does not refer to any kind of non-relativistic approximation it is natural to define the non-resonant and resonant contribution by minimal subtraction of Π (a) (q 2 ). This will be assumed in the following. To be precise, we write the expansion of the hard coefficient (2.6) in the form Then the expansion of Π (a) (q 2 ) to NLO in minimal subtraction is given by Π (a) (q 2 ) = Π (a) where the constant terms proportional to Γ t /m t arise from the divergent parts of G P (E) multiplying the order ǫ and ǫ 2 terms of the factor (d−2)/(d−1)×c 2 a , which relates G P (E) to Π (a) (q 2 ). 3 Note that we assume here that Γ t takes its numerical, four-dimensional, physical value, while for an analytic combination with the (yet unknown) non-resonant cross section, one eventually needs to use the analytic, d dimensional, leading-order expression for the top-quark width. The constant terms in (4.3) and (4.4) shift the cross section only by a tiny amount, so that the scheme-dependence related to the resonantnonresonant separation that can be resolved only once the non-resonant cross section is fully known, is not relevant for the following discussion.
In Fig. 2 we show the LO and NLO P-wave contributions to the R ratio in the PSS scheme, employing the expressions (4.3), (4.4), and including pole resummation up to the n = 6 bound-state pole. The overall size relative to the S-wave is below 1% in the threshold region, because in addition to the v 2 suppression the ratio of the couplings (v 2 e + a A comparison of the dependence on the two scales is shown in Fig. 3 for three values of the energy (above, at and below threshold from top to bottom). We again see that the dependence on the renormalization scale µ r is strongly reduced at NLO. This is not the case for the dependence on µ w , which remains almost the same as at LO. This is expected, since the finite-width scale dependence does not cancel by performing higherorder QCD calculations. Rather, it has to cancel only when the non-resonant corrections are added to the result. The fact that the finite-width scale dependence is dominant at NLO shows the importance of this cancellation. In particular, the lower-right plot shows that the finite-width scale dependence changes the cross section by a large factor below threshold, precisely where the non-resonant contributions are expected to be important.

Top-pair production in photon collisions
The photon collider option via the back-scattering method [43] was studied in the Technical Design Report for TESLA [44] and is also considered at the ILC (see Sec. 12.6 of [45]). We discuss here only γγ → tt collisions with opposite photon helicities, where the top pair is produced in a P-wave state. We define the inclusive normalized cross section R +− γ = σ +− γγ→ttX /σ 0 . Near the production threshold the cross section takes the  Figure 3: The scale dependence of the P-wave contribution to the R ratio. The left and right plots show the dependence on µ r for fixed µ w and vice versa. The dependence is shown above (upper plots), close to (middle plots), and below threshold (lower plots) corresponding to √ s = 347.733, 344.866, 342 GeV, respectively. In the left plots the finite width scale is fixed to µ w = m PS t . In the right plots the coupling-renormalization scale is fixed to µ r = 80 GeV.
form [18,19] with the hard matching coefficient In the absence of a d dimensional calculation of C +− h we define here the R ratio, to which non-resonant contributions should eventually be added, by minimal subtraction of the We show results within the PSS scheme in Fig. 4. The P-wave cross section is an order of magnitude larger than in e + e − collisions, due the different size of the electroweak couplings and numerical prefactors. It can be observed independently of the S-wave by adjusting the beam polarizations. Since the γγ induced cross section differs from the e + e − one only by the hard-matching coefficient and overall electroweak couplings, we observe essentially the same features as in the previous section. Most importantly, the theoretical uncertainty of the QCD contributions as measured by the residual µ r dependence is greatly reduced at NLO as seen from the (hardly visible) width of the dark-shaded band in Fig. 4.
We have compared our result for the Green function with a previous result from [18,19] and found agreement of the scheme-independent terms at LO. The scheme-dependent finite-width terms have been fixed in these papers by matching the non-relativistic, resonant computation to the full theory diagram with an off-shell top-quark self-energy. This procedure accounts for part of the non-resonant contributions, but does not eliminate the need for a complete calculation of the γγ → W + W − bb process with opposite photon helicities in the vicinity of √ s ≈ 2m t to achieve parametric LO accuracy. 5 Refs. [18,19] also present a calculation of the NLO P-wave Green function, but contrary to the closed expression (3.17) in dimensional regularization, the result is given in a sum representation that makes the comparison of even the scheme-independent terms difficult. We find, however, that we are able to reproduce the plots in [18,19] to a good approximation, if we choose the scale of the hard matching coefficient and µ w at the hard scale.
6 Stop-pair production at e + e − and hadron colliders The production of pairs of the scalar supersymmetric partner of the top quark is of interest at both hadron and e + e − colliders. The threshold cross section of stop-antistop pairs in e + e − collisions can be described in a fashion analogous to the previous sections [17].
Since the coupling of squarks to photons and Z bosons contains a derivative, stops are produced in a P wave. We focus here solely on production of the lighter mass eigenstatẽ t ≡t 1 =t L cos θt +t R sin θt. The ratio Rtt = σtt X /σ 0 is given by where zt is the coupling constant for the Ztt vertex, which depends on the mixing angle θt, and Π (∂) (q 2 ) is the two-point function of the derivative current, which matches to the NRQCD current Here ψ denotes the stop and χ the anti-stop field in the non-relativistic normalization ψ ∼ χ ∼ m 3/2 t . Including the hard matching coefficient, the two-point function takes the form where is the hard matching coefficient of the current j (∂)k [47] and the P-wave Green function. It is straightforward to evaluate this expression. Since (6.3) differs from the corresponding expression for top-antitop production only by the prefactor, we observe the same qualitative features. Most importantly, the theoretical uncertainty of the QCD contributions as measured by the residual µ r dependence is greatly reduced at NLO as compared to the LO calculation [17].
In the following we comment on the relevance of the P-wave for the production of stopantistop pairs at hadron colliders. In quark-antiquark annihilation the t-channel gluino exchange diagram that is dominant for pair production of light-flavour squarks is strongly suppressed for stops due to the negligible top parton distribution function in the proton. Thus, in quark-antiquark annihilation the s-channel, which produces stop pairs in a Pwave and colour-octet state, is the dominant contribution. Production of stop-antistop pairs in gluon-fusion can be described as for light-flavour squarks. In [26] a formalism for resummation of soft and Coulomb corrections in S-wave production of pairs of heavy coloured particles was derived at the next-to-next-to-leading logarithmic approximation (NNLL). This was generalized to stop pair production at NLL in [48,49]. The arguments presented there suggest that the factorization formula also holds at NNLL. The NLO Pwave Green function derived in this work constitutes the dominant part of the potential function J Rα accounting for NNLL terms beginning with α 2 s /β. To obtain the colouroctet Green function one only has to make the replacement −C F → D 8 = 1/(2N c ) in (2.12). Additionally, non-Coulomb potentials yield terms beginning with α 2 s log β [50], which also have to be included at NNLL order.

Conclusion
We have computed the P-wave Green function in dimensional regularization up to NLO. We further confirmed results for the NLO correction to energy levels and wave functions at the origin of P-wave bound states. The NLO correction reduces the renormalization scale uncertainty considerably. We have discussed the P-wave contribution to three different pair production processes. The NLO P-wave contribution to the top-quark pair production cross section near threshold is part of the complete NNNLO result. The P-wave production cross section turns out to be small relative to the dominant S-wave, below 1%. It is included in the forth-coming NNNLO result for the e + e − → ttX cross section [51]. The photon collider option further offers the possibility to produce tops in a pure P-wave and with a larger cross section.
In e + e − collisions squark-antisquark pairs are also produced in a P-wave. We have given the necessary formulas for the NLO cross section. If squarks that are sufficiently light for production at a future linear collider should be found, this will allow precision studies including a precise mass determination. The NLO P-wave Green function is also an important ingredient in the NNLL prediction of stop-antistop production in hadron collisions.

Acknowledgement
This work has been supported by the DFG Sonderforschungsbereich/Transregio 9 "Computergestützte Theoretische Teilchenphysik", the Gottfried Wilhelm Leibniz programme of the Deutsche Forschungsgemeinschaft (DFG), and the DFG cluster of excellence "Origin and Structure of the Universe."

A Computation of the diagrams
We compute the divergent part of the single-insertion function (3.16) in momentum space. The integrals contain only a single scale mE, which appears as a mass term in the non-relativistic heavy quark propagators after performing the integrations over the zero-component of the loop momenta. We rescale the integration momenta by this scale, p → √ −mE k, to make them dimensionless. We then use FIRE [35,36] to reduce the diagrams to a set of master integrals. In the master integrals solid lines denote rescaled massive propagators, which take the form 1/(k 2 + 1). Dashed lines denote potential gluons and wavy lines insertions of the NLO Coulomb potential, i.e. gluon propagators with an index 1 + u. Initially, we make the assignment µ r to all µ raised to powers of u and µ w to all µ raised to powers of ǫ, since the former arise from the running coupling contribution to the NLO Coulomb potential. There are however some subtleties associated with this scale separation, which will be discussed below. We obtain , (A.1) × (u + 4ǫ − 1) u 2u 2 + u − 1 + 8(4u + 1)ǫ 2 + 2(u(7u + 3) − 2)ǫ + 24ǫ 3 + 2(u + 3ǫ) 2u 2 + 6uǫ + u + 4ǫ(ǫ + 1) − 1 (A.4) The master diagrams with just two massive lines can be computed with standard methods. We find (A.8) The remaining master integral contains three massive lines and is therefore more complicated: (A.9) The solution for u = 0 can be found in [37] and agrees with our result obtained by using FIRE and the known master integrals. For u = ǫ we calculate the diagram as an expansion in ǫ up to order ǫ 2 . The quadratic term is required, because the coefficient of this integral in (A.3) and the potential each contain a factor ǫ −1 . We use the MB.m package [52] to perform an analytic continuation of the Mellin-Barnes integral in ǫ. We then close the integration contours to pick up single and double infinite sums over the residues of the integrand. These sums can be transformed into cyclotomic harmonic sums. They were reduced to a set of known basis sums using the Harmonic Sums package [53][54][55][56][57][58][59]. The result is given by This is due to the fact that the origin of poles cannot be unambiguously identified in dimensionally regulated multi-loop integrals. To obtain the correct scale assignment, we subtract from the results of (A.1)-(A.4) the scale dependent logarithms and then add the respective terms obtained without the identification u = ǫ, i.e. with the Coulomb potential expanded in ǫ. Furthermore, the scales µ r and µ w are set equal in pole terms ln(µ r /µ w )/ǫ to ensure that the pole terms are scale-independent.
The results for the individual diagrams are: (A.14)

B Evaluation of the hypergeometric function
The generalized hypergeometric function in the result for NLO Green function can be expressed in terms of harmonic sums. This is useful when our result is applied to particles with vanishing width, which means that λ has a large positive real part as one approaches the threshold from below. The necessary analytic continuation can be easily done for the harmonic sums, see for example [57,60]. We first use 4 F 3 (1, 1, 4, 4; 5, 5, 3 − λ; 1) = 4(λ − 2)(λ − 1)λ(λ + 1)