Sommerfeld-corrected relic abundance of wino dark matter with NLO electroweak potentials

Extending previous work, we calculate the electroweak potentials for all co-annihilation channels of wino dark matter at the one-loop order and obtain the wino relic abundance including the Sommerfeld effect at the next-to-leading order (NLO).


Introduction
A weakly interacting massive particle (WIMP) is one of the best motivated dark matter (DM) candidates. Despite their simplicity, WIMP extensions of the standard model (SM) exhibit a rich and interesting phenomenology. Hisano et al. [1][2][3] recognized that despite the fundamentally weak coupling, the annihilation cross section of DM particles with mass m χ above a TeV is substantially enhanced by attractive forces that become effectively strong between slowly moving DM particles, the so-called Sommerfeld effect. For minimal DM models and the minimal supersymmetric standard model, this nonperturbative effect is by now routinely included at leading order in the calculation of the forces in both, the prediction of signals of annihilating DM in cosmic ray fluxes [4][5][6][7][8][9], and computations of relic abundance [10][11][12][13]. The effect also appears in non-WIMP DM models, as long as there exists a suitable hierarchy between the DM mass and the mass of a light boson coupled to it [14].
Other loop effects are often important as well. In scenarios, where the DM particle originates from an electroweak multiplet, the Sommerfeld enhancement depends strongly on the mass differences [15] among the members of the multiplet. For the simplest multiplets, the mass splittings are known up to two loops [16][17][18]. Further, in annihilation to final states with identified particles with electroweak charge, perturbation theory also breaks down, because the radiative corrections to the Born cross section are further enhanced by large logarithms of the ratio of the DM particle mass to the electroweak gauge boson mass [19][20][21]. Consequently, the fixed-order computations must be complemented by all-order resummation of the dominant logarithmic corrections. This has been achieved for high-energy cosmic photons with the help of soft-collinear effective field theory for the electroweak fermionic triplet ("wino") [22][23][24][25][26] and doublet ("Higgsino") DM model [27]. Overall, the most advanced computations of the high-energy photon yields from DM annihilation reach one to few percent accuracy, depending on the DM mass and model. This motivates a closer scrutiny of the calculation of the Sommerfeld effect, which is usually calculatedwith the leading order (LO) potential generated by electroweak gauge boson and photon exchange.
In the previous letter [28], we discussed, for the first time, the Sommerfeld effect for wino DM with one-loop, next-to-leading order (NLO) corrections to the electroweak Yukawa potential. More precisely, we considered the pair annihilation of the DM particle χ 0 into γ + X, and found the NLO potential to give a sizeable correction to the LO Sommerfeld effect, including a shift of the Sommerfeld resonance positions by about 6%. For state-of-the-art theoretical predictions of high-energy photon signals for indirect detection experiments, such as the Cherenkov Telescope Array, the NLO computation of the Sommerfeld enhancement is therefore indispensable. The present paper serves two purposes. First, we extend the one-loop computation to the potentials in the coannihilation channels χ 0 χ ± , χ ± χ ± , which were not given in [28], and perform the first computation of a DM relic abundance with the NLO Sommerfeld effect. Second, we provide analytic results for the potentials in momentum space, and the technical details of the NLO computation. We provide all the one-loop integrals relevant for the evaluation of the NLO correction and discuss the properties of the NLO potential functions.
The outline of the paper is as follows. In Section 2, we discuss the construction of the EFT for the WIMPs, introduce the power-counting, and review the calculation of the tree-level potentials. Subsequently, in Section 3 we discuss the computation of the one-loop correction to the potential for the wino model in all channels including details on renormalization schemes, asymptotic behaviours, and the parameter dependence. In Section 4, we calculate the relic abundance and analyze the size of the correction. We conclude in Section 5. In a series of appendices, we collect additional technical details on the one-loop calculations in Feynman and general covariant R ξ -gauge, Fourier transforms between momentum and position space, and relevant expressions for the asymptotic behaviours.

