Heavy neutralino relic abundance with Sommerfeld enhancements - a study of pMSSM scenarios

We present a detailed discussion of Sommerfeld enhancements in neutralino dark matter relic abundance calculations for several popular benchmark scenarios in the general MSSM. Our analysis is focused on models with heavy wino- and higgsino-like neutralino LSP and models interpolating between these two scenarios. This work is the first phenomenological application of effective field theory methods that we have developed in earlier work and that allow for the consistent study of Sommerfeld enhancements in non-relativistic neutralino and chargino co-annihilation reactions within the general MSSM, away from the pure-wino and pure-higgsino limits.


Introduction
The experimental determination of the cold dark matter density in our Universe has reached percent level accuracy, Ω cdm h 2 = 0.1187 ± 0.0017, taking recent PLANCK data into account [1]. While the nature and origin of the cosmic dark matter component are still unknown it is intriguing that the observed abundance can be explained rather naturally as thermal relic of a TeV scale particle with weak interaction strength. A central ingredient in the relic abundance calculation of a particle dark matter (DM) candidate is its pair-annihilation rate. Consequently the increasing experimental precision on Ω cdm h 2 has triggered a particular interest in the calculation of radiative corrections to particle DM annihilation cross sections.
Probably the best motivated and certainly one of the most studied particle dark matter candidates is the neutralino LSP (χ 0 1 ) in the minimal supersymmetric standard model (MSSM) [2,3]. Several codes [4,5] allow for the calculation of the χ 0 1 relic abundance in the general MSSM, currently relying on co-annihilation rates calculated at tree-level. The calculation of radiative corrections to these rates follows two different directions. On the one hand, a great effort is undertaken in the calculation of nextto-leading order co-annihilation rates in fixed order perturbation theory in the MSSM. The complete next-to-leading order SUSY QCD corrections in χ 0 1 co-annihilations with potentially nearly mass degenerate charginos and sfermions has been performed [6][7][8][9][10] and the first steps in the calculation of the full one-loop electroweak corrections are undertaken [11][12][13]. On the other hand it has been noted some time ago that in nonrelativistic dark matter pair-annihilations a certain class of radiative corrections can be enhanced and requires a systematic resummation up to all loop-orders [14,15]. The resulting Sommerfeld enhancement arises naturally in theories with light mediator exchange between the co-annihilating non-relativistic dark matter particles prior to their actual annihilation reaction. For heavy χ 0 1 dark matter this effect should be addressed generically in the relic abundance calculation as well as in indirect detection rates: in both cases the annihilating particles move at non-relativistic velocities and the mutual exchange of electroweak gauge bosons -and to a lesser extent, Higgs bosons -prior to the actual annihilation gives rise to long-range potential interactions eventually requiring a systematic resummation of certain contributions up to all loop orders. In context of χ 0 1 pair annihilation reactions, Sommerfeld enhancements have been first addressed in the pure-wino and pure-higgsino χ 0 1 scenarios in [14,15] and were subsequently studied in [16,17]. The particular relevance in χ 0 1 indirect detection has been investigated for the pure-wino case in [18][19][20].
In [21][22][23] we have developed a formalism that allows to systematically address the calculation of enhanced radiative corrections in non-relativistic neutralino/chargino pairannihilation reactions in the general MSSM by means of a non-relativistic effective field theory approach, where our particular focus is the consistent calculation of Sommerfeld enhancements. By "general MSSM" we imply that the lightest neutralino can be an arbitrary admixture of wino, higgsino and bino. Analytic results for the short-distance coefficients encoding hard tree-level annihilation reactions of non-relativistic co-annihilating neutralinos and charginos including P -and next-to-next-to-leading order S-wave rates are given in [21,22]. Corresponding analytic expressions for the long-range potential interactions eventually causing enhanced annihilation rates are presented in [23]. In the latter work we also describe the technical details involved in a precise determination of Sommerfeld enhancements in the χ 0 1 relic abundance calculation. It is worth noting that in addition to covering the general case of χ 0 1 being an arbitrary admixture of the electroweak gaugino eigenstates, our approach extends previous investigations on the subject in several other aspects, such as the consistent treatment of off-diagonal annihilation rates, the separation into S-and P -wave components with their own, separate Sommerfeld factors, and the ability to deal with many nearly mass-degenerate states.
The purpose of this paper is a detailed investigation and discussion of Sommerfeld enhancements in the χ 0 1 relic abundance calculation in some popular MSSM scenarios. The underlying physics effects are analysed in detail in each step of the calculation. This allows to illustrate the general use of our method [21][22][23] applicable in the general MSSM and to address the question of viability of popular MSSM scenarios in light of a consistent treatment of the Sommerfeld effect. We choose to consider three scenarios taken from the set of Snowmass pMSSM benchmark models [24]. These models pass all constraints from so far unsuccessful SUSY searches at the LHC, additional collider, flavour and precision measurement bounds as well as constraints from dark matter direct detection experiments and indirect searches. The neutralino LSP relic abundance within these models, calculated from perturbative annihilation rates, is not larger than the WMAP bound, but can be smaller than the experimentally measured value. The latter allows for the case that neutralino dark matter does not make up all the cosmic cold dark matter. In addition to these benchmark scenarios we investigate the Sommerfeld enhancements in neutralino/chargino co-annihilations in a set of models interpolating between a scenario with almost pure-higgsino χ 0 1 to a wino-like χ 0 1 model. The MSSM spectra for the models on this "higgsino-to-wino" trajectory are generated with DarkSUSY [4]. As our work allows for the first time a consistent study of the Sommerfeld effect on the relic abundance calculation for models with mixed wino-higgsino neutralino LSP we provide an extensive discussion of the Sommerfeld effect in such a scenario.
Throughout this work we neither include thermal effects nor the effect of running couplings. As regards thermal effects in context of dark matter relic abundance calculations including Sommerfeld enhancements, the temperature dependence of the gauge boson masses has, for instance, been considered in [16,25]. Concerning the running of couplings, this effect can in principle be relevant to Sommerfeld-enhanced rates as well: the annihilation process involves the mass scale of the co-annihilating particles, that is associated with the hard annihilation reaction, as well as the much smaller scale of the non-relativistic kinetic energies of the co-annihilating particle pairs and the masses of the exchanged particles. The latter scales are connected with the physics that causes the Sommerfeld enhancement. Both the effect from thermal corrections and from running couplings will be investigated in future work.
The structure of this paper is as follows. In Sec. 2 we consider the case of a wino-like χ 0 1 benchmark scenario taken from [24], followed by the investigation of a higgsino-like χ 0 1 benchmark spectrum in Sec. 3. In both cases we compare to results obtained in the well-studied "pure" wino and higgsino scenarios where the χ 0 1 is assumed to be part of an unbroken SU(2) L triplet or two unbroken SU(2) L doublets. As Sommerfeld enhancements have been studied extensively in the particular case of a pure wino χ 0 1 in the literature, we address the question of the validity of conclusions inferred from these pure wino and higgsino scenarios to wino-and higgsino-like χ 0 1 spectra in the general MSSM. In Sec. 4 the effect of Sommerfeld enhancements in co-annihilations of wino-like neutralino and chargino states in a bino-like χ 0 1 benchmark scenario is considered. A "higgsino-towino" trajectory is defined in Sec. 5, by introducing 13 models that interpolate between a higgsino-and wino-like χ 0 1 spectrum while the relic density calculated from perturbative rates is kept fixed. Our discussion here is focused on the spectra and the obtained relic abundances omitting particular details on the Sommerfeld enhanced co-annihilation cross sections. The specific features of the Sommerfeld effect for a mixed wino-higgsino χ 0 1 are subsequently studied in detail in Sec. 6, where the selected spectrum is one of the trajectory models of the preceding section. We draw our conclusions and give an outlook to future work in Sec. 7.
The present paper is intentionally stripped of all technical details underlying the computation of the results and focuses on the nature of the Sommerfeld enhancement and its physics interpretation. The reader interested in technical aspects is referred to [21,22] for the computation of the annihilation cross sections and in particular to [23] for the computation of the Sommerfeld enhancements and the solution of the multichannel Schrödinger equation.
2 Wino-like χ 0 1 Wino-like χ 0 1 dark matter arranges into an approximate SU(2) L fermion triplet together with the two chargino states χ ± 1 . In the SU(2) L ×U(1) Y symmetric limit the triplet would be assigned zero hypercharge. All states χ 0 1 , χ ± 1 share the same O(TeV) mass scale, characterised by the wino mass parameter M 2 , m χ ∼ |M 2 |. Electroweak symmetry-breaking introduces a small mass splitting between the neutral and the charged components of the triplet. The tree-level mass splitting happens to be very small, O(m 4 W /m 3 SUSY ), and the one-loop radiative corrections dominate over the tree-level splitting.
A pMSSM scenario with wino-like χ 0 1 is provided by the SUSY spectrum with model ID 2392587 in [24]. A measure for the wino fraction of a given neutralino LSP state is the square of the modulus of the neutralino mixing-matrix entry Z N 21 . For pMSSM scenario 2392587 the χ 0 1 constitutes a rather pure wino, |Z N 21 | 2 = 0.999, with a mass m LSP ≡ m χ 0 1 = 1650.664 GeV. The mass of the chargino partner χ ± 1 is given by m χ + 1 = 1650.819 GeV, such that δm = m χ + 1 − m χ 0 1 turns out to be 0.155 GeV. Without any modification these values are taken from the spectrum card provided by [24] where the mass parameters refer to the DR-scheme. As the precise sub O(GeV)-scale χ 0 1 χ ± 1 mass splitting is an essential ingredient in the calculation of the Sommerfeld-enhanced coannihilation rates we have to assume an accuracy of the given mass spectrum at the level of 10 MeV for our analysis of the Sommerfeld enhancement in the pMSSM scenario to be meaningful. A rigorous analysis of Sommerfeld-enhanced co-annihilation processes in a given model should refer to the on-shell mass spectrum of the neutralino and chargino states instead of DR-parameters, where a sub-GeV scale precision of the mass parameters requires the consideration of one-loop renormalised quantities. For reference purposes, however, we do not modify the publicly available DR-spectra of [24] for all three pMSSM models discussed here.
In the context of minimal dark matter models [25], wino dark matter is realised as the neutral component of an approximate SU(2) L triplet state as well. In contrast to MSSM scenarios with wino-like χ 0 1 , the SU(2) L triplet minimal dark matter models (referred to as "pure-wino" models in the following) consider interactions of the dark matter states with the electroweak gauge bosons only. Two-particle final states in minimal dark matter pair-annihilation reactions are hence given by pairs of SM particles and the SM Higgs boson and all heavier states above the minimal dark matter mass scale are treated as completely decoupled. Such a scenario agrees with the decoupling limit in a MSSM scenario with wino-like χ 0 1 LSP. To the contrary, the wino-like pMSSM model that we consider here features non-decoupled sfermion states at the 2 − 3 TeV scale with non-vanishing couplings of the χ 0 1 and χ ± 1 to sfermions and to the (heavier) Higgs states, though the latter are suppressed with respect to the couplings to the gauge bosons, because any Higgs-χχ (tree-level) interaction takes place between the gauginocomponent of the one and the higgsino-component of the other χ. As the higgsino-like neutralino and chargino states in the pMSSM model under consideration reside at the O(3.9 TeV) scale any Higgs-χχ interaction plays a sub-dominant role in our analysis of pair-annihilation reactions of the wino-like χ 0 1 and χ ± 1 states. Due to the non-decoupled sfermion states though, some annihilation rates in the wino-like χ 0 1 pMSSM scenario are reduced with respect to the pure-wino dark matter case.
In the calculation of the relic abundance we have to take into account all possible two-particle co-annihilation reactions between the (approximate) SU(2) L triplet states χ 0 1 , χ ± 1 . In addition, in the pMSSM model 2392587, the bino-like χ 0 2 is only about 8% heavier than the χ 0 1 , m χ 0 2 = 1781.37 GeV. Hence the χ 0 2 is a potentially relevant coannihilating particle as well. It turns out though, that this state eventually plays no role for the relic abundance, as the corresponding cross sections are strongly suppressed with respect to those of the wino-like particles χ 0 1 and χ ± 1 due to the much weaker couplings of the bino-like χ 0 to gauge bosons and to the remaining χ 0 /χ ± states. All remaining heavier particles in the pMSSM scenario lie above the 2 TeV scale, so they are already Boltzmann suppressed and hence practically irrelevant during the χ 0 1 freeze-out. Sommerfeld enhancements on the co-annihilation rates are taken into account by including in the multi-state Schrödinger equation all χχ two-particle states with mass smaller than M max = 2 m χ 0 1 +m χ 0 1 v 2 max , where we set v max = 1/3. This choice is motivated by the fact that v max roughly corresponds to the χ 0 1 's mean velocity around freeze-out, hence these states are potentially relevant for co-annihilation processes, and can still be produced on-shell in a χ 0 1 χ 0 1 scattering process. The remaining heavier two-particle states with mass above M max are included in the computation of the Sommerfeld enhancement of the lighter states in the last loop before the annihilation, following the method developed and discussed in [23]. The χχ-channels, whose long-distance interactions are treated exactly, can be classified according to their total electric charge. The sector of neutral two-particle states comprises the χ 0 1 χ 0 1 and χ + 1 χ − 1 channels. In the pMSSM scenario considered here, this sector contains in addition the χ 0 1 χ 0 2 state. In the singlecharged and the double-charged sectors of a pure-wino dark matter scenario there is only one state present in each sector, χ 0 , whereas in the pMSSM scenario we have to add in addition a second state with χ 0 1 replaced by χ 0 2 , in agreement with the rule above that defines the channels which enter the Schrödinger equation. Since the bino-like neutralino essentially neither couples to the wino-like particles nor to gauge bosons, and because sfermion states are rather heavy, potential interactions as well as tree-level annihilation reactions involving the bino-like χ 0 2 are strongly suppressed with respect to the corresponding interactions with wino-like particles χ 0 1 , χ ± 1 . As a consequence, χ 0 2 plays essentially no role for Sommerfeld enhancements, and we focus the discussion that follows on the channels built from the wino-like χ 0 1 and χ ± 1 states only. In each of the charge sectors long-range interactions due to potential exchange of electroweak gauge bosons, photons and light Higgses are present. 1 Potential W -boson exchange leads to a Yukawa potential interaction that induces transitions between the χ 0 1 χ 0 1 and the χ + 1 χ − 1 state in the neutral sector. Hence the part of the neutral sector consisting of the channels χ 0 1 χ 0 1 and χ + 1 χ − 1 is characterised by a potential matrix with non-vanishing off-diagonals which are of the same strength as the diagonal entries. As the incoming χ 0 1 χ 0 1 pair cannot build a 3 S 1 or 1 P 1 state, potential interactions are responsible for transitions between the two neutral states χ 0 1 χ 0 1 and χ + 1 χ − 1 in a 1 S 0 or 3 P J configuration.
In Fig. 1 we plot the enhancement (σ SF v)/(σ pert v) of annihilation rates including long-range interactions, σ SF v, with respect to the perturbative tree-level result, σ pert v, for the two-particle states χ 0 1 χ 0 1 and χ + 1 χ − 1 in the neutral sector of the model as a function of the velocity v LSP of the incoming χ 0 1 's in their centre-of-mass frame. We define the velocity v LSP by √ s = 2m χ 0 1 + m χ 0 1 v 2 LSP with √ s the available centre-of-mass energy. The spin-averaged tree-level annihilation rates σ pert v are calculated in the non-relativistic approximation where v denotes the relative velocity of the annihilating particles. In case of the χ 0 1 χ 0 1 state the relation between the relative velocity v and v LSP is given by 1 Potentials from Higgs exchange are negligible compared to the leading contributions from gauge bosons in the pMSSM scenario with wino-like χ 0 1 , again because in any Higgs-χχ vertex the gaugino component of one χ is coupled to the higgsino component of the other χ. In the wino-like χ 0 1 Snowmass model the lowest-lying χ's relevant for the Sommerfeld effect are rather pure wino-like χ 0 and χ ± states with a very small higgsino component.  Figure 1: The enhancement of the χ 0 1 χ 0 1 and χ + 1 χ − 1 annihilation cross sections for Snowmass model 2392587 relative to the perturbative tree-level rate, (σ SF v)/(σ pert v). The solid lines refer to the calculation of the Sommerfeld-enhanced rates with off-diagonal entries in the annihilation matrices Γ properly included. The dashed curves show the enhancement with respect to the perturbative cross sections when off-diagonal annihilation rates are not considered. The dotted curve labelled "pure-Coulomb enh." shows the enhancement from photon exchange only in the χ + 1 χ − 1 channel.
The coefficients a and b are determined from the absorptive part of partial-wave decomposed Wilson coefficients given in [21,22]. In case of the Sommerfeld-enhanced rates σ SF v each partial wave contribution gets multiplied by an enhancement factor related to the two-particle wave-function of the respective incoming state, see [23] for the detailed expression. Unless otherwise stated, Sommerfeld-enhanced results include the one-loop corrections from heavy χχ-states in the last potential loop, following the approximation discussed in [23]. The results for the wino-like pMSSM scenario hence include perturbative corrections from heavy χχ-pairs involving the higgsino-like χ 0 3,4 and χ ± 2 particles. The effects of the latter nevertheless amount only to a negligible per mil level deviation on σ SF v. This can be traced back to the fact that the higgsino states lie at the rather high mass scale of around 3.9 TeV and thus are basically decoupled. The (σ SF v)/(σ pert v) curves in Fig. 1 show some characteristic features, which we describe next. As there is a small mass splitting between the χ 0 1 and the χ ± 1 , the threshold for the on-shell production of the heavier neutral state χ + 1 χ − 1 opens at v LSP /c ≃ 0.014. Well below this threshold, the enhancement for the χ 0 1 χ 0 1 system is velocity-independent and of O (10). This saturation effect is characteristic for Yukawa-type interactions in the kinematic regime where the relative momentum of the incoming state is well below the mass scale of the mediator: this is the case for the χ 0 1 χ 0 1 state at very small velocities, where off-diagonal Yukawa potentials are generated by W -boson exchange with m χ 0 1 v LSP ≪ m W . The actual strength of the enhancement is, however, a combined effect of the off-diagonal Yukawa potential from W -exchange that allows for χ 0 1 χ 0 1 → χ + 1 χ − 1 transitions and the QED Coulomb interaction in the (kinematically closed) χ + 1 χ − 1 channel. At velocities v LSP just below the χ + 1 χ − 1 threshold resonances in the χ 0 1 χ 0 1 channel can be observed. While the main plot in Fig. 1 displays a curve smoothed over this region, we show in the small sub-figure a closeup of the resonance pattern. The existence of resonance enhancements at the threshold of a heavier channel is well-known and has been described for instance in [26]. However, opposed to the pattern in the close-up in Fig. 1 no oscillating behaviour was found in [26], as only Yukawa potentials were considered. In fact the oscillatory pattern is related to the photon exchange in the χ + 1 χ − 1 subsystem. Going to even larger velocities, above the χ + 1 χ − 1 threshold, the enhancement in the χ 0 1 χ 0 1 channel decreases, approaching one as we depart from the non-relativistic regime. Turning to the enhancement in the χ + 1 χ − 1 channel, it shows quite a different behaviour right above its threshold compared to the χ 0 1 χ 0 1 system at small velocities: instead of approaching a constant value, the enhancement factor for χ + 1 χ − 1 rises increasingly as the velocities of the χ ± 1 get smaller. Such a behaviour is expected in the presence of long-range Coulomb-potential interactions, where the enhancement does not saturate because the mediator is massless. Indeed, the photon exchange between the charged constituents of the neutral χ + 1 χ − 1 pair dominates the potential interactions in the regime of very small velocities: the Yukawa potentials become very short-ranged and thus negligible compared to the Coulomb-interaction. The dotted (black) curve in Fig. 1 displays the enhancement factor in the χ + 1 χ − 1 system arising from Coulomb interactions due to photon exchange only. For small velocities the pure-Coulomb enhancement factor diverges as 1/v χ + 1 . The true enhancement curve, that involves all potential interactions affecting the χ + 1 χ − 1 system asymptotically reaches this Coulomb-like behaviour for velocities directly above the χ + 1 χ − 1 threshold. 2 For larger velocities in the χ + 1 χ − 1 system the presence of the Yukawa potentials leads to a larger enhancement than in case of Coulomb interactions only.
The dashed curves in Fig. 1 show the enhancements (σ SF v)/(σ pert v) for the χ 0 1 χ 0 1 and χ + 1 χ − 1 states when off-diagonal terms in the annihilation matrices are (incorrectly) left out. This can lead to a 30% underestimation of the actual enhancement in the χ 0 1 χ 0 1 channel. The effect is less pronounced for the χ + 1 χ − 1 channel, as in this case the cross section also gets significant contributions from 3 S 1 annihilations and not just from 1 S 0 ones. As the 3 S 1 sector is purely diagonal, the effect of off-diagonals, relevant in the case of 1 S 0 wave annihilations, becomes milder for the spin-averaged total cross section σ SF v. It is worth to stress that the overall order of magnitude of the enhancements is O(10), and becomes O(10 2 ) in the resonance region around the χ + 1 χ − 1 threshold. The quantity that enters the Boltzmann equation for the neutralino number density 2 Note that in spite of the ∝ 1/v χ + 1 divergence, the enhanced cross sections lead to a finite result in the average over the thermal velocity distribution due to the v 2 χ + 1 term in the integration measure, is the thermally averaged effective annihilation rate σ eff v . Fig. 2 shows σ eff v as defined in [23] as a function of the inverse scaled temperature x = m χ 0 1 /T . The lower solid (blue) curve represents the perturbative (tree-level) annihilation rates while the upper solid and the dashed (red) lines refer to Sommerfeld-enhanced cross sections including and neglecting off-diagonal annihilation rates, respectively. The plot can be divided into several regions with different characteristics. Let us first note that for x 10 the depicted behaviour of σ eff v is unphysical. The mean velocity of the annihilating particles in the plasma scales as 1/x and hence is no longer non-relativistic for x < ∼ 10 while the results of our framework strictly apply only to non-relativistic χχ pair-annihilations, i.e. for x 10. Around x ∼ 20 the annihilation rates of χ 0 1 and χ + 1 can no longer maintain chemical equilibrium and the particles start to decouple from the thermal plasma. Hence only the region above x ∼ 20 is important for the calculation of the relic abundance. Around x 10 4 the number densities of the χ ± 1 are so strongly Boltzmann suppressed with respect to the χ 0 1 number density despite the small mass splitting that the rates of the charginos basically play no role in the effective rate σ eff v , which is then essentially given by χ 0 1 χ 0 1 annihilations. Note that we can estimate the point of chargino decoupling between x ∼ 10 4 − 10 5 from the ratio of the Boltzmann distribu- , taking the O(10 −1 GeV) mass splitting into account. After χ ± 1 decoupling, σ eff v including the Sommerfeld enhancements becomes constant, which we can infer from the constant enhancement factor for the χ 0 1 χ 0 1 system for very low velocities shown in Fig. 1. Before χ ± 1 decoupling, σ eff v including the Sommerfeld enhancements rises with increasing x due to the contributions from the charginos but also due to the velocity-dependent enhancement on the χ 0 1 χ 0 1 system itself for larger relative velocities. On the contrary, the perturbatively determined σ eff v shows a constant behaviour before and after χ ± 1 decoupling with a rise only around the decoupling region; the contributions that dominate the perturbative cross sections in the non-relativistic regime are the velocity-independent leading-order S-wave terms. Fig. 3 compares the thermally averaged effective rates σ eff v as calculated from the wino-like pMSSM scenario and from a pure-wino SU(2) L triplet minimal dark matter model with the same χ 0 1 mass. In the pure-wino model the mass splitting between the χ 0 1 and χ ± 1 has to be kept in the Schrödinger equation as it is of the same order as the non-relativistic kinetic energy and the potentials. However in the hard annihilation rates the mass splitting is a subleading effect and is neglected; the annihilation matrices in the pure-wino model depend on the χ 0 1 mass only (the corresponding expressions can be found, for instance, in [23]). While the rates for χ 0 1 χ 0 1 annihilations agree at permille level, the cross sections involving χ ± 1 are generically larger by factors of O(1) in the pure-wino model as compared to the pMSSM wino-like model. This can be mainly traced back to the destructive interference between t-channel sfermion and s-channel Z (and Higgsboson) exchange amplitudes in χ + 1 χ − 1 → f f annihilations in the pMSSM scenario case, while the t-channel sfermion exchange amplitudes are absent in the pure-wino model. In addition the pure-wino case neglects all final state masses which in particular gives rise to larger annihilation rates into the tt and electroweak gauge boson final states as compared to the pMSSM scenario, where the non-vanishing masses of all SM particles are taken into account. This accounts for the deviation between the curves in Fig. 3 before χ ± 1 decoupling. Finally we consider the yield Y = n/s, defined as the ratio of the number density n of all co-annihilating particle species divided by the entropy density s in the cosmic co-moving frame. The dependence of the yield on the scaled inverse temperature x = m χ 0 1 /T is governed by a Boltzmann equation and the χ 0 1 relic abundance is obtained from the yield today. In Fig. 4 we show the ratio of the yield Y calculated from Sommerfeld-enhanced cross sections in both the pMSSM and the pure-wino model to the corresponding results using perturbative cross sections, Y pert , as a function of x.
First note, that the denominator Y pert in the ratio Y /Y pert differs for the pMSSM and the pure-wino model, which is a consequence of the different effective rates σ eff v , see Fig. 3. Further, in case of the pMSSM scenario we show results corresponding to a calculation of Y including and neglecting off-diagonal annihilation rates. Around x ∼ 20 the yields including Sommerfeld enhancements start to depart from the corresponding perturbative results; the enhanced rates delay the freeze-out of interactions, which leads to a reduction of the yield Y compared to the perturbative result Y pert . The most drastic reduction in Y /Y pert occurs between x ∼ 20 and x ∼ 10 3 . In this region the enhancement factors on the cross sections are of O(10) (and not yet O(10 2 ) as for very large x), leading to Y /Y pert values that deviate from 1 by a few 10%. For x 10 5 the fraction Y /Y pert stays constant, meaning that at these temperatures the particle abundances in both the perturbative and Sommerfeld-enhanced calculation are frozen in. In case of the wino-like model we find that the relic densities calculated from the yield today read Ω pert h 2 = 0.112 and Ω SF h 2 = 0.066. Hence taking into account the Sommerfeld effect leads to a reduction of the calculated relic abundance of around 40%. On the other hand, neglecting the offdiagonal annihilations in the calculation of Sommerfeld-enhanced rates overestimates the relic density by 15% compared to the correct Ω SF h 2 . Let us recall that the relic density calculated without corrections from heavy χχ-states in the last potential loop differs from the Ω SF h 2 value quoted above at most at the per mil level. Due to overall larger hard annihilation rates in the pure-wino model, the calculated relic density including Sommerfeld-enhanced rates turns out to be Ω SF pure-w h 2 = 0.034, while the corresponding perturbative result is Ω pert pure-w h 2 = 0.056. It is difficult to quantify the theoretical error on such numbers. Conventional treelevel calculations of annihilation cross sections and the ensuing relic densities neglect radiative corrections and are supposed to be accurate to O(5%) in the absence of enhanced corrections due to non-relativistic scattering, large Sudakov logarithms, or, potential strong-interaction effects for quark and gluon final states. When the Sommerfeld effect is included, the latter two restrictions still apply. The computation of the Sommerfeld effect itself neglects O(v 2 ) corrections to the scattering potentials as well as ordinary, non-enhanced corrections to the short-distance annihilation coefficients. Hence the accuracy of the Sommerfeld-corrected annihilation cross sections and relic densities is presumably again at the O(5%) accuracy level at best.
The higgsino-like neutralino χ 0 1 arises as the lightest out of four mass eigenstates χ 0 1,2 , χ ± 1 related to two SU(2) L fermion doublets. Note that the hypercharges of the two SU(2) L doublets are given by Y = ±1/2 respectively, which ensures the electric neutrality of the χ 0 1 . The common mass scale of the χ 0 1,2 , χ ± 1 states is set by the O(TeV) higgsino mass parameter, m χ ∼ |µ|. Electroweak symmetry breaking introduces a tree-level splitting between m χ 0 1 and the masses of the three heavier states of O(m 2 Z /m LSP ) ∼ O(1 GeV). This is considerably larger than the tree-level mass splitting in the wino-like χ 0 1 case; in particular loop corrections play a sub-dominant role in the mass splittings of higgsino-like neutralinos and charginos.
As in Sec. 2, it is instructive to compare the pMSSM scenario with higgsino-like χ 0 1 and co-annihilating χ 0 2 and χ ± 1 to a model with pure-higgsino χ 0 1,2 , χ ± 1 states and completely decoupled sfermions and heavy Higgses. We refer to the latter scenario as "pure-higgsino" model; such model is also discussed in the context of Minimal Dark Matter [25]. Pure-Higgsino states interact only with the SM gauge bosons W ± , Z, γ but not with the Higgs bosons. The accessible final states in 2 → 2 co-annihilation reactions of pure higgsinos are hence given by particle pairs formed out of SM gauge bosons and fermions as well as of the (SM-like) Higgs h 0 , where all these SM particles are taken to be massless, and only SM gauge bosons and higgsinos appear as intermediate states in tree-level annihilations. The co-annihilation rates of the higgsino-like χ 0 1,2 , χ ± 1 states in the pMSSM scenario 1627006 happen to be larger than the corresponding reactions in the pure-higgsino case. This can be traced back to the presence of non-decoupled sfermion and Higgs states in the higgsino-like χ 0 1 pMSSM model and in particular to non-decoupled wino-like states χ 0 3 , χ ± 2 at the scale of 1.6 TeV. In the determination of the χ 0 1 relic abundance for this pMSSM scenario including co-annihilations only the higgsino-like states are relevant. Other heavier states are already sufficiently Boltzmann-suppressed during χ 0 1 freeze-out. Hence we neglect the co-annihilations of the lightest sfermion statesτ 1 andν 3 , with masses around 1.44 TeV, although we include co-annihilation reactions of all heavier χ 0 /χ ± states. Yet the latter have basically no effect on the χ 0 1 relic density, as their abundances are already sufficiently suppressed at χ 0 1 decoupling. Obviously, in the pure-higgsino scenario only the co-annihilations between the higgsino-like species χ 0 1,2 , χ ± 1 are taken into account for the calculation of the relic abundance.
We consider Sommerfeld corrections to all co-annihilation rates between two higgsinolike particles in both the pMSSM scenario 1627006 and the pure-higgsino model by treating all channels built from the states χ 0 1,2 , χ ± 1 exactly in the corresponding Schrödinger equations. Moreover, the remaining heavier χ 0 /χ ± two-particle states in the higgsinolike pMSSM scenario are treated perturbatively in the last potential loop, see [23]. In case of the pure-higgsino model though, all heavier states are considered as completely decoupled. Dividing the co-annihilation reactions into sets corresponding to total electric charge, we identify a neutral sector with the four two-particle states χ 0 1 χ 0 1 , χ 0 1 χ 0 2 , χ 0 2 χ 0 2 and χ + 1 χ − 1 . The single-positive (negative) charged sector contains the two states χ 0 , whereas the double-positive (double-negative) charged sector features only one two-particle state relevant in co-annihilations with the higgsino-like χ 0 1 dark matter candidate: . Note that annihilations of the latter double-charged states The enhancement factor for the additionally relevant channel χ 0 1 χ − 1 agrees with the one for the χ 0 1 χ + 1 pair. Solid lines refer to the calculation of the Sommerfeld-enhanced rates with off-diagonal terms in the annihilation matrices properly included. Dashed curves show the enhancement when the off-diagonal annihilation rates are neglected. χ + 1 χ + 1 and χ − 1 χ − 1 are absent in the pure-higgsino model due to hypercharge conservation in this SU(2) L × U(1) Y symmetric limit, as they have a non-zero hypercharge, namely Y χ ± χ ± = ±1. In contrast, in the higgsino-like χ 0 1 pMSSM case with broken U(1) Y symmetry, annihilations of the double-charged channels into a W + W + or W − W − pair are possible, though the rates are suppressed by a factor ∼ m W /m χ 0 1 compared to the magnitude of the neutral sector's leading rates. Fig. 5 shows the enhancement (σ SF v)/(σ pert v) of the individual cross sections for those channels that have the most relevant contribution to the relic abundance calculation, that is χ 0 1 χ 0 1 , χ + 1 χ − 1 , χ 0 1 χ 0 2 in the neutral sector, and χ 0 1 χ + 1 in the single-charged sector (χ 0 1 χ − 1 gives the same contribution). First note that the enhancements are only of O(1), opposed to O(10 2 ) enhancements in case of the wino-like model in Sec. 2. This can be explained due to the larger mass splittings to the next-to-lightest states χ ± 1 , χ 0 2 in the higgsino-like χ 0 1 case and the fact that the couplings to SM gauge bosons and (light) Higgs particles are generically smaller for higgsinos than for winos. The enhancement of the χ 0 1 χ 0 1 rate as a function of the velocity v LSP shows again the saturated, velocityindependent behaviour typical for Yukawa type potentials in the low velocity regime well below the thresholds of the heavier two-particle states. As in the wino-model, both the off-diagonal Yukawa potential and the (diagonal) Coulomb potential in the kinematically closed χ + 1 χ − 1 channel contribute here to the actual size of the enhancement. At larger velocities, two resonance regions at the thresholds for χ + 1 χ − 1 and χ 0 2 χ 0 2 production are visible (the χ 0 2 χ 0 2 channel opens up at v LSP /c ≃ 0.127; the ratio (σ SF v)/(σ pert v) for this channel is very close to 1, and is not shown in Fig. 5). One might ask why no resonance at the χ 0 1 χ 0 2 threshold is visible in the χ 0 1 χ 0 1 channel: recall that Fermi-statistics forbids the χ 0 1 χ 0 1 -pair to build the totally symmetric partial-wave configurations 3 S 1 and 1 P 1 . In case of unbroken SU(2) L × U(1) Y symmetry it turns out, though, that the χ 0 1 χ 0 2 pair can build 3 S 1 and 1 P 1 configurations but not 1 S 0 and 3 P J states. Hence there are no off-diagonal entries in the neutral potential matrices encoding χ 0 1 χ 0 1 ⇌ χ 0 1 χ 0 2 interactions in the pure-higgsino limit. Departing from the SU(2) × U(1) Y symmetric limit gives rise to χ 0 1 χ 0 2 contributions to the enhancement (σ SF v)/(σ pert v) in the χ 0 1 χ 0 1 channel that are however suppressed by (m W /m χ 0 1 ) 3 with respect to the leading contributions; this explains why no χ 0 1 χ 0 2 threshold effect is visible in Fig. 5. Such restrictions due to non-accessible partial-wave configurations do not exist for the next-to-lightest neutral two-particle state χ + 1 χ − 1 , and resonances at the thresholds of all co-annihilating neutral χχ-pairs heavier than the χ + 1 χ − 1 are visible in the latter channel in Fig. 5. Furthermore, note the 1/v χ + 1 Coulomb-type enhancement in the χ + 1 χ − 1 channel directly above its threshold caused by potential photon-exchange between the χ + 1 and χ − 1 . The Coulomb potential surpasses the potentials from massive gauge boson and Higgs exchange at very small velocities in the χ + 1 χ − 1 channel, but for moderate velocities both the Coulomb and the (off-)diagonal Yukawa interactions are relevant. Turning to channel χ 0 1 χ 0 2 , the corresponding enhancement (σ SF v)/(σ pert v) increases as the velocity decreases. In particular, there is no saturation of the enhancement directly above threshold, because the lighter channels χ 0 1 χ 0 1 and especially χ + 1 χ − 1 are always kinematically open and accessible from an on-shell χ 0 1 χ 0 2 state via off-diagonal potential interactions. The ratio (σ SF v)/(σ pert v) for the charged state χ 0 1 χ + 1 that is additionally plotted in Fig. 5 (lowermost magenta line) shows that the Sommerfeld effect can also produce corrections that reduce the perturbative result. For the channel χ 0 1 χ + 1 the negative correction arises from the interference of amplitudes where, after multiple electroweak and Higgs boson exchanges, the state that annihilates into the light final state particles is the same as the incoming one, χ 0 1 χ + 1 , with amplitudes where the actual state that annihilates is χ 0 2 χ + 1 . In the EFT formalism such interferences arise from the off-diagonal annihilation terms χ 0 1 χ + 1 → χ 0 2 χ + 1 and χ 0 2 χ + 1 → χ 0 1 χ + 1 , combined with the off-diagonal potential term for χ 0 1 χ + 1 → χ 0 2 χ + 1 . The dashed magenta curve in Fig. 5 refers to the situation where off-diagonal short-distance rates are neglected in the calculation of the Sommerfeld enhanced χ 0 1 χ + 1 annihilation cross section. It is nicely seen that the destructive interference effect disappears in this case and the ratio (σ SF v)/(σ pert v) is always positive. The enhancement in the χ 0 1 χ + 1 channel also saturates as its on-shell production threshold is approached. This should be the case as the χ 0 1 χ + 1 channel is the lightest in the single positive-charged sector, and its behaviour should be similar to the one of the lightest neutral channel, χ 0 1 χ 0 1 , directly above threshold. However, such saturation is not visible in Fig. 5 because there we plot the χ 0 1 χ + 1 cross section as a function of v LSP and not as a function of the relative velocity of the channel, related to the latter by v 2 = 2(m χ 0 . 3 Let us also mention that the dip in the χ 0 1 χ + 1 cross section caused by interference effects is located at the velocity where the other state included in the Schödinger equation for this charge sector, χ 0 2 χ + 1 , opens up. As we have already noted in context of the χ 0 1 χ + 1 channel above, the dashed curves in Fig. 5 show the results for the corresponding enhancements of the pMSSM scenario 1627006 when off-diagonal annihilation rates are neglected. This disregard would lead to an underestimation of the actual enhancement due to the long-range potential interactions of around 30% in the χ 0 1 χ 0 1 channel. The effect is much milder for the χ 0 1 χ 0 2 and χ + 1 χ − 1 pairs and is explained by the contributions of 3 S 1 partial-wave annihilations to the cross sections (absent for the identical particle-pair channel χ 0 1 χ 0 1 ); off-diagonal 3 S 1 annihilation rates are suppressed relative to the leading (diagonal) rates by an order of magnitude, due to destructive interference effects between sfermion and gauge boson exchange amplitudes. Hence, as off-diagonals play a minor role in 3 S 1 annihilations, their effect in the spin-averaged cross sections σ SF v will also be less pronounced. As the conclusions on the enhancements in case of the pure-higgsino χ 0 1 model are similar to the results in Fig. 5 we do not show a corresponding plot here. Let us mention though again, that the hard co-annihilation rates in the pure-higgsino model are a few percent smaller than in the higgsino-like χ 0 1 model. Furthermore, the off-diagonal rates for 3 S 1 annihilations in the system of χ 0 1 χ 0 2 and χ + 1 χ − 1 states are of the same order of magnitude as the diagonal ones. Fig. 6 shows the thermally averaged effective annihilation rate σ eff v as a function of the inverse scaled temperature x. The lower solid (blue) curve represents the result using perturbatively calculated rates, while the upper two (red) curves with solid and dashed line style refer to computations with Sommerfeld-enhanced cross sections including and neglecting off-diagonal annihilation rates, respectively. Again the region for x 10 is unphysical, as the co-annihilating particles' mean velocities are outside the non-relativistic regime. Due to larger mass splittings between the higgsino-like neutralino and chargino states, the decoupling of the heavier states χ ± 1 and χ 0 2 takes place already around x ≃ 10 3 . As can be seen from Fig. 6, the Sommerfeld effect enhances the thermally averaged effective annihilation cross section by 3% up to 25% with respect to the perturbative result in the region of x around 10 − 10 3 which is most relevant in the relic abundance calculation. The effect of correctly treating off-diagonal annihilation rates is most essential for large values of x in the range 10 4 − 10 8 , where σ eff v would be underestimated by around 25% if off-diagonals were neglected in the hard annihilation rates. In the region x = 10 − 10 3 the effect of off-diagonal rates is also noticeable, leading to an overestimation of σ eff v that reaches 6% if off-diagonal rates are not taken into account. The latter difference with respect to the true result is traced back to the contribution to σ eff v of the charged χ 0 1 χ + 1 channel, which in the absence of off-diagonal annihilation terms does not get the negative interference term that lowers the Sommerfeld-corrected 3 If the χ 0 1 χ + 1 cross section behaves as σ SF v ≃ a + bv 2 close to threshold, the saturation is visible because of the zero slope of this function at v = 0; in terms of v LSP it reads σ SF v = a + b ′ (v 2 LSP − c), which does not have a zero slope at the threshold of the channel, v LSP = √ c. cross section when off-diagonal rates are consistently taken into account then explains why the correct σ eff v result crosses the dashed line for x > ∼ 10 3 in Fig. 6. Finally, Fig. 7 shows the ratio Y /Y pert . The solid (blue) and dashed (black) curves refer to calculations within the pMSSM Snowmass model 1627006 with off-diagonal annihilation reactions included and neglected, respectively. The dot-dashed (red) line applies to the pure-higgsino model. The relic abundances that we calculate within the pMSSM Snowmass model read Ω pert h 2 = 0.108 if perturbative annihilation reactions are considered and Ω SF h 2 = 0.100 taking Sommerfeld-enhanced rates into account. Accounting for the long-range potential interactions hence leads to a reduction of 8% on the predicted relic density for the pMSSM higgsino-like χ 0 1 model. Neglecting off-diagonal rates in the pMSSM Snowmass model calculation reduces the relic abundance to a value Ω SF, no-off h 2 = 0.096. This is because the effective thermal average cross section without the off-diagonal rates is larger in the region where chemical decoupling takes place, see Fig. 6. The error on Ω SF h 2 when disregarding off-diagonal rates therefore amounts to an underestimation of 4% in this case. The Sommerfeld-enhanced rates without the one-loop corrections from heavy χχ-states in the last potential loop before annihilation give a 1% deviation on the final Ω SF h 2 result. In contrast, the relic abundances in the pure-higgsino model, obtained using perturbative or Sommerfeld-enhanced rates, almost coincide, namely Ω pert pure-h h 2 = 0.127 and Ω SF pure-h h 2 = 0.126, where the latter result includes the off-diagonal rates. As can be expected, the overall smaller annihilation rates in the pure-higgsino scenario lead to a larger relic abundance than in the higgsino-like pMSSM scenario. The fact that the perturbative yield surpasses the Sommerfeld-corrected one right after chemical decoupling in the pure-higgsino model is explained by the slightly smaller σ eff v in the Sommerfeld-corrected result in that region of x, which is in turn produced by the Sommerfeld suppression in the charged channels χ 0 1 χ ± 1 . Overall, there is a strong cancellation between cross section enhancement in the neutral and suppression in the charged channels, leading to an almost vanishing net Sommerfeld correction.

