Relic density of wino-like dark matter in the MSSM

The relic density of TeV-scale wino-like neutralino dark matter in the MSSM is subject to potentially large corrections as a result of the Sommerfeld effect. A recently developed framework enables us to calculate the Sommerfeld-enhanced relic density in general MSSM scenarios, properly treating mixed states and multiple co-annihilating channels as well as including off-diagonal contributions. Using this framework, including on-shell one-loop mass splittings and running couplings and taking into account the latest experimental constraints, we perform a thorough study of the regions of parameter space surrounding the well known pure-wino scenario: namely the effect of sfermion masses being non-decoupled and of allowing non-negligible Higgsino or bino components in the lightest neutralino. We further perform an investigation into the effect of thermal corrections and show that these can safely be neglected. The results reveal a number of phenomenologically interesting but so far unexplored regions where the Sommerfeld effect is sizeable. We find, in particular, that the relic density can agree with experiment for dominantly wino neutralino dark matter with masses ranging from 1.7 to beyond 4 TeV. In light of these results the bounds from Indirect Detection on wino-like dark matter should be revisited.

The so-called "WIMP miracle" is the observation that a thermally produced, stable, massive particle χ with electroweak interactions (WIMP) naturally accounts for the observed dark matter relic density, if its mass is of order of the electroweak or TeV scale. Indeed, on adding a fermionic SU (2) triplet to the Standard Model (SM), its tree-level pair annihilation into electroweak gauge bosons yields Ω cdm h 2 = 0.1188 for m χ ≈ 2.2 TeV. Such models provide attractive dark matter (DM) candidates due to their minimal particle content [1], but further model building is required to explain why the mass of the χ particle should be close to the electroweak scale. The minimal supersymmetric standard model (MSSM) is a prime example of a model where the DM particle mass is tied to the electroweak scale by the desire to temper the quantum corrections to the Higgs. The underlying symmetry principle then leads to a proliferation of particles and interactions, allowing for different successful DM candidates. Within this context thermal scenarios with DM masses below 1 TeV require additional mechanisms or accidental degeneracies, e.g. resonant annihilation or coannihilation, in order to avoid overproduction in the early Universe. These are also becoming somewhat constrained by LHC and dark matter searches, for recent analyses see refs. [2][3][4][5]. Scenarios with heavier dark matter interpolate to minimal models, since the supersymmetric particles form approximate electroweak multiplets (except for degeneracies). In particular, when the lightest supersymmetric particle ("wino") is the partner of the electroweak gauge bosons, the model is similar to the minimal triplet model, but modifications arise due to the mixing with the Higgsino and bino states, as well as the interactions with sfermions. It is this "wino-like" region of the MSSM parameter space, which we focus on in this paper. The wino-like region deserves special attention, since the DM relic density cannot be calculated reliably from the tree-level annihilation cross section. Loop effects from electroweak gauge boson exchange are large in non-relativistic scattering before the annihilation of TeV-scale dark matter, and lead to the electroweak Sommerfeld effect [6,7], which is particularly strong in the wino-like region. This has been studied extensively in the purewino limit [6][7][8][9][10], which corresponds to the minimal triplet model. To be specific, in the analysis below we find that the observed relic density is attained at significantly larger mass m χ = 2.88 TeV when the Sommerfeld effect is accounted for, instead of 2.22 TeV at tree level, when M 1 = 3M 2 , µ = 2M 2 and the common sfermion mass M sf = 20 TeV, which corresponds effectively to the pure-wino limit. The Sommerfeld effect also displays a resonance at 2.33 TeV, where the relic density is reduced by a factor 3.9 relative to the computation based on the tree-level cross section. This highlights the importance of including the Sommerfeld effect in full MSSM calculations of the relic density in the wino-like region.
Away from the pure-wino limit, the lightest neutralino is a mixture of wino, Higgsino and bino eigenstates and interacts accordingly, which makes the computation of the Sommerfeld effect much more involved. This problem was first approached in refs. [9,11], however a framework that deals systematically with mixed states, multiple co-annihilating states and the corresponding off-diagonal reactions was only developed in refs. [12][13][14], which allowed the computation of the relic density including the Sommerfeld effect with a relative accuracy similar to state-of-the-art computations employing Born cross sections.