EFT of non-relativistic WIMPs
The low-energy effective field theory (EFT) of non-relativistic WIMPs was constructed in [29][30][31] in analogy with the respective non-relativistic EFTs of QED and QCD [32][33][34][35] for onium systems [36]. In this section, we review the structure of the potential nonrelativistic effective theory for wino DM, establish a consistent power-counting, and identify the leading corrections to the potential. The essence does not depend on the particular wino DM model and is applicable to general multiplets and cases that include hypercharge. Our starting point is the SM supplemented with the wino Lagrangian where χ denotes an SU(2)-triplet of Majorana fermions and D µ is the covariant derivative, D µ = ∂ µ − ig 2 W a µ T a . The modes relevant to construct the non-relativistic WIMP EFT are (i) hard (k 0 ∼ m χ , k ∼ m χ ), (ii) soft (k 0 ∼ k ∼ m W ), (iii) potential (k 0 ∼ m 2 W /m χ , k ∼ m W ), and (iv) ultrasoft (k 0 ∼ k ∼ m 2 W /m χ ). We introduce the power-counting parameter in terms of the W -boson mass, and assume that the DM mass is such that m W /m χ ∼ v ∼ α 2 , where v is the small velocity of the DM particles. A different relative size of v compared to α 2 or m W /m χ does not affect the construction of the non-relativistic EFT, but leads to different regimes (e.g., Coulombic if α 2 ∼ v m W /m χ ). In the first step, we match to the non-relativistic EFT, i.e., we integrate out the hard modes. This step is performed in the unbroken phase where the electroweak symmetry is still manifest. The theory resembles NRQCD with an SU(2) gauge group, hence, the Feynman rules are known [37] upon appropriately adjusting the group factors. The non-relativistic Lagrangian terms that we need in this paper are simply given by where only the soft, potential, and ultrasoft modes are the dynamical degrees of freedom.
Next, electroweak symmetry breaking is implemented. The resulting Lagrangian is of the same form as above; however, the relevant degrees of freedom change from weak eigenstates to mass eigenstates. The fields χ a , where a = 1, 2, 3 in the unbroken theory, are combined to mass eigenstates with χ ± = (χ 1 ∓iχ 2 )/ √ 2 and χ 3 = χ 0 . The electrically charged states χ ± acquire a radiatively induced mass splitting δm χ = 164.1 MeV [17] with respect to the neutral χ 0 state. 1 Finally, we integrate out the soft fields and potential gauge bosons. We obtain nonlocal (in space) four-fermion operators whose matching coefficients coincide at tree-level with the classical static potential. Loops of soft fields induce quantum corrections to the DM potential. We are left with a theory of potential fermions and ultrasoft gauge bosons, which is described by the potential non-relativistic DM (PNRDM) Lagrangian [29] L PNRDM = i=±,0 where e i is the electric charge of fermion i in units of the positron charge e. The electromagnetic covariant derivative is given by iD 0 i = i∂ 0 + e i eA 0 . In our convention, all DM fields are particle fields (cf. [29]). Structurally the above Lagrangian is the same as for QED or QCD. The phenomenology, however, is different as the gauge symmetry is broken. For example, the potentials can be off-diagonal for the mass eigenstates, unlike in QCD, where quarkonium states decompose into gauge eigenstates, which are simultaneously mass eigenstates (singlet/octet), and where the potentials are diagonal in the space of quarkonium mass eigenstates.
The ultrasoft fields in the gauge-covariant derivative D 0 i (t, 0) and the electric field are multipole-expanded and only include the photon field, as, in the broken theory, the W -and Z-bosons are too heavy to have ultrasoft scaling. 2 The term x·E originates from the ultrasoft interaction of charged fermions with the photon field and the additional coupling of the photon to the W -boson Yukawa potential after application of equation-ofmotion identities [38]. It makes the unbroken ultrasoft electromagnetic gauge symmetry manifest, and, relative to the leading kinetic term, it is suppressed by (m W /m χ ) 3/2 . Hence, it will play no role in the determination of the NLO correction. 3 1 In the non-relativistic theory the mass correction is obtained from the soft Z-and W -boson correction to the DM field propagator, see Appendix A.1.3.
2 For very large DM masses m χ ∼ O(100 TeV) there is a regime where m W /m χ ∼ α 2 2 . For this regime, the W -and Z-boson can contribute to the ultrasoft radiation. However, in such a case, the potentials would be effectively Coulombic unbroken-theory potentials, and the effective theory would be very similar to PNRQCD for scattering states, as α 2 ∼ v m W /m χ . 3 The leading ultrasoft correction comes from A 0 inside the covariant derivative, which is only suppressed by (m W /m χ ) 1/2 relative to the leading terms. For S-wave annihilation, the ultrasoft photons couple only to the total electric charge Q of the wino two-particle state. Although Q is non-vanishing for the χ 0 χ ± and χ ± χ ± states, there is nevertheless no contribution to the total annihilation cross section, The ultrasoft interactions, however, are relevant in determining the DM bound-state formation rates, which can modify indirect detection signals and the DM relic abundance. The effect is especially significant for large multiplets. However, for wino DM in the few TeV mass range it is important neither for indirect detection [39] nor the relic abundance [40].
The mass difference term δm i ∼ α 2 m W ∼ m 2 W /m χ is of the same order as the kinetic terms by power-counting and therefore contributes at leading order, even though δm i is a one-loop effect. Hence, to obtain NLO accuracy of the calculation, we include the two-loop result for the mass splitting [17].
The crucial new ingredient in obtaining NLO accuracy is the potential term in (2.3). NLO corrections to this term could be twofold. First, from potentials that are more singular than 1/r, but the structure of the non-relativistic EFT vertices implies that such potentials can appear only from the next-to-next-to-leading order in the powercounting (similar to QED/QCD). Second, from the one-loop correction to the tree-level Coulomb and Yukawa potentials, which are the subject of this paper and [28].
For completeness, we recall that the LO potential is obtained from tree diagrams involving the exchange of an electroweak gauge bosons between two wino particles. In the neutral two-particle sector we find (in momentum space) where the entries refer to the χ 0 χ 0 , χ + χ − , χ − χ + states, respectively, and T χχ→χχ denotes the T-matrix in the specific scattering channel. In the single-charged and double-charged sectors, one has where the entries refer to χ 0 χ + , χ + χ 0 and χ + χ + , respectively. The same expressions hold for the Q = −1 and Q = −2 charge sectors. For solving the Schrödinger equation, one transforms the potentials to coordinate space, using At tree level we need the Fourier transform which is related to the forward-scattering amplitude, that is, the matrix element of a local four-fermion operator whose net charge vanishes. We checked this statement by performing an explicit one-loop computation for the wino model. Therefore, at leading order we encounter Coulomb and Yukawa potentials only. In the above basis, which we refer to as method-I, following [31], the spin and angular momentum configuration of the initial states do not play a role. It is more conventional [3,31] to decompose the two-particle states into partial-wave configurations 2S+1 L J of definite total angular momentum L and spin S. The resulting potential is referred to as method-II and in the neutral sector reads (in coordinate space) Similar decompositions hold for all other two-particle states and the P -wave potentials.
The decomposition into the method-II two-particle states removes redundancies among the two-particle states. It automatically implements the symmetry properties of the underlying states, e.g., that the identical Majorana particles χ 0 χ 0 cannot exist in a 3 S 1 spin configuration. A detailed discussion of the correspondence between method-I and method-II can be found in [31].
The calculation of the one-loop correction to the potential proceeds in an analogous fashion. First, the momentum-space potential is calculated in the form of (2.4). Then the Fourier transformation to position space is performed. However, more complicated functions will lead to a wider variety of potentials at the one-loop order. The transition from method-I to method-II follows the same rules as the tree-level potential.

The NLO potential
In this section, we turn to the calculation of the NLO potentials. We describe in detail the case of the χ + χ − → χ + χ − scattering channel, for which all possible diagram topologies contribute, and the effects of EWSB play an important role (particularly through γ − Z-mixing). The results for the other channels are discussed subsequently, and the differences are highlighted. In doing so, we provide technical details and analytic expressions not supplied in [28]. We further discuss the renormalization scheme, the gauge invariance of the results, the large-and small-distance asymptotic behaviour, and the top-quark mass dependence of the NLO correction.

The
The one-loop correction to the potential in the χ + χ − → χ + χ − scattering channel originates from the non-relativistic EFT diagrams depicted in Fig. 1. The relevant loop momentum for matching the potentials has the soft scaling in the framework of the threshold expansion [41]. This amounts to replacing the DM propagators by static ones, as the soft momenta throw the non-relativistic propagators off-shell. In the calculation, it is important to ensure that the pinch poles at k 0 = ±iε are not picked up, as they are reproduced in the EFT by iterations of the LO potentials and belong to the potential momentum region. The results for the above diagrams are given in Appendix A. The box and crossed box diagram involving photons and Z-bosons cancel each other, similarly to the photon boxes in QED. However, there is a box contribution from Wboson exchange that has no crossed box counterpart. Furthermore, the self-energy diagrams are not diagonal but mix the photon and the Z-boson. The self-energies involve all SM particles and were evaluated in general covariant gauge using FeynArts [42], FORMCalc [42] and Package-X [43] and checked against the Feynman-gauge results (excluding tadpoles) from [44]. Lastly, let us comment on the second diagram in the second row of Fig. 1. This diagram vanishes in Feynman gauge, as the vertices on the external fermion lines project out the zero-component of the gauge-boson propagator. However, in general covariant R ξ -gauge this diagram is non-zero and required to obtain a gaugeparameter independent result. Before we assemble the full correction, we first discuss the on-shell renormalization scheme that we used in our computation.

