{'fraction_non_alphanumeric': 0.029848246289906535, 'fraction_numerical': 0.014371377843288333, 'mean_word_length': 4.345898997134671, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accurate relic densities with neutralino, chargino and sfermion coannihilations in mSUGRA

Neutralinos arise as natural dark matter candidates in many supersymmetric extensions of the Standard Model. We present a novel calculation of the neutralino relic abundance in which we include all so called coannihilation processes between neutralinos, charginos and sfermions, and, at the same time, we apply the state of the art technique to trace the freeze-out of a species in the early Universe. As a first application, we discuss here results valid in the mSUGRA framework; we enlight general trends as well as perform a detailed study of the neutralino relic densities in the mSUGRA parameter space. The emerging picture is fair agreement with previous analyses in the same framework, however we have the power to discuss it in many more details than previously done. E.g., we find that the cosmological bound on the neutralino mass is pushed up to ~565 GeV in the stau coannihilation region and to ~1500 GeV in the chargino coannihilation region.


Introduction
The latest years will be remembered in the history of Science as those that marked the entrance into the era of precision Cosmology. A number of experiments have been pinning down the values of cosmological parameters to a level of precision hardly foreseeable just a decade ago, with perspectives from upcoming measurements even more spectacular. Most remarkably, experiments with different focus, as well as exploiting complementary techniques, have all collected data which point to one single overall-consistent picture, a "concordance" model [1] in which the Universe is flat with about 30% of its present average energy density in matter and about 70% in some form with negative pressure (a cosmological constant or dark energy). More precisely, from the combined analysis of the latest data on the cosmic microwave background and large scale galaxy surveys, the cold dark matter (CDM) and baryonic contributions have been recently estimated [2] to be Ω CDM h 2 = 0.115 ± 0.009 and Ω b h 2 = 0.022 ± 0.002, respectively (here Ω is the ratio between mean density and critical density ρ c = 1.879 × 10 −29 h 2 g/cm 3 , and h is the Hubble constant in units of 100 km s −1 Mpc −1 ; the analysis in [2] gives h = 0.665 ± 0.047).
Among the ideas which have been put forward to account for the CDM term, the most natural solution is probably the scheme in which CDM appears as a thermal leftover from the early Universe: in this context, stable weakly interacting massive particles (WIMPs) are ideal CDM candidates, as their thermal relic abundance is naturally of the order of the measured one. Given the accuracy of the current and future measurements of Ω CDM , it is useful to have an equally accurate calculation of the relic abundance of WIMPs. This paper continues our program of accurate relic density computations.
Since the late seventies, when the idea of WIMP dark matter was formulated [3,4,5], the computation of the WIMP relic abundance has been constantly refined. One important step was to recognize the role of coannihilation effects [6,7]: if in the particle physics theory one considers, a stable WIMP appears together with a slightly heavier particle into which it can transform, when computing the present density of the lightest particle one needs to retrace the thermal history of both particles simultaneously. Such an effect is common even for the most popular WIMP dark matter candidate, the lightest supersymmetric particle (LSP) in supersymmetric extensions of the particle physics Standard Model. When the LSP is the lightest neutralino, the mass splitting between the LSP and the next-to-lightest supersymmetric particle may in some cases be small. Coannihilations are then important. Their inclusion in the computation of the neutralino relic density has been the subject of numerous studies, which vary in the degrees of refinement of the method implemented to compute the relic abundance and in the selection of which particles to include in the coannihilating set.
The novel calculation that we present in this paper is included in a new extended version of the DarkSUSY package [8], which will be publicly available in the near future. DarkSUSY now allows for the most generic coannihilation effect in the framework of the minimal supersymmetric extension of the Standard Model (MSSM), at the same time applying the state of the art technique to trace the freeze-out of a species in the early Universe [9]. In fact, on one side, we fully solve the density evolution equation (which determines the evolution of the number density of neutralinos) numerically, including all possible resonance and threshold effects, and avoiding approximations in the thermally averaged cross sections (such as the expansions in terms of the relative velocity that is often applied). On the other side, extending the work of two of us [10] who first applied such a technique to the case of coannihilations with charginos and neutralinos, we include here the possibility of having coannihilations with all sfermions as well. As a result, assuming masses, widths and couplings of particles in the MSSM are given with an adequate precision, we provide here a tool to compute neutralino relic abundances with an estimated precision of 1% or better.
Compared to other recent calculations, we believe this is the most accurate calculation available at present. The standard lore so far has been to calculate the thermal average of the annihilation cross section by expanding to first power in temperature over mass and implementing an approximate solution to the evolution equation which estimates the freeze out temperature without fully solving the equation (see, e.g., Kolb and Turner [11]). Sometimes this is refined by including resonances and threshold corrections [7]. Among recent studies, this approach is taken in e.g. Refs. [12,13]. Other refinements include, e.g., solving the density evolution equation numerically but still using an approximation to thermal effects in the cross section [14,15,16,17,18], or calculating the thermal average accurately but using an approximate solution to the density equation [19,20,21]. At the same time, only in a few of the quoted papers the full set of initial states has been included. As already mentioned, the present calculation includes all initial states, performs an accurate thermal average and gives a very accurate solution to the evolution equation. Though the inclusion of initial state sfermions in the DarkSUSY package is a new feature introduced in the present work, other groups [22,23,24] have earlier introduced some sfermion coannihilations in an interface with the old DarkSUSY version.
The DarkSUSY package has been written in a very general and flexible format, so that it can be used for any theory embedded in the MSSM. As a first example, we will present here results valid in the minimal supergravity (mSUGRA) framework. The mSUGRA framework is the framework considered in most of the previous analyses. Although rather restrictive in the way parameters are set, it is sufficient to enlight most coannihilation effects which can emerge in the MSSM.
The outline of this paper is the following: in the next section we discuss the supersymmetric model we work in; in section 3 we review the framework in which we calculate the relic density including coannihilations. In section 4 we examine the effects of coannihilation in detail, stressing the physical insights of the results.

Supersymmetric model
The particle physics model implemented in the DarkSUSY package is a minimal supersymmetric extension of the Standard Model, built with N = 1 generator of supersymmetry and containing the smallest possible number of fields. We restrict ourselves to the case in which the LSP is the lightest neutralino, i.e. the lightest of the four mass eigenstates obtained from the superposition of the supersymmetric partners of the neutral gauge and Higgs bosons. The supersymmetric part of the spectrum contains also two chargino mass eigenstates, the gluino and the scalar superpartners of leptons and quarks. The Higgs sector needs two Higgs doublets, which, after electroweak symmetry breaking, give five scalar fields, denoted as H ± , H 0 1 , H 0 2 and H 0 3 (for a more detailed description of the model and of our conventions see [25]). 1 As needed in the relic density calculation, we perform the computation of all two-body final state cross sections at tree level for all initial states involving neutralinos, charginos, sleptons and squarks. We do not neglect mass terms for fermions in final states (as sometimes done in the literature), nor implement any expansion to first order in relative momentum of the incoming particles (as most often done in the past). 2 A list of all processes included is given in Appendix A. 1 Apart from sfermion coannihilations, another major improvement in the DarkSUSY version that we use for this work is the inclusion of one-loop corrections to the Higgs widths, with formulas taken from Refs. [26,27]. 2 In the current version of DarkSUSY we do not include processes with the gluino in the initial state or any flavour changing processes in the sfermion coannihilations; these might be included in future work but we do not expect any result in the present analysis to be affected.
We compute analytic expressions for the amplitudes using standard Feynman rules with generic expressions for the vertex couplings (e.g. ig L φf f ′ P L + ig R φf f ′ P R for a vertex involving two fermions and a scalar). The amplitudes for all the different types of Feynman diagrams are obtained with the help of symbolic manipulation programs. For initial states involving neutralinos and charginos, helicity amplitudes for each type of diagram are computed analytically with Reduce [28] (this is described in [10] 3 ). For the diagrams including sfermions in the initial state, Form [29] is used to analytically calculate the amplitudes squared. The output from Form is then converted into Fortran with a Perl script. All of these analytic formulas can be converted into a compact form, but we do not consider it useful to reproduce those expressions here.
The actual values of the MSSM vertex couplings are introduced during the numerical calculation. We recomputed all vertex couplings from the MSSM lagrangian, and checked them against standard literature (e.g. [30,31]; DarkSUSY also contains couplings which do not appear in the literature, namely those involving generic intergenerational mixing in the sfermion sector).
When giving explicit examples of computations of the neutralino relic abundance we restrict ourselves in this paper to mSUGRA models. Under the assumption of universality at the grand unification scale, the mSUGRA action has five free parameters, m 1/2 , m 0 , sign(µ), A 0 and tan β. The parameters m 1/2 , m 0 and A 0 are the GUT unification values of the soft supersymmetry breaking fermionic mass parameters, scalar mass parameters and trilinear scalar coupling parameters, respectively (of the trilinear couplings, only A t , A b and A τ differ from zero). The absolute value of the Higgs superfield parameter µ follows from electroweak symmetry breaking, but its sign is free. Finally, in the Higgs sector, tan β denotes the ratio, v 2 /v 1 , of the vacuum expectation values of the two neutral components of the SU(2) Higgs doublets. Our convention on the sign of µ is that µ appears with a minus sign in the superpotential (i.e. following the convention used in e.g. [30]), while the definition of sign and dimension for the A 0 's are such that A t appears in the stop mass squared matrix as off-diagonal elements of the form (A t − µ cot β)m t (and analogously for A b and A τ ).
GUT scale values of the soft breaking parameters, as well as of the gauge and Yukawa couplings, have to be evolved down to the weak scale. To do that we make use of the ISASUGRA RGE package in the ISAJET 7.64 software [32] 4 . The interface to DarkSUSY is such that the whole SUSY spectrum of masses and mixings as given in the output at the weak scale in ISASUGRA is given as input in the DarkSUSY MSSM spectrum. Regarding Standard Model masses and couplings, we have implemented 1-loop renormalization group equations. We fix the top pole mass at 174.3 GeV, which is the central value stated by the Particle Data Group 2002 [33]. For the c and b quarks we input the running quark masses 3 Compared to the analysis in [10], we have here corrected a sign error in the chargino -(down-type) fermion -(up-type) sfermion vertex with clashing arrows. Consequently, the discrepancy between Dark-SUSY and micrOMEGAs (model B in [19]) has now disappeared. 4 The ISASUGRA output is rather unstable in regions where the convergence is slow (e.g. in the focus point region where the lightest chargino is nearly degenerate with the lightest neutralino); to improve the stability, we have converted ISASUGRA to double precision and increased the requirements on the convergence. m q (m q ) and choose the values of 1.26 GeV and 4.2 GeV respectively [33,34]. To a good approximation, the RGE scale at which (co)annihilations occur is twice the LSP mass, and therefore we evaluate the gauge and Yukawa couplings at this scale.
We have chosen to use the ISASUGRA RGE code since it is a widely used program that is kind of a standard in the field. However, there are other RGE codes on the market (e.g. SOFTSUSY [35], SPHENO [36] and SUSPECT [37]), all of which solve the RGEs with a different level of sophistication. We do not want to go into the details of these different packages here, but instead refer the interested reader to a recent comparison in Ref. [38]. The main message we want to convey here is that different RGE codes give slightly different results. E.g. the sparticle masses typically differ by a few GeV between the different codes.
When scanning over the parameter space, we compare with accelerator constraints, implementing limits recommended by the Particle Data Group 2002 (PDG) [33]. We have not included limits which are still preliminary as of this writing (most important among the preliminary bounds is the limit of about 104 GeV on the chargino mass; we have nevertheless indicated its effect in some of the figures). For the Higgs's H 0 2 and H 0 3 we have implemented the tan β dependent mass limits as provided by the PDG in the most conservative setup. This is particularly important for the mass of H 0 2 since for a substantial tan β interval this constraint is less restrictive than the lower mass limit (114.3 GeV) on the Standard Model Higgs. In this paper, which focusses on coannihilations, we do not emphasize the role of the b → sγ constraint (we refer the interested reader to recent dedicated analyses in the same mSUGRA context, e.g., Refs. [39]).

Relic density calculation
The importance of coannihilations in the computation of the density of a relic particle was recognized by Binetruy, Girardi and Salati [6], and independently by Griest and Seckel [7]. In Edsjö and Gondolo [10], this was further analysed and put in a form that allows for an accurate treatment of coannihilations in the same basic framework as without coannihilations. We do not repeat every step of that calculation here. Instead, we only briefly review the results we need for this paper. For more details, we refer the reader to Ref. [10].

The density evolution equation and thermal averaging
We consider the coannihilation of N species of particles (χ i , i = 1, . . . , N ) with masses m i and internal degrees of freedom (statistical weights) g i (see Appendix B for conventions on degrees of freedom adopted in this paper). We order the masses in increasing order, m 1 ≤ m 2 ≤ · · · ≤ m N , and use m 1 and m χ interchangeably for the mass of the lightest particle. If the lightest particle is stable and the others decay into it, instead of considering the thermal history of each particle separately, we can follow the evolution of the sum of the number densities. The problem is then formulated in terms of the density evolution equation [7,10] dn dt = −3Hn − σ eff v n 2 − n 2 eq (3.1) where σ eff v is the effective thermally-averaged annihilation cross section, H is the Hubble parameter, and n is the total number density summed over all coannihilating particles. The effective thermally-averaged annihilation cross section is where n eq is the equilibrium number density, which in the Maxwell-Boltzmann approximation (which is a good approximation for our case) reads and A is the annihilation rate per unit volume at temperature T Here K i (x), (i = 1, 2), are the modified Bessel functions of the second kind of order i; s = 4p 2 eff + 4m 2 χ is the usual Mandelstam variable giving the center-of-mass energy squared for the χ 1 χ 1 system; p eff is the center-of-mass momentum for the χ 1 χ 1 system; and W eff is the effective annihilation rate obtained by summing over all annihilation and coannihilation channels, For the coannihilation of particles i and j, W ij is the annihilation rate per unit volume and unit time given by 5 is the common magnitude of the 3-momentum of particle i and j in the center-of-mass frame of the pair χ i χ j . Defining W ij (s) = 0 for s ≤ (m i + m j ) 2 , the radicand in Eq. (3.5) is never negative. The effective thermally averaged annihilation cross section can then be written as This expression is very similar to the one in the case of no coannihilations given in Gondolo and Gelmini [9] (and correctly reduces to it in the absence of coannihilations). The only differences are in the denominator and in the replacement of the annihilation rate with the effective annihilation rate. The key feature in Eq. (3.8) is the definition of an effective annihilation rate independent of temperature, with p eff as integration variable. This gives a remarkable calculational advantage, as W eff can be tabulated in advance, before taking the thermal average and solving the density evolution equation. The steps we implement to compute the relic density are the following. We first rephrase the density evolution equation (3.1) as an equation for Y = n/s, with s being the entropy density. Then we tabulate W eff including thresholds, resonances and coannihilations, and spline it. As a third step we solve the density evolution equation starting from the boundary condition which sets particles in equilibrium at the temperature T = m χ /2 (a specially-devised implicit method is used in the numerical integration of the density equation). This implies actually a double integration since at each temperature step we need to calculate the thermal average σ eff v (this integration is relatively fast as we use the W eff tabulation). The integration of the density equation is performed until Y has reached a constant value (the freeze out value) which we finally convert into the value of the relic abundance.
As stated earlier, we estimate our calculation of the relic density to be accurate to about 1%. We base this estimate on the precision with which we do the tabulation of the effective annihilation rate W eff and on the precision with which we perform the numerical integrations. For example, we explicitly make sure that resonances and thresholds are tabulated and integrated with such a precision that the end result is accurate to at least 1%. One should keep in mind that this is the accuracy of the calculation of the neutralino relic density starting from given values of masses, widths and couplings of MSSM particles. In the mSUGRA framework considered here, the accuracy of the sparticle masses for a given set of input parameters is less than 1% (see e.g. [38] for a comparison of different RGEcodes) and the main uncertainty in evaluating the neutralino relic density in mSUGRA thus comes from the RGE code ISASUGRA.

A few examples of coannihilation effects
For illustrative purposes only, we rewrite Eq. (3.8) in the form: where E eff = p 2 eff + m 2 χ is the energy per particle in the center of mass for χ 1 χ 1 annihilations. The term we factor out, W eff /4E 2 eff , can be thought of as an effective σv term (compare with Eq. (3.6)). In the p eff → 0 limit, W eff /4E 2 eff reduces to the χ 1 χ 1 annihilation rate at zero temperature, which is the relevant quantity in indirect dark matter detection searches. The term κ we introduced in Eq. (3.9) contains the Boltzmann factor and the phase-space integrand term and can be regarded as a weight function, at the temperature T , that selects which range of p eff is important in the thermal average. As the phase-space integrand term dominates at small p eff and makes κ go to 0 in the p eff → 0 limit, κ shows a peak at an intermediate p eff and then rapidly decreases due to the Boltzmann suppression;  Table 2 in Appendix C). The solid line shows the effective annihilation cross section W eff /4E 2 eff as a function of momentum p eff , while the dashed line shows the thermal weight factor κ(p eff , T ). The thermally-averaged annihilation cross section is the integral over p eff of the product of the two. Note that when including coannihilations, not only new thresholds appear, but the freeze-out temperature is also changing, meaning that we sample a different region of the annihilation cross section. For this model, the relic density with coannihilations is Ω χ, coann h 2 = 0.135 and that without is Ω χ, no coann h 2 = 1.43.
the position and height of the peak depends on the temperature considered and on the particles involved.
We are now ready to show some examples of coannihilation effects. As already mentioned, the examples we display have the lightest neutralino as the LSP and are in the mSUGRA framework. In Fig. 1a we consider a case in which the neutralino, with mass of about 400 GeV, is nearly mass degenerate with the lightest stau. The lightest selectron, the lightest smuon and the lightest stop are relatively close in mass as well. (To fully specify the example models we present, the model parameters and some properties are given in Table 2 in Appendix C. The model in Fig. 1 is model A in that table.) The solid curve shows W eff /4E 2 eff , and one can nicely see coannihilations appearing as thresholds at √ s equal to the sum of the masses of the coannihilating particles (just as final state thresholds do). As usually happens when considering coannihilation effects with neutralinos as the LSP, the χ 0 1 -χ 0 1 contribution to W eff is small compared with the one provided by the coannihilating particles. The role of coannihilating particles can be quantified better with a look at the function κ (dashed curve, in units of GeV −1 , and with relative scale shown on the right-hand side of the figure). The factor κ is plotted at the freeze out temperature, defined as the temperature at which the abundance of the relic species is 50% higher than the equilibrium value 6 , in this case T = m χ /24.3. On the top of the panel, the tick mark labelled '0' indicates the position of the momentum p max eff corresponding to the maximum of κ, while the other tick marks indicate the momenta p (n) eff at which κ is 10 −n of its maximum value, κ(p (n) eff )/κ(p max eff ) = 10 −n . The tick marks provide a visual guide to the interval in p eff which is relevant in the thermal averaging. The integral of the product of W eff /4E 2 eff and κ gives σ eff v thermally averaged at the freeze out temperature (shown in the figure as a horizontal dotted line). This is the quantity which is sufficient to get a rough indication of the neutralino relic abundance through the rule of thumb [41] In Fig. 1b we consider the same model but ignore coannihilation effects. One can see that W eff /4E 2 eff is now, on average, much smaller, and therefore one can expect the relic abundance to be higher. This is indeed the case, with a shift from Ω χ h 2 = 0.135 including coannihilations (left panel) to Ω χ h 2 = 1.43 when coannihilations are neglected (right panel). Note, however, that the change in Ω χ h 2 is smaller than what one would naively expect from comparing the solid curves in the two panels. This is due to the fact that there is a significant change in the freeze out temperature as well, from T = m χ /24.3 to T = m χ /21.7. The weight function κ for this new temperature is shown in the figure, and comparing it to the one in the left panel, one clearly sees the change in normalization (partially due to the change in the number of degrees of freedom involved in the two cases, see the denominator in Eq. (3.8)) and in width (for the latter effect, note the shift in the scale shown on the top of the figure, while the displayed range in p eff has been kept fixed). The net result is that σ eff v at the freeze out temperature is lowered by just about an order of magnitude, and then the increase in the relic abundance is of the same order. From this discussion it is evident that a very accurate solution of the density evolution equation is needed to claim good accuracy on the estimate of the relic density.
We mentioned that usually, when considering neutralino dark matter, W eff increases sharply when coannihilating particles are included. The reason is that the coannihilating particles typically have non zero electric or colour charges, while the neutralino interacts only weakly. In the MSSM there are a few exceptions to this general trend of increasing W eff (see [10]), and we find one such exception in the mSUGRA framework as well. In Fig. 2a we show a model (model B in Table 2 in Appendix C) for which the neutralino mass is slightly below half of the H 0 3 mass and tan β is large. The χ 0 1 -χ 0 1 annihilation rate is dominated by the s-channel H 0 3 resonance and is quite large. (There is also a H 0 1 resonance at this mass, but it is subdominant with respect to the H 0 3 resonance). The effect of stau coannihilating particles comes on top of that, but the total contribution to W eff is just about of the same order as from the χ 0 1 -χ 0 1 term. If we now go to the case when coannihilations are neglected, Fig. 2b, we see that the freeze out temperature remains about the same, but still there is a shift in the normalization of the weight function due to different numbers of degrees of freedom in the two cases (see the denominator in Eq. (3.8)). We find that the net effect is a slight increase in σ eff v , which suggests that the relic abundance should be smaller if coannihilations are neglected. This is confirmed by the full solution of the density evolution equation: the calculation that includes coannihilations yields Ω χ h 2 = 0.155, while if stau coannihilation is neglected one obtains Ω χ h 2 = 0.137. Therefore one should keep in mind since we solve the density evolution equation numerically.  Table 2 in Appendix C), where coannihilations increase the relic density (in this case from 0.137 to 0.155).
that although the general trend is that coannihilation effects lower the neutralino relic abundances, in some cases they can increase it.
Note that due to the Boltzmann suppression of heavier particles, we do not need to include all supersymmetric particles in the calculation. By extensive scans of the mSUGRA parameter space and estimating the effects of including different particles in the calculation, we have found that a convenient criterion to have an accuracy of 1% or better on the relic density in cases that are cosmologically interesting is to include all supersymmetric particles with a mass below 1.5m χ . This is shown in Fig. 3 where we plot the relative difference in relic density without coannihilations and with coannihilations, versus the mass ratio of the lightest coannihilating particle and the neutralino in our sample of models: the coding in relic density shows that, for cosmologically interesting cases, only mass differences up to about 40% are important, hence the cut at the mass ratio 1.5 (vertical dashed line in the figure) is sufficiently conservative. If a 1% accuracy is required even for models that are cosmologically disfavoured, one would need to raise the cut to about 1.7: the coannihilation effects one picks in this way, however, correspond to the χ 0 1 − χ 0 1 annihilation cross section being very suppressed and, even adding coannihilation terms, the relic abundance still remains Ω χ h 2 ≫ 1. In Fig. 4 we show the effective annihilation cross section for such an example (model C in Table 2 in Appendix C). In this case stop coannihilations are important even though the stop mass is about 50% higher than the neutralino mass, but still the relic density, which is Ω χ h 2 = 73. 2   give is the easiest to implement in the numerical calculation, but clearly depends on the strength of the interactions we are dealing with; we can say that it is perfectly safe in the scheme we are dealing with, but need not be valid in other schemes.
A final check we perform is the following: sometimes in the literature only the next-tolightest (NLSP) sparticle has been included in the coannihilation calculation and one may question whether this can be considered a fair approximation. There is reason to believe that also heavier sparticles might change the relic LSP density by a non-negligible amount if they are either close in mass to the NLSP or if they have larger coupling strength than the NLSP (as e.g. χ ± 2 has compared to χ 0 2 ). Whether or not it is a good approximation to include only NLSP coannihilations should therefore be checked case by case. Even restricting to the cases of coannihilating particles with comparable couplings, such as staus, smuons and selectrons, it is not straightforward to provide a firm criterion telling when to include all of them or just the NLSP (i.e. theτ in this case) in the calculation. In practice, we find regions in the mSUGRA parameter space when smuon and selectron contributions  Table 2 in Appendix C) where stop coannihilations are important even though the lightest stop is about 50% heavier than the lightest neutralino.  Table 2 in Appendix C) where inclusion of all three lightest sleptons is important to get the correct relic density.
give corrections to the relic density of the order of 10% or higher. An example is given in Fig. 5 (model D in Table 2 in Appendix C). In this example, the lightest smuon and selectron are almost degenerate in mass with both the lightest stau and lightest neutralino. Including in the calculation the coannihilation with these three sfermions result in a relic density Ω χ h 2 = 0.109. If instead we only include χ 0τ andττ coannihilations, the relic density would be 0.128, i.e. we would be ∼ 18% off from the correct value. Hence, as this example has shown, it can be important to include the coannihilations of more than one particle. To be on the safe side, we always include all particles with masses up to 1.5m χ .

