Non-relativistic pair annihilation of nearly mass degenerate neutralinos and charginos III. Computation of the Sommerfeld enhancements

This paper concludes the presentation of the non-relativistic effective field theory formalism designed to calculate the radiative corrections that enhance the pair-annihilation cross sections of slowly moving neutralinos and charginos within the general minimal supersymmetric standard model (MSSM). While papers I and II focused on the computation of the tree-level annihilation rates that feed into the short-distance part, here we describe in detail the method to obtain the Sommerfeld factors that contain the enhanced long-distance corrections. This includes the computation of the potential interactions in the MSSM, which are provided in compact analytic form, and a novel solution of the multi-state Schr\"odinger equation that is free from the numerical instabilities generated by large mass splittings between the scattering states. Our results allow for a precise computation of the MSSM neutralino dark matter relic abundance and pair-annihilation rates in the present Universe, when Sommerfeld enhancements are important.


Introduction
For heavy neutralino dark matter (DM) there is class of radiative corrections to the tree-level pair annihilation cross section, which can become larger than naively expected by the electroweak nature of the interaction, and even exceed the lowest-order cross section by orders of magnitude. This so-called Sommerfeld effect arises when an attractive interaction between the non-relativistic DM particles significantly distorts their wave function, such that they have a larger probability to undergo annihilation. In terms of Feynman diagrams the effect arises from the exchange of the electroweak gauge bosons between the DM particles, which contributes a factor g 2 M DM /M W , such that an additional particle exchange is not suppressed by the electroweak coupling g 2 when the DM mass is much larger than the mediator mass. The Sommerfeld effect is non-perturbative in the sense that a resummation of diagrams to all orders in g 2 is needed in order to calculate the annihilation cross section. The relevance of the Sommerfeld effect was first pointed out for the annihilation cross section of (wino-or higgsino-like) neutralino DM into two photons [1], although it was not until 2008, when an anomalous positron excess was measured by PAMELA, that Sommerfeld enhanced DM models attracted more attention as a mechanism to boost the DM annihilation rates [2]. For these models, the larger velocities around DM freeze-out in the early Universe (v ∼ 0.3c) can yield an annihilation rate in accordance with the observed relic density, whereas the much smaller neutralino velocities today (v ∼ 10 −3 c) enhance the annihilation rate into electrons and positrons observed as cosmic rays. Since then the Sommerfeld enhancements in cosmic ray signatures and in the thermal relic abundance have been discussed extensively for relevant MSSM scenarios with neutralino DM [3][4][5][6][7][8][9][10], but also for generic multi-state dark matter models [11][12][13].
Here, and in our previous papers [14] (from now on referred to as paper I) and [15] (paper II), we develop a formalism that improves the calculation of the DM annihilation cross section in the general minimal supersymmetric standard model (MSSM) with neutralino DM by including the Sommerfeld corrections. By allowing the lightest supersymmetric particle (LSP) to be an arbitrary admixture of the electroweak gauginos and higgsinos, the framework can be used to study regions of the MSSM parameter space where the Sommerfeld correction may no longer be an order one effect, but still constitute the dominant radiative correction. Our method builds upon the non-relativistic nature of the pair of annihilating particles and separates the short-distance annihilation process (taking place at distances O(1/m LSP )) from the long-distance interactions characterized by the Bohr radius of order 1/(m LSP g 2 ) or 1/m W , responsible for the Sommerfeld effect, in analogy to the NRQCD treatment of quarkonium annihilation. However, in the MSSM co-annihilation effects of the LSP with heavier neutralino and chargino species have to be accounted for in regions where mass degeneracies are generic. Dealing with many nearly mass-degenerate scattering states (channels) requires an extension of the conventional NRQCD setup, which is provided in this work to accommodate any number of neutralino and chargino species. Previous works on the Sommerfeld enhancement focused on the wino and higgsino limits of the MSSM with at most three neutralino-• Our multi-channel approach implies that an incoming state χ i χ j can scatter to heavier two-particle states that are kinematically closed for the available centerof-mass (cms) energy. Eventually, when some mass-splittings become larger than M 2 W /M DM , numerical instabilities in the solution of the matrix-Schrödinger equation for the scattering wave functions prevents one from obtaining accurate results for the Sommerfeld factors. We circumvent this problem by reformulating the Schrödinger problem directly for the entries of the relevant matrix that provides the Sommerfeld factors, instead of solving for the wave functions. Leaving aside limitations related to the CPU time needed to solve a system of many coupled differential equations, this method allows to compute the Sommerfeld factors reliably also when many non-degenerate two-particle channels are present.
• The short-distance annihilation rates are obtained in the non-relativistic limit including corrections of O(v 2 rel ), where v rel = | v 1 − v 2 | denotes the relative velocity of the annihilating particles in their center-of-mass frame, i.e.
While the leading order describes only S-wave annihilations, the subleading term b contains both S-and P -wave contributions, that get multiplied by different Sommerfeld factors when building the full annihilation cross section from the rates above. A separation of the two partial-wave components of b is thus needed, which is automatically achieved in our EFT framework, since the Sommerfeld factors in question arise from different non-relativistic operator matrix elements.
In the EFT approach the short-distance contributions to the neutralino and chargino pair-annihilation processes are encoded into the Wilson coefficients of local four-fermion operators, whereas the absorptive part of the matrix element of these four-fermion operators gives the full annihilation rates, including the Sommerfeld corrections. In papers I and II we have written down the dimension-6 and dimension-8 four-fermion operators that mediate the short-distance annihilation rates at leading and next-to-next-to leading order in v rel , respectively, and presented the complete results for the corresponding Wilson coefficients in the general MSSM. In the present paper we turn to the longdistance part of the annihilation process and provide all the technical details required for the computation of the matrix elements of the four-fermion operators. Apart from the annihilation rates, the other model-specific ingredient, the non-relativistic potentials generated by electroweak gauge boson and Higgs boson exchange in the MSSM, are also provided in compact analytic form in this work. The contents of papers I-III therefore allow us to compute the full Sommerfeld-enhanced neutralino and chargino co-annihilation rates for an arbitrary MSSM parameter-space point and the corresponding relic abundance. A detailed investigation and discussion of Sommerfeld enhancements in the relic abundance calculation in some popular MSSM scenarios is the subject of an accompanying paper [16]. The general phenomenological study of Sommerfeld enhancements in the MSSM parameter space is postponed for a future publication.
A number of limitations needs to be mentioned. The formalism in its present form cannot be applied to DM annihilation through a narrow-width s-channel resonance, since in this case the annihilation process is no longer short-distance. Furthermore, we computed the pair annihilation and potentials of neutralinos and charginos, but not of neutralinos and sfermions, which excludes MSSM scenarios where neutralino-sfermion co-annihilation is important. Both, resonant annihilation and neutralino-sfermion coannihilation scenarios require the accidental coincidence of some MSSM parameters and are in this sense less generic than neutralino-chargino co-annihilations, especially when the LSP is heavy (m LSP ≫ m W ), and therefore a member of an approximate electroweak multiplet. On the technical side, we presently neglect the scale dependence of the electroweak couplings as well as electroweak double logarithms of the Sudakov type. The latter, of course, are only important when the formalism is applied to exclusive final states, as is relevant to indirect detection signals of DM. EFT methods have been proposed recently for summing up the Sudakov logarithms in multi-TeV DM annihilation [17][18][19]. More important to relic-density computations is the fact that for heavy dark matter the freeze-out occurs at temperatures where the temperature-dependence of the electroweak gauge-bosons masses can be relevant, as can be the temperature-dependence of the gaugino masses. While in general the thermal corrections to DM freeze-out are tiny [20], the case of Sommerfeld enhancements is special, since they depend sensitively on the range of the potential and the small mass splittings between the co-annihilating particles. In the wino and Higgsino limits some of the above effects have already been studied in the context of Sommerfeld enhancement [4,11]. For the general case, this is left for future work.
The contents of the paper are the following. The theoretical framework is summarized at some length in Sec. 2. We discuss the structure of the EFT Lagrangian and derive the master formula for the Sommerfeld-corrected cross section in Subsec. 2.1. The standard formalism to obtain the thermal relic abundance is reviewed in 2.2. The computation of the potentials for neutralino and chargino scattering in the MSSM is one of the main results of this paper, and we have devoted Sec. 3 to outline the details. The computation of the Sommerfeld factors is the subject of Sec. 4. In the first part of this section we relate the four-fermion matrix elements with the scattering wave functions that are determined by solving a multi-channel Schrödinger equation with MSSM potentials. Subsec. 4.2 describes the standard method to solve the Schrödinger equation and points out to the numerical instabilities that are caused by kinematically closed channels. The improved method that solves this problem is then explained in Subsec. 4.3. As in NRQCD, the matrix element of second-derivative operators can be related to the leading order ones, though here it becomes more involved due to the presence of Yukawa-like interactions generated by massive gauge bosons; the exact relation is derived in Subsec. 4.4. In the final part of Sec. 4 we give an approximation to the treatment of heavier channels which can be used to reduce the number of channels to be treated exactly in the multi-channel Schrödinger equation, and thus the time required for its numerical solution. A summary of our main results is given in Sec. 5. Finally, we collect the analytic expressions for the potentials in the MSSM in Appendix A, and we show the equivalence between the two different basis of two-particle states that can be used to evaluate the Sommerfeld factors in Appendix B. As mentioned above, the present papers concentrates on technical aspects of the framework, which is rather general and applicable to any model of heavy dark matter with electroweak interactions. A non-technical discussion of results is presented in an accompanying paper [16].

Theoretical framework
Together with previous results given in [14,15], the formalism that we present in this work allows to describe Sommerfeld-enhanced neutralino and chargino pair-annihilation rates within generic R-parity conserving MSSM scenarios, including the most general form of flavor mixing in the squark and slepton sector. After electroweak symmetry breaking the four neutralino states χ 0 i , i = 1, . . . , 4 result as combinations of the four SU(2) L × U(1) Y eigenstates bino, the electrically neutral wino and the two electrically neutral higgsinos. The chargino states χ ± j , j = 1, 2 are a mixture of the charged wino and higgsino electroweak eigenstates. While the bino and the winos are related to the soft SUSY breaking mass-parameters M 1 and M 2 , respectively, the higgsino states are associated with the mass-parameter µ in the underlying SUSY Lagrangian.
A MSSM scenario can be obtained with publicly available MSSM spectrum generators, for example [21][22][23], where the parameters M 1 , M 2 and µ among other required SUSY parameter inputs have to be specified. In constrained MSSM scenarios, as for instance models with grand unification of gauge couplings, certain relations among the input SUSY parameters are assumed. Our setup is not restricted to such cases, but allows to analyze Sommerfeld enhancements in χ 0 , χ ± co-annihilations in the most general MSSM models. Generically we require for our calculations a (SLHA formatted) MSSM spectrum including mass parameters, mixing matrices and angles, typically provided as output of a MSSM spectrum calculator. From this spectrum we determine the corresponding Sommerfeld-enhanced χ 0 , χ ± co-annihilation rates as well as the χ 0 1 relic abundance. Our formalism requires positive mass parameters of the neutralino and chargino states, which we automatically account for by an appropriate rotation of the neutralino and chargino mixing matrices, as explained in Appendix A of paper I [14].
A rigorous analysis of Sommerfeld-enhanced χ 0 , χ ± co-annihilation processes in a given model should refer to the on-shell mass spectrum of the neutralino and chargino states, instead of to the DR-parameters that are provided by most spectrum calculators. Furthermore, the mass splittings between the co-annihilating states play an essential role in the precise calculation of Sommerfeld enhancements, requiring one-loop on-shell renormalized masses in some cases. Results on the one-loop on-shell renormalized χ 0 , χ ± sector of the MSSM are available [24][25][26][27], but have not yet been implemented in public spectrum generators.
We perform our calculation in the framework of an effective field theory (EFT) of non-relativistic neutralinos and charginos, which generically covers the case of models with n 0 ≤ 4 neutralino-states χ 0 i , i = 1, . . . , n 0 and n + ≤ 2 chargino-states χ ± j , j = 1, . . . , n + , that are nearly mass-degenerate with the χ 0 1 . The description of our EFT framework is the purpose of Sec. 2.1, starting with the discussion of the Lagrangian in the effective theory in Sec. 2.1.1. In Sec. 2.1.2 we recall the notation for the four-fermion operators that reproduce the short-distance annihilation reactions and that were already introduced in [14,15]. The relevant formulae for the Sommerfeld-enhanced annihilation matrix elements, which determine the co-annihilation cross sections entering the χ 0 1 relic abundance calculation, are given in Sec. 2.1.3. Finally the χ 0 1 relic abundance calculation is summarized in Sec. 2.2. Our EFT framework does not cover scenarios where coannihilations with nearly mass-degenerate sfermion such as thet 1 orτ 1 are important in the relic abundance calculation. These cases require a straightforward extension of our EFT setup but are beyond the scope of this work. For the time being we exclude MSSM scenarios from our analysis, where other than χ 0 /χ ± co-annihilations are important in the χ 0 1 relic abundance calculation.