The renormalization scheme
We choose to renormalize the ultraviolet (UV) divergences in the on-shell scheme following [44]. 4 As the matching between the non-relativistic and the potential non-relativistic theory is performed at the scale m Z it is natural to choose the input parameters at this scale. As the input parameter set, we use and set the CKM-matrix to the unit matrix. The counterterms in this scheme are for the gauge-boson masses and for the electromagnetic coupling at the Z-resonance. The self-energies on the right-hand side are evaluated in dimensional regularization with all fermions other than the top quark taken to be massless, and scaleless integrals are dropped. For later convenience it is also useful to define the counterterms for the Weinberg angle Equipped with these definitions, we can now assemble the potential correction. Before doing so, let us comment on the treatment of tadpole diagrams. Tadpoles in the electroweak theory are ubiquitous, and several treatments are possible. In the end, regardless of the scheme adopted, their contribution cancels in physical observables [44], and hence, in principle, we do not need to consider them. However, keeping the tadpoles has its merits as the gauge boson self-energies, including tadpoles, are gauge-invariant on-shell [45], and so are the mass and coupling counterterms, which helps to demonstrate the gauge-invariance of observables [45,46]. Therefore, we will keep the tadpole contributions to the self-energies as this will make the discussion of gauge-invariance more transparent.

The correction to the momentum-space potential
The previous considerations allow us to assemble the full one-loop correction in the on-shell renormalization scheme. We obtain in terms of box, vertex and self-energy functions, and the counterterms. The explicit results are lengthy and provided in Appendix A. The first large square bracket corresponds to the correction to the Coulomb potential, namely the vertex corrections (A. 6 16). These are the only terms that are not directly associated with one of the tree terms (though the γ-Z-mixing contribution could be partial-fractioned and grouped with the tree terms). We checked that the UV and IR poles cancel in (3.5) and that the expression is gauge-invariant, as discussed in detail below.

Gauge-invariance of the potential
The potential (3.5) is explicitly gauge-invariant. As discussed above, the inclusion of tadpoles ensures that the on-shell self-energies, charge counterterm δZ e , Weinberg angle counterterm δs W , and the Z-mass counterterm δm 2 Z are separately gauge-invariant. The further cancellations between box, vertex, and self-energy diagrams are analogous to those for SM processes (see, e.g., [46] for an extensive discussion of this issue), but with diagrams expanded in the non-relativistic/soft region. Contrary to e.g., HQET, the cancellation between wave-function renormalization and the vertex correction (lower left in topology in Fig. 1) is only partial, because we work in the mass eigenbasis and not in the weak eigenbasis. The remnant pieces are precisely the ones needed to complete the cancellation with the other diagrams.
Since the fermion self-energies are separately gauge-invariant, we can split the potential into three separately meaningful corrections: (1) A pure electroweak correction, which includes all contributions of gauge and Higgs bosons, except the fermionic selfenergies, and also incorporates the parts of the counterterms that originate from the respective topologies. (2) The light fermionic contribution incorporates all massless fermion loops except (3) the third generation quarks, which are separated for illustrative purposes (again also including the respective parts of the counterterms). Note that although the bottom quark is taken to be massless, its contribution is not separable from the top quark, e.g., in the W -self-energy, as they form an SU(2) doublet, hence it belongs to (3).

The remaining channels
The calculation of the potential correction in the other channels follows a similar logic. In the χ 0 χ 0 → χ 0 χ 0 scattering channel, where the tree-level potential is vanishing, only the W -boson box and crossed box topologies are possible. They exactly cancel each other, such that The potential in the off-diagonal channel χ 0 χ 0 → χ + χ − is also easily assembled. Except for the on-shell counterterms associated with the tree-level contributions, the other topologies are similar to the χ + χ − → χ + χ − channel (adapted to the exchanged Wboson). The only topologies that again deserve a special comment are the boxes. Crossed boxes are not possible, as the χ 0 couples only to W -bosons. The boxes are comprised of a W -boson and either a photon or a Z. The complete correction to the off-diagonal potential reads The individual terms are-as before equation These channels are sufficient for the calculation of the indirect detection cross section from χ 0 χ 0 annihilation as outlined in [28]. For the DM relic abundance computation, the charged co-annihilation channels are also needed. For the singly charged channel χ 0 χ ± → χ 0 χ ± , the topologies are similar to the χ 0 χ 0 → χ + χ − channel as the tree-level potentials are the same up to a minus sign. For the box diagrams, only a crossed box is possible due to the charge flow, which in turn leads to the conclusion that the correction in this channel is exactly the negative one of the χ 0 χ 0 → χ + χ − channel: Finally, we consider the doubly charged channels. Similar to before the correction is identical up to a minus sign to the χ + χ − → χ + χ − channel, as the same tree-level structures are involved, i.e.
The only difference stems from the fact that only the crossed W -boson box is possible, compared to the normal W -box topology in the charge-0 channel. This leads to the overall minus sign. Gauge invariance for all these channels can be checked as for the χ + χ − → χ + χ − channel above. Finally, we also checked for all channels that in the limit m W → m Z (i.e., s W → 0, c W → 1) we reproduce previously known results for the Higgsed SU(2) theory [47,48]. More precisely, the unrenormalized potential (Eq. 16 of [47]) was compared analytically. The renormalized result was not compared, as the renormalization scheme was not fully specified in [47].

Analysis of the channels
For an investigation of the Sommerfeld effect and other applications, one solves the position-space Schrödinger equation. The analytic and numerical Fourier transforms required to obtain the position-space NLO potential are given in Appendix B. Here we discuss the charge-neutral channels, as the remaining channels have equal corrections up to a minus sign. The results are shown in Figs. 2 and 3. In the following, we discuss the asymptotic behaviours and relative importance of the NLO correction. As input parameters for the numerics and plots, we use the following: the on-shell electromagnetic coupling α = α os (m Z ) = 1/128.943 at the Z-mass, and the gauge boson masses m W = 80.385 GeV and m Z = 91.1876 GeV. The SU(2) coupling and the Weinberg angle are determined by the on-shell relations α 2 = α os (m Z )/s 2 W and c W = m W /m Z . Furthermore, we need the Higgs-boson and top-quark mass, for which we take the onshell masses m H = 125 GeV and m t = 173.1 GeV. The uncertainty of these parameters is small enough to be ignored, except for the top-quark mass.