The role of coannihilations in the mSUGRA framework
The mSUGRA framework [42] has been extensively discussed in the literature. We summarize here the main features that are relevant for understanding the results of the relic abundance computation. In line with most previous analyses, we sample the 5-dimensional mSUGRA parameter space choosing a few values of tan β and A 0 , and slices along the m 1/2 , m 0 planes for both sign(µ).
Consider first the case A 0 = 0. Two regimes with cosmologically interesting relic abundances have been identified: The region m 0 ∼ < m 1/2 and the region m 0 ≫ m 1/2 . In the first region, the lightest neutralino is a quasi pure bino with mass set essentially by m 1/2 alone; the parameter m 0 sets the sfermion mass scale, with the slepton sector lighter than the squark sector and with the lightest stau always being the lightest sfermion, possibly lighter than the lightest neutralino if m 0 ≪ m 1/2 . In the second regime, the relevant region is a narrow band, sometimes dubbed the "focus point" region [43], close to the region where there is no radiative electro-weak symmetry breaking: in such a band the parameter µ is driven to small values and forces a mixing between the gaugino and Higgsino sectors. As a consequence the lightest neutralino may contain a large Higgsino fraction and the nextto-lightest supersymmetric particle is a chargino. Large values of A 0 can introduce a third regime in the intermediate m 0 range: off-diagonal entries in the stop mass matrix can become sufficiently large to drive the lightest stop to masses smaller than the mass of the lightest stau or even the mass of the lightest neutralino. On top of these generic trends, details are sensitive to the value of tan β and the sign(µ).
In each of the regimes above, the mass splitting between the LSP, which we require to be the lightest neutralino, and the next-to-lightest SUSY particle can be small enough for coannihilation effects to become important. We label the three regions as the slepton coannihilation region, the chargino coannihilation region and the stop coannihilation region, and describe each of them separately.