JHEP03(2016)119
This was studied in a number of models that interpolate from a pure-wino to a pure-Higgsino DM particle [15], but a detailed investigation of the MSSM parameter space was left for the future.
We report on this investigation in the present work, focusing on the wino-like region of the full MSSM. We note that in this region the Sommerfeld effect is not a small correction and should be included in any reliable relic density computation and in particular when the relic density is correlated with other observational constraints. The most important is from indirect dark matter searches. For instance, the thermal pure-wino scenario is often said to be excluded (barring some astrophysical uncertainties, see ref. [16]) by the non-observation of a photon line signal from the Galactic Centre [16][17][18]. Other search channels, especially the cosmic ray antiprotons and the diffuse gamma rays from dwarf spheroidal galaxies, also start to give competitive limits [16,[19][20][21].
This conclusion need not hold in the full MSSM, when the mixed nature of wino-like dark matter is taken into account. The framework adopted here follows refs. [12][13][14], with several improvements applied relative to ref. [15]. We now include the running of the electroweak couplings from the electroweak to the dark matter scale, and use the exact one-loop neutralino and chargino on-shell masses to compute the mass splitting, which is important in the resonance region. We justify neglecting thermal effects due to the fact that the freeze-out happens at temperatures close to the electroweak scale. On the practical side, a considerable speed-up of the numerical evaluation has been achieved, which now allows a systematic investigation of the relevant MSSM parameter space in the wino-like region.
The outline of the paper is as follows: in section 2 we define the ranges of the parameters of the phenomenological MSSM, and discuss the theoretical and observational constraints we apply to select viable models. We further briefly summarize the computation of the Sommerfeld correction, the implementation of mass splittings and the running of the electroweak coupling. The set-up is rather general, but the present version does not include sfermion-neutralino/chargino potentials and hence excludes models with sfermion co-annihilation, as well as s-channel resonant annihilation, in which case the annihilation process is not short-distance. Section 3 contains our main results. Here we show and discuss, in order, the dependence of the relic density and the relative importance of the Sommerfeld effect on the sfermion masses M sf (all assumed degenerate for simplicity), on the heavy MSSM Higgs bosons, further on the Higgsino admixture via the difference µ−M 2 of the Higgsino and wino mass parameters of the MSSM, and similarly on the bino admixture. We shall see that away from the pure-wino limit, the observed relic density is obtained for a wide range of wino-like dark matter particle masses and we quantify and explain the parameter dependence. We further study the dependence on other MSSM parameters, which generally turns out to be minor, except in the vicinity of the Sommerfeld resonance. We summarize in section 4.
The investigation of thermal effects is contained in appendix A. We consider the temperature dependence of the electroweak gauge boson masses, which in turn affects the range of the electroweak Yukawa potential, and of the neutralino-chargino mass difference. The dependence arises from the temperature-dependent Higgs vacuum expectation value and the one-loop self-energies. Despite the fact that freeze-out may begin in the symmetric JHEP03(2016)119 phase of the electroweak interactions, where the Higgs field has no expectation value and the thermal effects are large, we find that the impact on the relic density is negligible within other uncertainties. We explain why previous work [8,9] overemphasised the effect.

MSSM definition and parameter ranges
We are interested in exploring the parameter space of the CP-conserving, minimal flavour violating MSSM defined at the electroweak scale. Within this space we focus primarily on the calculation of the relic abundance of neutralino dark matter in the close-to-wino region, which only depends strongly on a subset of parameters. Clearly a central role is played by those parameters describing the chargino and neutralino sector: the bino mass M 1 , the wino mass M 2 and the Higgsino parameter µ. 1 The tree level mass matrix for the charginos is given by where s β /c β ≡ sin β/ cos β, tan β being the ratio of the vevs of the two MSSM Higgs doublets, and m W is the mass of the W boson. The mass matrix for the neutralinos is given by where s W ≡ sin θ W , c W ≡ cos θ W for the Weinberg angle θ W and m Z is the mass of the Z boson. On diagonalising the hermitian squares of these matrices one obtains the values of the masses of the charginos and neutralinos mχ+ i and mχ0 j respectively, numbered i = 1, 2, j = 1, . . . , 4 in increasing order.
We concentrate on the region where the lightest supersymmetric particle (LSP) mass is at the TeV scale, as this is where the wino-like neutralino can provide the correct thermal relic density and the electroweak Sommerfeld effect is non-negligible. Here we assume that either the bino, Higgsino or both are much heavier than the wino. We can study the mixing angles of the wino with the bino or the Higgsinos and the resulting mass eigenstates by expanding in m Z /µ etc. The mixing as well as the mass difference between the lightest chargino and neutralino can play an important role in determining the size of the Sommerfeld enhancement. In the region where the bino is decoupled, provided that m W |µ| − M 2 , the splitting δmχ+ 1 ≡ mχ+ 1 − mχ0 1 is given by

JHEP03(2016)119
is close to −1, the Higgs mass constraint does not have much impact on our results. For M sf > 6 − 7 TeV the effect of the neglected corrections to the Higgs mass could therefore be compensated by a change in X t , leaving the relic density unaltered.
ρ parameter. We require that the value of ∆ ρ computed in the MSSM [31] does not exceed two standard deviations from the SM expectation [32]: Since the SUSY contribution can only be large when the mass splitting in the sfermion SU(2) doublets is large, and in the scenario we consider all the sfermion doublets are nearly degenerate, it does not have a significant effect on our parameter space.
b → sγ. In general MSSM scenarios this branching ratio provides a strong constraint, as the contribution from broken SUSY is generically large, while the SM prediction is compatible with measurement. The experimental [33] and SM theory [34] values, with the corresponding uncertainties, we use are The SUSY contribution ∆B B → X s γ is computed with FeynHiggs and the implemented criterion reads There are three classes of diagrams which contribute to b → sγ in the MSSM: these are diagrams involving either charged Higgs bosons, charginos or gluinos. The first always interfere constructively with the SM contribution, and decouple as the Higgs mass increases beyond the TeV scale. The chargino contribution can take either sign, depending on the sign of µ and A t , but also decouples with increasing |µ| and M 2 . At the scales that are relevant to this study, i.e. above 1 TeV, in general the MSSM contribution lies within the uncertainties. B s → µ + µ − . The correction to B s → µ + µ − from SUSY should also lie within the errors from the experimental measurement and the SM calculation. To this end, we check whether the result of the calculation in the MSSM [31] is consistent with the combined CMS and LHCb result, (2.9 ± 0.7) × 10 −9 [35]. The 3 sigma error on the experimental result is added to the uncertainty on the theoretical result in quadrature, where the updated SM prediction is (3.56 ± 0.30) × 10 −9 , using latest values on the B 0 s lifetime and relative B 0 s decay width difference [35,36]. We note that as we consider the wino-like region with masses of the LSP of O(TeV), and masses of the heavy Higgs bosons also of O(TeV), this constraint does not have much influence on our parameter space. Another related constraint is of course the branching ratio of B → τ ν, measured precisely at the B-factories [37,38]. However, we do not consider this constraint as the parameter space of interest in our analysis, in particular the large masses of the charged Higgs bosons and values of tan β, do not result in MSSM contributions beyond the combined experimental and theoretical uncertainty [39].

JHEP03(2016)119
g µ − 2. The experimental and SM theory values adopted for the muon anomalous magnetic moment, a µ = gµ−2 2 , are given by [32], a exp µ = (1165920.91 ± 0.63) × 10 −9 , a SM µ = (1165918.03 ± 0.48) × 10 −9 , and we require that ∆a µ , the MSSM contribution, satisfies where for the error on the difference between experimental and SM values we take (2.14) This means that we do not insist that the MSSM contribution explains the deviation between the experimental and SM theory values. Note that as the SUSY contribution is proportional to tan β and inversely proportional to the square of the masses of the sparticles, it is typically strongly suppressed in the region of interest where M sf lies at the TeV scale.

Theoretical constraints
Higgs potential. Theoretical consistency demands that the scalar potential is free from charge and/or colour breaking minima (CCB). For the tree-level scalar potential in the MSSM, the corresponding criteria read [40,41] One can always choose the trilinear couplings low enough such that the CCB constraint is satisfied without altering the nature of the neutralino.
s-channel resonances. Our calculation relies on the factorisation of the annihilation cross section into the short-range tree-level annihilation and the long-range potential interaction. However, this factorization does not hold in the case that the final light particles are produced through an s-channel propagator which is resonant, as such a contribution cannot be attributed to the short-distance part of the annihilation. Therefore, we need to exclude regions of parameter space where this may occur. In the MSSM this means that we need to avoid s-channel resonances through the Higgs bosons, and to be conservative we assume that the masses of the heavy Higgses lie outside the interval It follows that in this work we are not in a position to study the H-and A-funnel regions [42,43].

Cosmological and direct DM detection constraints
In choosing suitable points to calculate the Sommerfeld effect on the relic density, we insist that certain basic constraints are fulfilled. First we require that the lightest neutralinoχ 0 1 is the LSP. We further insist on compatibility with Direct Detection bounds. The details of how these conditions are imposed is described in this subsection. We choose not to include any limits coming from Indirect Detection experiments or measurements of the CMB, as although these may be relevant they are subject to large systematic uncertainties and their discussion goes beyond the scope of this work; we plan to address such constraints in the future.
Direct detection. We require that the DM-nucleon spin-independent cross section σ SI is less than twice the LUX limit [44]. The theoretical prediction of this cross section within the MSSM is obtained using micrOMEGAs. The spin-independent cross section is sensitive to the Higgs exchange between the LSP and the quarks of the nucleon. The interaction with the Higgs relies on the LSP containing both gaugino and Higgsino components, and therefore this constraint is most relevant for the scenarios we study where |µ| ∼ M 2 .
Note that the limits of the spin-dependent cross section coming from Direct Detection experiments and neutrino signals from the Sun are always much weaker than those coming from spin-independent results for the scenarios we are interested in here.

One-loop mass splittings
The differences in mass between the LSP and the heavier neutralinos and charginos can have an effect on the relic density. The most relevant case is the small mass difference between the lightest chargino and neutralino state,χ + 1 andχ 0 1 , respectively. In order to be consistent with the accuracy of the rest of the calculation we calculate these masses at one-loop. In doing so we adopt an on-shell renormalisation scheme, which is described here in brief. For further details we refer the reader to refs. [45][46][47][48][49].
The mass matrix in the chargino sector is renormalised via X → X + δX, where δX is defined by containing the renormalisation constants (RCs) for the wino parameter M 2 and Higgsino parameter µ, i.e. δM 2 and δµ. In addition, the matrix δX contains the RCs of c β and s β , i.e. δc β and δs β (which can be expressed in terms of δ tan β), and of the W boson mass m W , δm W . Definitions of and expressions for δ tan β and δm W can be found in ref. [49]. The neutralino mass matrix, Y , is renormalised in a similar manner via Y → Y + δY , where δY is defined in analogy to δX in eq. (2.17) and further contains RC of the bino parameter M 1 , δM 1 . In the on-shell scheme, we must fix the RCs δM 1 , δM 2 and δµ (as in e.g. ref. [46]) by requiring that three out of the total six physical masses of the charginos and neutralinos satisfy on-shell conditions, i.e. that the tree-level masses, mχ i , coincide with the one-loop renormalised masses, Mχ i = mχ i + ∆mχ i ,

JHEP03(2016)119
Scenario Particles on shell Table 2. Choice of particles whose masses are required to be on shell for the various scenarios corresponding to the possible orderings of M 1 , M 2 and µ that we consider.
Note that we define the coefficientsΣ The left-and right-handed vector and scalar coefficients,Σ L/R ij (p 2 ) andΣ SL/SR ij (p 2 ) of the renormalised self-energy are defined analogously. Expressions for the renormalised selfenergies can be found in e.g. ref. [49]. The mass shifts for the remaining three chargino and neutralino masses are therefore given by ∆mχ± i and ∆mχ0 j in eq. (2.18). For the calculation of these mass shifts we used the program FeynArts [50,51], together with the packages FormCalc [52] and LoopTools [52], using the model files presented in ref. [53].
The choice of which masses should be chosen on shell is non-trivial, as certain choices can lead to unphysical divergences when e.g. |M 1 | = M 2 or |µ| = M 2 , and we follow the prescription discussed in refs. [48,54] as follows to avoid this situation as far as possible. We therefore employ the NNC scheme, that is, two neutralinos and one chargino are chosen on-shell, of which the chargino should be wino-like, and the neutralinos should be bino and Higgsino-like. Note however that there is an ambiguity here given that there are two Higgsino-like neutralinos. In this work we are particularly interested in the region where the neutralino has a large wino component, i.e. M 2 < |M 1 |, |µ|, and may in addition contain a sizeable bino or Higgsino component. We therefore find that in order to obtain results free from scheme-dependent divergences, the choice of particles whose masses are required to be on shell should be made as in table 2. This corresponds to the Higgsino closer in mass to the wino being on shell. Note that when all three parameters are very close (< 0.1% splittings) the situation may arise that the ordering of the neutralinos changes, and one should exercise caution in these regions. This has been accounted for in the code.

Running couplings
Due to the multi-scale nature of the considered problem, the running of the coupling constants has to be treated consistently. In different parts of the calculation the couplings should be taken at a different energy scale Q, in particular Q = m Z for the potential interactions, Q = m LSP for the mass splittings in the neutralino/chargino sector and Q = 2 m LSP for the short-range annihilations.
We perform the running in the unbroken SU(2) L × U(1) Y theory, since most of the running occurs above the electroweak scale. The starting values of the SU(2) L and U(1) couplings at Q = m Z are taken as α 2 (m Z ) = 0.034723 and α 1 (m Z ) = 0.009986, respectively. Since the short-range annihilation is evaluated at tree-level, we run the couplings to JHEP03(2016)119 Q = 2 m LSP with the one-loop renormalisation group equation. In the computation of the one-loop mass splittings discussed above, the couplings are evaluated at Q = m LSP . The energy range from m Z to 2m LSP that we are interested in can be divided into five regions 4 where the beta functions are constant, delimited by the scales M A , |M 1 |, M 2 , |µ|, at which we decouple respectively the heavy Higgs doublet Φ 2 , the bino, the wino, and the Higgsinos. At the required level of accuracy there are no threshold effects to be considered and the leading order beta function β 0,i at a scale Q is given by where T

(i)
R are the generators of the group i in the representation R and the three terms correspond respectively to gauge bosons (always in the adjoint representation A), fermions, and scalars. The sums extend only to particles with mass smaller than Q, and the contributions are listed in table 3.

Annihilation matrix implementation
The rate at which neutralinos and charginos annihilate into the (light) standard model particles in the early Universe is a necessary input for the calculation of the present-day amount of dark matter. For a given two-particle stateχ iχj ≡ [χχ] a formed out of two neutralino or chargino species, the annihilation rate including long-distance Sommerfeld corrections can be parametrised as [14] 4 We neglect the sfermions' contribution to the beta functions, when the sfermion mass lies between mLSP and 2mLSP. The error introduced in this way is small: for a 2.5 TeV LSP the total running of α2 from mZ up to 2mLSP is around 5-6%, and the maximum contribution from sfermions (when they are all decoupled at their smallest allowed mass 1.25 mLSP) is only 0.4%.

JHEP03(2016)119
. . , the relative momentum of the annihilating particles in their centre-of-mass frame, with M a , µ a the total and reduced mass, respectively, of the two-particle state. The quantitiesf ab ( 2S+1 L J ),ĝ ab ( 2S+1 L J ), . . . are the absorptive part of the Wilson coefficients of local four-fermion operators which reproduce the short-distance annihilation of the chargino and neutralino pairs into SM and light Higgs final states in the non-relativistic EFT framework [12][13][14]. They were determined by matching the tree-level MSSM amplitudes for the process [χχ] a → X A X B → [χχ] b with SM and Higgs intermediate states X A X B in refs. [12,13]. 5 The definition of the various Wilson coefficients appearing in eq. (2.21) can be found in ref. [14]. The Sommerfeld factors S a [. . . ] in eq. (2.21) account for the long-distance interactions of the two-particle states prior to the short-distance annihilation. Details on the computation of these factors are given below. The tree-level annihilation rate with no long-distance corrections is readily recovered by setting all the Sommerfeld factors in eq. (2.21) to one. The tree-level annihilation cross section thus obtained depends only on the diagonal entry of the Wilson coefficients corresponding to channel [χχ] a , i.e.f aa ( 2S+1 L J ),ĝ aa ( 2S+1 L J ), . . . . As shown in eq. (2.28) below, the computation of the Sommerfeld factors also requires knowledge of the off-diagonal terms,f ab ( 2S+1 L J ),ĝ ab ( 2S+1 L J ), . . . , with a = b, since the interference of loop diagrams where the two-particle states that undergo short-distance annihilation are different are accounted for in the Sommerfeld-corrected cross section. A word on the notation for labelling the two-particle states is relevant here. The two-particle statesχ iχj formed out of charginos and neutralinos are denoted by a single label a = 1, . . . N |Q| , where N |Q| is the total number of states (channels) for each electriccharge sector, |Q| = 0, 1, 2, corresponding to neutral (χ 0χ0 ,χ +χ− ), single-charged (χ 0χ± ) and double-charged (χ ±χ∓ ) sectors. If all four neutralinos and the two charginos are considered, in the charge-0 sector the single label runs over the 14 different states whereas in the charge ±1 sectors we have 8 channels each, and just three each in the charge ±2 sectors, The coefficientsf ab ( 2S+1 L J ) for each partial wave can then be considered as the entries of a matrix whose dimension is equal to the number of channels in each sector. Since the JHEP03(2016)119 such annihilation matrices turn out to be hermitian. The computation of each of the annihilation matrices appearing in the annihilation cross section formula (2.21), requires the evaluation of 105, 2×36 and 2×6 distinct entries for neutral, single-and double-charged sectors, respectively. Ten of such matrices are needed for a complete calculation of the Sommerfeld-corrected annihilation cross section including O(v 2 ) corrections (see ref. [14] for details on this), making up a total number of 1890 independent entries. In the CP-conserving case, the annihilation cross sections of the charged-conjugated sectors,χ 0χ+ andχ 0χ− ,χ +χ+ andχ −χ− , become equal, and the number of independent annihilation matrix entries is reduced to 1470. A code to obtain the analytic results for the entries of the annihilation matrices at O(α 2 2 ) in the MSSM has been developed following the conventions and recipes of refs. [12,13,55]. The expressions account for the sum of all possible X A X B exclusive states with X A/B being a SM particle (including the light Higgs) or heavy MSSM Higgs (the mass of the state X A X B must however be smaller than 2m LSP ). 6 For the neutral, single-and doublecharged sectors, the number of exclusive final states is 31, 16 and 3, respectively, including the possible heavy Higgs final states; a complete list can be found in appendix A of ref. [12]. Despite coming from the product of tree-level amplitudes, the analytic expressions for the Wilson coefficients are very large, which is traced back to the fact that there are several diagrams with different topologies and/or virtual intermediate particles contributing to a given exclusive state, and because of the non-relativistic expansion performed. Recall as well that we keep the general dependence on all MSSM parameters in the coefficients. The numerical evaluation of all matrix entries for a given MSSM parameter set is done using precompiled functions within Mathematica, taking on average approximately 300 sec of CPU time. If only the annihilation matrices necessary for the leading-order cross section,f ab ( 1 S 0 ) andf ab ( 3 S 1 ), are evaluated, the cost in CPU time reduces to less than 40 sec per model.
We should mention here a modification of a part of the analytic expressions for the Wilson coefficients given in refs. [12,13] that we have implemented in the present code. The Wilson coefficients obtained in refs. [12,13] is the average of the masses of the two-particle states taking place in the short-distance part of the annihilation process (for diagonal reactions, a = b, this is just an expansion around the [χχ] a threshold). When the annihilation proceeds through s-channel boson exchange, such an expansion implies for the boson propagator (with generic mass m φ ) that up to linear terms in √ s − M . The first term on the right-hand side of eq. (2.25) contributes to the leading-order Wilson coefficients, whereas the second goes to the S-wave v 2 -suppressed ones. For heavy Higgs exchange, the following problem may arise: once JHEP03(2016)119 radiative corrections are included, any (virtual) states a, b can participate in the shortdistance part, such that we can find a situation where M gets very close to the Higgs mass m φ ≈ M A , producing arbitrarily large contributions in the right-hand side of eq. (2.25). Those resonance contributions are spurious, since the annihilating cross section of the external state [χχ] I in the non-relativistic regime should be expanded for energies close to the mass of that state, i.e. around √ s = M I , which produces terms from s-channel contributions proportional to 1/(M 2 I −m 2 φ ) instead of those in eq. (2.25). For the relevant co-annihilation channels, the latter terms cannot become resonant in our analysis because we have explicitly excluded Higgs masses inside the range [1.7 mχ0 1 , 2.3 mχ0 1 ], see eq. (2.16). Therefore, the problem of spurious resonances is absent if we have a set of annihilation matrices for each co-annihilation channel I where the s-channel propagators have been expanded around √ s = M I . In practice, that solution is unfeasible, since the number of co-annihilations channels in a mixed scenario can be rather large and evaluating several annihilation matrices would increase the required CPU time beyond reasonable limits. We can adopt, however, another solution that avoids the occurrence of spurious resonances in s-channel propagators that only requires minimal changes in the Wilson coefficients obtained in refs. [12,13]. It amounts to modifying the expanded s-channel propagators from the Wilson coefficients such that they correspond to their expansion around √ s = 2mχ0 1 , regardless of which is the external co-annihilating state. We note that since the relevant channels that are included in the long-distance radiative corrections are very close in mass (see next section), the differences between the annihilation amplitudes expanded around 2mχ0 1 or around any of the other masses of the co-annihilating states are in any case negligible, and the suggested prescription is a very good approximation. The necessary modifications can be immediately read off by rewriting the right-hand side of eq. (2.25) using M = 2mχ0 where we have dropped terms of second order in the small quantities (M − 2mχ0 1 ) and ( √ s−M ). We notice that the dependence on M cancels out in the second line of eq. (2.26), and the resulting expression matches the expansion of the Higgs The replacements that have to be performed in the Wilson coefficients of refs. [12,13] thus read: (2.27)

JHEP03(2016)119
The factor of 2M in front of the square of a scalar propagator gets replaced by 4mχ0 1 in the S-wave v 2 -suppressed Wilson coefficients to get exactly the form in the right-hand side of eq. (2.26), though the difference between both expressions is formally of higher order. In self-energy contributions, the replacement (2.27) in the scalar propagator of LO Wilson coefficients produces an O(v 4 ) term from the product of the right and left s-channel propagators in the diagram, which is consistently dropped in our code in order to keep the expansion of Wilson coefficients to O(v 2 ) everywhere.

Sommerfeld-corrected cross section
The annihilation cross sections for the processes [χχ] a = χ i χ j → X, eq. (2.21), are computed by multiplying every term in the partial wave expansion of the Born cross section by its specific Sommerfeld factor When the Sommerfeld factors are neglected, eq. (2.21) reproduces the Born annihilation cross section including O(v 2 ) terms. The Sommerfeld factors are computed by solving the Schrödinger equation for a system of coupled two-particle states with the leading-order Yukawa and Coulomb potentials generated by the exchange of electroweak gauge bosons, Higgs bosons 7 and the photon. For further details including notation, we refer to ref. [14]. The calculation can be done separately in the sectors of two-particle states with different electric charge 0, ±1, ±2. Since we restrict ourselves to the CP-conserving MSSM, the annihilation cross sections for the negatively charged two-particle states are identical to the corresponding positively charged ones, and do not have to be calculated explicitly.
In every charge-sector, the Sommerfeld factors are computed for all two-particle states with mass less than 1.2 × 2m LSP unless the number of such states is larger than four, in which case the four lightest two-particle states are selected. For the other, heavier twoparticle states, we employ the Born cross sections. Furthermore, in the computation of the Sommerfeld factor we include the light states (at most four) exactly in the solution of the Schrödinger equation and the others approximately in the last loop near the annihilation vertex as described in [14]. The mass cut at 1.2×2m LSP is motivated by the fact that heavier states are either strongly Boltzmann-suppressed and irrelevant for freeze-out or they are sufficiently off-shell within the ladder diagrams to not contribute substantially to the Sommerfeld effect of the lighter states. The restriction to at most four light states is motivated by CPU considerations, since the time needed for the matrix Schrödinger equation solution increases rapidly with the number of stated treated exactly. The restriction is certainly sufficient for models close to the pure-wino case, when the degenerate states areχ 0 1χ 0 1 , χ + 1χ − 1 in the neutral sector, andχ 0 1χ + 1 ,χ + 1χ + 1 in the charge-1 and and charge-2 sectors, respectively. When the LSP acquires a substantial Higgsino or bino component, the number of degenerate states increases and may exceed four in the neutral and charge-1 sector. An example of a strongly mixed wino-Higgsino LSP model has been analysed in ref. [15], which

JHEP03(2016)119
demonstrated that in this case the effect of the additional states is accurately reproduced by the approximate treatment in the last loop before the annihilation. In the analysis of strongly mixed wino-Higgsino LSP models with a nearly decoupled bino discussed below, all possible 10 neutral states fall below the mass cut 1.2 × 2m LSP in much of the interesting region. We checked on a subset of 1575 analysed model points that the relic density is always accurately reproduced by the approximate treatment. The largest difference we find is 4%, but it is below 1% in 96% of these points, and most of the times closer to the permille level. In any case, this is not a restriction, since the code can always be run with the full set of states treated exactly, at the expense of an increase in CPU time of about a factor of ten.
The Sommerfeld factors are computed from the asymptotic behaviour at r → ∞ of radial solutions of the Schrödinger equation with boundary conditions near the origin. In practice, evolution of the differential equation system to large r is costly, and a finite value of r ∞ must be chosen. We determine this value by requiring that the Sommerfeld factor changes by less than 0.3%, when r ∞ is doubled. This accuracy is often difficult to achieve for very small velocities v, defined by E = m LSP v 2 = √ s − 2m LSP or near values, where new two-particle channels with mass above 2m LSP open, especially forχ +χ− states which experience the long-range Coulomb interaction. Hence we fix (50), when v < 0.03 (within 0.0002 of a threshold). This can lead to local inaccuracies of several percent. However, we find that the deviation from the exact result is oscillatory, and mostly averages out in the thermal average. Once again, this treatment is not necessitated by a limitation of the code but a convenience, since one can set always x ∞ to larger values if needed.
We generate tables of annihilation cross sections (σv rel ) a of two-particle states with on average around 50 velocity points chosen adaptively from 10 −4 . . . 1 with more sampling points near thresholds and the characteristic velocities near the freeze-out temperature. We interpolate these functions and compute the thermally averaged effective cross section σ eff v , summed over all co-annihilating two-particle states for around 60 suitably chosen values of x = m LSP /T between 1 and 10 8 . This table is interpolated and the interpolating function is employed in the Boltzmann equation for Y = n/s. Here G is the gravitational constant, and the parameter g 1/2 * is defined in the standard way as in terms of the effective degrees of freedom g eff and h eff of the energy and entropy densities: For g 1/2 * (T ) and h eff (T ) we use the values derived in ref. [56], which can be found conveniently tabulated as a function of temperature among the package files of the automated JHEP03(2016)119 programs DarkSUSY [57] and micrOMEGAs [25,26]. 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 ref. [58].
Given the annihilation matrices, the calculation of all Sommerfeld factors, cross section tables, thermal averages and, finally, the evolution of the Boltzmann equation through freeze-out takes about 400 sec of CPU time, leading to a total computation time (including the evaulation of the annihilation matrices) of somewhat above 10 min per MSSM parameter point.

Analysis
The departure from the pure-wino limit can be obtained by lowering the sfermion masses and/or introducing non-negligible Higgsino or bino fractions of the lightest neutralino. Therefore we organise the analysis and results in three parts: effect of the sfermion masses (section 3.1), Higgsino admixture (section 3.2) and bino admixture (section 3.4). The residual dependence on remaining parameters is discussed in section 3.5.

Impact of sfermions
The role played by the sfermions in the production of the thermal neutralino relic density is threefold: i) they appear in the t-and u-channel annihilation into SM fermions, ii) they introduce corrections to the neutralino and chargino masses indirectly, via loop effects, and iii) if light enough, they can contribute to the effective annihilation cross section through additional co-annihilation channels. The last of these is beyond the scope of this work and we leave a detailed analysis of general sfermion co-annihilation regions with the inclusion of the Sommerfeld enhancement for future work. Therefore in our results for the perturbative and Sommerfeld corrected relic density, shown in figure 1 and the ratio of these results shown in figure 2, we require that all the sfermions are at least 25% heavier than the LSP. 8 From points i) and ii), the indirect effect of changing the spectrum is sub-dominant, even for the regions of parameter space where the Sommerfeld effect exhibits a resonance and where the resulting cross section is extremely sensitive to the mass difference betweeñ χ ± 1 andχ 0 1 . The reason is that the main contribution to this quantity comes from loops involving gauge bosons and the ones with sfermions are suppressed by their large masses.
The direct impact is on the other hand quite important. Decoupled sfermions mean that the only contributions to the effective co-annihilation cross section from the processes with SM fermion final states arise due to the s-channel annihilation through gauge or Higgs bosons. When the sfermions become lighter the t-and u-channel processes start to be non-negligible. This is especially relevant for the co-annihilation channels, while the direct LSP annihilation to SM fermions is helicity or p-wave suppressed. These t-and u-channel diagrams involving sfermions interfere destructively with the s-channel gauge boson exchange effectively lowering the co-annihilation cross section [59].