Light scenario
Light neutralino dark matter with a relic abundance of the order of the observed value is realised for a χ 0 1 with a sizable bino component. The bino is a SU(2) L singlet with zero hypercharge. As for a pure bino there are no interactions with electroweak gauge bosons nor photons we can already expect that there will be essentially no long-range potential interactions for the bino-like χ 0 1 and hence no Sommerfeld enhancements in χ 0 1 χ 0 1 annihilations. Yet it is interesting to confirm this expectation and to investigate the relevance of Sommerfeld enhancements in possible co-annihilations with (slightly) heavier neutralino and chargino states. As an example for such a bino-like χ 0 1 we chose to study the pMSSM Snowmass model with ID 2178683 that features wino-like NLSP states with masses around 6% heavier than the χ 0 1 state: m χ 0 1 = 488.8 GeV, m χ 0 2 = 516.0 GeV and m χ + 1 = 516.2 GeV. In the calculation of the χ 0 1 relic abundance we consider co-annihilation reactions among all χ 0 /χ ± two-particle states, although only the two-particle annihilations between the states χ 0 1,2 , χ ± 1 are relevant since the higgsino-like states χ 0 3,4 , χ ± 2 lie at the 2 TeV scale and their abundances are strongly Boltzmann-suppressed at χ 0 1 freeze-out. The lightest sfermions are theτ 1 andν τ with masses around 770 GeV and we neglect their effect in the relic abundance.
Sommerfeld corrections on the co-annihilation cross sections from all two-particle states built from χ 0 1,2 and χ ± 1 are determined exactly through the solution of the corresponding Schrödinger equations in each charge sector. The outcome for the enhancement (σ SF v)/(σ pert v) in the neutral sector, which entails the two-particle states χ 0 1 χ 0 1 , χ 0 1 χ 0 2 , χ 0 2 χ 0 2 and χ + 1 χ − 1 , is shown in Fig. 8. Solid (dashed) curves correspond to a calculation with (without) off-diagonal annihilation rates in the Sommerfeld-enhanced reactions. Due to the absence of interactions with the electroweak gauge bosons in case of a purebino state, the χ 0 1 of the pMSSM Snowmass model 217868 also experiences basically no long-range potential interactions and there is essentially no coupling between the binolike χ 0 1 and the NLSP χ 0 2 . As a consequence, both the absolute (perturbative as well as Sommerfeld-enhanced) χ 0 1 χ 0 1 and χ 0 1 χ 0 2 annihilation rates are strongly suppressed and there is no enhancement in these reactions; the ratio (σ SF v)/(σ pert v) is equal to one in both cases. As it cannot be inferred from Fig. 8, let us note in addition that the absolute χ 0 1 χ 0 1 (χ 0 1 χ 0 2 ) annihilation cross section is suppressed with respect to the dominant χ 0 2 χ 0 2 and χ + 1 χ − 1 rates by four (two) orders of magnitude. In the subsystem of the neutral wino-like two-particle channels χ 0 2 χ 0 2 and χ + 1 χ − 1 , the Sommerfeld enhancement due to long-range potential interactions is effective, see the corresponding curves in Fig. 8. Note that χ 0 2 and χ ± 1 co-annihilations should still be relevant in the χ 0 1 relic abundance calculation within the pMSSM scenario 2178683, as the threshold velocities for χ 0 2 χ 0 2 and χ + 1 χ − 1 on-shell production are v χ 0 1 < ∼ 0.34 c and thus of the order of typical χ 0 1 velocities during thermal freeze-out. This scenario provides an example showing that the criterion established before for including long-distance effects among two-particle states with masses smaller than M max = 2m χ 0 1 + m χ 0 1 v 2 max and v max = 1/3 should not be considered rigidly. Rather it has to be reassessed according to the given MSSM spectra to avoid overlooking interesting effects. Consequently, in order to account for the wino-like subsystem formed by the states χ 0 2 χ 0 2 and χ + 1 χ − 1 we have set v max = 0.34 in the light scenario. At very small velocities the enhancements in the χ 0 2 χ 0 2 and χ + 1 χ − 1 channels show the characteristics discussed already for the wino model in  Sec. 2: In the χ 0 2 χ 0 2 system we find resonances just below the χ + 1 χ − 1 threshold, smoothed out in Fig. 8. The strength of the enhancement below and above this resonance region is a combined effect of the (off-diagonal) Yukawa and the diagonal Coulomb potential interactions in the χ + 1 χ − 1 system. In particular the enhancement is finite below the χ + 1 χ − 1 threshold. To the contrary, the χ + 1 χ − 1 channel shows the typical Coulomb-like 1/v χ + 1 enhancement from the dominating photon-exchange potential at velocities directly above its on-shell production threshold. Opposed to the O(10 2 ) enhancements found in Sec. 2, the overall enhancements of the neutral wino-like two-particle channels here reach factors of O(1) only. These less pronounced enhancements result from the lower masses of the wino-like states, since as m χ 0 1 decreases the Yukawa potentials from electroweak gauge boson exchange eventually become short-ranged as compared to the Bohr radius of the system proportional to (m χ 0 1 α EW ) −1 , where α EW = g 2 2 /(4π) and g 2 denotes the SU(2) L gauge coupling. Fig. 9 displays the effective annihilation cross section σ eff v (x). The dominance of the wino-like χ 0 2 , χ ± 1 particle annihilation rates by more than three orders of magnitude before their decoupling near x ∼ 100 is clearly visible. The Sommerfeld enhancement affects only the annihilation of the wino-like particles and therefore disappears for x > 100. Although the Sommerfeld factors for these channels lead to O(1) enhancements of the cross sections above the threshold near v LSP ∼ 1/3, similar in magnitude to the model with wino-like LSP for the same velocities, the thermal average over v LSP dilutes the enhancement, since the cross section for the heavy channels vanishes below the threshold. pert. Figure 9: The thermally averaged effective rate σ eff v (x) within the pMSSM Snowmass model 2178683 with Sommerfeld enhancements (upper red curve) and in the perturbative computation (lower blue curve). The result from disregarding off-diagonal rates in the Sommerfeld-enhanced processes is plotted by the dashed line. However the latter curve basically overlays with the upper (red) curve in this plot. This is because the Sommerfeld-enhanced σ eff v (x) is dominated by the χ 0 2 χ 0 2 and χ + 1 χ − 1 rates (before χ 0 2 and χ ± 1 decoupling), and the effect of disregarding off-diagonals in the latter gives a correction of around 10% only, see Fig. 8.
Nevertheless, the small enhancement visible in Fig. 9 occurs precisely in the x range most relevant for freeze-out. The effect of co-annihilations with the wino-like NLSP states therefore leads to a reduction of the yield when taking into account Sommerfeld enhancements with respect to the perturbative case, as is shown in Fig. 10. The relic density with perturbative annihilation rates is found to be Ω pert h 2 = 0.120. There is a ∼ 15% reduction of this result when considering the Sommerfeld-enhanced rates, Ω SF h 2 = 0.102. The latter sizable reduction of the relic density is attributed purely to the co-annihilating heavier wino states. Note that in the sector of wino-like states the potentials from massive gauge boson and photon exchange are equally important for the Sommerfeld enhancement, while in the χ + 1 χ − 1 system the Coulomb potential dominates over the Yukawa potentials only for very small velocities of the charginos. Neglecting the perturbative correction from the heavier χχ-states not included in the Schrödinger equation leads essentially to no difference (below per mil level) in the relic density, as the heavy higgsino-like χ 0 3,4 , χ ± 2 species lie at the scale of around 2 TeV. If no off-diagonals in the calculation of Sommerfeld-enhanced rates were considered, the relic abundance  would be overestimated by 3.5%.