Lagrangian
In [14] we have introduced an effective field theory (EFT), the non-relativistic MSSM (NRMSSM), designed to describe the dynamics of charginos and neutralinos which are off-shell by an amount of the order of (m LSP v) 2 , where m LSP is the mass of the lightest neutralino and v its small velocity. The framework allows us to compute the inclusive annihilation rates of pairs of charginos and neutralinos moving at small velocities in a systematic expansion in the coupling constant and the velocity. The NRMSSM shares many similarities with the non-relativistic EFT of QCD [29], that has been employed for computing heavy quarkonium properties and radiative corrections to heavy quark production with remarkable accuracy, and can be considered as an extension of the latter in two aspects. First, the NRMSSM can account for several non-relativistic particle species, namely those neutralinos and charginos whose masses are nearly degenerate with m LSP , and second, it includes potential interactions generated by massive gauge bosons (i.e. Yukawa potentials) and not just Coulomb potentials. The structure of the EFT Lagrangian has already been discussed in [14]. It consists of three parts: where the dots stand for terms of higher order in the non-relativistic expansion, that are not required for the present calculation of the Sommerfeld-enhanced annihilation rates. L kin contains the bilinear terms in the two-component spinor fields ξ i and ψ j = η j , ζ j that represent the non-relativistic neutralinos (χ 0 i ) and charginos (χ − j and χ + j ), respectively. For n 0 ≤ 4 non-relativistic neutralino species and n + ≤ 2 non-relativistic chargino species, L kin is given by In order to have a consistent power-counting in the amplitudes describing transitions between two-particle states formed from the neutralino and chargino species included in the EFT we need that the mass differences (m i − m LSP ) are formally considered of order m LSP v 2 [14], 1 of the same order as the time-derivative and kinetic-energy term in the Lagrangian. This implies that heavier neutralinos and charginos (as well as further heavy SUSY particles and higher mass Higgs) are not among the degrees of freedom of the effective theory, and their virtual effects can only appear as short-distance corrections to the operators in L NRMSSM . In the same way, the hard modes associated to the SM and light Higgs-particle produced in neutralino and chargino pair-annihilations are encoded in the Wilson coefficients of four-fermion operators in δL ann , which are local in spacetime because the annihilation takes place at short-distances of O(1/m LSP ), as compared to the characteristic range O(1/m LSP v) or O(1/m W ) (whichever is smaller) of the nonrelativistic interactions between the charginos and neutralinos. The explicit form of the operators in δL ann , which are relevant for this work, and the details on the associated Wilson coefficients are given below in Sec. 2.1.2. The term L pot accounts for the exchange of SM gauge bosons and Higgs particles between the two-particle states χ e 1 χ e 2 and χ e 4 χ e 3 with non-relativistic relative velocity. In the non-relativistic limit, such interactions become instantaneous but spatially nonlocal, and are described in the EFT by four-fermion operators whose matching coefficients neutral reactions single-charged reactions double-charged reactions Table 1: Possible χ e 1 χ e 2 → χ e 4 χ e 3 scattering reactions classified according to the total charge. The labels e i on the fields χ e i are suppressed in the above table. If χ e i represents a field χ 0 e i , the label e i can range over e i = 1, . . . , n 0 , whereas e i = 1, . . . , n + for the case of a χ ± e i field. Redundant reactions where the type of particle in position 1 and 2 and/or 3 and 4 are interchanged (e.g. the reaction χ + χ − → χ 0 χ 0 ) are not explicitly written.
are Yukawa-and Coulomb potentials depending on the relative distance r = x ′ − x (r ≡ | r |) in the two-body system: The sum in (4) is taken over all χ e 1 χ e 2 → χ e 4 χ e 3 neutral, single-charged and doublecharged reactions with χ e i = χ 0 e i , χ ± e i , which have been summarized in Table 1. The particular particle species (χ 0 or χ ± ) participating in the reaction is indicated by the label χχ → χχ in the potentials. The explicit form of the terms in (4) are generated by replacing the generic fields χ e i by the field symbols ξ e i , η e i or ζ e i , corresponding to particle species χ 0 e i , χ − e i and χ + e i in all possible ways compatible with charge-conservation in the reaction. At leading order in the non-relativistic expansion, the potentials depend only on the total spin of the two-particle states, which is thus conserved. An individual potential contribution from the exchange of a gauge boson or Higgs particle with mass m φ has, at leading order, the generic form where we have written explicitly the spin indices α i omitted in (4), which are contracted with the corresponding indices in the field operators χ e i . The total spin operator S is built from the spin operators acting on the particles interacting at points x 1 and Since we shall decompose the two-particle states χ e 1 χ e 2 and χ e 4 χ e 3 undergoing potential interactions into 2S+1 L J partial-wave states with defined spin S = 0, 1, we can drop the spin indices in the potentials in what follows and replace the S 2 operator acting on χ e 1 χ e 2 ,χ e 4 χ e 3 by its eigenvalue S(S + 1) = 2S for S = 0, 1.
In this work we account for O(v 2 ) effects in the (co-)annihilation of neutralino and chargino pairs coming from the short-distance part of the annihilation but ignore O(v 2 ) contributions from the long-range part. Consequently, only the leading-order Coulomb-and Yukawa potential interactions need to be considered in L pot . Details on the calculation of the potentials at O(g 2 2 ) in the MSSM, where g i are the generic SU(2) L ⊗ U(1) Y gauge couplings, will be given in Sec. 3.

Annihilation matrices
The short-distance annihilation of the chargino and neutralino pairs into SM and light Higgs final states is reproduced in the EFT by local four-fermion operators contained in δL ann . The Wilson coefficients of these operators can be determined by matching the MSSM amplitudes for the process χ e1 χ e 2 → X → χ e 4 χ e 3 with SM and light Higgs intermediate states to the tree-level matrix element of the four-fermion operators for the same incoming and outgoing states. For the computation of the neutralino and chargino inclusive annihilation rates, only the absorptive part of these Wilson coefficients are required, see (13), and consequently the matching can be done for the absorptive part of the amplitude only. At lowest order in the electroweak gauge couplings g i , the contributions to the Wilson coefficients arise from the absorptive part of χ e1 χ e 2 → χ e 4 χ e 3 one-loop scattering diagrams with two SM or light Higgs particles in the intermediate state, Note that the annihilation rates include the absorptive parts of off-diagonal amplitudes in the space of two particle states χ e χ e ′ , since the exchange of gauge and Higgs bosons before the short-distance annihilation may transform one two-particle state into another.
The leading-order contributions to δL ann are given by dimension-6 four-fermion operators, that describe leading-order S-wave neutralino and chargino scattering processes χ e1 χ e 2 → χ e 4 χ e 3 . They read [14] δL d=6 ann = χχ→χχ S=0,1 where f χχ→χχ {e 1 e 2 }{e 4 e 3 } 2S+1 S J are the corresponding Wilson coefficients, which will be often abbreviated as f 2S+1 S J . The explicit form of the dimension-6 S-wave operators with O χχ→χχ The first sum in (6) is taken over all neutralino and chargino scattering reactions χ e 1 χ e 2 → χ e 4 χ e 3 specified in Tab. 1, including redundant ones where the particle species at the first and second and/or third and fourth position are interchanged. This redundancy implies that several operators describe one specific process with a χ e 1 and χ e 2 (χ e 3 and χ e 4 ) particle in the initial (final) state, and as a consequence their respective Wilson coefficients obey certain symmetry relations under the exchange of the labels e 1 ↔ e 2 and/or e 3 ↔ e 4 [14], see (12) below. The absorptive parts of the Wilson coefficients, f 2S+1 S J , have been calculated in [14] in the MSSM at O(α 2 i ). A master formula and necessary ingredients to obtain the contributions tof 2S+1 S J from individual states X A X B in analytic form can be found therein.
At O(v 2 ) in the non-relativistic expansion in momenta and mass differences, dimension-8 four-fermion operators contribute to δL ann : 2 The operators O ( 1 P 1 ) and O ( 3 P J ), J = 0, 1, 2, contain one derivative acting on each of the initial and final bilinear operators (χ c † e 2 χ e 1 and χ † e 4 χ c e 3 , respectively), while in the P 2S+1 S S operators the two derivatives act either on the initial or in the final state. The explicit form of these operators can be read off from Tab. 1 in [15]. The remaining next-to-next-to leading S-wave operators Q i 2S+1 S S , i = 1, 2, are the same as the dimension-6 operators O 2S+1 S S written in (8) with m e k the masses of the χ e k particles involved in the reaction χ e 1 χ e 2 → χ e 4 χ e 3 , and The mass scale M and mass differences δm, δm are process-specific quantities, their value being determined by the curly bracket labels of the short-distance coefficients f , g, h i multiplying them. Since δm = δm = 0 for diagonal annihilation reactions χ e 1 χ e 2 → χ e 1 χ e 2 (where the absorptive part of the amplitudes is related to the corresponding annihilation cross section), the Q i 2S+1 S S are only relevant for the computation of the off-diagonal rates. The convergence of the non-relativistic expansion for these offdiagonal terms requires that the mass differences are considered as O(v 2 ) effects [14]. 2 In a parity-violating theory such as the MSSM, there exist also O(v) contributions to δL ann , corresponding to dimension-7 four-fermion operators which describe 1 S 0 − 3 P 0 , 3 S 1 − 1 P 1 and 3 S 1 − 3 P 1 transitions where the spin and/or the orbital-angular momentum of the χχ pair is changed. They are not considered here because they contribute to the annihilation rates only in conjunction with O(v)-suppressed potential interactions in L pot , which we neglect. For the same reason, dimension-8 four-fermion operators in δL ann causing 3 P 1 → 1 P 1 and 3 S 1 → 3 D 1 transitions are ignored. Therefore, in scenarios where the off-diagonal annihilation terms can be relevant, the nondegeneracies among the particle species whose long-distance interactions are described within the EFT formalism are limited to δm ≪ m LSP . Particles with masses such that δm ∼ m LSP or larger should be decoupled explicitly and integrated out -they cause small modifications of short-distance coefficients in the effective Lagrangian, which are not relevant to the relic density computation.
Analytic results for the absorptive parts of the Wilson coefficients appearing in δL d=8 ann in the general MSSM can be extracted from the expressions given in [15]. For the Pand O(v 2 ) S-wave Wilson coefficients in (9) there are symmetry relations under the exchange of the particle labels analogous to those for the leading-order S-wave Wilson coefficients [15]. We reproduce them here for later reference: where k refers to any of the Wilson coefficients, f, g and h i in (6) and (9).