Slepton coannihilations
We consider first the case A 0 = 0 and m 0 ∼ < m 1/2 . The slepton coannihilation region was recognised in Ref. [14] and has been the focus of several recent studies, including [15,13,21,20,44]. As an example, we take tan β = 10 and study the m 1/2 − m 0 plane for both signs of µ. In the top panels of Fig. 6, we plot isolevel curves for the neutralino relic abundance, including coannihilation effects, for 8 different values of Ω χ h 2 starting from Ω χ h 2 = 0.3 and decreasing down to Ω χ h 2 = 0.025. The latter value corresponds to the case in which neutralinos would be a subdominant dark matter component in the Universe but could still account for a major part of the dark matter in galaxies. The shaded area on the left in each panel is excluded by accelerator constraints, while the shaded area towards the bottom right corner in each panel is removed because in this region the LSP is the lightest stau rather than the lightest neutralino: its upper bound marks the line along which the (bino-like) neutralino and the lightest stau have equal mass. We can give a schematic interpretation of the results displayed starting with the isolevel curves on the top left corner of each panel, where all isolevel curves converge to a narrow band. There, the model has a relatively heavy sfermion sector, and the lightest neutralino mass is just a few GeV larger than half the Z 0 boson mass. The bino pair annihilation rate into fermions is dominated by the diagram with Z 0 in the s-channel at energies just slightly displaced from the Z 0 resonance: this resonant annihilation leads to acceptable values of Ω χ h 2 in a narrow band. Thus Ω χ h 2 is very sensitive to the parameter m 1/2 as the value of the bino mass has to be fine-tuned in such a way that the thermally averaged cross section picks the right portion of the resonance; for larger m χ , σv drops rapidly and Ω χ h 2 becomes large, closer to the resonance σv becomes exceedingly large and Ω χ h 2 very small. When we follow the isolevel curves from the top left down to smaller values of m 0 , sfermion masses decrease and the amplitude for neutralino annihilations into fermions mediated by sfermion exchange in t-and u-channels increase, eventually becoming dominant. The largest contribution to σv is given by the τ − τ + final state which is mediated by the lightest sfermion, in this case the lightest stau. Such an increase in the cross section when moving to lower m 0 can be compensated by increasing the bino mass, i.e. by shifting to larger m 1/2 . This explains the trend of the isolevel curves starting to align along the diagonal in the figure. Further down this diagonal the mass of the lightest stau is more or less constant, but the mass of the heaviest stau increases. This increased mass splitting between the staus causes the cross section to increase, but this increase is again compensated by a larger neutralino mass. Before reaching the lower right corner where the lightest stau is the LSP, coannihilation effects take over: σ eff v becomes rapidly dominated by neutralinoslepton and slepton-slepton contributions. The coannihilation cross section decreases with increasing initial state masses, an effect which can be compensated by decreasing the mass splitting between the LSP and the stau. (This effect is further discussed in connection with Fig. 8.) Hence, the isolevel curves bend almost parallel to the bound of the excluded region, confined to a band which gets progressively narrower towards larger values of m 0 and m 1/2 . As can be seen by comparing the top left with the top right panel, the flip in sign(µ) does not alter this overall picture but just slightly displaces the position of the isolevel curves in the m 1/2 − m 0 plane. The region where coannihilation effects become important is highlighted in the bottom panels of Fig. 6. There dashed lines show the isolevel curves for the relative difference in relic abundance computed neglecting and including coannihilations, ∆Ω/Ω ≡ (Ω χ, no coann − Ω χ, coann )/Ω χ, coann . The values we display span from 1% to 500%. Isolevel curves for the relative mass splitting between the lightest stau and the lightest neutralino are shown as well and, as expected, the correlation with the isolevel curves of ∆Ω/Ω is evident.
In Fig. 7 we repeat the exercise for tan β = 30 and the picture we see is analogous to what we found in the tan β = 10 case. The shift of the isolevel curves to slightly larger m 0 values is mainly due to the fact that, comparing corresponding points in the plane, sleptons are driven to slightly lighter masses at larger tan β.
From Figs. 6 and 7, as well as from other sample checks performed for other values of tan β, we can infer that, as a rule of thumb, in the A 0 = 0 case a 25% mass splitting between the bino and the lightest stau is approximately the borderline below which slepton coannihilations have to be included for an estimate of Ω χ with a 1% accuracy level.
We have briefly mentioned the reason why the stau coannihilation region takes the form of a tail extending to large values of m 1/2 . When m 1/2 , and therefore m χ , is increased, it is necessary to increase the effect of coannihilations, and therefore to lower the relative mass splitting between the LSP and the NLSP, if the relic density should not increase. Fig. 8 gives a detailed illustration of this effect. We select here two models with A 0 = 0, tan β = 10, sign(µ) positive, and the m 1/2 , m 0 pair chosen in such a way that the relic abundance for both models is Ω χ h 2 = 0.115, but ∆Ω/Ω is 100% for one model (left panel, model E in Table 2) and 1000% for the other (right panel, model F in Table 2). Analogously to some of the figures in Section 3.2, we plot W eff /4E 2 eff , as well as the weight function κ computed at the freeze out temperature. Going from the left to the right panel, m χ increases from 138.5 GeV to 371.1 GeV, the lightest stau from 148.0 GeV to 371.8 GeV. The mass splitting consequently goes from 6.8% to 0.21%. The increase in neutralino mass  Table 2 in Appendix C) with tan β = 10 and Ω χ h 2 = 0.115. The figure illustrates how coannihilation effects become more important the further out in the coannhilation tail we get. See text for a discussion. causes a sharp drop in the χ 0 1 -χ 0 1 annihilation cross section 7 ; on the other hand such a drop can be compensated in the thermally averaged cross section by increasing the role of coannihilating particles. Even though their annihilation cross sections have also been slightly reduced by the mass shift, we have forced their contribution to the thermal average to be larger by shifting them deeply in the region where the weight function κ is large. More precisely, we have required the coannihilation threshold to move to lower p eff , i.e. to have a lower mass splitting. 8 With these changes in the mass splitting, we have found two models with about the same thermally averaged cross sections, and then, with the rule of thumb Ω χ h 2 ≃ 10 −27 cm 3 s −1 / σ eff v one can expect similar relic densities. In fact, for the two cases discussed here we have chosen the parameters to give exactly the same neutralino relic density Ω χ h 2 = 0.115 from the full calculation. 7 As W eff /4E 2 eff at p eff = 0 is the neutralino annihilation rate at zero temperature, which is the relevant quantity for indirect dark matter detection, we find that, roughly speaking, it is a factor of 140 harder to detect indirectly the model on the right, a factor of 20 due to the reduced value of the cross section, plus a factor of 7 due to the square of the ratio of mχ's (which takes into account the neutralino number density scaling for a fixed dark matter mass density). and Ωh 2 less than: Figure 9: For a set of values of Ω χ h 2 , we show the importance, ∆Ω/Ω, of including slepton coannihilations, as well as the mass splitting between lightest neutralino and lightest stau, both as functions of the neutralino mass. We have here chosen tan β = 30, µ > 0, A 0 = 0 and low m 0 (the latter to be in the stau coannihilation region). In the right-hand panel, we also indicate by arrows where the upper limit on the mass would be if coannihilations were not included (for these limits the vertical axis is meaningless). The shift to higher masses when coannihilations are included is clearly seen.
From the above discussion it should also be evident that we cannot play the same game to arbitrarily large neutralino masses. At some point one reaches the edge at which the role of coannihilating particles cannot be further strengthened and one finds an upper bound on the neutralino mass (for that specific choice of tan β, sign(µ) and Ω χ h 2 ). The introduction of slepton coannihilations have extended the cosmologically allowed mass interval of the bino-like LSP in the mSUGRA context. This fact has already been pointed out by several authors. In the work of some authors [20,44] it seems that the neutralino mass might be unbounded from above, but the majority [21,14,15] found that there was still an upper limit to the allowed mass. Given the high precision in the relic abundance calculation we present in this analysis, we are able to make a qualitative analysis of the high m χ region, and indeed we find upper limits to the neutralino mass. We will show our results in the case tan β = 30 and positive µ considered in Fig. 7. Keeping in mind that we are in the m 0 ∼ < m 1/2 regime, there is a one to one correspondence between the pair (m 1/2 ,m 0 ) and the pair (m χ , (mτ − m χ )/m χ ) or (m χ ,∆Ω/Ω). We therefore show first in Fig. 9a the difference in relic density, (Ω χ, no coann − Ω χ, coann )/Ω χ, coann , versus the neutralino mass for a few values of Ω χ h 2 . We clearly see the importance of stau coannihilations in the high mass region. In Fig. 9b  much larger values, when coannihilation effects are included, is evident, as well as the fact that we do find a new maximum value of m χ . The results for negative µ are very similar, while those for tan β = 10 are analogous, but show slightly more stringent upper bounds on the neutralino mass.