JHEP03(2016)119
At the TeV scale the degeneracy between the charginos and neutralinos is more pronounced resulting in co-annihilation channels not being Boltzmann suppressed. Therefore the effective annihilation cross section is strongly affected by the co-annihilation channels e.g. from the processesχ + 1χ − 1 →f f ,χ 0 1χ ± 1 →f f . Due to interference between the sfermion t-channel and W boson s-channel diagrams, the lower the sfermion masses, the smaller the contribution from these processes, leading to lower total annihilation cross section and higher thermal relic density. In other words, the contours of constant relic density move towards lower m LSP values as the sfermion masses decrease. This is indeed what is observed in the left panel of figure 1, where the contours of constant perturbative relic density are plotted in the M 2 -M sf plane for the case of a wino-like LSP (µ = 2M 2 , M 1 = 3M 2 ). In particular note that, by varying the sfermion masses one can obtain the perturbative thermal relic density in agreement with the observed abundance over a large range (∼ 800 GeV) of LSP masses. It is also worth pointing out that the fact that the contours become denser as M 2 increases is a simple result of the approximate quadratic dependence of the relic density on the wino mass.
The situation becomes more involved at the non-perturbative level, as shown in the right panel of figure 1. The main features that were previously discussed for the perturbative case are still present, but two important modifications arise. First, the contours are seen to be shifted towards larger values of M 2 . This is simply the effect of the Sommerfeld enhancement on the annihilation cross section, such that one requires a wino mass of around 2.9 TeV rather than 2.2 TeV in order to obtain the correct thermal relic density in the decoupled sfermion case. The size of the shift however depends on the masses of the sfermions, in particular the lowest wino mass giving the correct relic density Ωh 2 | exp = 0.1188 ± 0.0010 [60] without sfermion co-annihilations is around 2.3 TeV. This is related to the second effect, namely the resonance in the Sommerfeld enhancement, which is also responsible for lowering the constant relic density contours in the sfermion mass at m LSP of around 2.3-2.4 TeV. The presence of the resonance is most clearly seen in figure 2 where the impact of the Sommerfeld effect on the relic density is shown. It can be seen that, as expected, the Sommerfeld effect gets stronger for larger values of M 2 until the resonance region is reached, and that in the resonance region the relic density can be suppressed by nearly an order of magnitude. What is worth stressing is that the Sommerfeld effect is also approximately independent of the value of the sfermion masses. This can be easily understood by noting that the largest impact of the Sommerfeld effect comes from its contribution on theχ 0 1χ 0 1 annihilation, which does not depend in any significant way on the nature of the sfermions. We also note that the Sommerfeld effect changes the relic density by almost a factor of two in the region where the observed relic density is attained (green bands in the figures), and by an even larger factor for smaller M sf , when the observed relic density is produced near the Sommerfeld resonance.
The generic behaviour of the results for the relic density as a function of sfermion masses shown and discussed above holds when one departs from wino-like neutralino as well, but with the details depending on the precise neutralino composition and the spectrum of the sfermions. The latter comes from the fact, that while the coupling of the sfermions with the wino is purely gauge, the one with the Higgsino is Yukawa-type, and therefore