Sommerfeld-corrected cross section
The spin-averaged center-of-mass frame χ i χ j annihilation cross section summed over all accessible light final states is given by the imaginary part of the forward-scattering amplitude χ i χ j → χ i χ j by virtue of unitarity, see Fig. 1. In the non-relativistic effective theory, including up to O( p 2 ) corrections, this observable is obtained as with v rel = | v i − v j | the relative velocity of the annihilating particles in the center-of-mass frame, 3 and assuming the non-relativistic normalization p| p ′ = (2π) 3 δ (3) ( p − p ′ ) for the incoming chargino and neutralino one-particle states. In order to make the notation simpler we have omitted in (13) the sum symbol over the intermediate states χ e 1 χ e 2 and χ e 4 χ e 3 that undergo annihilation and the labels in the Wilson coefficients and operators (as well as in the quantities M, δm and δm). This means that, for instance, the first term in (13) in full form reads and the sum extends over all χ e 1 χ e 2 → χ e 4 χ e 3 annihilation reactions where the states χ e 1 χ e 2 , χ e 4 χ e 3 have the same charge as the incoming χ i χ j pair. In what follows we will always omit the symbol χχ→χχ , and when repeated indices e i appear in an expression a summation over the particle species will be implied. To obtain the last equality in (13) we have used that following from the definition of the absorptive part of the Wilson coefficientsf [14], which implies 4 Note also that we have related the spin-average of the matrix element of the operators O( 3 P J ), J = 1, 2, with that of O( 3 P 0 ) in the last line in (13). The matrix-elements of four-fermion operators in (13) account for the long-distance interactions between the annihilating pair, while the short-distance annihilation into light particles is described by the Wilson coefficients. The matching calculation of the absorptive part of the Wilson coefficients can be performed for exclusive two-particle final states X A X B at the tree-level [14] (i.e. at O(α 2 i )), since infrared divergences are absent at that order. Therefore (13) also applies separately for every final state X A X B to yield the exclusive annihilation rates σ χ i χ j →X A X B v rel with O(α 2 i ) short-distance corrections. It is well-known from quarkonium physics that matrix-elements of four-fermion operators analogous to those in (13) can be expressed in terms of non-relativistic wave functions and their derivatives evaluated at the origin. This relation becomes clear, if we explicitly insert the operator |0 0| that projects onto the Fock space with no neutralino and chargino states into the four-fermion operators, which is exact at the level corresponding forward scattering amplitudes; the former contain an additional factor (2π) 4 δ 4 (p f − p i ) which should not be included in the relation (13). 4 It follows directly from the definition off thatf χχ→χχ In papers I and II this relation was incorrectly written for the Wilson coefficients themselves and not for their absorptive parts; we take the opportunity to correct it here. of terms included here in L kin + L pot . For instance, the matrix element of the operator O χe 1 χe 2 →χe 4 χe 3 ( 1 S 0 ) can be written as where ψ (L,S) e 1 e 2 , ij is the χ e 1 χ e 2 -component of the scattering wave function for an incoming χ i χ j state with center-of-mass energy √ s, orbital quantum number L and total spin S, evaluated for zero relative distance and normalized to the free scattering solution.
The symbols ξ i , ξ j in the second line of (18) denote the Pauli spinor of the incoming particles χ i and χ j , and . . . stands for the trace in spin space. The multi-component wave function ψ (L,S) ij accounts for the potential interactions of the incoming χ i χ j state with all possible intermediate two-body chargino and neutralino states with the same charge and identical spin and partial-wave configuration. Recall from Section 2.1.1 that we consider only leading-order potential interactions in L pot , which cannot change the spin and orbital angular momentum of the two-particle states. Both wave-function components e 1 e 2 and e 2 e 1 are generated by the matrix-element of operator χ c † e 2 χ e 1 ; for an operator with quantum numbers L and S, there is a relative sign (−1) L+S between the two components, see (61). The precise definition of ψ (L,S) ij is postponed to Section 4. Let us just mention here that the lowest-order perturbative result for the matrix-elements of four-fermion operators is obtained by replacing ψ (L,S) eae b , ij → δ eai δ e b j , as can be easily checked by explicit computation.
In terms of the non-relativistic wave functions the spin-averaged annihilation cross section (13) takes the form where we have used the symmetry relations (12) for the Wilson coefficients, and the spin sums For a given two-particle state, the relative momentum of particles χ i , χ j in their centerof-mass frame is related to the center-of-mass energy of the collision by where M ij and µ ij are the total and reduced mass, respectively, of the two-particle system. In addition we have used the relation between the matrix-elements of operators O( 2S+1 S S ) and P( 2S+1 S S ), wherê and In (22), (23) the value of M 2 is determined by the particle labels of the short-distance coefficient in the numerator of the fraction. The sum in the second term in (23) extends over all potential interactions χ e 1 χ e 1 → χ e ′ 1 χ e ′ 2 arising from φ a -boson exchange. The coefficients of the potentials, c , are given in Tab. 2 of Appendix A. The derivation of (21) is postponed to Sec. 4. It is interesting to note that when we reduce the spectrum of two-particle states to just one state and the exchanged bosons are massless, the relation (21) between the leading and the v 2 -suppressed S-wave operators acquires the simpler form familiar from the NRQCD applications to heavy quarkonium, We define for an incoming state χ i χ j with center-of-mass energy √ s the Sommerfeld enhancement factor associated to a generic Wilson coefficient or combination of Wilson coefficients,f , describing the short-distance annihilation of χχ states with spin S and orbital-momentum L as the ratio The subscript "LO" in the denominator of (24) means that only the leading order in α 2 of the Wilson coefficientf ( 2S+1 L J ) should be kept in the denominator. For our purposes, 5 this is only relevant for the case of S ij [ĝ κ ( 2S+1 S S )], where we have to set the α 2 term in κ to zero, so thatĝ χχ→χχ The Sommerfeld factors are functions of √ s or, equivalently, of the relative velocity v rel of the incoming state. They allow us to parametrize the long-distance corrections to the annihilation rate of the state χ i χ j in a convenient way. Indeed, in terms of the Sommerfeld factors, the spin-averaged annihilation cross section (13) acquires the simple form where we have introduced the combinations of Wilson coefficientŝ The last relation in (26) reminds us that the knowledge of the spin-weighted sum over the three different 3 P 0 , 3 P 1 and 3 P 2 partial-wave Wilson coefficients is sufficient, since the leading-order potential interactions while being spin-dependent, do not discriminate among the three spin-1 P -wave states 3 P J with different total angular momentum J = 0, 1, 2. The pure tree-level annihilation rate with no long-distance corrections is readily recovered by setting all the Sommerfeld factors in (25) to one. The tree-level annihilation cross section thus obtained depends only on the diagonal entry of the Wilson coefficients corresponding to channel χ i χ j . Note that we have written in (25) f

Relic abundance
The thermal relic abundance of neutralino dark matter can be obtained by solving the Boltzmann equation that describes the time evolution of n = N i=1 n i , the sum of the particle number densities of all supersymmetric particles χ i taking part in the relevant annihilation reactions χ i χ j → X, which change the lightest neutralino number density n 1 , either directly or indirectly through the later decay into χ 1 of the species involved in those reactions. Note that one can solve (27) for n to obtain n 1 in the present Universe because R-parity conservation implies that all other χ i must decay into the lightest one by today. In this work the supersymmetric particles χ i considered in the relic density calculation are the neutralinos and charginos. For SUSY models where other particle species (staus for example) may have important co-annihilation effects with the neutralino, an extension of the EFT framework presented here to a larger set of two-particle states would be needed. The co-annihilation rates enter the Boltzmann equation through the thermal average of the effective cross section, σ eff v , whose specific form is given below. The other quantities entering (27) are the Hubble parameter H and n eq , the sum of the equilibrium number densities of each particle n i, eq . Eq. (27) is derived using Maxwell-Boltzmann statistics for all species in thermal equilibrium instead of Fermi-Dirac for fermions and Bose-Einstein for bosons, which is a very good approximation for the temperatures relevant for the relic density calculation of heavy dark matter (T T f with T f ∼ m χ 0 /20 the typical temperature where the departure from equilibrium takes place). Furthermore, the evolution equation (27) is valid under the assumption that the particle's phase-space distributions are proportional to those in equilibrium, f i , with a factor of proportionality that depends on T but not on the energy. As argued in [30], this is true if the χ i particles are maintained in kinetic equilibrium through elastic scattering with the much more abundant standard model particles in the thermal plasma during all of their evolution, even after they leave chemical equilibrium. In this work we assume that dark matter is kept in kinetic equilibrium after chemical decoupling sufficiently long until freeze-out is completed. The possibility that the elastic scattering processes stop to be effective and the consequences of early kinetic decoupling in the context of Sommerfeld-enhancement have been investigated in [31][32][33][34].
The quantity σ eff v that feeds into the equation for the dark matter relic density reads where σ ij v ij is the thermal average of the co-annihilation cross section of χ i χ j into standard model particles and v ij is the so-called Møller velocity of particles χ i and χ j . 6 In the cosmic comoving frame, i.e. the frame where the plasma is at rest as a whole, the thermal average σ ij v ij can be calculated using f i = e −E i /T for the Maxwell-Boltzmann equilibrium distributions. A general expression for σ eff v written in terms of Lorentz invariants was first derived in [35] and further generalized to include co-annihilations in [36]. It reads with g i the internal degrees of freedom of particle species χ i , which for neutralinos and charginos is g i = 2, and K n the modified Bessel function of the second kind of order n. The integrand in (29) evaluated in the center-of-mass frame allow us to use the expressions obtained in Sec.
Note that the thermal average of a heavy co-annihilation channel is typically suppressed by a factor e −(m i +m j −2m 1 )/T with respect to σ 11 v 11 , which arises from the asymptotic expansion of the Bessel function The efficiency of dark matter production and annihilation processes to maintain chemical equilibrium with the plasma in an expanding Universe can be better understood when the Boltzmann equation is written in terms of yield Y ≡ n/s, defined as the ratio of the particle density to the entropy density in the comoving cosmic frame. This is because, assuming that the entropy per comoving volume is conserved, the change of n and s due to the expansion of the Universe is the same, namely ds/dt = −3Hs, and gets scaled out from the Boltzmann equation. Since σ eff v is computed as a function of the temperature rather than time, it is also convenient to trade the independent variable t in the evolution equation by T . In terms of the yield, and defining x = m 1 /T , the Boltzmann equation (27) takes the form where Y eq = n eq /s. The final step requires the Friedmann equation that relates the expansion of the Universe (Hubble rate H) to its energy density. In a radiation-dominated universe at early times, the equation for the evolution of Y can be finally written as [35] with G the gravitational constant. The parameter g 1/2 * is defined as in terms of the effective degrees of freedom g eff and h eff of the energy and entropy densities: The effective degrees of freedom are slowly-varying functions of the temperature except at the QCD quark-hadron phase transition (T ∼ 150 − 400 MeV), where there is a sharp increase [35,37]. The equilibrium yield Y eq can be obtained using the formula for n eq given in (30) and the parametrization of the entropy density in (34) above.
In the form (32), we see that the change in the yield tries to compensate for deviations from the equilibrium, such that if Y was initially close to Y eq at high temperatures (low x) it will continue to track Y eq as the temperature decreases. At a given temperature the prefactor σ eff v /x 2 (take g 1/2 * constant) is not large enough to provide the change in Y needed to keep up with the rapid exponential drop of Y eq , and chemical decoupling of the massive particle from the plasma takes place. The yield then remains constant as the universe cools down (phenomenon known as "sudden freeze-out") or, if the annihilation cross section is further increased for lower temperatures the yield continues to be reduced; that is the case when the annihilation rates get enhanced by Sommerfeld corrections.
In order to obtain the relic density we solve (32) numerically from x = 1, where one can safely assume that the yield still tracks the equilibrium value and take Y = Y eq as initial condition, up to x 0 = m 1 /T 0 , where T 0 is the current photon temperature in the Universe. Once the present value of the yield, Y 0 ≡ Y (x 0 ), is known, the neutralino relic density today in units of the critical density is determined by the relation with ρ crit = 3H 2 /8πG the critical density and s 0 the entropy density of the present universe, that can be calculated from (34). It is customary to provide the relic density as Ωh 2 , where h is the value of the Hubble parameter determined at present in units of 100 km s −1 Mpc −1 .
In the presence of Sommerfeld corrections the effective annihilation cross section σ eff v acquires a complicated velocity dependence that can no longer be parameterized in the simple form a + bv 2 . For the results that we present in [16], Y 0 has been obtained by solving (27) using numerical routines for differential equations built in Mathematica. Due to the stiffness of the evolution equation, the numerical integration of the equation is done by splitting the region from x = 1 to x = 10 8 into several pieces and adapting the starting and maximum step sizes in each of them. For g 1/2 * (T ) and h eff (T ) we have used the values derived in [35], which can be found conveniently tabulated as a function of temperature among the package files of the automated programs DarkSUSY [38] and micrOMEGAs [39,40]. Other numerical values needed for the computation of the relic density are T 0 = 2.7255 K and ρ crit = 1.05368 × 10 −5 h 2 GeV cm −3 , both taken from [41].