Higgsino-to-wino trajectory
In case of the wino-like χ 0 1 model of Sec. 2 we have seen that the relic abundance including Sommerfeld enhancements on the co-annihilation rates is reduced by about 40% with respect to the result calculated from tree-level annihilation rates. In contrast, the model with higgsino-like χ 0 1 in Sec. 3 shows a less strong reduction, which is however still of the order of Ω SF h 2 /Ω pert h 2 ≈ 0.9. The difference in the reduction factor Ω SF h 2 /Ω pert h 2 between the wino-and the higgsino-like χ 0 1 model was explained by the smaller Sommerfeld enhancements in the latter case due to larger mass splittings between all co-annihilating particles and the fact that the potential interactions happen to be generically weaker for higgsino-like compared to the wino-like χ 0 1 models. In addition, we observed a Sommerfeld suppression effect in the single-charged sector of the pure higgsino scenario as well as the higgsino-like Snowmass model. Departing from the scenarios with rather pure wino, higgsino or bino χ 0 1 , we may ask ourselves about the features of a model with χ 0 1 LSP that contains both significant wino and higgsino contributions. It is worth to mention here that previous work in the literature focused on the wino-or higgsino-like χ 0 1 cases only, due to the lack of expression for potentials and annihilation matrices for a generically composed χ 0 1 state. Our results allow for the first time to perform a rigorous study of Sommerfeld enhancements in χχ pair-annihilations within models with mixed gaugino and higgsino composition of the co-annihilating neutralinos and charginos. We find it particular instructive to consider a series of models in the MSSM parameter space that describes the transition from a model with higgsino-like χ 0 1 to a model with primarily wino-χ 0 1 . In the following we will refer to this series of models as models on a "higgsino-to-wino" trajectory. We are interested in the case of reductions of Ω SF h 2 relative to Ω pert h 2 by 10% here and hence will not consider a significant bino-admixture to the χ 0 1 ; as we have seen in Sec. 4 the bino-like χ 0 1 itself does not experience any Sommerfeld enhancement. In such a situation a reduction of Ω SF h 2 can only arise due to co-annihilating particles with Sommerfeld-enhanced rates, see for example the model discussed in Sec. 4 with co-annihilating wino-like NLSPs.
In order to define the models for the higgsino-to-wino trajectory, we should note first that the proper choice of the two SUSY parameters µ and M 2 controls the higgsino and wino content of the mass eigenstate χ 0 1 . In order to avoid a bino-admixture to the χ 0 1 state we will choose the parameter M 1 , that controls the neutralinos' bino-content, to be sufficiently larger than both µ and M 2 throughout this section. Our setup excludes accidental mass degeneracies of the MSSM sfermions with the χ 0 1 , which implies that the actual parameters of the sfermion sector play a minor role in the choice of adequate models on the trajectory. Let us recall that the sfermion sector is irrelevant for Sommerfeld enhancements in our setup, as the latter are caused by potential gauge boson and light Higgs exchange between neutralino and chargino two-particle states prior to the hard annihilation reactions. The sfermion sector parameters only affect the precise value of the hard (tree-level) annihilation rates. The sfermion -basically the stop -sector however controls the value of the Higgs h 0 mass and we will adjust its parameters such that the experimental value for m h 0 is reproduced within 2.5% accuracy. Yet matching the precise experimental Higgs mass value is in fact not important to us here, as potential exchange from the h 0 gives always a sub-leading contribution to the potentials compared to the effects from SM gauge boson exchange.
In order to generate MSSM scenarios on a higgsino-to-wino trajectory we hence make the following choice for MSSM input parameters in the spectrum generation: • fix a common sfermion mass scale of 9 TeV, • set the trilinear couplings to A t = A b = 9 TeV, • fix m A 0 = 500 GeV and • choose tan β = 15.
All other trilinear couplings are assumed to vanish. The gluino mass parameter M 3 is fixed by M 3 = α s /(sin(θ w ) α e ) M 2 , but this choice is completely irrelevant to our discussion. To avoid a significant bino-admixture to the χ 0 1 we further restrict to models with M 1 = 10 M 2 . This leaves us with yet-to-choose parameter pairs in the µ − M 2 plane.
We require that the trajectory models allow for an explanation of the observed cosmic cold dark matter in terms of the neutralino relic abundance without including radiative corrections: in order to do so we employ the program DarkSUSY [4] and identify (µ, M 2 ) pairs such that the DarkSUSY calculated relic density Ω DS h 2 matches the most accurate determination obtained from the combination of PLANCK, WMAP, BAO and high resolution CMB data, Ω cdm h 2 = 0.1187 ± 0.0017 [1] 4 . In such a way we define 13 models on the higgsino-to-wino trajectory. The position of these models in the µ − M 2 plane is shown in Fig. 11. For each of the 13 models, given the pairs (µ, M 2 ) as well as the remaining input parameters defined above, we run our code and determine the corresponding relic densities including and neglecting Sommerfeld effects. The comparison between our perturbative results Ω pert h 2 with the corresponding DarkSUSY expressions Ω DS h 2 provides a cross-check of our perturbative calculation.
There is one important point to note concerning the MSSM spectrum generation from the SUSY input parameters. The DarkSUSY spectrum calculated from the inputs refers to tree-level DR-parameters. It is well-known that the mass splitting between a wino-like neutralino and its chargino partner is dominated by radiative corrections; the leading one-loop contribution to the splitting is of O(160 MeV) and dominates over the O(1 MeV) tree-level contribution. Both for the calculation of the Sommerfeld enhancements and in the determination of the relic abundance including co-annihilations a precise knowledge of the mass splitting between the χ 0 1 LSP and the NLSP particles is crucial and in a rigorous analysis we should therefore consider the spectra determined with one-loop accuracy. To this end we have been provided by one-loop on-shell renormalised SUSY spectra for all 13 models on the trajectory by a member of the collaboration [27,28]. The values of the input parameters µ, M 2 , . . . are the same as for the corresponding calculation within DarkSUSY with the difference that for the one-loop on-shell spectrum generation these inputs are considered as on-shell parameters and no renormalisation group running of the mass parameters is performed. Hence there are small differences in the values for the masses and mixing-matrix entries between the spectra that we use in our code and the corresponding DarkSUSY spectra. In particular the mass splittings between the χ 0 1 LSP and the NLSPs obtained from the on-shell masses renormalised at one-loop can be significantly different from the splittings derived using tree-level DRparameters. There exist different renormalisation schemes for on-shell renormalisation in the neutralino/chargino sector [27][28][29][30]: for all trajectory models apart from model 8 the on-shell renormalisation has been performed requiring that the values of the two chargino masses as well as the heaviest (in all our models bino-like) neutralino mass at one-loop are given by their tree-level values ("CCN-scheme"). Such a scheme works well as long as the two charginos are rather pure wino-and higgsino-like states. As soon as the charginos are (strongly) mixed wino-higgsino states -as in case of our model 8, where the input parameters µ and M 2 happen to be very close to each other -a more suitable scheme is obtained when only one chargino, one lighter neutralino and the heaviest bino-like neutralino mass are fixed to their tree-level value ("CNN scheme"). in Tab. 1, together with the one-loop renormalised LSP mass m χ 0 1 as well as the oneloop on-shell mass splitting δm χ + 1 = m χ + 1 − m χ 0 1 . The χ ± 1 is the NLSP in all models considered in this section. As additional information we give the χ 0 1 's wino fraction |Z N 21 | 2 and collect the results for Ω SF h 2 including Sommerfeld effects as well as for the suppression Ω SF h 2 /Ω pert h 2 of the former relic density with respect to the perturbative result. Both Ω SF h 2 and Ω pert h 2 are calculated from our programs, and the latter shows small deviations of the order of a few percent from the DarkSUSY value Ω DS h 2 = 0.1187. As can be read off Tab. 1 we can categorise the models on the trajectory to feature either a higgsino-like χ 0 1 with wino fraction below 10% but a higgsino fraction |Z N 31 | 2 + |Z N 41 | 2 above 0.9 (models 1−6), a mixed wino-higgsino χ 0 1 where both the wino and the higgsino their code.  Table 1: Information on the models on the higgsino-to-wino trajectory. The first column is the model ID while the second and third column contain the input parameter values for µ and M 2 . The one-loop on-shell renormalised χ 0 1 LSP mass is given in the fourth column and we provide the one-loop mass splitting to the lighter chargino, δm χ + 1 = m χ + 1 − m χ 0 1 in the fifth column. The χ ± 1 are the NLSP states in all models considered here. In the sixth column the wino fraction, |Z N 21 | 2 , of the χ 0 1 is specified. The second-to-last and the last columns give the relic density including Sommerfeld-enhanced cross sections as well as the suppression factor of the Ω SF h 2 with respect to the perturbative result Ω pert h 2 . The results including the Sommerfeld enhancements involve corrections from heavier χχ-pairs in the last potential loop. fraction lie within 0.1 − 0.9 (models 7 − 9) or a predominantly wino-like χ 0 1 with wino fraction above 0.9 (models 10 − 13). For all models we collect the relic density results Ω pert h 2 and Ω SF h 2 in Fig.12. The bars with dotted (black) hatching indicate Ω pert h 2 . Bars with solid-line (red) and dashed (blue) hatching give the corresponding results including Sommerfeld enhancements with and without off-diagonal rates, respectively. In particular for the higgsino-like models 1−6 but also for models 7−9 our relic densities Ω pert h 2 agree very well with the relic density Ω DS h 2 = 0.1187 calculated with DarkSUSY for the same set of input parameters. The latter relic density value is indicated by the black horizontal line and the grey horizontal band comprises all values deviating at most by 5% from the Ω DS h 2 value. For the wino-like models our relic density results deviate by < ∼ 8% from the corresponding DarkSUSY value. Let us discuss the characteristics of the models in the three different classes corresponding to their wino and higgsino admixture in turn. The models 1 − 6, with Figure 12: Relic densities Ωh 2 for models 1 − 13 on the higgsino-to-wino trajectory calculated with our code. The charts with dotted (black) hatching are the perturbative results Ω pert h 2 . Bars with dashed (blue) and solid-line (red) hatching refer to a calculation with Sommerfeld-enhanced cross sections neglecting and properly including off-diagonal rates, respectively. The grey shaded band comprises Ωh 2 values within 5% around the mean experimental value Ω cdm h 2 = 0.1187 [1]. The latter value is indicated by the black horizontal line and agrees with the DarkSUSY result for all 13 MSSM models on the trajectory. predominant higgsino composition, resemble the higgsino model of Sec. 3. This applies also to the corresponding shapes of the Sommerfeld-enhanced rates σ SF v, σ eff v , as well as to the yields Y /Y pert , that we do not show here. The reduction in the relic density when taking the Sommerfeld effect into account ranges from 3% to 14% for trajectory models 1 − 6. Models 1 − 3, with a 3% to 4% reduction are close to a pure-higgsino limit behaviour, whereas models 4 − 6 yield a similar outcome as for the Sec. 3 higgsinolike χ 0 1 Snowmass model. The potential interactions among all two-particle states built from the higgsino-like particles χ 0 1,2 , χ ± 1 have been accounted for exactly by solving the corresponding multi-state Schrödinger equation in models 1 − 6. This is in agreement with the criterion introduced in Sec. 2 that considers the long-distance effects among all χχ-states with mass smaller than M max = 2 m χ 0 1 + m χ 0 1 v 2 max , where v max = 1/3 is of the order of the χ 0 1 's mean-velocity during freeze-out. Heavier χχ channels enter the calculation through the perturbative corrections to the annihilation rates of the lighter channels treated exactly, and their tree-level co-annihilation rates are also included in the calculation of the χ 0 1 relic density, as done in the previous sections. The effect of neglecting off-diagonal annihilation rates in the determination of Ω SF h 2 yields an error of about 9% to 3% for models 1 − 5, underestimating the true result. In case of model 6 the Ω SF h 2 results obtained when neglecting or correctly including off-diagonal annihilation rates happen to agree. This can be understood from the Sommerfeld suppressions in the two single-charged sectors that arise when correctly accounting for off-diagonal annihilation rates and that can lead to a partial compensation of enhancements encountered in the neutral sector. While there is no suppression effect if off-diagonal annihilation rates are neglected, also the Sommerfeld enhancements in the charge-neutral sector are milder in that case, see for instance Fig. 5. Relic density results with and without off-diagonal annihilation rates can therefore accidentally agree, as it happens for model 6. If corrections from heavier states in the last potential loop were not included in the calculation of the relic abundance, the corresponding result would be larger by 2% for model 1 to 6% for model 6 as compared to the Ω SF h 2 values quoted in Tab. 1. As expected, the latter effect gains importance as the mass splitting of the heavier states to the higgsino-like χ 0 1,2 and χ ± 1 becomes smaller; while the wino-like states χ 0 3 , χ ± 2 in model 1 are rather heavy (m ∼ 3.3 TeV), these states have a mass of about 1.6 TeV in case of model 6.