The asymptotic behaviour of the NLO potentials
We start with the discussion of the small-and large-distance asymptotics of the various channels. For technical reasons (not all Fourier transforms are analytically available), we discuss the results in |k| space, separately for the various gauge-invariant pieces to outline the origin of the corrections and the dominant contributions.

Full correction
Light fermions 3rd gen. quarks Electroweak Figure 2: The NLO correction to the potential in the channel χ + χ − → χ + χ − . The upper panel shows the modulus of the potential |r · V (r)| for the LO and NLO potential, the NLO contribution only, and the small and large-distance asymptotic behaviour. In the lower panel, we show the ratio of the full NLO potential to the LO potential (blue solid), and separately for the three gauge-invariant pieces identified in the text (other curves).
This behaviour is similar to QED/QCD, namely, the prefactor of the logarithmic term is proportional to the light-fermion contribution to the SU(2) beta function. For the contribution of the third generation quarks we find

Full correction
Light fermions 3rd gen. quarks Electroweak Figure 3: The NLO correction to the potential in the channel χ 0 χ 0 → χ + χ − . The upper panel shows the modulus of the potential |r · V (r)| for the LO and NLO potential, the NLO contribution only and the asymptotic behaviours. The change from solid to dashed for the blue δV (r) curve marks its change of sign. In the lower panel, we show the ratio of the NLO to the LO potential for the full correction, and for the three gauge-invariant pieces, which illustrates their different behaviour.
and Z mass given in Appendix C. To per mille accuracy in the interval of ±10 GeV around the on-shell top mass m t,os = 173.1 GeV it can be approximated by . The most complicated contribution originates from the electroweak corrections, i.e., the gauge and Higgs bosons. Contrary to the two former contributions, the non-self-energy diagrams also contribute here. We find where the logarithmic term is proportional to the non-fermionic part of the SU(2) beta function. B(m W , m Z , m H ) is a function of the Higgs, W -and Z-mass given in Appendix C and evaluates for on-shell parameters to −1.03577. The analytic result displays an interesting manifestation of the screening theorem [49]. Even though individual terms are Higgs-mass dependent up to m 6 H , B itself is only logarithmically dependent on the Higgs mass for large m H .
Adding all three separately gauge-invariant pieces, and defining where γ E is the Euler-Mascheroni constant. The last line provides the numerical value for the on-shell parameters as was already given in [28]. Of the numerical coefficient, the light-fermion term makes up −1.3188, the third generation quarks 8.59038, the electroweak terms −0.51788, and Euler-Mascheroni constant associated with the logarithm −1.82785. The identical short-distance behaviour of the diagonal and off-diagonal channels is also visible by comparing the lower panels of Figs. 2 and 3. The logarithmic behaviour implies a breakdown of perturbation theory, as ln(m Z r) grows arbitrarily large for small r. This is a consequence of renormalizing the parameters on-shell. The logarithmic behaviour can be absorbed by using MS running couplings as will be discussed in Sec. 3.3.3. However, let us note that using the on-shell renormalized potentials is sufficient in the calculation of the Sommerfeld effect, where the dominant contribution comes from the region m W r ∼ 1, where the difference between various renormalization schemes is of higher-order (which we also checked numerically).
The r → ∞ / k 2 → 0 limit In the opposite limit r → ∞, we have to distinguish between the χ 0 χ 0 → χ + χ − and χ + χ − → χ + χ − scattering potentials. We begin with the latter and the light-fermion contribution that scales according to the U(1) em beta function contribution of the massless fermions. At large distances, respectively, small momenta, the potential is dominated by photon exchange, which explains the transition to the electromagnetic beta function. This also holds for the third generation quarks, where the coefficient is determined by the massless bottom-quark contribution to the U(1) em beta function: The top-quark contribution is cut off due to the finite mass and therefore does not contribute to the asymptotic behaviour. The electroweak contribution does not play a role in the large-distance behaviour of the χ + χ − → χ + χ − potential, as it is cut off by the boson masses. It starts with a constant term, where the function C(m W , m Z , m H ) is given in Appendix C and for on-shell parameters evaluates to 3.67219. Even though it does not contribute to the asymptotic behaviour, the result is a good check of the calculation through its m H dependence. As required by the screening theorem [49], it is logarithmic in m H , even though the individual terms depend on the Higgs mass with up to m 6 H . Overall we find that the dominant r → ∞ behaviour in the χ + χ − → χ + χ − channel is given by the purely abelian correction to the Coulomb potential due to massless fermions, where β 0,em = −80/9 is the electromagnetic beta function coefficient for all SM fermions except the top quark. This asymptotic behaviour dominates the correction to the potential for m W r ≥ 5 as can be seen in Fig. 2.
In the channel χ 0 χ 0 → χ + χ − , the large-distance asymptotics also originates from the light-fermion terms. The relevant terms are which scales as k 2 ln(k 2 /m 2 Z ) for k 2 → 0. n ld denotes the number of massless fermion doublets, in our case n ld = 5. The Fourier transforms for the individual terms are discussed in detail in Appendix B. After expanding for large r, we find   This power-like long-range behaviour is a consequence of taking the SM fermions to be massless. 5 It is also a manifestation of a breakdown of perturbation theory, as for r 1/m W , the correction exceeds the exponentially decreasing tree-level potential. For the later physics applications this does not pose a problem, since the Sommerfeld effect is governed by distances m W r ∼ 1. We checked this numerically and confirmed that the region where the power-like long-range potential dominates does not affect the calculation of the Sommerfeld factors in any significant way.
In reality, the SM fermions are, of course, not massless. A formal treatment of the r → ∞ limit would require a further matching procedure, where the W -mass scale is integrated out. The resulting theory predicts the same r −5 asymptotics as above. In the next step, one would successively match on theories where the individual fermions acquire mass m f , which will then cut off the contributions at distances r ∼ 1/m f , similar to the top-quark contribution discussed below.
The third-generation quarks and the electroweak piece start with a constant in the Taylor expansion around k 2 = 0 and are therefore exponentially suppressed for large r. Explicitly, where the functions D, E are given in Appendix C and evaluate for on-shell values to 14.6515 and 2.76239, respectively. Again the screening theorem is fulfilled by these expressions. A breakdown of perturbation theory at large r manifests itself also in these channels, as the tree-level potential is exponentially suppressed. Terms such as from the gauge-boson mass renormalization behave as compared to the tree-level exp(−m W r)/r. However, these terms are subdominant compared to the light-fermion tail and therefore contribute even less to the Sommerfeld factor. Dyson resummation would cure this behaviour and result in potentials of the form exp(−(m 2 W + δm 2 W ) 1/2 r)/r. Overall the contribution in the χ 0 χ 0 → χ + χ − channel for large r is given by the light-fermion contribution as discussed above and reads δV r→∞ which explains the steep increase in Fig. 3 around m W r ≈ 10. As also seen in this figure, the third generation quark and electroweak contributions are exponentially suppressed for large r.