Chargino coannihilations
In the focus point region, the value of the soft mass parameter (at the electro-weak scale) for the Higgs doublet that couples to up-type quarks, m Hu , is naturally of the electro-weak scale, regardless of m 0 [43]. As a consequence, the parameter µ is forced to be light, and can be at the level of the gaugino mass parameter m 1/2 or even lower. This implies that the neutralino LSP may have a large Higgsino fraction and be nearly degenerate in mass with the lightest chargino and the next-to-lightest neutralino. Especially at higher m 0 -values, the Higgsino fraction can be very large, close to one. Hence, in this high m 0 focus point region, chargino (and neutralino) coannihilations are expected to be important. Chargino coannihilations have been extensively studied in the generic MSSM context [45,46,10], but have been rarely stressed in the mSUGRA framework (although they are included in some recent analyses, e.g. [20,44,47]). In Fig. 10 we show the lower part of the focus point region for tan β = 30 and A 0 = 0. The top-left corners of these figures are excluded due to no radiative electro-weak symmetry breaking, but close to that region, we see the focus point region emerge. In this region, the Higgsino fraction is usually small, but non-negligible, and the same is true for the effect of chargino coannihilations. This is the part of the focus point region most often discussed in the literature. However, if we continue to higher masses, we get a band of  cosmologically interesting relic densities where the Higgsino fraction increases as we go up in mass (at the highest masses, it is close to 1). In this case, coannihilations with the lightest chargino (and the next-to-lightest neutralino(s)) occur and are important. In Fig. 11a we show the relic density isolevel curves without including coannihilations and in b) with coannihilations included in this high-mass region. (Note the scale m 0 −2m 1/2 on the y-axis, which is chosen to clearly show the coannihilation region.) We clearly see the effect of coannihilations pushing the cosmologically allowed region to higher masses. Note that in the focus point region of the parameter space, there is no longer a simple relation between the neutralino mass and m 1/2 . Some neutralino mass contours are therefore indicated in the figures. As shown in Fig. 12, we find that (just as for the MSSM case [10]) chargino coannihilations are non-negligible below the W mass where the dominant neutralino-neutralino annihilation channel is annihilation into fermion anti-fermion pairs via s-channel Z-boson exchange. The coupling for a Higgsino-like neutralino to the Z-boson is very suppressed though, whereas that of charginos is not and coannihilations could then give a big effect. However, in the mSUGRA framework, the neutralinos below the W mass are rather mixed than Higgsino-like and the effect is not as dramatic as in the MSSM. Above the W mass, where annihilation into W + W − dominates for neutralino-neutralino annihilation (annihilation into Z 0 Z 0 is also significant above the Z mass), the effect of coannihilations is smaller since the annihilation rate into W + W − is comparable for neutralinos and charginos. Still, annihilation into fermion anti-fermion pairs is not suppressed for chargino-chargino annihilation and thus the chargino-chargino annihilation cross section is typically slightly higher than the neutralino-neutralino one. Hence, coannihilations do change the relic density even for higher masses, and the further up in mass (and down in mass splitting) we go, the more important they are. In Fig. 12a we plot ∆Ω/Ω versus the neutralino mass. We see that ∆Ω/Ω is of the order of 1% for models with m χ just below the W mass and drops when we get above m W . Then, ∆Ω/Ω slowly increases as we go up in mass with a small drop at m χ = 175 GeV, where the neutralino annihilation into tt becomes significant, thereby reducing the importance of coannihilations somewhat. In Fig. 12b we show the relative mass splitting between the lightest chargino and the lightest neutralino versus the neutralino mass. We here clearly see the same effect as we have seen with slepton coannihilations, i.e. as we go up in mass, we need to make the mass splitting smaller to maintain the same value of the relic density. As for slepton coannihilations, we cannot continue this game to arbitrarily high masses and as we see in Fig. 12b, we get upper limits on the neutralino mass for a given relic density. The corresponding curves for µ < 0 are very similar and we do not show them separately.