JHEP03(2016)119
discriminates the three generations, as well as squarks from sleptons. The analysis of such scenarios, with Higgsino-and bino-like neutralinos, will be provided in the future.

Higgsino admixture
The Higgsino-wino mixing predominantly depends inversely on the difference between µ and M 2 , as discussed in section 2.1. Increasing the Higgsino component of the predominantly wino-like LSP has several effects on the relic density: i) it modifies the LSP annihilation cross section due to different couplings of the wino and Higgsino components, ii) it changes the relevant number and weights of the co-annihilation channels and finally iii) it significantly alters the Sommerfeld effect. The first two effects are very well known, we therefore concentrate on the non-perturbative effects. We choose to parametrise the Higgsino admixture via the difference between the input parameters µ and M 2 . For definiteness, we restrict ourselves to positive µ in the following analysis. In the m LSP range considered, values of µ − M 2 500 GeV lead to nearly decoupled Higgsinos, and the LSP is practically purely wino-like, while values of around 300-500 GeV correspond to a Higgsino fraction of around a few %, growing up to 50% for µ = M 2 .
The results of the analysis are displayed on figures 3 to 5. In figure 3 the contours of constant relic density are shown in the M 2 vs. (µ − M 2 ) plane for the perturbative (left plot, blue lines) and Sommerfeld enhanced (right plot, red lines) cases. In the upper region of the plot the contour lines flatten as we recover the pure wino scenario, while in the lower region they tend to lower values of m LSP because the large Higgsino fraction suppresses the annihilation cross section. Equivalently, for a fixed LSP mass, increasing the Higgsino admixture increases the relic density. At the perturbative level this is mainly a consequence of the lower value of the coupling to gauge bosons, while in the case of the Sommerfeld effect it also is a result of the larger mass splitting between the LSP and the lightest chargino. In figure 4 we show the ratio of the above plots in order to display the impact of the Sommerfeld enhancement over the M 2 vs. (µ − M 2 ) plane.
In figure 5 we show those contours giving the correct thermal relic density, for three different values of the sfermion mass parameter M sf , and show both the perturbative (left three lines, blue) and Sommerfeld-corrected (right three lines, red) results in one plot. The M sf = 12 TeV lines correspond to the (green) bands in the previous figures 3 and 4. We change here from displaying contours of constant relic density to the correct relic density in order to highlight the effect of the sfermion mass parameter. Note that, in agreement with what was discussed in the previous section, the lower the sfermion masses, the larger the relic density and hence the lower the LSP mass at correct relic density -both for the perturbative and non-perturbative results.
Comparing the perturbative result to the full one in figure 5, one observes the following: i) the contours shift to higher masses, indicating the decrease of the relic density with respect to the perturbative result; in particular the mass of the lightest neutralino (when µ > M 2 ) giving correct thermal relic density is around 1.    resonance leads to much more efficient annihilation, strongly suppressing the relic density; it appears at different positions in m LSP depending on the neutralino composition. In particular, on increasing the Higgsino fraction, the resonance occurs for heavier LSPs, which is mainly due to the increasing of the mass splitting between the lightest chargino and neutralino, the decreasing coupling, and the fact that the resonance depends on the splitting through (mass splitting)/(m LSP α 2 2 ).
The position of the peak of the resonance was clearly visible in figure 4. It is therefore evident that in figure 5 the contours of the correct relic density cluster around this peak for higher values of M 2 and lower values of µ − M 2 . This is easily understood when one recalls that in this region the neutralino at a perturbative level has a thermal abundance larger than that observed by a factor of a few. The proximity to the resonance enhances the cross section, reducing the relic density to agree with the measured value. In particular, it follows that the largest value of M 2 giving the correct thermal relic density is close to 3.3 TeV, approximately 20% higher than that for the pure-wino scenario.
Note also that, in contrast to the pure-wino scenario with decoupled sfermions, a region of parameter space exists where the thermal relic density is obtained in very close vicinity to the resonance, leading to strong bounds on such scenarios coming from dark matter indirect searches. Previously only limiting wino and Higgsino cases have been studied from this perspective with the inclusion of the Sommerfeld effect [8,61,62], and even slightly mixed scenarios remain unexplored. 9 It also follows that some regions of the pMSSM JHEP03(2016)119 parameter space can be effectively constrained by non-observation of any dark matter signal in cosmic or γ-rays. The precise analysis of such phenomenologically interesting regions will be presented in upcoming work [66].