The complete NLO corrections
The exact NLO potential interpolates between the r → 0 and r → ∞ asymptotics. The most significant deviations from the asymptotics are observed around m W r ∼ 1, which is the crucial region to determine the Sommerfeld effect accurately. Therefore it is not sufficient to simply glue the asymptotics together. For an accurate determination, either the full numerically calculated potential or the fitting functions provided in [28] have to be used. They read for the off-diagonal potentials in (2.8) and (2.9) and for the diagonal ones where x = m W r, x 0 = 555/94 and L = ln x. The fitting functions provide per mille accuracy for the Sommerfeld factors [28]. In general, the correction to the Coulomb and Z-Yukawa potential is closer to the full numerical result. The reason is the sign change for the W -Yukawa potential at x 0 = 555/94. The position of this sign change is set by the distance where the light-fermion contribution starts to dominate the correction. Although the correction for very small (large) r is significant due to the logarithmic (power-like) behaviour, these regions contribute little to the Sommerfeld factors. In the relevant region m W r ∼ 1, the NLO correction to the potentials is in the few percent range. The complete NLO result is determined by the interplay of the various corrections. For example, for r → 0 in both, the χ + χ − → χ + χ − and χ 0 χ 0 → χ + χ − channels, the correction due to light fermions is of opposite sign to the electroweak contribution.

Scheme conversion to MS-couplings
The on-shell scheme employed for the calculations so far exhibits large short-distance logarithms related to the beta function. This behaviour originates from on-shell renormalization at the scale m Z , which is suitable for the calculation of the Sommerfeld effect, but leads to logarithms of the form ln(m Z r).
It is more appropriate to use running couplings at the scale µ = e −γ E /r or µ 2 = k 2 in position or momentum space, respectively, if one is interested in the potential at short distances. To this end, we convert the on-shell coupling to the MS-scheme using where δs W was defined in (3.4) and on-shell parameters were used for all terms involved. For numerics in the MS scheme, we use the MS top mass m t (m t ) = 163.35 GeV. To keep notation short, from here on couplings in the MS-scheme are denoted by a hat. With these ingredients, the issue of large logarithms in the r → 0 asymptotics can be further investigated. To see the cancellation of the large logarithms, the MS-coupling at m Z is converted to the coupling at an arbitrary scale µ by expanding the running couplings to fixed order, which, e.g., for the tree-level W -Yukawa potential in momentum space leads to For |k| → ∞, this exactly cancels the logarithmic contribution in the asymptotic behaviour (3.14). Splitting the beta function contribution into β 0, SU(2) = 19/6 = 43/6 − 1 − 3, where terms correspond to the electroweak, third-generation quark, and lightfermion contributions, respectively, this can also be done for each of the separately gauge-invariant pieces. In position space, a similar expansion can be performed using which cancels the logarithms for r → 0. The MS scheme presented here applies to momenta and distances of k 2 > m 2 W and 1/r > m W , respectively. 6 For the tree-level Coulomb and Z-Yukawa potential, the logarithms for r → 0 can be eliminated in a similar fashion using where we have used the beta function for the hypercharge β 0,Y = −41/6 = −1/6−11/9− 49/9 (which can be split into electroweak, third-generation quarks and light fermions, respectively). As expected from the tree-level potential, sinceα +α 2ĉ 2 W =α 2 , the hypercharge contribution drops out in the limit k 2 → ∞ when the two above potentials are summed, and the logarithmic contribution cancels in the asymptotic behaviour (3.14). A similar argument works in position space.
In Fig. 4, we show the absolute value of the potential in the MS-scheme using the position-space conversion (3.30) and one-loop running couplings at the scale µ = e −γ E /r. While the NLO and LO potentials diverge for small r in the on-shell scheme due to the breakdown of perturbation theory, the NLO correction remains always small in the MSscheme and the correct short-distance behaviour is already attained at tree-level, due to the use of the running scale. This behaviour is also exemplified by the inset of Fig. 4, which shows the ratio of the NLO to the LO potential. It also shows that it does not matter whether the running coupling is implemented in position or in momentum space, as it should be.
The MS-scheme is clearly the better scheme for large momenta. Solving the Schrödinger equation to obtain the Sommerfeld effect technically probes all momentum regions, but we find that the changes for the Sommerfeld factor for various MS approximations and the on-shell result are compatible with differences of the size of well-behaved electroweak corrections beyond the one-loop order considered here. We therefore stick limit. The reason is hidden in the fact, that the Fourier transform ofα 2 (k 2 ) is proportional to r −3 . For example, the Fourier transform of (3.29) together with the NLO terms (3.7) produces the large-r asymptotics while in position space the asymptotics is the same as for the on-shell potential (3.20) (exchanging the couplings). A more detailed discussion of the Fourier transform that leads to this behaviour is found in Appendix B. The r −3 dependence is not a conceptual problem for two reasons. First, similar to the on-shell case (3.20) that shows an r −5 behaviour, another EFT would need to be constructed that integrates out the massive bosons and keeps only light fermions dynamical. Secondly, the MS-scheme is designed to absorb the logarithms that grow large for r → 0 / k 2 → ∞ and is therefore not expected to work in the opposite limit anyway.  with the on-shell scheme for the computation of the relic abundance. However, let us note that the conceptual control over the r → 0 logarithms is an essential check of the calculation and demonstrates perturbative control over the potential correction.