Stop coannihilations
The new version of DarkSUSY offers the possibility to include all squark coannihilations in the relic density calculations, but only stop coannihilations have proven to affect the result for the mSUGRA framework. In this context, the lightest stop is always the lightest squark; its mass is usually much larger than the LSP mass, unless off-diagonal entries in the stop mass matrix become sufficiently large to drive the lightest stop mass to small values: this can happen when |A 0 | is large. If the mass of the lightest stop is close to the neutralino mass, stop coannihilation effects become important in the neutralino relic density calculation [12,16,48]. In Fig. 13  In the hatched region to the left,t 1 is the LSP (except for a narrow band at the lower part of the right edge where the neutralino is the LSP, but where the H 0 2 mass constraint is not fulfilled). In the hatched region at the bottom right,τ 1 is the LSP. In the almost vertical band of interesting relic densities, a lightt 1 is important in two respects, both by boosting the χ 0 1 − χ 0 1 annihilation above the tt threshold (as seen in a)) and by coannihilations (as seen in b)). The most prominent effect of thet 1 coannihilations is the narrow band for neutralino masses less than m(t) ∼ 174 GeV.
for tan β = 10, A 0 = −2750 GeV and µ > 0. The hatched region is excluded mainly due to the neutralino not being the LSP here. The lower right part has theτ 1 as the LSP, whilet 1 is the LSP in most of the 'bump' to the left. A narrow band close to the lower part of the right edge of this 'bump', i.e. to the left of the Ωh 2 curves, is excluded because of the mass limit on H 0 2 . In other words, the H 0 2 mass limit cuts away some part of the parameter space with the smallest mass splitting betweent 1 and the LSP. In Fig. 13a we show the relic density isolevel curves without including coannihilation processes. We have also included the isolevel curves for the neutralino mass. It is clearly seen that the neutralino-neutralino annihilation cross section dramatically increases at the tt threshold. The relic density then decreases to cosmologically interesting values in a band were the mass splitting between thet 1 and the neutralino is small. It is the t-channel exchange of the light stop that boosts the neutralino annihilation into tt. In Fig. 13b we show the results of including coannihilation processes. Three effects are seen. First of all, the stop coannihilations open a narrow band of cosmologically interesting density for neutralinos below the top mass. Secondly, in the almost vertical band above m(χ 0 1 ) = m(t), that was already present in Fig. 13a, the relic density is decreased further due tot 1 coannihilations. Finally, we also find interesting relic densities for high m 1/2 and low m 0 whereτ 1 -coannihilations lower the relic density. The same three features were also found by Ellis, Olive and Santoso,  Fig. 13, we plot ∆Ω/Ω ≡ (Ω no coann − Ω coann )/Ω coann versus the neutralino mass. We clearly see the effect that below m t , we only get interesting relic densities if coannihilations witht are included. When we get close to m t , though, even χ 0 1 χ 0 1annihilation gets efficient and the relative effect of including coannihilations go down. When coannihilations withτ get important at higher masses, the relative difference in Ω χ goes up again. The small dip at m χ ≃ 280 GeV corresponds to the low-m 0 corner where the curves in Fig. 13 bend sharply and arises because χ 0 1 χ 0 1 -annihilation to τ + τ − gets significant here.
[16] Fig. 6d. Beside the trivial sign difference in the convention used for A 0 , we also have chosen the numerical value of A 0 slightly lower than done in [16], in order to obtain a figure as similar to theirs as possible. The possible explanation for the required change of |A 0 | is that we use different RGE codes, and a given |A 0 | therefore results in different low energy masses. This is not only true for the neutralino and sfermion masses but also for the Higgs masses. As mentioned, part of the edge of the excluded region is set by the H 0 2 mass limit. The excluded region would be even larger if we applied the Standard Model Higgs limit (∼ 114 GeV) as done in [16] instead of the value m(H 0 2 ) ∼ 92 GeV that applies for tan β = 10 [33].
To see more clearly the effect of coannihilations, in Fig. 14 we show the difference in relic density, (Ω χ, no coann − Ω χ, coann )/Ω χ, coann , versus the neutralino mass along the isolevel curves shown in Fig. 13. We see that in the coannihilation region in the upper left corner of Fig. 13b,t 1 -coannihilations are indeed very important changing the relic density by a factor of several 100. When we get to the tt threshold, the effect of coannihilations decreases drastically since the χ 0 1 χ 0 1 annihilation cross section is high. When the neutralino mass is increased further, the importance of coannihilations again increases and above ∼ 275 GeV theτ -coannihilations start getting important.