Effect of the heavy Higgs bosons
In the MSSM, the only particles beyond the SM having positive R-parity are the additional Higgs bosons. These can therefore act as an s-channel mediator and, if light enough, as end-products of the (co-)annihilation. As the effect of these Higgs bosons is greatest when the Higgsino mass parameter is close to M 2 , we discuss this first before moving on to the wino-bino mixed case. The additional Higgs bosons can affect the relic density in the following two ways: • by contributing to the (co-)annihilation rate via s-channel diagrams, particularly if M A lies in the vicinity of 2 m LSP , thereby typically reducing the relic density • by providing additional final states with one heavy Higgs plus one gauge or light Higgs boson, or with two heavy Higgs bosons, if the combined mass of the final state lies below 2m LSP , which leads to a reduction in the relic density.
The former is only relevant when the coupling of the annihilating particles to the Higgs bosons is non-negligible. This requires one of the two annihilating neutralinos or charginos to contain a considerable gaugino component and the other a considerable Higgsino component. Forχ 0 annihilation this implies that the LSP is mixed. We remind the reader that we do not consider the resonant annihilation region when M A is inside the interval 1.7-2.3 m LSP as explained in section 2.2. As for the latter point, the heavy Higgs and gauge boson final state is obtained via a s-channel gauge boson, or a t-channel neutralino or chargino. This is again more relevant when the LSP contains some Higgsino admixture, as this also allows the coupling of neutralinos to Z bosons. However, in contrast to the case of the heavy Higgs boson in the s-channel, this contribution does not vanish when the Higgsino decouples, as a coannihilating chargino and neutralino can annihilate into a heavy Higgs and gauge boson via a s-channel W boson even in the pure-wino limit. We explore these issues in figure 6 where we show contours of constant relic density in the M 2 vs. M A plane both at the perturbative level (left) and on taking account of the Sommerfeld effect (right). The region corresponding to the measured relic density is shown by the green band. As the s-channel resonance cannot be accurately calculated in our framework we do not provide results near M A = 2 M 2 . It is seen that for the perturbative case above the excluded region the lines are approximately vertical, just bending slightly towards higher values of M 2 on approaching this region. Below the excluded area we find that there is a slight shift to the right as the heavy Higgs bosons are accessible in the final state. The difference in M 2 giving the correct relic density is approximately 150 GeV when M A changes from 10 TeV to 500 GeV. For the Sommerfeld-enhanced case the result is qualitatively similar, however the difference in M 2 giving the correct relic density is around 250 GeV for the same change of M A . approximate way and/or without inclusion of recent developments [12][13][14].