Top mass dependence
The input parameter uncertainties have a negligible impact on the accuracy of the potential except for the top-quark mass. The top-quark mass first enters the potential at NLO, and the dependence on it is not only logarithmic, but also quadratic. At this point it would be possible to use the pole mass m t = 173.1 GeV or the corresponding MS-mass (at four loops) m t (m t ) = 163.35 GeV. Since the scheme ambiguity is not fixed at NLO accuracy for the potential, both choices are legitimate input values. The difference of 10 GeV causes by far the largest uncertainty of all input parameters. For the other parameters, the dependence is negligible, as they either already enter at leading order (W/Z-mass and couplings) thus the scheme dependence is reduced by the NLO correction, or the dependence is only logarithmic due to screening [49] (Higgs mass), or they are known precisely anyway.
To estimate the top-mass dependence, we investigate the function A(m W , m Z , m t ) given in Appendix C that controls the size of the Coulomb term in the r → 0 asymptotics. which is a 17 % decrease of the coefficient of the Coulombic behaviour for r → 0 for the full correction to the potential. However, in this region also the logarithmic term contributes, which is of similar size, decreasing the effect of the correction to roughly 10 % of the NLO correction.
In Fig. 5 we show the ratio of the full NLO potential to the LO potential for a range of r centred around 1/m W for different values for the top mass (left panel). The figure shows that the top mass dependence is more relevant in some regions than in others, but it does not change the gross features of the NLO correction. In the right panel of the figure, the ratio to the default value m t = 173.1 GeV is depicted. In both diagonal (+−) → (+−) and off-diagonal (00) → (+−) channels, the top mass affects the NLO correction by up to 10 % with the largest change around m W r ∼ 1. At large r, the precise top-mass value employed does not matter, as the top contribution becomes negligible in comparison to the light fermions. In the χ 0 χ 0 → χ + χ − channel, the singularity of the ratio around m W r ∼ 5 in the right panel is an artefact of showing the ratio, as δV (r, m t ) changes sign at slightly different m t -dependent values of r.
The top-mass uncertainty of the potential translates into a small effect on the Sommerfeld factor, which has already been investigated [28]. For example, the location of the first Sommerfeld resonance is shifted due to the NLO potential correction from 2.283 TeV to 2.408 TeV instead of 2.419 TeV when the MS-mass m t (m t ) = 163.35 GeV instead of m t = 173.1 GeV is adopted. This effect is small enough to be ignored for the present, hence in the following we will stick with the pole mass value m t = 173.1 GeV.

Wino relic abundance
In this section, we compute the WIMP relic abundance under the thermal freeze-out assumption and discuss the importance of the new NLO correction to the Sommerfeld potential.

Technical details on the DM abundance calculation
The computation divides into the calculation of the Sommerfeld factors for various partial-wave cross sections in all co-annihilation channels, followed by the thermal average and solution of the Boltzmann equation. Our implementation follows [31], to which we refer for a detailed description and notation employed here.
As the freeze-out process starts for DM relative velocities v ∼ 0.2, the inclusion of O(v 2 ) corrections to the annihilation cross sections is necessary to obtain percentlevel accuracy. We therefore include P -wave and O(v 2 )-suppressed S-wave annihilation. All required annihilation matricesf ( 2S+1 L J ) are conveniently tabulated in Appendix C of [30]. In these short-distance quantities, we use the MS couplings evolved with one-loop accuracy to the scale µ = 2m χ . The partial-wave Sommerfeld factors for the annihilation cross section of the two-particle DM state χ i χ j are given by where ψ (L,S) e 1 e 2 , ij denotes the wave function at r = 0 for the initial state ij in partial-wave configuration with angular momentum L and spin S to scatter into the state e 1 e 2 under the influence of the potential. The annihilation cross section in the channel ij that enters the thermal average is then obtained to O(v 2 ) accuracy by weighting each Born partial-wave term by its respective Sommerfeld factor, resulting in 7 Technically, we determine the Sommerfeld factors using the variable-phase method to solve the Schrödinger equation developed in [31]. This requires a fast and numerically stable evaluation of the NLO Sommerfeld potential in coordinate space, which we obtain by precalculating and interpolating the numerical Fourier transform where necessary. The Schrödinger equation is then solved from an initial value x 0 = m χ vr 0 = 10 −7 to some large x ∞ = m χ vr ∞ , which is determined using an adaptive procedure, which terminates when doubling an already large initial x ∞ changes the Sommerfeld factor by less than three per mille. For points near the χ + χ − -threshold this convergence criterion is sometimes hard to reach, and we abort the above procedure if x ∞ > 10 4 . The behaviour around the true value is oscillating and since we scan the threshold accurately, the 1% inaccuracies incurred by the abortion tend to average out.
We tabulate the Sommerfeld factor as a function of velocity in all relevant channels using 100 velocity points distributed logarithmically between v = 10 −4 and 1 and additional points around the two-particle thresholds v = 2δm χ /m χ and δm χ /m χ resulting in around 150 points for every partial-wave Sommerfeld factor. The resulting /(k 2 + m 2 φa ). The second term in (4.2) arises from a linearly divergent integral, and is finite but scheme-dependent in dimensional regularization, which has been used in obtaining (4.2). The scheme-dependence cancels with a one-loop correction to the short-distance annihilation matrix, but this is not available here. The generalization of the above identity to NLO potentials is not straightforward, since it would require a treatment of the singular short-distance behaviour in dimensional regularization. We therefore use the LO potentials here. This is justified, since the issue of the uncancelled scheme dependence for the O(v 2 ) suppressed S-wave terms is already present at LO and would not be improved by adding the NLO correction, but more importantly since in practice, the term in question represents a small correction to the cross section, as will be discussed at the end of this section.
cross-section tables are then monotonically interpolated for use in the velocity integration to obtain the temperature-dependent thermally-averaged effective annihilation cross section including co-annihilation. The thermal average is calculated in the variable x = m χ /T for 160 logarithmically distributed points between x = 1 and x = 10 8 .
These points are again monotonically interpolated. The resulting function forms the input to the Boltzmann equation solver. The differential equation is solved numerically with different methods, one of them simply Mathematica's built-in NDSolve, with an implicit solver to determine the yield Y (x) between x = 1 and x = 10 8 with initial condition Y (1) = Y eq (1). The relic abundance is obtained from Y (10 8 ). As input for the effective number of degrees of freedom in the Boltzmann equation we adopt the implementation from [54], extracted from the plots and tables therein, supplemented by results of [55] in the regions above T = 280 GeV and below T = 1 MeV. The critical energy density value equals ρ crit. = 1.05368 · 10 −5 GeV cm −3 h 2 .

NLO relic density for the wino model
The first zero-energy bound-state resonance for the χ 0 χ 0 ( 1 S 0 ) total annihilation cross section is located at m χ = 2.282 TeV for the LO potential and at 2.419 TeV for the NLO potential. The Sommerfeld factor for the total cross section is slightly different from the one for γ + X considered in [28], since in the latter case only the χ + χ − → χ + χ − component of the annihilation matrix enters. Therefore in that case only the wavefunction components ψ (00)(+−) and ψ (+−)(00) are probed. On the other hand, for the relic abundance calculation the annihilation matrix is non-zero in all entries χ 0 χ 0 /χ + χ − → χ 0 χ 0 /χ + χ − and the Sommerfeld calculation is sensitive to all components of the wave function. Nevertheless, the resonance masses are the same within sub-GeV accuracy as the ones found in [28] for the annihilation to γ + X.