MSSM potentials
Sommerfeld enhancement occurs when the interaction between the incoming DM particles significantly distorts their two-body wave function away from the initial plane wave. In terms of Feynman diagrams the effect arises from the exchange of electroweak gauge (and to a lesser extent, Higgs) bosons between the neutralinos and charginos. In the non-relativistic limit the dominant contribution arises from the potential loop momentum region in ladder box graphs, which have to be summed to all orders in α EW , as is shown in Sec. 4 below. A necessary input to perform such resummation are the potential interactions generated by electroweak-and Higgs-boson exchange, which are given in this section.
In the center-of-mass system of the incoming χ i χ j pair, p i + p j = P = ( √ s, 0 ), the potential loop momentum p ′ running inside a box graph with internal fermions χ e 4 χ e 3 can be chosen to be the relative momentum of the χ e 4 χ e 3 pair, i.e. p ′ = (p 4 − p 3 )/2; the momentum carried by the exchanged boson is then equal to the difference between the relative momenta in the χχ pair before and after the interaction (see Fig. 2). The potential loop-momentum routed in this way is characterized by the scaling p ′ 0 ∼ p ′ 2 /m LSP ≪ m LSP (idem for p), and fermion and boson propagators can be expanded accordingly. To leading order in the potential-region expansion, the denominator of the boson propagator, and thus represents an instantaneous interaction between the χ e 1 χ e 2 and χ e 3 χ e 4 pairs, that in the non-relativistic EFT is taken into account by the potentials in L pot . We note that the (p ′ 0 − p 0 ) 2 term dropped in the boson propagator has a term proportional to the mass difference squared (m 4 − m 1 ) 2 , which we neglect consistently given that we are dealing with a set of nearly mass-degenerate with the neutralino LSP and mass differences are assumed to be at most of order of The leading-order potential interactions in the EFT can be obtained by expanding the full-theory tree-level scattering amplitude χ e 1 χ e 2 → χ e 4 χ e 3 for on-shell charginos or neutralinos exchanging a gauge or Higgs boson. Three tree diagrams which differ in the direction of the fermion flow have to be considered depending on the particle nature of the scattering particles, see Fig. 3. By means of an example, Z-boson exchange, we illustrate in the following the essential steps to perform this tree-level matching.
Let us denote the vector (v), axial-vector (a) interaction vertices of charginos and neutralinos with gauge bosons in the MSSM generically as where the spin-1 field A V µ stands for either the Z-boson (then χ e i χ e j = χ 0 e i χ 0 e j , χ + e i χ + e j ), the photon field (χ e i χ e j = χ + e i χ + e j ), or the W + -field (which sets χ e i χ e j = χ + e i χ 0 e j ) 7 . In the latter case, the hermitian conjugate of (36) has to be added to describe the chargeconjugated interaction, i.e.
Likewise, the interaction vertices of charginos/neutralinos with Higgs particles has the generic form where the scalar field φ can be any of the CP -even (H 0 Note, however, that our results apply as well to the CP -violating MSSM where the neutral Higgs particles are no longer CP eigenstates. For charged Higgs-boson exchange the hermitian conjugate of (38) has to be considered as well, which introduces the hermitian conjugates of the coupling matrices, (s H + m † ) ij and (p H + m † ) ij . Explicit expressions for all the v ij , a ij , s ij , p ij couplings are given in the Appendix A.
The amplitude of Fig. 3a with Z-boson exchange is obtained at leading order in the non-relativistic expansion by taking the limit of small relative momenta, p, p ′ ∼ m LSP v rel ∼ M Z and p 0 , p ′ 0 ∼ m LSP v 2 rel , in the particle spinors and Z-boson propagator, here defined in the R ξ -gauge. At leading-order, the non-relativistic expansion of the product of Dirac bilinears in the amplitude of diagram 3a,ū( Figure 3: Tree-level diagrams with t-channel boson exchange that generate the leadingorder potential among non-relativistic neutralinos and charginos. The arrows in the neutralino/chargino lines indicate the fermion flow, such that each diagram contributes only to the scattering processes written below. in the full-theory amplitude, where the right-hand-side of these relations should be understood as the matrices acting on the two-component Pauli spinors of the non-relativistic neutralinos and charginos at the upper and lower interaction vertices of Fig. 3a. Written in components, the spin operators above read σ i ⊗σ i ≡ σ i α 4 α 1 σ i α 3 α 2 and 1⊗1 ≡ δ α 4 α 1 δ α 3 α 2 . The use of the non-relativistic normalization for the relativistic spinors, u † (p)u(p) = 1, is implied in the replacements (40). The relation in the third line of (40) can be obtained using the equation of motion for the relativistic spinors; we comment on the non-relativistic scaling of this term below. Using (40) we thus obtain for the leading-order term in the expansion in p, Fig. 3a when the exchanged particle is a Z-boson. For simplicity we have considered the case where the couplings in both vertices are equal, which applies to χ 0 χ 0 → χ 0 χ 0 and χ + χ + → χ + χ + processes.
The term in the second line of (41) arises from the ξ-dependent term of the Zboson propagator (39). Gauge invariance requires that this term is cancelled against the contribution from the exchange of the Goldstone boson A 0 2 ≡ G 0 , which has a mass ξM 2 Z . The cancellation holds because certain conditions among the vector and scalar couplings of the neutralinos and charginos are fulfilled in the MSSM. To see this, write the leading-order contribution corresponding to diagram 3a with scalar-boson φ exchange by expanding the scalar propagator, i/(q 2 − m 2 φ ), and the Dirac bilinear in the non-relativistic limit; the result reads The pseudoscalar interactions do not survive in (42) (41) one obtains the condition for the cancellation of the ξ-dependent terms. Eq.
ij , a relation which can be proven using the explicit definition of the couplings in terms of the mixing matrices given in Appendix A and their diagonalization properties. We can therefore drop the second line in (41) together with the pseudo-Goldstone contribution to the potential.
The term proportional to the mass differences in the first line of the potential (41) typically yields a very small contribution, since in the EFT we treat only those species for which δm ∼ O(m LSP v 2 ). For a very heavy LSP, where we could nevertheless have m LSP v 2 ≫ M Z , the size of this term would still be at most of O(1) due to suppressed vector couplings. To see this, consider the decoupling limit m LSP → ∞. If the mass difference δm ij refers to particles within the same electroweak multiplet then δm ij ∼ M 2 EW /m LSP and the mass-difference terms in (41) are suppressed as δm ij /M Z ∼ M EW /m LSP . If particles i and j belong to different multiplets, δm ij /M Z can be large, but is multiplied by v Z ij ∼ M EW /δm ij , as follows (43). Since the axial couplings go to zero in the decoupling limit, it follows that the gauge-independent, off-diagonal terms in square brackets in the first line of (41) are given in this limit by the pseudo-Goldstone couplings −s G 0 e 4 e 1 s G 0 e 3 e 2 . The corresponding results for diagrams b and c in Fig. 3 can then be easily obtained from (41): for every line where the fermion flow has been reversed we need to interchange the labels in the vertex couplings (for instance, X e 3 e 2 → X e 2 e 3 to go from a to b), change the sign of the scalar, pseudo-scalar and axial-vector couplings (which arises when writing the antiparticle spinors as the charge-conjugate of particle ones), and add a global sign due to Wick ordering. For the sake of clarity, let us write the results for Z-and neutral Higgs-boson exchange for diagram b: The potentials generated by photon exchange can be obtained from the result for the Z-boson case by plugging the corresponding vector coupling and keeping only the g µν part of (39) with M Z → 0 for the photon propagator. W and charged Higgs-boson exchange are relevant for the processes χ 0 χ 0 → χ ± χ ∓ and χ 0 χ ± → χ ± χ 0 . For these amplitudes, we have to hermitian conjugate the coupling matrices at vertices involving either an initial χ + or a final χ − . The potential from charged boson exchange can be obtained from the neutral boson amplitudes with trivial replacements in the couplings and the boson mass.
The amplitudes above provide the potentials in momentum space. However, the Schrödinger equation which sums the all-order exchange of bosons among the chargino and neutralino pairs acquires a simpler form in coordinate space. The coordinate-space potentials are obtained by taking the Fourier transform where r ≡ | x |, and T χχ→χχ e 1 e 2 e 4 e 3 stands for the momentum-space amplitude as given above in (41,42,44). From the identity we immediately obtain the well-known Yukawa-like potential for amplitudes with exchange of a force carrier of mass m. Applied to the massless case, (46) gives the Coulomb potential. 8 Before we write the result for the potentials in coordinate space, let us write the spin operator σ i ⊗ σ i in terms of the total spin operator as where S = s 1 + s 2 ≡ 1/2 ( σ ⊗ 1 + 1 ⊗ σ), with s 1,2 the spin operators acting on the particles at the upper and lower vertices in the scattering diagram, and we have replaced s 2 1,2 by s(s + 1) (1 ⊗ 1) = 3/4 (1 ⊗ 1) for charginos and neutralinos. Working in the basis of eigenstates of total spin for the neutralino and chargino pairs, we can as well replace S 2 by S(S + 1) (1 ⊗ 1) = 2S (1 ⊗ 1) for S = 0, 1 in the potentials. Since the operator 1 ⊗ 1 is the identity operator, we find that a leading-order potential contribution in the basis of total spin and in coordinate space thus has the form e −m X r r , for the case of the exchange of a vector boson with mass m X , among the incoming and outgoing χχ pairs. For leading-order scalar boson and photon exchange, the coefficient b e 1 e 2 e 4 e 3 vanishes. For Z-and W -boson exchange the spin-dependent part of the potential arises from the axial-vector coupling, see (41), (44).
We show in Sec. 4 by analysing the perturbative expansion of the annihilation amplitude χ i χ j → X A X B with ladder-like exchanges in the non-relativistic limit, that the amplitude can be described as a product of potential interactions like (48) times the two-particle propagator of the internal pairs χ m χ n in the diagrams, integrated over corresponding loop momenta, and finally multiplied by the short-distance coefficientf with the appropriate quantum numbers. No additional combinatoric factors arise as long as for the internal states we consider as different the states χ e 1 χ e 2 and χ e 2 χ e 1 and sum over all of them. Diagrammatically this means for instance that we have to consider a graph like Fig. 2 and the one where particles χ e 4 , χ e 3 are interchanged. At this point, it is convenient to switch from the two indices {e 1 , e 2 } that denote a two-particle state χ e 1 χ e 2 to a single label m = 1, . . . N |Q| , where N |Q| is the total number of channels for each total charge sector, |Q| = 0, 1, 2. For the case of n 0 neutralinos and n + charginos, we have n 2 0 χ 0 χ 0 and 2n 2 + χ ± χ ∓ neutral states, 2n 0 n + different states with Q = ±1 (i.e. of the type χ 0 χ ± or χ ± χ 0 ), and n 2 + χ ± χ ± states of charge Q = ±2. Therefore, the total number of neutral states is N 0 = n 2 0 + 2n 2 + , whereas N 1 = 2n 0 n + and N 2 = n 2 + . If all four neutralinos and the two charginos are relevant for the long-distance part of the annihilation (n 0 = 4 and n + = 2), then we have to consider the interactions among N 0 = 24 states in the neutral sector, as well as those among N 1 = 16 and N 2 = 4 states in the singly-charged and doubly-charged sectors, respectively. For instance, in the charge-0 sector with all four neutralinos and two charginos, the single label runs over states with 24 different states in total, whereas in the charge-1 sector we have 16 channels, and just 4 in the charge-2 sector, The ordering of the states in each sector is of course a matter of convention. An alternative basis of two-particle states with a smaller number of states can be used, which takes into account the fact that states χ e 1 χ e 2 and χ e 2 χ e 1 have identical particle content. The new basis, that we shall refer to as "method-2 basis" in order to distinguish it from the basis just described above, is built only by states χ e 1 χ e 2 with e 1 ≤ e 2 , provided we identify χ i = χ 0 i for i = 1, . . . 4 and χ 5 = χ + 1 , The two potential contributions of method-1 corresponding to χ e 1 χ e 2 → χ e 4 χ e 3 and χ e 1 χ e 2 → χ e 3 χ e 4 (the so-called "crossed" contribution, different from the former if χ e 3 = χ e 4 ) are accounted for in method-2 by a single potential entry describing the scattering of state (χχ) e 1 e 2 into (χχ) e 4 e 3 , where we have introduced the notation (χχ) ... to denote a state in the new basis, and assumed that e 1 ≤ e 2 and e 4 ≤ e 3 . The potential entries for method-2, V The (−1) L+S in front of the crossed diagram amplitude in the first line of (52) arises from the product (−1) × (−1) L × (−1) S+1 , where the (−1) comes from Wick ordering, and the exchange e 3 ↔ e 4 produces a factor (−1) S+1 in the spin wave function and a change of sign in the relative momentum which translates into the factor (−1) L . An additional rule has to be accounted for when building the potential matrix for method-2: since identical spin-1/2 particles cannot form a 2-particle state with odd L + S (i.e. with quantum numbers 3 S 1 and 1 P 1 in the case at hand), the entries involving such states for the potential when L + S is odd have to be set to zero. This prevents that method-2 yields a non-zero annihilation amplitude for a forbidden state of two identical neutralinos or charginos through an intermediate transition to an allowed state, such as (χ 0 χ 0 ) 11 → (χ + χ − ) 11 → light particles. In method-1 this transition is automatically zero because there are two annihilation terms, χ + 1 χ − 1 → light and χ − 1 χ + 1 → light, which cancel each other by virtue of the symmetry properties of the Wilson coefficientsf under exchange of the particle labels, see (8) and (6) in [14] and [15], respectively. We note that the the relations (52) and the selection rule just mentioned imply that the potential in method-2 depends on the orbital angular momentum L even at the leading-order in the non-relativistic expansion.
Likewise, the annihilation matrices built from the Wilson coefficients of the fourfermion operators must be supplemented with additional factors in the method-2 basis. Namely, the absorptive part of the one-loop annihilation amplitude (χχ) e 1 e 2 → (χχ) e 4 e 3 involving method-2 states in a 2S+1 L J partial-wave configuration is given by the Wilson coefficientf χχ→χχ {e 1 e 2 }{e 4 e 3 } ( 2S+1 L J ) multiplied by (1/ √ 2) n id , where n id = 1, 2 if the twoparticle states (χχ) e 1 e 2 or/and (χχ) e 4 e 3 are formed of identical particles, and n id = 0 otherwise.
The basis of states in method-2 for the case of n 0 neutralinos and n + charginos, containsÑ 0 = n 0 (n 0 + 1)/2 + n 2 + neutral states,Ñ 1 = n 0 n + states with Q = ±1, and N 2 = n + (n + + 1)/2 states with charge Q = ±2. The reduction in the number of states in method-2 as compared to method-1 is more significant as more particle species are considered in the non-relativistic EFT framework. For instance, if all neutralino and chargino species are considered then (Ñ 0 ,Ñ 1 ,Ñ 2 ) = (14,8,3), whereas for method-1 (N 0 , N 1 , N 2 ) = (24,16,4). The list of method-2 states can be read off from (49,50,51) by dropping states that differ only in χ i χ j where i > j. In what follows, quantities with indices referring to two-particle states, such as the Wilson coefficients that built the annihilation matrices and the potentials, will be written using the single-label notation with lower-case letters, i.e.f mn ( 2S+1 L J ) and V mn (r), with m, n = 1, . . . N |Q| in the case of method-1, and m, n = 1, . . .Ñ |Q| for method-2. We shall also adopt the convention that the potential matrix element V mn describes the scattering m → n, whereas for the annihilation matrices, the entryf mn ( 2S+1 L J ) is the absorptive part of the one-loop amplitude for m → n. The formula for the Sommerfeld enhancement factor (24) using this notation reads for an incoming two-particle state a, and the same equation holds in both methods.
In Appendix B we demonstrate explicitly for a simplified example that method-1 and method-2 yield the same Sommerfeld factors.
To provide an example, let us give both the potential and the annihilation matrices in the well-known wino limit of the MSSM. The relevant particles in this limit are the purewino lightest neutralino χ 0 1 , with mass m χ , and its mass-degenerate chargino partners χ ± 1 . In the neutral sector, the potential matrix of method-1 for both S = 0, 1 derives from the formulae given in Appendix A after plugging in the values of the neutralino and chargino mixing matrix entries relevant to the calculation in the pure-wino, i.e. Z N i1 = δ i2 and Z ± i1 = δ i1 : where the matrix indices (m, n = 1, 2, 3) correspond to channels χ 0 1 χ 0 1 , χ + 1 χ − 1 , χ − 1 χ + 1 . The potential for method-2 then follows from the rules given in (52) together with the rule that sets to zero entries involving a forbidden 3 S 1 or 1 P 1 χ 0 1 χ 0 1 state: where the matrix indices m, n = 1, 2 refer now to channels χ 0 1 χ 0 1 , χ + 1 χ − 1 . The annihilation matrices for any MSSM model to O(v 2 rel ) can be obtained from the analytic formulae given in [14,15]. In the pure-wino limit, the explicit expressions for the absorptive parts of the Wilson coefficients have been worked out in detail in Appendix C of [15], and can be read off from (147) and Tables 5-7 therein. In the neutral sector, the leading order S-wave annihilation matrices read in method-1, and in method-2. The results (58) have been given before in [3]. Let us also write the P -wave annihilation matrices in both methods: Recall that the entries in the annihilation matrices of method-1 which differ only in the replacement χ + 1 χ − 1 ↔ χ − 1 χ + 1 as in or out state are equal up to a factor (−1) L+S , which arises as a consequence of the exchange of labels in the Wilson coefficients (see Eqs. (2.7) and (2.4) of [14,15], respectively). The potential and annihilation matrices for the other charge sectors can be obtained similarly.

Sommerfeld enhancement
In this section we discuss the computation of the Sommerfeld factors S ij [f ( 2S+1 L J )] defined in (24). Within the non-relativistic MSSM the Sommerfeld enhancement arises from the matrix elements of annihilation operators such as O χχ→χχ Here Γ = 1 2×2 , σ for S = 0, 1, respectively, and K is a polynomial in relative momentum (derivatives) corresponding to a given angular momentum L. p denotes the relative momentum of the χ i and χ j particle as defined below in (63). With this definition ψ (L,S) e 1 e 2 , ij is normalized to one, when the Sommerfeld effect is neglected, and the matrix element is evaluated in the tree approximation.
In the first part of this section, we establish the relation between the non-relativistic (NR) matrix element above, diagrammatic resummation and the solution of a multichannel Schrödinger equation. We then describe the actual solution of the Schrödinger equation. We shall find that the standard method fails to provide an accurate result in many relevant cases with kinematically closed channels and derive an alternative method that solves this problem. This is our main result, since it allows to compute the Sommerfeld factors for general MSSM parameter points. In practice, the solution with the full number of channels is time-consuming and often not relevant, since the heavier channels have only a small effect on the final result. In the final part of this section we therefore describe an approximation to the treatment of heavier channels, which is accurate and fast for most practical purposes.

NR matrix elements and the Schrödinger equation
The enhancement of non-relativistic scattering originates from the potential loop momentum region, when the momentum k of the exchanged virtual particle satisfies k 0 ≪ | k | ≪ µ, where µ is the reduced mass of the two non-relativistic particles. It can be shown by power-counting arguments, see for instance [42], that no other region requires resummation, and that the enhancement from the potential region is only present in ladder diagrams. Relating the sum of ladder diagrams shown in Fig. 4 to the solution of a Schrödinger equation involves a number of steps, which we now sketch. The first consists of kinematic simplifications of the propagators in the potential region. We already mentioned that the (k 0 ) 2 term can be dropped in the propagator of the exchanged particle, so that the interaction corresponds to a potential between the heavy particles. The pair of heavy-particle propagators in each ladder rung can also be simplified. Let be the momenta of the on-shell external state χ i χ j , and µ ij the reduced mass. Since | P |, | p | ≪ µ ij , the center-of-mass energy of the annihilation process is up to higher-order terms in the non-relativistic expansion. For later purposes it proves convenient to introduce the variable E, which measures energy from the common ref-erence value 2m LSP . Since mass differences m i + m j − 2m LSP scale as m LSP v 2 by assumption, both energy variables have the non-relativistic scaling E ∼ E ∼ m LSP v 2 . The masses of the two-particle state a in a given ladder rung are denoted by m a 1 , m a 2 , and M a ≡ m a 1 + m a 2 . The loop momentum energy-component integrals in each ladder rung can then be performed employing To arrive at this result we systematically neglected higher-order terms in the nonrelativistic expansion. As discussed before, we must require that the mass differences along a given heavy particle line are small, that is m a 1 − m i and m a 2 − m j should be of order m LSP v 2 , but we do not have to assume that m j −m i and m a 2 −m a 1 are small. Note than within these approximations the denominator of the kinetic energy term, µ a , could be equivalently substituted by µ ij or m LSP /2, or any other two-particle state reduced mass.
With these simplifications to the ladder-diagram sum, the annihilation matrix element, including the tree diagram with no exchange, can be written as where i refers to the initial two-particle state ij and e (ē) to the state e 1 e 2 (e 2 e 1 ) which annihilates into light particles. Below we employ this compact notation to label twoparticle states by compound indices a, b, . . . . The functionG is defined through and (For the term n = 0 we set a 0 = a and k 0 = 0.) The first term on the right-hand side of (66) accounts for the tree diagram, a term with given n in (67) for the (n + 1)loop ladder diagram with n + 1 exchanges. 9 The factorsV ab (k) contain the propagator of the exchanged particle and coupling factors from the two vertices. These are the momentum-space potentials discussed in Sec. 3. Note that using the on-shell condition for the external particles and non-relativistic approximations, the Dirac matrices from the numerator can be reduced to the tree structure ξ c † j Γξ i . Depending on whether Γ = 1 2×2 or σ, the potentials for S = 0 or S = 1 must be used in the above equations, which causes the dependence of ψ (L,S) e 1 e 2 , ij on the total spin S. The angular-momentum dependence arises through K[ q ] in (65), which reads K = 1 for L = 0 and K = q for L = 1. The limiting procedure in (65) is required since the factor [Ê − p 2 /(2µ ij )] vanishes forÊ = E, whilẽ G ab ( p, q; E) develops a singularity that corresponds to scattering states with relativemomentum kinetic-energy E.
The second step consists of verifying thatG ab ( p, q; E) is in fact the momentum-space Green function of a certain Schrödinger operator. Indeed, one easily checks that it satisfies the Lippmann-Schwinger equation hence the Fourier-transform G ab ( r, r ′ ; E) is the Green function for the Schrödinger equation with the coordinate-space potential Note that this result could be obtained directly as an equation of motion from the effective Lagrangian given in Sec. 2 by including the potential in the unperturbed Lagrangian. The third and final step uses elementary scattering theory to express (65) in terms of eigenfunctions with energy E rather than the full Green function. Let us define the Green operators for the interacting and free Hamiltonians of the Schrödinger problem and the corresponding momentum eigenstates Note that G, G 0 are matrix-valued operators with dimensionality related to the number of two-particle states (χχ) a = χ a 1 χ a 2 . Accordingly, the eigenstates carry a compound index a that refers to the component of the wave-function proportional to the twoparticle state a. The states | p+, a are the exact stationary scattering states of H with eigenvalue E = p 2 /(2µ a ), while | p , a are plane waves. We may represent the former in the plane-wave basis by [ψ E ( q )] ab = q , a| p+, b . The Lippmann-Schwinger equation can be used to derive | p+, a = GG −1 0 | p , a . With these preparations we use (61), (65) to compute The same result can also be derived from the spectral representation in terms of coordinate-space scattering wave-functions at the origin. Note that the result carries two compound indices referring to two-particle states. The second, ij, refers to the incoming two-particle state | p , i with kinetic energy E − (M i − 2m LSP ) in the cms frame of the annihilation. The first, e = e 1 e 2 , specifies that only the component proportional to the two-particle state e of the wave-function at the origin for this incoming state is picked out by the annihilation operator χ c † e 2 χ e 1 that defines ψ (L,S) e 1 e 2 , ij . The wave-function [ψ E ( r )] a,ij can be obtained directly from the matrix-Schrödinger equation with the potential (70). The label ij refers to the fact that this equation should be solved with the initial condition corresponding to a particular incoming two-particle state ij.
Let us finally mention that the complex conjugated scattering wave-function appears in (73) because of the convention used for the left and right states in the definition of the Green function (66). Using the opposite convention we would end up with [ψ E ( q )] ei in (73) but with V ba instead of V ab in the Schrödinger equation (76). 10

Solution of the multi-channel Schrödinger equation
We first solve the matrix-Schrödinger equation for the Sommerfeld factors by following closely the method described in [43]. Although, eventually, we will not use this method, we discuss it to set up the framework for the subsequent improvement. Under our assumption for the mass splittings, we can approximate µ a ≈ m LSP /2 in the kinetic term, since the difference is a v 2 correction, which we consistently neglect in the longdistance part. The Schrödinger equation now reads and we recall the definitions We also use the velocity variable v defined by This implies that the threshold for co-annihilation channels occurs at finite velocities v = ((m i + m j − 2m LSP )/m LSP ) 1/2 . With the above redefinitions the dependence on the initial scattering state appears only in the initial condition for the solution as indicated by the subscript i of [ψ E ( r )] bi , but not in the equation itself. We note that the MSSM matrix potential is hermitian, but in general not symmetric, since the entries of the potential matrix can be complex numbers. The potential is spherically symmetric. For large r, the potential approaches the diagonal matrix V inf with diagonal entries V inf,a = M a − 2m LSP . We therefore define the wave numbers which determine the asymptotic behaviour of the solutions. When E > V inf,a , channel a is kinematically open and k a is real. Otherwise the channel is closed and k a is imaginary. The scattering solution describing an incoming plane wave propagating along the z-direction and an outgoing scattered spherical wave takes the asymptotic form where, due to the azimuthal symmetry of the problem f ai (θ) does not depend on the azimuthal angle. 11 Since the potential is non-diagonal, the scattered spherical wave is not proportional to δ ai . The asymptotic behaviour (81) should be matched to the behaviour of the general solution where P L (cos θ) denotes the Legendre polynomials. This expresses [ψ E ( r )] ai as a superposition A bi of basis solutions [u L (r)] ab of the radial Schrödinger equation for each partial wave. When V is a N × N matrix, there exist 2N linearly independent solutions to (83) corresponding to the N independent initial conditions i, each of which is an N-component vector with index a. N solutions are irregular at the origin, hence restricting us to the set of N regular ones. The asymptotic behaviour of the ath component of the regular linear independent solution for initial state i is given by with constant coefficients n ai and scattering phases δ ai . 12 Matching the asymptotic behaviours of (81) and (82), we obtain potentials arise exclusively in the diagonal entries of the potential matrix V (r), when written in the two-particle mass-eigenstate basis. For large values of r, V (r) then always approaches a diagonal matrix (the non-diagonal Yukawa potentials being exponentially suppressed) with entries V inf,a + c a α/r, containing the Coulomb potential contributions as well as constant mass splitting terms. Accounting for the presence of long-range Coulomb potentials, the exp (ik a z) factor in (81) should be replaced by the incoming wave-function in presence of the Coulomb potential in the aa component of V (r), and in the outgoing scattered wave, one has to replace exp (ik a r) by exp (i(k a r + m LSP c a α/(2k a ) ln(2k a r))). Since it can be shown that the subsequent derivation of the Sommerfeld factor including Coulomb potentials on the diagonals of V (r) is completely analogous to the short-range potential case, and leads to the same result for the enhancement factor, we will for brevity refer to the asymptotic behaviour (81) in the following. 12 Both depend on L, but we do not indicate this explicitly. Eq. (84) assumes short-range potentials. In the presence of Coulomb potentials in the diagonal entries of the potential matrix V (r) into account, the asymptotic behaviour of the regular basis solutions reads [u L (r)] ai r→∞ = n ai sin k a r − Lπ 2 + m LSP c a α 2k a ln(2k a r)δ ai .
The modifications match those of the scattering solution, see footnote 11.
The basic ingredients (75) for the Sommerfeld factors now require the wave function and its derivatives at the origin. We suppose that the annihilation operator is constructed such that it overlaps only with a single partial wave L, the relevant cases in practice being L = 0 (S-wave annihilation) and L = 1 (P -wave annihilation). The leading term in the Taylor expansion of [u L (r)] ai around the origin is given by with [u (L+1) L (0)] ai denoting the (L + 1)-th derivative at the origin. Hence using (85) in (82), we have The quantity ψ (L,S) e 1 e 2 , ij is defined in (61) as the ratio of the annihilation matrix element to the same matrix element evaluated in the tree approximation, which corresponds to replacing V by V inf in the Schrödinger equation. The free Schrödinger equation can be solved exactly, and one finds Comparison with (88) gives We can avoid the computation of the phase-shift matrix M bi by making use of the fact that the Wronskian matrix is r-independent. Let [v L (r)] ai be the N linearly independent singular basis solutions with asymptotic behaviours which defines the matrix T . The Wronskian is defined as where the prime denotes the derivative. 13 Inserting the asymptotic behaviour of the regular and singular solutions, we obtain Both expressions must be equal due to the constancy of the Wronskian, which implies Plugging this into (90) gives which expresses ψ (L,S) e 1 e 2 , ij in terms of the large-r behaviour of the singular solutions. This can be obtained from the regular solutions as described below. The Sommerfeld factor is defined through the square of the annihilation amplitude. The definition (24) can now be evaluated in the form which is equivalent to the result first derived in [43]. Note, in compound-index notation, In practice, the matrix-Schrödinger equation must be solved numerically. We obtain the matrix T from the regular solutions by the following three steps: (1) Determine the N linearly independent regular solution vectors [u L (r)] ai for i = 1, . . . , N for given L by solving the radial Schrödinger equation with initial condi- wherer 0 is close to zero. In practice, we work with the dimensionless, scaled variabler = m LSP vr and user 0 = 10 −7 . The normalization is chosen such that [u (2) Pick a large value r ∞ , such that the asymptotic behaviour (91) of the irregular solution applies, and the Wronskian evaluates to and [u L (r)] aj known from step 1. Hence the matrix T appearing in (96) follows from matrix inversion, (3) Since in practice the Schrödinger equation can only be solved up to some finite value of r, check the stability of the result by varying (and increasing) r ∞ until T is independent of r ∞ within a certain target accuracy.
The procedure described here has been implemented in Mathematica and works well, when all N states included in the multi-channel Schrödinger equation are degenerate to a high degree. This is the case in MSSM parameter regions where the Sommerfeld enhancement is most effective, such as the wino or Higgsino limit for the neutralino, and when the other states not related to the wino or Higgsino electroweak multiplet are decoupled and ignored. However, our aim is to compute the Sommerfeld-enhanced radiative corrections in a larger part of the MSSM parameter space, when the mass splittings become larger than in the wino or Higgsino limit. In this case the method outlined above encounters severe numerical problems, as we now describe, and fails to provide the desired solution.
For the purpose of illustration we consider a MSSM parameter point where the LSP is wino-like. The relevant masses are m LSP = m χ 0 1 = 2749.4 GeV, m χ + 1 = 2749.61 GeV and m χ 0 2 = 2950.25 GeV. Further details of the models are not relevant for the following discussion. We first compute the Sommerfeld factor S ≡ S χ 0 ] as function of x ∞ = m LSP vr ∞ by keeping only the χ 0 1 χ 0 1 and χ + 1 χ − 1 two-particle states in the Schrödinger equation and the annihilation process. The velocity is chosen v = 0.012, slightly below the threshold for the χ + 1 χ − 1 state. The result S(x ∞ ) is shown as the light-grey (red), dashed curve in the upper panel of Fig. 5. After a rapid variation with a peak structure, the Sommerfeld factor reaches a plateau and for x ∞ > 50 stays at the constant value S(∞) ≈ 199.59. Next, we add the χ 0 1 χ 0 2 state to the problem, so that the Schrödinger system is now for a 3 × 3 matrix. Since the new state is 200 GeV heavier and moreover rather weakly coupled to the two lowest, nearly degenerate wino states, we expect that it should have little effect on the value of the Sommerfeld factor. However, now the numerical solution fails when x ∞ is slightly large than 1, as seen from the dark-grey (blue) curve in Fig. 5, which drops to 0 after a few spikes. It is not possible to reach the plateau, where S(x ∞ ) stabilizes. A similar behaviour impedes the calculation of the P -wave Sommerfeld factor, as shown in the lower panel of Fig. 5.
The numerical instability originates from the presence of kinematically closed twoparticle state channels, here the χ 0 1 χ 0 2 state. When the mass splittings become larger or v small, the situation 2m LSP +m LSP v 2 < M b is rather generic. To quantify the issue suppose that b in the Schrödinger equation (83) is a closed channel, while a refers to the two-LSP χ 0 1 χ 0 1 channel. The solution [u l (x)] bi for the closed channel involves an exponentially growing component proportional to . In the absence of channel mixing, the solutions for the kinematically open channels would be oscillating. The off-diagonal potentials V ab (r) that couple the different channels decrease as e −M EW r /r, where M EW is the mass scale set by the electroweak gauge boson exchange that mediates the channel mixing. Hence, the open-channel solutions [u l (x)] ai inherit the exponential growth from the closed channels. For the two-LSP χ 0 1 χ 0 1 channel, exponential growth occurs when at least one of the included kinematically closed channels b satisfies Since typically m LSP ≫ M EW for the dark matter scenarios of interest, this condition The light-grey (red) dashed curve shows S(x ∞ ), when only the χ 0 1 χ 0 1 and χ + 1 χ − 1 twoparticle states in the Schrödinger equation and the annihilation process are kept and the asymptotic regime is reached for x ∞ > 50. The dark-grey (blue) solid curve shows the result when the χ 0 1 χ 0 2 state is included. In this case, the evaluation fails for x ∞ > 2 and no reliable result is obtained.
is easily satisfied unless all two-particle states included in the computation are very degenerate within a few GeV or less. In consequence the formally linearly independent solutions [u l ] ai degenerate and the matrix U ai becomes ill-conditioned with exponentially growing entries in the rows corresponding to open channels a. The matrix inversion (100) can no longer be done in practice for r ∞ large enough such that the asymptotic regime is reached, which causes the instability seen in Fig. 5.

Improved method
To avoid the matrix inversion (100) we seek a method where the inverse of U defined in (99) follows directly from the solution of a differential equation system. This can indeed be found by an adaptation of the reformulation of the Schrödinger problem described in [44].
We work with the dimensionless radial variable x = m LSP vr and dimensionless wave We separate the asymptotic value of the potential by defining V (x) = V inf +V (x), such that V inf contains the mass splittings, and is diagonal, whileV (x) x→∞ = 0. The method proposed in [44] starts with the variable phase ansatz 14 Here are linearly independent Bessel function solutions of the free Schrödinger equations The Wronskian of the free solutions is The ansatz doubles the set of unknown functions by introducing α ai (x) and β ai (x) and the additional freedom can be eliminated by imposing the condition for every ai [44]. The condition reduces the N second-order differential equations (102) to a system of 2N coupled differential equations of first order for α ai (x) and β ai (x). The important quantities in the present approach are the matrix-functions N,α, defined by with O ab defined through β ai = O ab α bi , and From the original Schrödinger equation for [u L ] ai and the above definitions, one derives that N andα satisfy the first-order differential equations The initial conditions (97) translate into where x 0 is a number close to zero, which we set to 10 −7 . Up to this point we essentially reviewed the modification of the variable phase method proposed in [44]. The crucial observation is now that which implies that the matrix U ai from (99) is asymptotically trivially related toα ai , since The last equality follows, because g ′ a /g a − ik a vanishes for large x. Specifically, for the relevant cases L = 0, 1, We recall that our task is to compute the inverse U −1 for large enough x, such that the result is practically x-independent. This is now easily accomplished, since using d(α −1α )/dx = 0, (111) can be turned into the differential equatioñ with Z given in (111), and then from the solution to (116) for sufficiently large x ∞ . To summarize, the matrix T and hence the Sommerfeld factors (96) are computed by the following three steps: (1) Solve for N ab the first-order differential equations (110) with initial conditions (112) for every b = 1, . . . , N.
(3) Evaluate (117) for several x ∞ and check (by varying and increasing x ∞ ) that T is independent of x ∞ within a certain target accuracy.
To illustrate the improved performance of the method developed in this section, we show in the upper plot of Fig. 6 the same result as in Fig. 5, that is, the S-wave Sommerfeld factor S χ 0 in the wino-LSP model described at the end of Sec. 4.2, but now using the new method. The solution for the three-state case that was impossible to obtain with the standard method can now be evolved to sufficiently large x ∞ , such that the Sommerfeld factor can be extracted without difficulty. In fact, the two curves referring to the two-state solution and the one including the heavier χ 0 1 χ 0 2 state cannot be distinguished on the scale of the plot. As expected, the heavier state has no effect on the Sommerfeld enhancement in this particular model -the enhancement factor changes only by a tiny amount from S(∞) ≈ 199.59 to S(∞) ≈ 199.72.
The lower plot of Fig. 6 shows the corresponding results for a P -wave Sommerfeld factor, here for S χ 0 Once again the two-state solution (light-grey/red curve [invisible]) and the one including the heavier χ 0 1 χ 0 2 state (dark grey/blue curve) cannot be distinguished on the scale of the plot. The plateau is reached for x ∞ > 50 with S(∞) ≈ 4.31. We also show the two-state solution with the standard method (light-grey/red dashed curve). All results agree for large enough x ∞ , but contrary to the S-wave case there is a visible difference for x ∞ < 10. The difference arises from the approximation in the second line of (114), which implies that the U matrix in the standard and new method agree only for large x ∞ in general. The "approximation" is exact only for L = 0. For the P -wave case the difference vanishes as 1/x 2 ∞ . Since there is no numerical restriction on x ∞ with the new method, the difference can always be made sufficiently small as seen in the lower plot. Note that it would not be possible to compute the P -wave Sommerfeld factor with the standard method for the three-state problem, since the matrix inversion becomes unstable at x ∞ ≈ 2, as in the S-wave case.

Second-derivative operators
The above method allows us to calculate the matrix elements of the S-and P -wave operators O χχ→χχ  ) for v = 0.012 in a wino-like model as described at the end of Sec. 4.2 using the improved solution method. The light-grey (red) curve shows S(x ∞ ), when only the χ 0 1 χ 0 1 and χ + 1 χ − 1 two-particle states in the Schrödinger equation and the annihilation process are kept. The dark-grey (blue) curve shows the result when the χ 0 1 χ 0 2 state is included. In both plots the two cases are now indistinguishable. The lower plot shows in addition (light-grey (red) dashed) the P -wave Sommerfeld factor with the "old" method for the two-state problem.
differences. In order to fully describe the Sommerfeld enhancement of the short-distance annihilation cross section at order v 2 , we further require the matrix element of the second-derivative operators P χχ→χχ {e 4 e 3 }{e 2 e 1 } . For spin-0 it is defined as [15] P χχ→χχ {e 4 e 3 }{e 2 e 1 } ( 1 S 0 ) = with obvious generalization to spin-1. We show in this section that the matrix element of P χχ→χχ {e 4 e 3 }{e 2 e 1 } can be obtained from the non-derivative S-wave operators by the relations (21), (23). To this end, we show that where in compound index notation (23) reads We start from (65) and (73), which for the specific case at hand imply The momentum-space Schrödinger equation (76) leads to The momentum-space potential is a sum of terms from the exchange of gauge and Higgs bosons,V where m φ is the mass of the exchanged gauge boson, α 2 = g 2 2 /(4π) the SU(2) gauge coupling and c (a) ee ′ the coefficients of the potentials given in Tab. 2 of Appendix A. Shifting q → q + k in (122) factorizes the two integrations. Employing dimensional regularization gives the finite result for the linearly divergent integral. Inserting this into (121), we obtain × 2µ e (2m LSP − M e + E)δ e ′ e + 2µ e α 2 a m φa c (a) * ee ′ where the last equality follows from the definition (120) of the matrix κ and the expression for the relative momentum of the two-particle state e, p 2 e = 2µ e (2m LSP −M e +E)+O( p 4 e ). We can rename the dummy index e ′ →ē ′ in the second term of (125), and use that κēē′ = κ ee ′ (since the potentials for χ e 1 χ e 2 → χ e ′ 1 χ e ′ 2 and χ e 2 χ e 1 → χ e ′ 2 χ e ′ 1 involve the same coupling structure), to write (125) as which proves (119). Note that the second term in the definition (120) of κ is only present for the exchange of massive particles, and not for the Coulomb potential. The factor κ can be absorbed into an effective annihilation matrixĝ κ as has been done in (23).
In parameter-space regions where the Sommerfeld effect is important and the exchange of a particle with mass m φ ≪ m LSP leads to an enhancement of the radiative correction, all terms in κ are parametrically of the same order, since M e − m LSP ∼ E ∼ m LSP v 2 , and α 2 m φ ∼ vm EW ∼ m LSP v 2 . A general MSSM parameter point, however, may also lead to heavy Higgs bosons, which though strongly suppressed in the potential, may give a large contribution to the last term in square brackets. The origin of this unphysical power-counting breaking contribution is the linearly divergent integral (124). Since heavy Higgs-boson exchange causes a suppressed, local (χ † χ) 2 potential interaction, the simplest solution is to decouple the heavy Higgs bosons by not including them into the MSSM potential of Section 3. After this decoupling, the generated local interaction is a O(v 2 ) correction to the long-distance part, which we consistently neglect. In practice, we simply eliminate Higgs-exchange from the potential for Higgs masses m H /m LSP > 0.5 unless m H < 100 GeV.

Approximate treatment of heavy channels
For a generic MSSM parameter space point the two-particle states will span a certain mass range and not be degenerate. While the improved method for the Schrödingerequation solution allows us to cover the non-degenerate case without numerical instabilities there are still practical limitations related to the increasing CPU time, as the number of states, and hence the dimensionality of the matrix increases. For example, if for a given MSSM model and scattering energy E, the computation of the 1 S 0 Sommerfeld factor in the charge-neutral channel for a velocity close to the threshold of a nearly degenerate state takes 0.1 s when only two out of the total of 14 channels are included in the computation, the CPU time increases to 14 s for four channels, 5 min for 8 channels and reaches nearly three hours for the solution of the full 14 × 14 matrix problem. 15 Since even for a single MSSM parameter set, many Sommerfeld factors must be computed for a single scattering energy, and further the thermal average requires an integral over the scattering energy, CPU considerations restrict the number of channels, for which the Schrödinger equation is solved exactly.
However, an exact solution of the 14 ×14 problem is only necessary if all 14 states are nearly degenerate, which is rarely the case. As soon as the mass splitting becomes large, the heavier two-particle states should have little effect on the Sommerfeld enhancement of the lighter states, since the propagator of the heavier state in one of the loops of the ladder diagram series is off-shell and does not lead to an enhancement. Yet there is the possibility that a heavy state couples more strongly to the annihilation process than the lighter states and hence effectively enhances the annihilation rate. To cover this case, we allow the heavy channels to appear in the last loop before the annihilation vertex, but not elsewhere in the ladder diagrams. Moreover, the non-relativistic power-counting shows that there is a suppression factor of order [E/(M h −2m LSP )] a when a light channel is replaced by a heavy one, with a = 1/2 for the contribution in the last loop before the annihilation, but a = 3/2 when the heavy state appears inside the ladder, which justifies the omission. In the following, we describe how the effect of the heavy channels in the last loop can be absorbed into an effective annihilation matrix for the lighter channels.
We divide the total number N of two-particle states (14 in method-2 in the chargeneutral sector) into n light states, which will be treated exactly, and N − n heavy states, which we include only in the last loop. 16 Suppose that h = {h 1 h 2 } refers to one of the heavy states, then by extension of (65) we obtain where l refers to the light channels summed over, andG il is the Green function of the Schödinger operator for the n × n problem for the light-states. The potential interaction V lh converts the light channel into the annihilating heavy channel h. The annihilation cross sections matrix elements such as above are multiplied by annihilation matricesf ee ′ ≡f χχ→χχ {e 1 e 2 }{e 4 e 3 } and summed over all two-particle channels e, see (14). Dividing the sum over e into the sum over light and the sum over heavy channels and proceeding for the heavy channels from (127) to the integral involving the wave-function similar to (73), we find where the sum over l (h) runs only over the light (heavy) channels, and the equality in (128) holds when heavy channels contribute only in the last loop before the annihilation of the perturbative expansion of the matrix element.
Our aim is to write the square bracket as an effective annihilation matrixf eff le ′ K[ q ] for the light channels. We write the potential as in (123) and perform the integration over loop momentum k to obtain with for S-wave matrix elements (L = 0, K[ q ] = 1) and for the P -wave case (L = 1, K[ q ] = q ) The superscript "right" reminds us that (127) represents only "one half" of the fourfermion operator that accounts for the square of the annihilation amplitude. Accounting for the fermion bilinear matrix element that multiplies the annihilation matrix from the left then givesf eff ll ′ =f ll ′ + I lhfhl ′ + I * l ′ hf lh + I lh I * for the effective annihilation matrix in the space of light states including the heavy states in the last loop of the ladder. We wish to note that the result (129) taking channel index h to run over all relevant two-particle states regardless of their mass provides the leading-order non-relativistic approximation to the 1-loop annihilation amplitude of the state l, which has also been obtained by direct expansion in [45]. It can be checked that the loop integrals I S,P defined in the latter reference, are equal to (2π| q |)I 0,1 . 17 17 The non-relativistic expansion performed in [45] keeps in addition terms proportional to mass differences between the incoming χ's and the virtual χ's inside the loop that originate from the numerator of the full 1-loop amplitude. In our approach, the latter can correspond to corrections of O(v 2 ) to the long-distance part, which we do not consider, as well as to the short-distance annihilation, where we have kept all (m i − m LSP )/m LSP ∼ O(v 2 ) terms in the Wilson coefficients of the annihilation operators. We note that similar (m i − m LSP )/m LSP terms arise as well from other sources, for instance from the anti-particle pole contribution in the q 0 -integration of the full amplitude, which have been neglected in [45] and also in our approach, since they are a class of O(v 2 ) long-distance corrections.
The expressions for I L above are not yet useful, since they define functions of q 2 , hence the evaluation of (128) would require the full momentum-space wave-function [ψ E ] li ( q ), which is not available, and not just the value and derivative at r = 0 in coordinate space. However, we would like to make use of these expressions only for heavy channels when the mass splitting M h − 2m LSP is large compared to the scattering energy E. Since the typical relative momentum scales as q 2 ∼ m LSP E, we can expand the integrals I L in q 2 /(M h − 2m LSP ). Keeping only the leading term in this expansion, we approximate Inserting this into (130) gives , from which the effective annihilation matrix (133) follows.
The case of the second-derivative S-wave operator is somewhat more complicated, since one has to apply the equation-of-motion relation discussed in the previous section to the factor K[ k ] = k 2 in (128). This is done by writing which leads to the same factor as in (125). The result is that the effective annihilation matrix for the second-derivative operators readŝ where The factor M in each term in (138) is built from the particle species that enter the twoparticle states specified by the indices of the accompanying g matrix, as defined in (10).
To illustrate the quality of approximation to the heavy channels described in this section, we employ the same wino-LSP model as discussed at the end of Section 4.2.
The spectrum of neutralino masses is m χ 0 = {2749.4, 2950.25, 3062.98, 3083.54} GeV, the chargino masses are m χ + = {2749.61, 3074.26} GeV. We consider the charge-neutral two-particle state sector with method 2, in which case there are 14 two-particle states of increasing mass. We study the approach to the full treatment of the 14 × 14 matrix problem by plotting in Fig. 7 the dependence of the Sommerfeld factor S ≡ S χ 0 of the lightest χ 0 1 χ 0 1 state in the S-wave, spin-0 configuration on the number of states included in the computation. As n increases from 1 to 14, the dashed curve shows S(x ∞ = 200), when only the n lightest states are included and the remaining ones are ignored. The solid curve shows the result, when the n lightest states are included in the solution of the Schrödinger equation, and the remaining 14 − n heavier channels are treated approximately by including them in the last loop before the annihilation. For the purpose of this comparison we keep the δm e 4 e 1 δm e 4 e 1 /M 2 Z,W terms in the potential for all 14 two-particle states, not only for those n treated exactly, such that the potential is independent of the number of states n treated exactly. 18 We show the result for two velocity points, v = 0.012 (upper plot), slightly below the threshold for the nearly degenerate χ + 1 χ − 1 state, and v = 0.15 (lower plot). The figure shows that n = 1 is always a poor approximation. This says that it is necessary to treat the nearly degenerate χ + 1 χ − 1 state exactly, as is of course expected, since the bulk of the enhancement comes from the mixing of the χ 0 1 χ 0 1 state with the charged state in models where the LSP is almost pure wino. Once n > 1, the dependence of the Sommerfeld enhancement on the remaining states is no longer large, resulting in only another 5% increase. We observe that the results that include the heavy channels approximately have a smoother approach to the full result attained for n = 14 than the ones that leave out the heavier channels completely. This demonstrates that it is a good approximation to include the heavy channels only in last loop before the annihilation process, but to neglect them in the Schrödinger equation. In the two cases shown, it is essentially sufficient to treat only two states exactly and all other states approximately, as the solid curve for n = 2 is already very close to the n = 14 value.
While this shows that the one-loop approximation for heavy states is often a good approximation to the full, resummed result, it does not prove that either of the two, the one-loop or full resummed treatment of very heavy channels in the non-relativistic approximation is a reasonable approximation to the true loop corrections, since the non-relativistic approximation in the potential region is not expected to be a good approximation to the exact loop integral for very heavy channels. Since such states have little influence on the Sommerfeld factor anyway, in practice, we simply include them by 18 This is explains why S(x ∞ = 200) for v = 0.012 shown in the upper plot approaches the value 212.18 rather than the previously quoted value S(∞) ≈ 199.59, which corresponds to n = 2 on the dashed curve. If, one the other hand the mass-difference terms in the potential are neglected entirely, we obtain S(x ∞ = 200) ≈ 206.06, hence their effect is small. This is expected, since the Sommerfeld enhancement is mainly generated by the two lowest, nearly degenerate states, for which the mass difference terms are completely negligible. The rise of S in Fig. 7 when adding states 11, 13 and 14, which suggests that these states are comparatively more important than the lower ones (except the first and second), should be taken with a grain of salt, since it originates from the (δm/M EW ) 2 terms of the potential, which here are included for all states.  the approximate procedure described in this section.

Summary
With this paper we conclude the presentation of the non-relativistic MSSM effective theory formalism, designed to calculate the enhanced radiative corrections to pair-annihilation of neutralinos and charginos at small relative velocity. The construction of the NRMSSM is similar to the non-relativistic EFT for heavy quarkonium annihilation. Characteristic features of the NRMSSM are the presence of more than one heavy particle species; massive mediators exchanged among the heavy particles, which generate electroweak Yukawa (rather than colour Coulomb) potentials; and weak couplings that allow one to compute the long-distance matrix elements; these contains the so-called Sommerfeld enhancements. The framework developed here applies to general MSSM scenarios, where the LSP is an arbitrary admixture of the electroweak gauginos, and accounts for the co-annihilations with all the neutralino and charginos, which are nearly mass degenerate with the LSP. Papers I [14] and II [15] focused on the short-distance part of the annihilation process, providing analytic results for the Wilson coefficients that reproduce the tree-level annihilation rates of neutralinos and charginos including up to P -and next-to-next-toleading order S-wave effects, i.e. up corrections of O(v 2 rel ). The present paper describes all the necessary ingredients to account for the long-range interactions that produce the Sommerfeld enhancements in the annihilation rates at the leading order in the nonrelativistic power-counting, which treats α 2 ∼ v rel ∼ M Z,W /m LSP . The calculated enhancements thus sum up large O(α 2 /v rel ) n , n = 1, 2, . . . , quantum corrections to the tree-level annihilation rates to all orders. The main results of this this work are: • A general formula for the Sommerfeld-corrected annihilation cross-section for a given channel χ i χ j that includes O(v 2 rel ) corrections in the short-distance part has been written down in (25). The formula is a generalization of the tree-level approximation σ χ i χ j v rel = a + b v 2 rel . Since the Sommerfeld factors depend on the partial-wave configuration of the annihilating state, σ χ i χ j v rel has been decomposed in (25) as a sum over partial waves, where the annihilation coefficients for each partial wave were computed in papers I and II.
• Compact analytic results for the potential interactions among the chargino and neutralino pairs for a general MSSM parameter point, which feed into the multichannel Schrödinger equation that has to be solved to obtain the Sommerfeld correction factors. Transitions among different states are accounted for by the off-diagonal entries in the potential.
• A novel method to solve the multi-channel Schrödinger equation. The standard method does provide an accurate result for the Sommerfeld factors when kinematically closed channels are included due to numerical instabilities caused by the exponential growth of some of the solutions. The new method solves this problem by reformulating the Schrödinger equation for the quantity relevant for the Sommerfeld factor instead of solving for the scattering wave functions. This important improvement yields accurate numerical solutions when many non-degenerate twoparticle channels are included, up to CPU time limitations.
• An approximate treatment of heavy channels, which can includes them in an effective annihilation matrix for the Schrödinger equation of a lower-dimensional problem. We demonstrated that including heavy channels only in the last loop before annihilation provides a good approximation to the full resummed result, and thus reduces the number of channels to be treated exactly.
While this paper, as well as papers I and II, have been dedicated to the technical details describing the computation of the Sommerfeld-enhanced annihilation cross sections, the investigation of the effect on the relic abundance in well-motivated MSSM scenarios with heavy neutralino LSP is the subject of an accompanying publication [16], which further illustrates the general use and potential of the framework.  Table 2: Potentials describing the non-relativistic interactions among chargino and neutralino pairs in the MSSM. The potential from neutral scalar exchange φ has to be summed over the "physical" neutral Higgs bosons, φ = H 0 , h 0 , A 0 . The expressions obtained from the table correspond to the potential entries in method-1. The potentials for the channels not shown are obtained trivially by interchanging indices (like χ − χ + → χ − χ + or χ + χ 0 → χ + χ 0 ) or are vanishing (like χ − χ + → χ + χ − ). We have introduced λ S ≡ (3 − 4S) and λ Z/W = 1 + δm e 4 e 1 δm e 3 e 2 /M 2 where H 0 1 ≡ H 0 , H 0 2 ≡ h 0 and A 0 1 ≡ A 0 , A 0 2 ≡ G 0 . We have used the abbreviations s W = sin θ W and c W = cos θ W , with θ W denoting the Weinberg angle. The coupling factors for the three-point interaction of an incoming neutralino χ 0 j , an outgoing chargino χ + i and either an incoming charged Higgs particle G + or H + or an incoming W + -boson are given by where H + 1 ≡ H + and H + 2 ≡ G + , and the mixing matrices are defined as in Ref. [48]. Finally, three-point interactions of an incoming χ 0 j and an outgoing χ 0 i with a (pseudo-) scalar Higgs particle or the Z-boson involve the following (axial-) vector or (pseudo-) scalar coupling factors: The (axial-) vector and (pseudo-) scalar coupling factors in (140) and (142), which are all related to interactions with neutral SM and Higgs particles, satisfy v * ij = v ji , a * ij = a ji , s * ij = s ji , p * ij = −p ji .
as a consequence of the hermiticity of the underlying SUSY Lagrangian.

B Equivalence between method-1 and method-2
We wish to show in this appendix by means of an example, that using any of the two basis of two-particle states introduced in Sec. 3, named as method-1 and method-2, gives the same outcome for the Sommerfeld enhancement factors. Let us consider the simple case of a system which in method-1 consist of three states labelled with a = 11, 12, 21, corresponding to χ 1 χ 1 , χ 1 χ 2 and χ 2 χ 1 , respectively. We can think of χ 1 and χ 2 as two different neutralinos (or two different charginos), and we note that only state 11 is composed of identical particles. Conversely, the basis of method-2 is built from only two states, (χχ) 11 and (χχ) 12 .
The Schrödinger equation in method-1 is a 3 × 3 matrix where the generic potential matrix V (1) ab (a, b = 1, 2, 3) reads and we have used that the potential for scattering of state 11 to either 12 or 21 is the same, since the exchange of state 12 by 21, which corresponds to crossing the lines of χ 1 and χ 2 in the diagram, gives the same amplitude due to the identical particle nature of the state 11. In addition we have taken into account that V 32 = V 23 (as well as V 33 = V 22 ), since the amplitudes that produce these potential terms only differ by a relabeling of the internal vertices of the associated diagram. Note that hermiticity requires as well that V ab = V * ba . The Schrödinger equation (83) in method-1 when applied to a generic multi-component wave function Ψ = (u 1 , u 2 , u 3 ) reads D u 1 + V 11 u 1 + V 12 (u 2 + u 3 ) = 0 , D u 2 + V 21 u 1 + V 22 u 2 + V 23 u 3 = 0 , D stands for the differential operator multiplying δ ab in (83), whose specific form is irrelevant for the current discussion, and we also omit writing explicitly the dependence of the quantities u a and V ab on the independent variable r. We have to seek for three independent solutions of the system (145). Let us first choose two solutions of the form Ψ i = ( √ 2 u 1i , u 2i , u 2i ), with i = 1, 2. Then the set of equations (145) applied to Ψ 1,2 reduce to two independent ones (since eqs. (145) are symmetric in u 2 and u 3 ), and acquire the form where we have used the symmetry relations under the exchange of the particle labels e 1 ↔ e 2 derived in [14,15], namelyf {e 2 e 1 }{e 4 e 3 } 2S+1 L J = (−1) S+Lf {e 1 e 2 }{e 4 e 3 } 2S+1 L J and likewise for the exchange e 3 ↔ e 4 , to write the annihilation matrix entries involving state 21 in terms of those of 12. These same relations also tell us that the annihilation rates involving channel 11 vanish for odd L + S, since such an state cannot be formed out of two identical spin-1/2 particles. Note also thatf 21 =f * 12 . On the other hand, the matrix T is determined following (94) from the matrices [u The Sommerfeld factors for the state 21 are found to be equal to those of state 12; this should be the case since both states are physically equivalent. We have used that f 21 =f * 12 to write (150) in a more compact form. Let us now turn to method-2. As mentioned above in this case the basis is given by the two states a = 11, 12. Applying the rules described in Sec. 2 (see (52) and the discussion following that equation), we can build the potential matrix entries corresponding to method-2 from those of method-1. Distinguishing even and odd S + L, we find In particular, note that for odd L + S the identical-particle state 11 does not exist, and the entries of the potential matrix involving that state have to be set to zero, to avoid that a potential transition 11 → 12 could give a non-vanishing annihilation amplitude for an incoming state 11, as argued in Sec. 2. In method-1 this is taken into account automatically and one gets the result S odd 11 = 0, as found in (151), at the expense of using two different states (12 and 21) to describe the physical state built from χ 1 and χ 2 . Coming back to method-2, we then see that in the odd L + S sector we deal with just a single annihilating state. We can therefore switch to a one-dimensional problem with V (2), odd (r) = V 22 − V 23 .
The Schrödinger equations which derive from the potential V (2), even for the wave function Ψ = (ũ 1 ,ũ 2 ) of method-2 read Eqs. (153) agree with (146) obtained previously for the components of the solutions Ψ 1,2 of method-1. Therefore we can choose Ψ 1 = (u 11 , u 21 ) and Ψ 2 = (u 12 , u 22 ) as linear independent solutions for method-2 in the case L + S = even. The matrices [u For the case of odd L + S, there is just one equation, which is the same differential equation as (147) Finally, we require the annihilation matrices in method-2 to compute the Sommerfeld factors. These are obtained from those of method-1 by adding a factor of 1/ √ 2 for each identical-particle state involved in the entry, as explained in Sec. 2: From the expressions for [u (L+1) L (0)] and M obtained for the even and the odd cases in (154) and (156), respectively, it is immediate to compute the corresponding T matrices. The Sommerfeld factors for method-2 then follow from (96) by inserting the annihilation matrices (157). This straightforward computation yields results for the Sommerfeld factor of states 11 and 12 for even and odd L+S that match the formulae (150) and (151) derived above using method-1.