JHEP03(2016)119
In figure 7 we further investigate the effect of the heavy Higgs bosons on the contours showing the correct relic density in the M 2 vs. µ − M 2 plane. The blue lines show the perturbative result while the red lines include the Sommerfeld enhancement. We see that on decreasing M A from 10 TeV to 800 GeV, the shift in the value of M 2 giving the correct density is indeed dependent on the proximity of µ to M 2 , increasing from 50 to 150 GeV in the perturbative case and 100 to 250 GeV in the Sommerfeld-enhanced case. As mentioned earlier, an increased Higgsino admixture allows a stronger coupling to the Higgs and Z bosons in the s-channel (where Z bosons can give rise to heavy Higgs bosons in the final state), increasing the effect of the heavy Higgs boson. Nevertheless when the Higgsino is decoupled a dependence on M A persists; for the M A = 800 GeV contours coannihilation via a W boson to a final state containing a heavy Higgs boson and a gauge boson is allowed but not for the M A = 10 TeV contours.

Bino admixture
The bino only mixes with the wino via the off-diagonal terms in the Higgsino block of the neutralino mass matrix. It follows that the mixing is weak, depending of course on the Higgsino parameter µ, and is further sensitive to tan β and the sign of M 1 and µ as seen in eq. (2.5). In order that the wino-like neutralino contains a substantial bino component, either µ should be of the same order as M 1 and M 2 or the M 1 and M 2 parameters should be highly degenerate. For example, when δM 1 = M 2 −M 1 = 10 GeV, µ = 2 M 2 and tan β = 15 the mixing is about 1%, decreasing to 0.1% when δM 1 = 100 GeV. Such situations may arise and are worth studying as the resulting features are of phenomenological interest. In this section we focus on the second case M 1 ∼ M 2 µ, since the first (M 1 ∼ M 2 ∼ µ) falls into the category of a mixed wino-Higgsino state and shares the gross features with the case of a decoupled bino analysed in the previous section. We further assume M 1 > 0, since M 1 < 0 entails an essentially decoupled bino.
When M 1 is close to M 2 , the perturbative relic density is affected both by the modification of the LSP annihilation cross section due to the change in composition and by the co-annihilation with the bino-like NLSP, with mass close to M 1 . On top of that the Sommerfeld effect is modified, analogously to the Higgsino-wino mixed scenario, by the weakened coupling and larger mass splitting betweenχ 0 1 andχ ± 1 as given in eqs. (2.6), (2.7). Qualitatively the behaviour observed on increasing the bino component is largely the same as in the Higgsino case, with an important quantitative difference: a larger sensitivity of the results to the remaining parameters.
Contours of constant relic density in the M 2 vs. (M 1 − M 2 ) plane for M 1 > 0, and both the perturbative and Sommerfeld enhanced case are displayed in figure 8. Note that the logarithmic scale for the vertical axis, chosen due to the weak mixing between the bino and the wino, changes the appearance of the resonance with respect to the Higgsino case. As one increases the bino component the mass of the LSP resulting in the correct relic density is approximately 1500 GeV rather than 1800 GeV for the perturbative case. This changes rather dramatically when the Sommerfeld enhancement is taken into account, notably for strong mixing, i.e.  relic density. One can interpret the larger two of these values as a result of the resonance in the Sommerfeld enhancement.
In figure 9 we study the ratio between the relic densities shown in figure 8, in terms of a density plot with contour lines overlaid. The correct relic density for the full calculation including the Sommerfeld enhancement is highlighted by the (green) band. We observe the maximal effect of the Sommerfeld enhancement is in fact in the region where the relic density agrees with observation, in particular when the difference between M 1 and M 2 is below approximately 10 GeV, the Sommerfeld enhanced relic density in agreement with that observed is three times smaller than the perturbative result at the same parameter values. Note that over the entire region covered by the plot the effect of the Sommerfeld enhancement is greater than 30%.
In figure 10 contours with the correct relic abundance for three choices of the sfermion mass parameter are shown. The results again resemble the Higgsino admixture case, up to differences already commented on. Note that in the region around the resonance, the effect of the sfermion masses is less pronounced than elsewhere. One observes that the lowest mixed wino-bino neutralino mass giving the observed relic density is around 1.8 TeV for M sf = 4 TeV, marginally higher than the wino-Higgsino case. The highest value is 2.9 TeV (for M sf = 12 TeV) compared to 3.3 TeV in figure 5. However, as can be seen in figure 11, the highest value of M 2 resulting in the correct relic density is strongly dependent on the value of the µ parameter, as this mediates the mixing. This dependence is demonstrated via contours for five different choices of µ. The contour for µ = 1.1M 1 bears a closer resemblance to the wino-Higgsino case, as suggested earlier. Note that as µ decreases, JHEP03(2016)119 the lightest chargino-neutralino mass splitting for a given point in the plane increases, 10 resulting in the resonance moving to higher values of m LSP . Due to the presence of the resonance, it appears that by making an appropriate choice in µ and M 1 the entire region could be covered, at least for values of M 2 from 2100 to 4200 GeV if not even higher. All these points would be on or around the resonance, having implications for Indirect Detection. Moreover, interestingly adding a bino component to the LSP can extend the possible neutralino masses giving observed dark matter abundance up to and even beyond 4.1 TeV.