Sommerfeld factors in individual channels
We compare the Sommerfeld-enhanced cross section to the Born cross section for the most important charge-neutral annihilation channels χ 0 χ 0 and χ + χ − in Fig. 6. The two panels refer to two mass values, one below and one above the first resonance. As the velocity decreases from right to left, the Sommerfeld factor increases, crosses the χ + χ − threshold at v = 2δm χ /m χ and reaches saturation for v < 10 −3 . The spikes slightly below the χ + χ − threshold are due to the Coulomb bound states [12], enlarged in the inset of Fig. 6. Formally, at the threshold the χ + χ − Sommerfeld factor becomes infinite as the relative velocity in this channel where v LSP is the velocity of the lightest DM particle. However, this will not spoil the abundance calculation, as in the thermal average, the 1/v +− divergence is removed by the integration measure d 3 v +− = v 2 +− dv +− dΩ. Below each panel, the effect of the NLO correction to the potential is highlighted by showing the ratio of the NLO to the LO Sommerfeld factor. For small velocities of the lightest two-particle state χ 0 χ 0 , the NLO potential correction is very important, causing O(1) changes in the Sommerfeld factor, as already observed in [26] for the γ + X final state. For velocities above 0.01−0.1, the modification is closer to the few percent of a typical electroweak one-loop correction. Both velocity regimes are important for the relic abundance calculation, and their precise weight depends on the temperature at which freeze-out ends. In particular, for DM masses near the Sommerfeld resonances, the small velocity region is important and freeze-out is delayed.

Thermally-averaged annihilation cross section
The thermally-averaged cross section σ eff v is shown in Fig. 7 for mass values chosen below and above, and at the resonance masses of the LO and NLO potential. At x 10, the non-relativistic approximation is not accurate, and in principle, a relativistic treatment is necessary. This, however, has a negligible impact on the final abundance, as freeze-out starts for x ∼ 20, where the non-relativistic expansion is already very precise, given that O(v 2 ) terms have been included. Around x ∼ m χ /δm χ , a small drop in the Born cross section due to the decoupling of the heavier χ + χ − channel can be seen. Except at the resonance, where the 1/v 2 ∼ x enhancement persists, the cross section For mass values below the first resonance, the NLO Sommerfeld-enhanced cross section is always smaller than the one computed with LO potential. After the first resonance of the NLO potential, the late-time thermally averaged cross section can be larger than for the LO potential, but even near the first resonance around m χ = 2.42 TeV, the NLO potential cross section exceeds the LO one only around x > ∼ 100. The reason is that the NLO correction weakens the potential, see Sec. 3. Therefore, the Sommerfeld effect is generally reduced, as seen in the small-x regime. Due to the weaker NLO potential, the first Sommerfeld resonance is shifted to a larger mass value. In a significant mass range above the resonance, the resonant behaviour can outweigh the suppression due to the NLO potential and the NLO cross section can be larger than at LO (see, in particular, the lower panel of Fig. 4 in [28]). Fig. 8 shows the suppression of the dark matter density due to the Sommerfeld effect at LO and NLO relative to the Born treatment. Until x ∼ 20 all yields are the same, i.e., the ratio is unity, as the freeze-out process has not yet started. Afterwards, there is a steep drop when the Sommerfeld-enhancement of the annihilation cross section leads to a faster depletion of the DM density, and simultaneously delays the end of the freeze-out process. Around x ∼ 10 4 , the late-time annihilations cease, and the yield ratio again becomes constant. An exception occurs for DM masses near the resonance values, where the cross sections are enhanced by 1/v 2 ∼ x in the low-velocity regime, leading to an additional O(1) late-time change of the yield, as seen in the second (resonance of the LO potential at 2.28 TeV) and third panel (resonance of the NLO potential at 2.42 TeV).

The dark-matter yield Y
It follows from the previous discussion of the thermally-averaged cross section that at early times (low x), the NLO yield is always larger than the LO one. Only at late times, the 1/v 2 (1/T ) enhancement of the (thermally-averaged) cross section for masses around the resonance value can lead to substantial late-time annihilations, and the NLO result drops below the LO result (see lower left panel of Fig. 8).

The wino relic abundance
Our main result, the wino relic density as a function of wino mass with NLO potential is shown in Fig. 9 and compared to the previous LO result, and the Born approximation.
The importance of accounting for the Sommerfeld effect is well-known and clearly seen in the figure. In comparison, the NLO correction is moderate, except near the location of the resonance, whose shift is very visible. Aside from the resonance mass region, the NLO relic abundance is larger by about (2 − 5) % than the LO one. The observed relic abundance Ω DM h 2 = 0.1205 [56] is attained for wino mass 2.842 TeV with NLO potentials compared to 2.886 TeV at LO (and 2.207 TeV for Born cross-section calculations). In comparison, the omission of the O(v 2 ) terms in the cross sections (4.3) would decrease the relic abundance by about (1 − 3) %.
It is difficult to quantify the accuracy of the NLO relic abundance. Probably the largest uncertainty arises from the missing one-loop radiative corrections to the hard annihilation process. The use of running couplings in the tree-level process has accounted for the dominant logarithmically enhanced correction, which has a considerable effect. For example, for m χ = 3 TeV, close to the value for the observed relic abundance, the relevant ratio of couplings isα 2 2 (2 × 3 TeV)/α 2 2 (m Z ) = 0.867. In Sudakov resummation studies of wino annihilation into γ + X [25,26], the non-logarithmic one-loop correction to the hard annihilation process was found to be at the 2 % level. Together with other sources of uncertainties, such as in the effective number of degrees of freedom, we estimate that the above result for the relic abundance is accurate to a few percent.

Conclusion
We provided a detailed description of our computation [28] of the NLO corrections to the electroweak Yukawa potential relevant to the Sommerfeld effect for wino dark matter. The calculation was performed in a general covariant gauge, and all necessary loop integrals were provided. The previous result was extended to include NLO effects in all co-annihilation channels. The main result of this paper is the first computation of a relic abundance with NLO accuracy for the Sommerfeld effect. Since the NLO correction weakens the potential, the value of the relic abundance increases by a few percent, except near the Sommerfeld resonances, where large effects are present and which are shifted towards larger masses. Our numerical investigations accentuate the importance of the NLO potential corrections for WIMP searches, and we advocate using the NLO Yukawa potential, which can easily be implemented using the fitting functions provided in Sec. 3.3.2.
The computations presented here might be extended to more general models, such as the minimal models with general SU(2) electroweak multiplets. Further, it would be interesting to investigate the Higgsino model and DM particles with hypercharge, and combine the NLO potential with the state-of-the-art computations of the high-energy photon spectrum [27], for which the missing NLO correction to the potential is likely to be the current largest theoretical uncertainty.
Acknowledgements We thank M. Drees and M. Laine for correspondence. This work was supported in part by the DFG Collaborative Research Centre "Neutrinos and Dark Matter in Astro-and Particle Physics" (SFB 1258).