ID µ/GeV
For models 7 − 9 with mixed wino-higgsino χ 0 1 , where the wino content increases with higher model ID, Fig. 12 shows a reduction of Ω SF h 2 the larger the wino admixture of the χ 0 1 . The ratio Ω SF h 2 /Ω pert h 2 ranges from ∼ 0.78 for model 7 over ∼ 0.66 for model 8 and gives ∼ 0.55 in case of model 9. In the region of mixed wino-higgsino χ 0 1 , where the masses of the states χ 0 1,2,3 , χ ± 1,2 lie close to each other, more two-particle states have been considered exactly in the multi-state Schrödinger equation. Precisely, the set of neutral χχ-states considered in the Schrödinger equations for model 7 comprises the seven states χ 0 while for model 8 the state χ 0 2 χ 0 3 is included in addition, and for model 9 only the six states χ 0 1 χ 0 1 , χ + 1 χ − 1 , χ 0 1 χ 0 2 , χ 0 1 χ 0 3 , χ ± 1 χ ∓ 2 are treated exactly in the neutral sector. While in the three models 7 − 9 (particularly in the neutral sector), the mutual interaction among a large number of channels is solved through the Schrödinger equations, it is mainly the larger wino fraction of the χ 0 1 that controls the increasing relevance of the Sommerfeld enhancements on the final relic abundance. While the wino fraction of the χ 0 1 in model 7 is 20% it becomes 46% for model 8 and finally reaches 83% in case of model 9. The larger wino admixture of both the χ 0 1 and χ ± 1 states also manifests itself in the decreasing mass splitting δm χ + 1 between these two states, ranging from 0.971 GeV (model 7) over 0.601 GeV (model 8) to only 0.266 GeV (model 9). A larger wino component of the χ 0 1 implies stronger potential interactions between the co-annihilating channels, in particular the χ 0 1 χ 0 1 and χ + 1 χ − 1 , where the latter is composed of χ ± 1 states with similar wino fraction as the χ 0 1 . The stronger potential interactions finally lead to a more pronounced Sommerfeld enhancement effect for models with larger wino admixture to the χ 0 1 state. Neglecting off-diagonal annihilation rates would lead to a result enhanced by 5% (model 7), 10% (model 8) and 14% (model 9) with respect to the actual Ω SF h 2 values given in Tab. 1. On the other hand, corrections to the Sommerfeld-enhanced rates from heavy χχ-states in the last potential loop reduce the final relic abundances Ω SF h 2 for models 7 − 9 by around 2 − 4%. The latter reduction is not as large as for model 6, despite the fact that the mass differences in models 7 − 9 are smaller. This is simply because there are less heavy channels contributing perturbatively now, as more χχ-states have been considered exactly in the Schrödinger equation.
Finally let us consider the subclass of wino-like χ 0 1 models with IDs 10 − 13. Here we account for Sommerfeld effects on the annihilation rates for χχ-states built from the wino-like χ 0 1 and χ ± 1 particles. The Schrödinger equations in the neutral sector for models 10 − 13 hence contain the two states χ 0 1 χ 0 1 and χ + 1 χ − 1 only. The models can be further subdivided into two groups with different impact of Sommerfeld enhancements: in case of models 10 and 11, Ω SF h 2 is significantly reduced by around 60% with respect to the result from a perturbative calculation. This happens to be the strongest reduction we find along the trajectory. The reason for the especially pronounced Sommerfeld-enhanced annihilation rates in case of models 10 and 11 can be attributed to the presence of a so called zero-energy resonance [14] in the χ 0 1 χ 0 1 annihilation channel: as already discussed, for velocities well below the χ + 1 χ − 1 threshold the enhancement in the χ 0 1 χ 0 1 system is controlled by the Yukawa potential due to electroweak W -exchange. As any short-ranged potential, a Yukawa-potential features a finite number of bound states. By varying the potential's strength and range it is possible to arrange for the presence of a bound state with (almost) zero binding energy [14] (see also [26]). In the presence of such a (loosely) bound state, the scattering cross section for incoming particles with very low velocities is strongly enhanced. This effect leads to O(10 4 ) enhancements in the χ 0 1 χ 0 1 channel for velocities below the χ + 1 χ − 1 threshold and eventually translates into the pronounced reduction of about 60% of the relic density. If off-diagonal annihilation rates were not taken into account, the Ω SF h 2 result would be larger by about 25% (model 10) and 23% (model 11), thus representing a rather large effect for both models: Off-diagonal annihilation rates are particularly important if the corresponding off-diagonal potential interactions are sufficiently strong. In wino-like χ 0 1 models, the only sector with relevant off-diagonal potential interactions is given by the two neutral states χ 0 1 χ 0 1 and χ + 1 χ − 1 in a 1 S 0 wave configuration. 5 For models 10 and 11, where the neutral χ 0 1 χ 0 1 channel experiences particularly large enhancements due to the presence of a (loosely) bound state resonance related to the off-diagonal W -exchange potential, also the impact of offdiagonal annihilation rates is therefore found to be significant. Regarding the corrections from heavier χχ-states treated perturbatively in the last potential loop, they are rather mild: Ω SF h 2 would be smaller by around 3% without this effect. Compared to model 6, where we found a corresponding 6% reduction in Ω SF h 2 , this suggests that the effect from heavier χχ-states in the last potential loop is most significant if these states are built from wino-like particles. The latter have in overall stronger (off-) diagonal annihilation rates compared to higgsino-like states with similar mass. Let us recall that the effect from heavier χχ-states in the last potential loop was at the per mil level in case of the pMSSM scenarios in Secs. 2 and 4 and around 1% for the higgsino-like scenario in Sec. 3, because heavier states were essentially decoupled in these models, opposed to the case for the models on the higgsino-to-wino trajectory.
At last, for models 12 and 13 we find a reduction of Ω SF h 2 relative to Ω pert h 2 of roughly 50% in both cases. This is still larger than the 40% reduction arising in case of the wino-like χ 0 1 pMSSM Snowmass model discussed in Sec. 2. To explain this effect note first that although the input value µ differs for models 12 and 13, this does not affect the parameters of the corresponding wino-like sectors. The masses of both χ 0 1 and χ ± 1 as well as their wino fractions are essentially the same in model 12 and 13, see Tab. 1. We can hence expect that the results for the χ 0 1 relic abundance calculation are very similar for both models. The presence of a zero-energy resonance in the χ 0 1 χ 0 1 annihilation channel is still noticeable for models 12, 13 -although it is less pronounced, as increasing the χ 0 1 mass moves us away from the exact resonance region. To conclude with the comparison to the wino-like χ 0 1 pMSSM Snowmass model in Sec. 2, recall that the mass of the wino-like χ 0 1 there was m χ 0 1 = 1650.664 GeV; in that case the Yukawa potential does not exhibit (almost) zero-energy bound states. Consequently no additional strong resonant enhancement takes place, such that in comparison to the wino-like models on the trajectory the Sommerfeld effect on the relic density is less prominent in Sec. 2, though still around 40%. Finally the calculated relic density Ω SF h 2 for both models 12 and 13 is increased by 17% and 16%, respectively, if off-diagonal annihilations are neglected. Not including the one-loop effects from heavy χχ-states increases the corresponding results for Ω SF h 2 in Tab. 1 by 2% in both cases.
6 Mixed wino-higgsino χ 0 1 As our framework allows for the first time to investigate Sommerfeld enhancements of χχ co-annihilations in scenarios with a χ 0 1 in an arbitrary wino-higgsino admixture, let us discuss here in more detail the mixed wino-higgsino χ 0 1 trajectory model with ID 8 considered in the previous section. Recall from section 5 that the neutral sector of the Schrödinger equation for this model is composed of the eight states χ 0 Fig. 13 shows the enhancements (σ SF v)/(σ pert v) in the two neutral channels χ 0 1 χ 0 1 and χ + 1 χ − 1 with (solid lines) and without (dashed lines) off-diagonal annihilation rates. The characteristic velocity-independent enhancement from the W -exchange Yukawa potential in the low velocity regime of the χ 0 1 χ 0 1 channel is visible, as well as the Coulomb-type 1/v χ + 1 enhancement for the χ + 1 χ − 1 system at very low velocities. Long-range potential interactions, although stronger than in case of higgsino-like χ 0 1 models are still weaker than in case of a wino-like set of states χ 0 1 , χ ± 1 ; as a consequence enhancement factors of O(1 − 10) result. We do not show (σ SF v)/(σ pert v) for the remaining six neutral twoparticle states in Fig. 13, but the resonance regions below their corresponding on-shell production thresholds can be seen as small enhancements in the χ 0 1 χ 0 1 and χ + 1 χ − 1 channels. The threshold for χ 0 1 χ 0 2 production opens at v LSP /c ≃ 0.18 but is hardly visible in the curves for channels χ 0 1 χ 0 1 and χ + 1 χ − 1 in Fig. 13. We can notice a broader (smoothedout) resonance region around v LSP /c ≃ 0.25, which comprises the thresholds for the four channels χ 0 2 χ 0 2 , χ 0 1 χ 0 3 and χ ± 1 χ ∓ 2 . Finally, the χ 0 2 χ 0 3 threshold shows up at v LSP /c ≃ 0.30. The enhancements for these channels, not shown in Fig. 13, are somewhat smaller than for the cases of χ 0 1 χ 0 1 and χ + 1 χ − 1 . Eventually, at v LSP /c ≃ 0.35 the threshold for on-shell production of the χ 0 3 χ 0 3 state is visible in the χ + 1 χ − 1 channel. The χ 0 3 χ 0 3 state is among the heavy states considered perturbatively in the last potential loop for the calculation  of the annihilation rates of the channels treated exactly in the neutral sector. Note that apart from the bino-like χ 0 4 state, which is very heavy (m χ 0 4 ∼ 19 TeV) and -being bino-like -couples very weakly to the gauge bosons and the other χ 0 /χ ± species, all χ states in the neutralino/chargino sector are relevant in co-annihilation reactions for the χ 0 1 relic abundance calculation of model 8. The thermally averaged effective annihilation rates σ eff v (x) including (upper solid (red) line) and neglecting (dashed red line) off-diagonal rates in the Sommerfeld-enhanced cross sections are depicted in the upper panel of Fig. 14. The corresponding perturbative result is given by the lower solid (blue) curve. The perturbative annihilation rates of two-particle states χχ heavier than the χ 0 1 χ 0 1 pair are larger than the perturbative rate of the latter, leading to a drop in the perturbative σ eff v (x) curve after decoupling of the heavier co-annihilating χχ states. As can be already inferred from Fig. 13, the effective rate including Sommerfeld enhancements turns out to be larger than the corresponding perturbative result by factors of at most O(1−3) in the x range x = 10 . . . 10 3 relevant to the relic abundance calculation. These enhancements finally give rise to the behaviour of the ratio of yields Y /Y pert shown in Fig. 14, lower panel. Including Sommerfeld corrections on the co-annihilation rates leads to a reduction of the relic density by 34%. For this model the effect of neglecting off-diagonal rates in the relic abundance calculation turns out to be milder than in the wino-like χ 0 1 models: with the off-diagonal entries we get Ω SF h 2 = 0.0791 while neglecting these would lead to a value larger by 10%.
It is interesting to analyse the impact on the calculated relic abundance Ω SF h 2 when  neutral χχ-states χ 0 1 χ 0 Table 2: χχ-states and corresponding masses M χχ in model 8, ordered according to their electric charge, that are relevant in the calculation of the χ 0 1 relic abundance Ω SF h 2 . Two-particle states involving the bino-like neutralino χ 0 4 are not shown. As their masses M χχ lie above the scale of 20 TeV, they are irrelevant in the calculation of Sommerfeld enhancements to the lighter χχ-channels and in the determination of the χ 0 1 relic abundance. The vertical double lines separate the states with masses below 3762 GeV and above 3893 GeV.
the number of channels included in the multi-state Schrödinger equation is changed, or the number of heavier states contributing to corrections from the last potential loop is varied. Let us recall that the results presented so far in this section correspond to calculations where all χχ-states with masses below M max = 3893 GeV are treated exactly in the Schrödinger equation, 6 while the remaining heavier states are included only at treelevel and in the last loop near the annihilation vertex in the Sommerfeld-corrected rates of the lighter states. Further we have considered δm 2 corrections in the potentials for the channels included in the Schrödinger equation but not in the approximate treatment of the heavier states (see [23] for details on these corrections). In order to compare the cases where the number of channels treated in the Schrödinger equation is changed, we neglect these δm 2 corrections in the potentials throughout in the following, so that all cases are computed with the same potential. We calculate Ω SF h 2 for the cases of M max = 3762 GeV and M max = 3893 GeV, corresponding to v max = 0.2 and 1/3, as well as for M max = ∞.
In the latter case all χχ-channels are taken into account in the Schrödinger equation. To investigate the accuracy of the approximate treatment of heavier states in the last potential loop compared to the case where these states are accounted for exactly in the Schrödinger equation, we introduce the variable M cut ≥ M max . χχ-states with a mass larger than M cut are ignored completely. States with mass below M max are included in the Schrödinger equation exactly, while those with mass between M max and M cut are treated approximately through the one-loop corrections in the last potential loop. The relevant χχ-states together with their masses are given in Tab. 2, from which the number of exactly and approximately treated states in each charge sector for each of the cases covered in Tab. 3 can be read off. The results on Ω SF h 2 that we obtain for our three  Table 3: Relic abundances Ω SF h 2 in trajectory model 8 with a different number of channels accounted for in the Schrödinger equation and with a different number of heavy χχ-states treated approximately in the last potential loop. Two-particle channels χχ with masses below M max are included in the Schrödinger equations. One-loop corrections of heavier χχ-channels with masses between M max and M cut are accounted for, while all χχ-channels heavier than M cut are ignored. All results are derived neglecting δm 2 corrections in the potentials.
choices for M max and for M cut set to M cut = 3762 GeV, 3893 GeV and M cut = ∞ are collected in Tab. 3. Let us first discuss the Ω SF h 2 values on the diagonal of Tab. 3, which display the effect of increasing the number of states in the Schrödinger equation while ignoring one-loop corrections from heavier states. Expectedly Ω SF h 2 decreases the larger M max . There are more χχ-channels for which Sommerfeld enhancements on their individual annihilation cross sections are taken into account. This leads to an increase of the thermally averaged effective rate σ eff v entering the Boltzmann equation, which in turn decreases the relic abundance. By increasing M max by the steps indicated in Tab. 3 the resulting Ω SF h 2 is reduced by 5% and 2% respectively. The effect on Ω SF h 2 from more channels in the Schrödinger equations is rather mild as compared to the 33% reduction with respect to the tree-level relic density. 7 The milder reduction mainly derives from the fact that the Sommerfeld enhancement of the heavier channels' cross sections is less pronounced than in case of the most relevant lighter channels χ 0 1 χ 0 1 , χ + 1 χ − 1 and χ 0 1 χ ± 1 . Further, as noted previously, the heavier χχ-channels enter the thermally averaged rate σ eff v with a Boltzmann suppression factor such that their contribution is generically sub-dominant, unless the individual rates are particularly enhanced. The main effect that leads to the respective 5% and 2% change of Ω SF h 2 comes from the slight increase of the Sommerfeldenhanced cross sections of the dominant light channels χ 0 1 χ 0 1 , χ + 1 χ − 1 and χ 0 1 χ ± 1 when more states appear in the potentials of the Schrödinger equations.
Let us now consider the reduction of Ω SF h 2 for fixed M max and increasing M cut . This happens because the effect of heavier channels amounts to a positive correction to the Sommerfeld-enhanced cross sections: the dominant potential interactions are attractive, such that the heavier states in the last potential loop typically give an additional positive contribution. For instance we find a significant reduction of Ω SF h 2 by 5% from 0.858 to 0.817, when for M max = 3762 GeV the value of M cut is increased from 3762 GeV to 3893 GeV. This indicates that the newly added heavier states in the last loop give a large positive contribution to the Sommerfeld-enhanced cross sections of the χχ-states in the Schrödinger equation. When CPU considerations make the restriction to fewer states treated in the Schrödinger equation necessary, the approximate treatment of heavy channels should give a reasonable approximation to the case where these heavy channels are included fully in the Schrödinger equation. This is nicely confirmed by the numbers shown in Tab.3: when the states with mass between 3762 GeV and 3893 GeV are treated approximately, the reduction of Ω SF h 2 from 0.0858 to 0.0817 is very close to the value 0.0816 obtained from the exact treatment of all states with mass below 3893 GeV. The same observation holds for the comparison between the approximate treatment of all states with masses above 3762 GeV, Ω SF h 2 = 0.0804, and the exact result Ω SF h 2 = 0.0801. The agreement becomes even better when the the perturbative treatment involves only the heavier channels with mass above 3893 GeV.