Residual dependence on other parameters
In the previous sections 3.1 to 3.4 certain parameters were fixed in order to obtain a clearer understanding of the dependence of the results on the central parameters M 1 , M 2 , µ, M sf , as well as on M A . However, it is important to confirm whether these are indeed the most relevant, and to investigate the effect of the other parameters, e.g. tan β, which was so far neglected.
One case in which additional parameters may play a significant role is when the lightest neutralino is a wino-bino mixture as in section 3.4. This is because, as seen in eqs. (2.5) to (2.7), the mixing of the wino with the bino, and the splitting between the lightest neutralino and chargino is sensitive to |µ|, the sign of µ and tan β. The dependence on µ was already examined in section 3.4, and the results can be found in figure 11. Here we further consider the effect of the sign of µ and the choice of tan β. Our results can be found in figure 12. The benchmark choice for the results presented in previous subsections was tan β = 15 and µ > 0; in addition here we consider µ > 0 with tan β = 5, 30 and µ < 0 with tan β = 15. Large deviations from the benchmark scenario are seen in the resonant region.  This can be understood by examining the expressions for the mass splitting between the lightest neutralino and chargino in eqs. (2.6) and (2.7). The splitting is seen to increase when tan β decreases, and also when µ > 0 compared to µ < 0, resulting in the position of the resonance moving towards higher values of m LSP , i.e. the correct relic density is observed for higher M 2 .
We further examine the sensitivity to the remaining parameters for both the cases of The perturbative results are shown by the dark-grey/blue histograms. We see that in all cases the distribution is strongly peaked near the central value with the variation of order of at most a few per cent, and that the distribution widens with the departure from the pure-wino limit. The situation changes when considering the full result (light-grey/red histograms), as all the distributions become broader and asymmetric. The position of the resonance in the Sommerfeld effect is greatly sensitive to values of the neutralino mass,

JHEP03(2016)119
couplings and the chargino. Therefore, slight changes in these values caused by different choice of remaining less relevant MSSM parameters, especially close to the resonance, can lead to observable differences in the relic density. Indeed, the broadening of the full result with respect to the perturbative one is strongest in the middle panel (due to the vicinity of the circle benchmark point to the resonance) and in the upper left-hand plot being not far from the resonance as well. The asymmetry in the distributions is caused by the fact that deviations around central parameters may go towards or away from the resonance, leading to larger or smaller Sommerfeld effects respectively. The bottom line of this analysis is that away from the resonance the residual MSSM parameters have a very mild impact, justifying our choice of central parameters, while in the vicinity of the resonance regions the variation is very significant.
To study the dependence on the residual parameters even further, we have generated a large number of points (50000 and 90000 for the Higgsino and bino case, respectively), where we considered the wino mass in the range M 2 ∈ {1, 3.5} TeV and (different from table 1 and the analyses in the previous sections) fixed the gluino mass parameter via M 3 = 2M 2 . The sfermion masses were fixed to the values given below, but we varied all other parameters in the following ranges: where f in A f includes all fermions except the top. In addition, for the Higgsino case we chose: M sf = 6 TeV, |M 1 | = 2.01M 2 , M 2 ∈ {1, 3.5} TeV, µ ∈ {M 2 , M 2 + 0.5 TeV}, (3.2) and for the bino case: From the generated points, we selected those where either the perturbative or the Sommerfeld enhanced relic density was found to be between 0.1168 and 0.1208, i.e. within two sigma of the central value. In figure 14 we overlay these points on the relevant plots shown earlier, figures 5 and 10 for the Higgsino (upper plot) and bino (lower plot) case, respectively. We observe that those points for which the perturbative relic density lies within 2σ of the central value are located very close to the respective sfermion mass contours, the spread of the points being comparable to twice the width of the 1σ contours. We conclude that the dependence on the residual parameters is very mild in this case. As could be expected, when including the Sommerfeld enhancement, the residual parameters can have a larger effect, especially close to the resonance. In order to investigate this effect further, for the wino-Higgsino case we have divided these points according to whether 2.3m LSP < M A < 2.5m LSP , M A > 2.5m LSP and M 1 > 0, or M A > 2.5m LSP and M 1 < 0. These are indicated in figure 14 by the cross, filled circle and open circle respectively. The former division is made in order to isolate those points in proximity to the heavy Higgs funnel region and the latter due to the effect of the sign on M 1 on the lightest neutralino-chargino mass splitting.

JHEP03(2016)119
We do not plot the points with M A < 1.7m LSP in this case. For the wino-bino case we separated the points according to whether µ > 0 or µ < 0, indicated in figure 14 by the filled circle and open circle respectively, as the sign of µ also plays a role in the size of the lightest neutralino-chargino mass splitting.
The effect of the residual parameters is sub-dominant with respect to e.g. that of the sfermion masses, but both the sign of µ and M 1 are seen to play a role in the resonance region for the wino-bino and wino-Higgsino cases, respectively. This can be understood in terms of the expressions for the mass difference between the lightest chargino and neutralino in eqs. (2.3), (2.4) and (2.6), (2.7) to which the resonance is sensitive. As the splitting increases the position of the resonance moves towards higher values of m LSP . Whether the heavy Higgs is below, above, or, in particular, close to the excluded window also has a noticeable effect for the case of wino-Higgsino mixing, and this extends beyond the resonance region and holds for the perturbative case as well. This is because for states with larger mixing the coupling to the heavy Higgs is enhanced, and therefore when M A decreases the s-channel annihilation cross section increases, and one has to go to higher values of M 2 to obtain the correct relic density. This is not relevant for the wino-bino case, where the dependence on the value of M A is negligible.
To summarise, we find that the assumption that our results of the previous sections were more or less independent of certain parameters was largely justified. Only for the wino-Higgsino case there is some dependence on the value of M A , and in the resonance region the sensitivity to these parameters is somewhat enhanced, particularly to the values of tan β and the sign of µ for the case of wino-bino mixing and M 1 for the case of wino-Higgsino mixing.

Summary
We have studied the Sommerfeld effect on the relic density of neutralino dark matter beyond the pure-wino limit. This involved a scan of parameter space for three scenarios where the lightest neutralino contained a large wino component: one with non-decoupled sfermions, and the remaining with either a Higgsino or bino admixture. We aimed to determine how in these scenarios (a) the mass of the LSP where the relic density constraint is satisfied and (b) the size of the Sommerfeld enhancement is altered in comparison to the pure-wino case.
The calculation of the Sommerfeld enhancement for the scenarios in question required a consistent treatment of mixed neutralinos including multiple co-annihilation channels and off-diagonal contributions as well as of O(v 2 ) contributions. As the Sommerfeld effect, in particular the position of the resonance, depends strongly on the mass splittings between the neutralinos and charginos, we used a dedicated on-shell renormalisation scheme scheme for one-loop masses. As the relic density depends strongly on the precise value of the gauge coupling, we adopted running couplings. Further we have argued that the size of thermal corrections is sufficiently below the uncertainty of our calculation, i.e. the percent level, that they can be neglected. Finally, all points in MSSM parameter space considered were checked for consistency with current experimental measurements. Our calculation was JHEP03(2016)119 carried out by a code which will be made available to the public: this will be presented in more detail in a separate publication.
Due to t-channel interference, we found that when the sfermions are non-decoupled they reduce the annihilation cross section such that the relic density constraint is satisfied at lower values of m LSP , from 2.9 TeV for decoupled sfermions down to 2.4 TeV. For the mixed neutralino scenarios we found a much larger dependence than expected on those parameters affecting the mass splitting between the lightest chargino and neutralino. This was particularly evident in the mixed bino-wino region, where the position of the resonance was seen to be sensitive to µ, tan β and the sign of µ. As the splitting increases we observed that the resonance lies at higher values of m LSP . This led to a large range of neutralino masses from 1.8 to beyond 4 TeV satisfying the relic density constraint. For the Higgsinowino mixed region the position of the resonance depends primarily on µ − M 2 , but whether the mass of the heavy Higgs boson lies below or above 2 m LSP also plays a significant role. Here we found values of m LSP ranging from 1.7 to 3.3 TeV.
This is the first time that the sensitivity of the Sommerfeld enhancement to MSSM parameters for mixed neutralinos has been studied systematically and to such accuracy, and the large range of possible m LSP masses providing the correct relic density was previously unknown. In most of the cases the Sommerfeld effect changes the relic density by a factor of two or even higher relative to the tree-level computation. This underscores the fact that the relic density of TeV scale MSSM dark matter can usually not be predicted correctly without accounting for this effect.
In light of these results a re-investigation of the bounds on the Sommerfeld enhanced scenarios coming from Indirect Detection experiments is imperative. It is likely that so far unexplored regions exist, sufficiently far away from the resonance that they are not excluded, but with Sommerfeld enhanced annihilation cross sections which could be probed by upcoming experiments, e.g. the Cherenkov Telescope Array (CTA). This will be the subject of a dedicated study in the near future.