A Expressions for the loop diagram topologies
In Sec. 3, we expressed the NLO correction to the potential in terms of functions representing diagram topologies. In this appendix, we provide the explicit results for these functions.

A.1 Feynman gauge results
We first collect the results for the diagram topologies depicted in Fig. 1 in Feynman gauge. In the following, we use the vector v µ = (1, 0) and the abbreviation where d = 4 − 2 and γ E is the Euler-Mascheroni constant. The calculation of the diagrams can be done with standard methods. Note that pinch poles at v · l = ±iε must not be taken into account in the integration. The residues of these poles are part of the Sommerfeld factor calculated with the leading-order potential.

A.1.1 Box topologies
We first consider the box topologies in the first line of Fig. 1. The crossed-box diagrams are related trivially to the box topologies by a minus sign in the case of the wino. We define λ(a, b, c) = (a − b − c) 2 − 4bc. In the case of two unequal-mass bosons, we find The box integral with masses is finite in d = 4. q 2 ≈ −q 2 refers to the exchanged potential three-momentum. If the masses of the exchanged bosons are equal, this simplifies to Finally there are the cases where one or both of the exchanged bosons are massless. In these cases the box diagrams, expanded in including the finite terms, result in

A.1.2 Vertex corrections
The soft correction to the vertex (first diagram in the second row of Fig. 1) is known as it appears also in soft corrections to exclusive DM annihilation [26]. We choose the convention The soft vertex correction involving the triple gauge-vertex (cf. Fig. 1) vanishes in Feynman gauge, as the vertices project on the zero-component of the propagator

A.1.3 DM field renormalization
The renormalization of the heavy DM field is analogous to HQET. The required integral is From the term proportional to v · p, we extract the on-shell static fermion field renormalization constants The v·p independent term generates the mass difference between the neutral and charged DM fermions in the non-relativistic theory. For the wino we find [57] δm χ 0 = −α 2 m W , (A.11)

A.1.4 Gauge boson self-energies
The self-energies in in Feynman gauge are not given explicitly for brevity, as we give the self-energies in R ξ -gauge in the following section. They can e.g. be found in [26,44].
where the gauge parameter with index i is associated with the boson of mass m i . For one vanishing mass (where ∆ j , ξ j are associated to the massless boson coupling with α j ) The massless box is given by (∆ i , ξ i and ∆ j , ξ j are associated to massless bosons coupling with α i and α j , respectively)

A.2.2 Vertex topologies
For the vertex diagram without the triple gauge-boson vertex, we find The vertex diagram with the triple gauge interaction vanishes in Feynman gauge. Which particles circle in the loop depends on the tree potential. For the Coulomb and Z-boson Yukawa potential, the loop is comprised of two W -bosons and, therefore, only depends on one gauge parameter ξ W . In the case of the W -boson Yukawa potential, the loop is a mixture of W -boson and either photon or Z-boson and hence depends on two gauge parameters and one, respectively two masses. The integral with W/Z bosons in the loop is given by If the Z-boson is replaced by the massless photon, we find Finally, if both bosons involved are W -bosons, we have

A.2.3 DM field renormalization
The heavy DM field renormalization constants in R ξ -gauge are

A.2.4 Gauge boson self-energies
The transverse self-energies (including tadpole diagrams) are split into the three separately gauge-invariant pieces discussed in the main text according to The photon self-energy The fermionic self-energy contributions are The electroweak part of the self-energy reads The photon-Z self-energy The Z-boson self-energy The fermionic contributions are given by The Z-and W -boson self-energies receive contributions from tadpole diagrams. In terms of the gauge-invariant parts, the tadpoles belong to the electroweak contributions and are related by a simple prefactor such that The electroweak part is then given by  Performing the |k|-integral analytically for as many terms as possible is crucial, as numerical instabilities can be avoided this way. For example, in the transforms of terms arising from Coulomb propagators, the 1/k 2 behaviour for k → 0 is not always manifest, but only in the linear combination with other terms.

B.1 Analytic transforms
In Table 1 we collect Fourier transforms that are helpful for the NLO potential. All NLO potential terms except logarithms involving the triangle function (these will be discussed in Section B.2) can be transformed using these results. In some cases partial-fractioning identities such as are necessary to bring the expressions into manageable forms. Some logarithmic terms can be obtained using Furthermore, some Fourier transforms are related by taking derivatives, e.g., (B.4) Lastly, the derivative of the Bessel K ν function for ν = ±1/2 is (the case relevant here) is ∂K ν (z) ∂ν ν=±1/2 = ± π 2z Γ(0, 2z) e z , where Γ(s, z) =    is the incomplete Gamma function.
As an example, let us apply the above tricks to (3.19). Using the second to last result of Table 1  which explains why the asymptotics is proportional to r −5 . In the case of the MS-scheme in momentum space (see footnote 6) the second term on the left-hand side of (B.8) is not present. Therefore in this case the asymptotic behaviour goes as r −3 (cf. (3.31)).

B.2 Numerical transforms
For certain terms involving the triangle function λ(a, b, c), we were not able to perform the Fourier transformation analytically. In the following, we discuss an example of how we obtain the numerical transform in such cases. The example arises from the self-energies, namely from the Higgs-Z loops, and is given by The numerical Fourier transform of this function is not always stable, depending on the value of r in (B.1) that determines the scales probed in the integrand. The leading behaviour for k → 0 is 1/k 2 . In the final expression, however, it may happen that the corresponding 1/r behaviour in position space cancels against another term in the complete expression for the loop integral. The subleading terms are exponentially suppressed for large r, which is hard to resolve numerically. We solve this issue by undoing the Feynman-parameter integration that led to the logarithm, writing .

(B.10)
We then subtract the Coulombic behavior in the limit k → 0 by adding the term 8 The Fourier integral over |k| can now be performed analytically, resulting in . (B.12) The above expression looks complicated, however, the numerical integration that has to be performed is now only an integral from x = 0 to 1 instead of an integral from |k| = 0 to ∞. This stabilizes the large-r tail as the exponential suppression of the final result is already captured in the integrand before the Feynman-parameter integration, and the numerical difficulty of accurately determining the exponential tail is thus circumvented.

C Expressions for the asymptotic behaviour
In Sec. 3.3.1 we discussed the asymptotic behaviour of the NLO correction to the potential in the limits r → 0 and r → ∞. In this appendix, we provide the mass-dependent functions appearing in the main text. The arctan terms in the expressions below stem from simplifying the real parts of self-energies.