Conclusions
We have presented a novel tool for calculating the neutralino relic abundances with an estimated precision of 1% or better, assuming masses, widths and couplings of particles in the MSSM are given. We allow for the most generic coannihilation effect in the framework of the MSSM, applying at the same time the state of the art technique to trace the freezeout of a species in the early Universe. The code for numerical computations will be publicly available in the near future together with a new extended version of the DarkSUSY package [8].
As a first example, in this paper we have discussed results valid in the mSUGRA setup, the rather restrictive framework most recent analyses of coannihilations have focussed on, but, at the same time, general enough to illustrate most of the effects that can arise in the MSSM. For a selection of mSUGRA models we have given the reader a visual guide into general trends and their exceptions, such as the tendency of coannihilation effects to lower the neutralino relic abundances and the cases in which the opposite can happen. We have then performed a broad study of the neutralino relic density in the mSUGRA parameter space. The features we have spotted are in agreement with the picture emerging from previous analyses; in particular we have discussed the cases of neutralino coannihilations with sleptons (most notably stau leptons), with stop squarks, and with charginos. Especially the chargino coannihilations have been treated in more detail than before.
The accuracy of the calculation we performed allowed us to go into great details in the cases presented. Novel features we have discussed include, e.g.: (i) a rule of thumb to discriminate the case when slepton coannihilations are relevant: 25% mass splitting between lightest stau and lightest neutralino for a 1% accuracy on the cosmologically allowed neutralino relic density; (ii) a shift to heavier masses of the cosmological bound on the neutralino mass in both the slepton and chargino coannihilation cases: in the tan β = 30 case we studied, we found 565 GeV as an upper limit on the neutralino mass from the cosmological bound Ω χ h 2 < 0.2 in the regime m 0 ∼ < m 1/2 , where slepton coannihilations are important, and 1500 GeV in the regime m 0 ≫ m 1/2 , where chargino coannihilations occur; (iii) a correlation between mass splittings, differences in Ω χ h 2 and neutralino masses in all of the cases presented.
In an upcoming paper, we will investigate the effects of the coannihilation channels in a more general MSSM context.