JHEP03(2016)119
Furthermore, large thermal masses and mass splittings may be generated. While thermal effects on the short-distance annihilation process are small [67], the gauge-boson mass determines the range of the potential, which is an important quantity for the Sommerfeld effect. Furthermore, the mass splitting of the lightest chargino and neutralino influences the location of the Sommerfeld resonance. In the following, we investigate the thermal modification of the gauge boson mass and neutralino-chargino mass splitting through a combination of estimates, analytical calculations and numerical checks.
The relevant temperature range for this investigation is limited from above by T m χ /20, when freeze-out begins, which allows us to treat T /m χ as small. The temperature of the Universe together with the Boltzmann distribution sets the characteristic scale of the three-momentum of the scattering dark matter particles to | p | ∼ (m χ T ) 1/2 . The Sommerfeld effect is caused by ladder diagrams with loop momentum k satisfying k 0 | k| m χ where k 0 is determined by the pole of the dark matter particle propagator. The characteristic scale of k is also (m χ T ) 1/2 until m χ T ∼ m 2 W , where m W is the mass of the exchanged electroweak gauge boson (zero for the photon), at which point the Sommerfeld enhancement saturates. When the Universe cools below T s m 2 W /m χ , the external momentum p continues to decrease while k ∼ m W remains constant, and the thermal modification of the Sommerfeld effect fades out. Hence the temperatures of interest are limited from below by T s in the range from 1 to 4 GeV. An exception is the Sommerfeld enhancement due to photon exchange between the charginos, whose effect on the relic density turns off only when the charginos decouple from the thermal plasma at a temperature set by the neutralino-chargino mass difference.

A.1 Higgs vacuum expectation value
We approximate the temperature dependence of the Higgs field vacuum expectation value by and zero above T c . The critical temperature is taken to be T c = 165 GeV, as follows from the effective potential given in ref. [68]. The expansion of the Universe proceeds adiabatically such that the particle masses are given by the standard expressions with the instantaneous value v(T ).

A.1.2 Lightest neutralino-chargino mass difference
The temperature dependence of the neutralino and chargino masses is not by itself of interest, since always m χ T . However, the temperature dependence of small mass splittings must be considered, since the mass splitting determines, for instance, the location of the Sommerfeld resonance, and further appears in the nearly on-shell propagator of the two-neutralino/chargino state in the ladder diagrams.
Close to the pure-wino limit the neutralino-chargino mass difference is dominated by the radiatively induced splitting, which in the pure-wino limit is given by The expression refers to the approximation M Z m χ and the numerical value employs the SU(2) coupling α 2 (m χ ) = 0.032810 at the scale m χ = 2.5 TeV. Whenever the radiative mass splitting dominates over the tree-level splitting, it changes very little compared to the pure-wino value. We therefore assume that it is proportional to v, and implement the temperature-dependence by multiplying the zero-temperature radiative contribution to the mass difference with v(T )/v.
In the presence of a Higgsino-or bino-component of the wino-like neutralino, there is an additional tree-level mass splitting, which can be computed from the mass matrices X, Y in section 2.1. The dependence on the Higgs vacuum expectation value changes from quartic towards the pure-wino limit to quadratic when the Higgsino or bino admixture increases. The cross-over to quadratic dependence occurs in about the same region in the parameter space shown in figures 3 (Higgsino admixture) and 8 (bino admixture),

JHEP03(2016)119
This is the saturation regime for the Sommerfeld enhancement of the Yukawa potential. (3) The neutralino-chargino mass difference provides an O(1) modification of the Sommerfeld enhancement factor I whenever m χ δm +0 ∼ max (| p| 2 , m 2 W ) and reduces or cuts off the enhancement when m χ δm +0 is larger than the right-hand side of this relation.
In the following we use the above expression to estimate the impact of the thermal modifications of m W (T ) and δm +0 (T ) discussed in the previous subsection. In doing so, we correlate the external neutralino or chargino momentum with the temperature of the Universe according to p 2 ∼ m χ T = m 2 χ /x. We further support these estimates by implementing the thermal effects into our Sommerfeld code as described below.

A.3.1 Yukawa potential
We first study the modification of the Yukawa potential generated by W exchange. As discussed above, the leading effect is the temperature dependence of the Higgs vacuum expectation value. The thermal self-energy correction has the same temperature dependence, but is smaller. The modification of I is due to the W propagator 1 max (m χ T, m 2 W ) + m 2 W − T 2 /4 , (A.11) where we assume T < T c and approximate m 2 W /T 2 c → 1/4. At the beginning of freeze-out when m χ T > m 2 W , the relative size of the thermal correction is 1/(4x) which for x f ∼ 20 does not exceed 1.2%. When saturation is reached at T s ∼ m 2 W /m χ , the relative correction is only of order m 2 W /(4m 2 χ ) ∼ 0.02% for a reference dark matter mass m χ = 2.5 TeV. An analytic expression for the Sommerfeld effect is available in a one-state model, when the Yukawa potential is replaced by the so-called Hulthén potential, which provides a good approximation [73]. Using this expression we find a maximal change of the Sommerfeld factor of 0.2% at the beginning of freeze-out and decreasing afterwards, confirming the above simple estimate.
The modification of the relic density is expected to be even smaller, since the suppression due to the Sommerfeld effect builds up from the beginning of freeze-out, where it is least significant, until about x ∼ 10 4 , where annihilations terminate, see for instance figure 4 of ref. [15]. We have implemented the thermal modification of the potential in our code to check this explicitly. The modification is CPU expensive, since the Sommerfeldcorrected cross section, which is thermally averaged for given T , must now be recomputed for every value of T . We created two-dimensional cross section tables in velocity and temperature, adopting 59 temperature points. We then compute the Sommerfeld effect and relic density including the temperature-dependent potentials for the Higgsino-to-wino trajectory in MSSM parameter space considered in ref. [15], to which we refer for details on these models. For the mostly wino models 8 to 13 of the trajectory, we find that the relic density change is below one permille in all cases, in good agreement with the above estimates, and not visible within the numerical accuracy of the code.

A.3.2 Neutralino-chargino mass splitting
It is evident from the temperature dependence of the two contributions (vacuum expectation value and thermal self-energy) that the largest relative effect of the thermal correction

JHEP03(2016)119
to the neutralino-chargino mass splitting again arises at the beginning of freeze-out. This has two immediate consequences. (1) While the Sommerfeld resonance depends sensitively on the mass splitting (see, for instance the two-state model in ref. [74], which shares the essential features regarding the mass splitting with wino-like MSSM models), the resonance effect develops sufficiently late after the beginning of freeze-out. We therefore conclude that the thermal effect on the resonance region is negligible. (2) At the beginning of freeze-out λ ∼ 1 in (A.10), hence the relative modification of I due to the mass splitting is of order as the thermal correction to the mass splitting is as large as the mass splitting itself. The numerical estimate is based on the assumption that the entire zero-temperature mass difference of up to 0.5 GeV vanishes due to the vanishing of the Higgs expectation value at T = T c , which gives the largest possible effect. We studied the impact of the temperature-dependent neutralino-chargino mass splitting on the relic density with the extended numerical code described above for the wino-like trajectory models of ref. [15]. The thermal modification is again in the permille range, in agreement with the analytic estimates, reaching 0.7% at maximum. Once again we find that the observed thermal effect is of the same order at the numerical uncertainties due to sampling and the choice of x ∞ , hence we can only state that the thermal effect is well below 1% in all cases studied.

A.3.3 Summary
We conclude that thermal modifications of the Sommerfeld effect change the relic density at most in the upper permille range, which is negligible for all practical purposes. We point out that our investigation of thermal effects is not complete. For example, we did not discuss the direct modification of the neutralino and chargino two-particle wave function due to interactions with gauge bosons in the thermal plasma, an effect that would be referred to as "dissociation" by soft gauge bosons in the case of bound states. Power counting suggests that this effect is of the same order as the ones investigated here. However, since all these are far smaller than the theoretical uncertainty from perturbative higher-order corrections, which is probably a few percent, we do not attempt a complete analysis in this work.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.