Summary
In this work we presented a detailed investigation of Sommerfeld enhancements in the χ 0 1 relic abundance calculation for several popular models with heavy neutralino LSP in the general MSSM. Our analysis is based on the effective field theory formalism that we developed and described in [21][22][23]. This framework allows us to calculate the χ 0 1 relic abundance consistently including Sommerfeld-enhanced neutralino/chargino co-annihilation rates, taking off-diagonal rates into account and accounting for many nearly mass degenerate co-annihilating two-particle states. We focused on three benchmark models with wino-, higgsino-and bino-like χ 0 1 taken from [24] as well as on a set of DarkSUSY generated spectra interpolating between the cases of a higgsino-to a wino-like χ 0 1 spectrum. With the latter set we defined a "higgsino-to-wino" trajectory in the parameter space of the general MSSM. It is worth to stress that our work allows for the first time to investigate the Sommerfeld enhancement in neutralino/chargino co-annihilations for a mixed wino-higgsino χ 0 1 . In scenarios with wino-like χ 0 1 we find a pronounced effect from Sommerfeld enhancements on the calculated χ 0 1 relic abundances, whereas for higgsinolike χ 0 1 , the effect becomes milder. This is in agreement with previous investigations in the literature in the pure-wino and pure-higgsino limits. In general the relic abundance obtained including the Sommerfeld effect is reduced the more the stronger the wino admixture to the χ 0 1 . In addition, there are cases of particular pronounced effects related to the existence of loosely or zero-energy bound states in the spectrum of the model. We show that the precise value of the calculated relic density depends on the particular details of the spectrum, such that results from a study in the pure-wino or pure-higgsino χ 0 1 scenarios do not apply directly. It is interesting to note that Sommerfeld enhancements in the co-annihilating sector of a bino-like χ 0 1 can affect the result on Ω SF h 2 at the 10% level. This is found for a benchmark model with bino-like χ 0 1 and slightly heavier winolike χ ± /χ 0 states. The knowledge of precise mass splittings between the co-annihilating neutralinos and charginos is essential in the calculation of Sommerfeld-enhanced rates and will typically require the knowledge of spectra with a one-loop on-shell renormalised neutralino/chargino sector.
We used three pMSSM benchmark models as well as the set of models on our "higgsino-to-wino" trajectory in order to show the general features of Sommerfeldenhanced rates and their effect on the relic abundance calculation. The results demonstrate that it will be necessary to systematically include the Sommerfeld effect when MSSM parameter space constraints on heavy neutralino dark matter from direct and indirect searches as well as from collider physics are combined with the requirement to reproduce, or at least not exceed, the observed abundance of dark matter. A future project is the investigation of the parameter space of the general MSSM as regards the relevance of Sommerfeld enhancements in the relic abundance calculation. Our aim is to identify regions where the Sommerfeld effect is not necessarily as pronounced as in the previously studied wino limit but constitutes the dominant radiative correction. To this end a scan of the MSSM parameter space is prepared and our findings will be reported in future work.