Acknowledgements
We thank L. Bergström for discussions regarding the calculation of the cross section of squarks annihilating into two gluons. We also thank K. Matchev for comments. J.E. thanks the Swedish Research Council for support. P.U. was supported in part by the RTN project under grant HPRN-CT-2000-00152 and by the Italian INFN under the project "Fisica Astroparticellare". P.U., J.E., and P.G. thank the Kavli Institute for Theoretical Physics, where parts of this work were completed, for its hospitality. This research was supported in part by the National Science Foundation under grant No. PHY99-07949.

A. Included coannihilations
Here we list all 2 → 2 tree-level coannihilation processes with sfermions, charginos and neutralinos. All the processes are included in the DarkSUSY code, but not all of them have been used in the calculation we report on. As mentioned earlier, we have found that an accuracy of 1% on the relic neutralino density in the cosmologically interesting interval is obtained by including coannihilations for particles lighter than 1.5m χ . Consequently, some sparticles never get included in the initial state. For the sparticles that satisfy the mass difference criterium, we have included all coannihilation processes and for each of these, all the exchange channels. Furthermore, the DarkSUSY relic density code always includes the one loop neutralino annihilation into γγ, γZ and gg.
It should be noted that we have not included all flavour-changing charged current diagrams. The DarkSUSY vertex code for the charged current couplings is written in a general form that includes all possible flavour-changing (and flavour-conserving) vertices. The flavour-conserving couplings are much larger than the flavour-changing. For the sfermion coannihilations with charged currents we only take the flavour-conserving contributions, while for the chargino coannihilations we include the flavour-changing contributions as well. In a future version of DarkSUSY, we may as well include the flavour-changing processes for the sfermion coannihilations, even if they are not expected to be important.
We have used the notationf for sfermions and f for fermions. Whenever the isospin of the sfermion/fermion is important, it is indicated by an index u (T 3 = 1/2) or d (T 3 = −1/2). The sfermions have an additional mass eigenstate index, that can take the values 1 and 2 (except for the sneutrinos which only have one mass eigenstate). A further complication to the notation is when the sfermions and fermions in initial, final and exchange state can belong to different families. Primes will be used to indicate when we have this freedom to choose the flavour. So, e.g.f u and f u will belong to the same family whilef u and f ′ u can belong to the same or to different families. Note that the colour index of (s)quarks as well as gluons (g) and gluinos (g) is suppressed.
Besides the sfermions we also have neutralinos and charginos in the initial states. The notation used for these are the following. The neutralinos are denoted by χ 0 j with the index running from 1 to 4. The charginos are similarly denoted χ ± j with the index taking the values 1 and 2.
In the table in this appendix, a common notation is introduced for gauge and Higgs bosons in the final state. We denote these with B with an upper index indicating the electric charge. So B 0 means H 0 1 , H 0 2 , H 0 3 , Z, γ and g while B ± is H ± and W ± . We will use additional lower indices m and n when we have more than one boson in the final state. Thus indicating that the bosons can be either different or identical. Note that the case of two different bosons also includes final states with one gauge boson and one higgs boson.
The table has been made very general. This means that when a set of initial and final state (s)particles have been specified, the given process might not run through all the exchange channels listed for the generic process. Exceptions occur whenever an exchange (s)particle does not couple to the specific choice of initial and/or final state. As an example we see that since the photon does not couple to neutral (s)particles, none of the exchange

C. Example models
Here we list the parameters and some properties of the example models considered in some of the figures.