Full one-loop corrections to the relic density in the MSSM: A few examples

We show the impact of the electroweak, and in one instance the QCD, one-loop corrections on the relic density of dark matter in the MSSM which is provided by the lightest neutralino. We cover here some of the most important scenarii: annihilation into fermions for a bino-like neutralino, annihilation involving gauge bosons in the case of a mixed neutralino, the neutralino-stau co-annihilation region and annihilation into a bottom quark pair. The corrections can be large and should be taken into account in view of the present and forthcoming increasing precision on the relic density measurements. Our calculations are made possible thanks to a newly developed automatic tool for the calculation at one-loop of any process in the MSSM. We have implemented a complete on-shell gauge invariant renormalisation scheme, with the possibility of switching to other schemes. In particular we will report on the impact of different renormalisation schemes for tan beta.


Introduction
The last few years have witnessed spectacular advances in cosmology and astrophysics confirming with an unprecedented level of accuracy that ordinary matter is a minute part of what constitutes the Universe at large. At the same time as the LHC will be gathering data, a host of non collider experiments will be carried out in search of Dark Matter, DM, (either direct or indirect) as well as through ever more precise determination of the cosmological parameters. In this new paradigm, the search for DM at the LHC is high on the agenda, as is of course the search for the Higgs. In fact these may be two facets of the New Physics that provides a resolution to the hierarchy problem posed by the Higgs in the Standard Model, SM. The epitome of this New Physics is supersymmetry which among many advantages furnishes a good DM candidate through the lightest neutralino,χ 0 1 . If future colliders discover supersymmetric particles and probe their properties, one could predict the dark matter density of the Universe and would constrain cosmology with the help of precision data [1,2] provided by WMAP [3] and PLANCK [4]. It would be highly exciting if the precision reconstruction of the relic density from observables at the colliders does not match PLANCK's determination, this would mean that the post-inflation era is most probably not entirely radiation dominated [5]. Already now the accuracy on the relic dark matter density is about 10% from WMAP and will soon be improved to about 2% from PLANCK. Such level of accuracy must be matched by precise theoretical calculations. From the particle physics point of view this means precision calculations of the annihilation and co-annihilation cross sections at least at one-loop. Quite sophisticated codes now exist [6,7] for the calculation of the relic density, however they are essentially based on tree-level cross sections with the inclusion of some higher order effects essentially through some running couplings, masses or some effective couplings. Some of these corrections [6] have already been shown to be essential like the corrections to the Higgs couplings that can completely change the picture in the so-called Higgs funnel (annihilation mainly through the pseudo-scalar Higgs, A). The use of other approximations needs to be justified by complete higher order calculations which contain more than just the effect of effective couplings. In a word, the level of accuracy that will soon be reached requires that one is ready to tackle in a general way a full calculation at one-loop for any annihilation (or co-annihilation) of the neutralinos in supersymmetry, just as what one has been doing for the cross sections at the colliders. The aim of this letter is to report on the progress towards automatisation of these calculations and to show and discuss some results on the one-loop corrected annihilation and co-annihilation cross sections of the LSP neutralino in the MSSM. In particular, we study here some of the most important scenarii: i) annihilation in the so-called bulk region into fermions for a bino-like neutralino, ii) co-annihilation involving the neutralino and the lightest stauτ 1 , iii) annihilation into a pair of massive gauge bosons in the case of a mixed neutralino and iv) annihilation into bb where the pseudo-scalar Higgs pole can play a role. We concentrate on the electroweak corrections, although iv) is an excuse to show how we handle some classes of QCD corrections. The couple of very recent calculations of loop corrections to the relic density tackled either QCD corrections in extreme, though highly interesting, scenarii such as annihilation into top quarks at threshold [8] and the nice study of stop-neutralino co-annihilation [9] 1 . Some important non-perturbative electroweak effects of the Coulomb-Sommerfeld type that occur for TeV winos or higgsinos with a relative mass splitting between the lightest supersymmetric particle (LSP) and the next-to-lightest supersymmetric particle (NLSP) of O(10 −4 ) have been reported in [15,16]. Let us also note that, though not to be seen as radiative corrections to the annihilation cross sections, the temperature corrections to the relic density have been considered and found to be totally negligible at the level of 10 −4 [17]. A better simulation of the cosmological equation of state to derive the effective number of relativistic degrees of freedom has been done giving corrections ranging from 1.5% to 3.5% [18] compared to the usual treatment as done in DarkSUSY [7] or micrOMEGAs [6].
2 General set-up and details of the calculation 2.1 Set-up of the automatic calculation of the cross sections Even in the SM, one-loop calculations of 2 → 2 processes involve hundreds of diagrams and a hand calculation is practically impracticable. Efficient automatic codes for any generic 2 → 2 process, that have now been exploited for many 2 → 3 [19,20] and even some 2 → 4 [21,22] processes, are almost unavoidable for such calculations. For the electroweak theory these are the GRACE-loop [23] code and the bundle of packages based on FeynArts [24], FormCalc [25] and LoopTools [26], that we will refer to as FFL for short. With its much larger particle content, far greater number of parameters and more complex structure, the need for an automatic code at one-loop for the minimal supersymmetric standard model is even more of a must. A few parts that are needed for such a code have been developed based on an extension of [27] but, as far as we know, no complete code exists or is, at least publicly, available. Grace-susy [28] is now also being developed at one-loop and many results exist [29]. One of the main difficulties that has to be tackled is the implementation of the model file, since this requires that one enters the thousands of vertices that define the Feynman rules. On the theory side a proper renormalisation scheme needs to be set up, which then means extending many of these rules to include counterterms. When this is done one can just use, or hope to use, the machinery developed for the SM, in particular the symbolic manipulation part and most importantly the loop integral routines including tensor reduction algorithms or any other efficient set of basis integrals. The results we will report are based on the development of a new automatic tool that uses and adapts modules, many of which, but not all, are part of other codes like FFL. This is the package SloopS whose main components and architecture we briefly sketch.
In this application we combine LANHEP [30] (originally part of the package COMPHEP [31]) with the FFL bundle but with an extended and adapted LoopTools [11]. LANHEP is a very powerful routine that automatically generates all the sets of Feynman rules of a 1 We do not list here loop induced annihilation processes suchχ 0 1χ 0 1 → gg [10][11][12][13]. A very recent paper discusses the QCD correction to annihilation into bb in the funnel [14], however the bulk of all contributions has been known for sometime and implemented in micrOMEGAs already.
given model, the latter being defined in a simple and compact format very similar to the canonical coordinate representation. Use of multiplets and the superpotential is builtin to minimize human error. The ghost Lagrangian is derived directly from the BRST transformations. The LANHEP module also allows to shift fields and parameters and thus generates counterterms most efficiently. Understandably the LANHEP output file must be in the format of the model file of the code it is interfaced with. In the case of FeynArts both the generic (Lorentz structure) and classes (particle content) files had to be given. Moreover because we use a non-linear gauge fixing condition [23], see below, the FeynArts default generic file had to be extended.

Renormalisation and renormalisation schemes
In the last half decade there has been an upsurge and flurry of activity constraining models of supersymmetry and other New Physics with the limit on the relic density delimiting most of the parameter space of these models. All these investigations are based on treelevel, sometimes with improved effective couplings, estimates of the relic density. Only in the last few months have some investigations [32], within mSUGRA, added a theoretical error estimate of O(10%), i.e, of the same order as the current experimental error. In these analyses based on renormalisation group running, a substantial uncertainty is due to the impact of the running itself [1,33]. Even if the weak scale spectrum is known, loop corrections to the cross sections are needed. In fact the precision one-loop calculations we are carrying will be most useful when confronting a measurement of the relic density once the microscopic properties of dark matter would have been pinned down at the collider and in direct/indirect detection. Henceforth we rely on the physical masses of the SUSY particles with the addition of some physical observables to fully reconstruct the model. We therefore work, as far as possible, within an on-shell scheme generalising what is done for the electroweak standard model. i) The Standard Model parameters: the fermion masses as well as the mass of the W and the Z are taken as input physical masses. The electric charge is defined in the Thomson limit, see for example [23]. Because we are calculating corrections to processes at a scale 2mχ0 1 ∼ 2M Z , the effect from the running electromagnetic coupling due to the light fermion masses will, alone, rescale the tree-level cross section leading to a correction of about 15% to the cross sections. The light quark (effective) masses, are chosen such as to reproduce the SM value of α −1 (M 2 Z ) = 127.77 including the light fermions contribution, which give the ∼ 15% corrections compared to the use of α(0). For the SM input masses see the last papers of Ref. [19] with the exception of m top = 170.9 GeV. We will keep this rescaling in mind. This effect can be reabsorbed by using a scheme where the effective α(M 2 Z ) eff. is used as input. ii) The Higgs sector: The Higgs sector is conceptually the trickiest. First we take M A the pseudoscalar Higgs mass as an input parameter and require vanishing tadpoles. The extraction and definition of the ubiquitous tan β, which at tree-level is identified as the ratio of the vev of the two Higgs doublet is the tricky part. Most schemes define the tan β counterterm at one-loop from a non physical quantity, such as the A 0 Z transition two-point function at q 2 = M 2 A . It has become customary to take a DR definition, by only taking into account the "universal" ultraviolet part from such quantities, leaving out all finite parts. These prescriptions are however not gauge invariant, see for exam-ple [34]. Moreover the "universal" part is only universal in the usual linear gauge. With the non-linear gauge fixing we implement, see section 2.3, our results would not be gauge invariant and one has to be very careful with the Ward identities. We leave this important issue to a forthcoming paper [35]. Nonetheless to conform with this widespread general usage, we also implement a DR scheme defined from a physical quantity, to be discussed shortly, which reproduces the usual counterterms defined from other quantities in the linear gauge. As known, the others Higgs masses m h 0 (for the lightest CP-even), m H 0 (for the heaviest CP-even) and m H ± (for the charged) receive corrections that can be very important. To be able to stick with an on-shell definition and in order to weigh the effect of the tan β scheme dependence, we also define two other schemes. One is based on the use of A 0 → τ + τ − partial width to which the QED corrections have been extracted, we will refer to this scheme as the A τ τ scheme. For the third one, we take m H 0 as an input parameter and trade it for "tan β" hence loosing one prediction, we will call this scheme MH. This scheme is also used in [29]. With tan β fixed, we can turn to the other sectors. iii) Neutralino and charginos: For the neutralino and chargino sector, we implement an on-shell scheme taking as input parameters the masses of the two charginos (this defines the counterterms to the SU(2) gaugino, winow, mass M 2 and to the higgsino,h, parameter µ) and the mass of the LSP mχ0 1 (which completes the extraction of the U(1) gaugino, binob, mass M 1 ). The other neutralino masses m χ 0 2 , m χ 0 3 and m χ 0 4 receive corrections to their tree-level value. Obtaining finite corrections for the masses and decays is a not trivial test of the procedure. Here our implementation is quite similar to the one adopted in [36] when one sticks to the DR tan β . iv) Sfermions: For the slepton sector we use as input parameters the physical masses of the two charged sleptons which in the case of no-mixing define the R-slepton soft breaking mass,Ml R and the SU(2) mass,Ml L , giving a correction to the sneutrino mass at oneloop. In the case of mixing one needs to fix the counterterm to the tri-linear coupling. The best option would have been to define this from a decay such asτ 2 →τ 1 Z. In the present letter we take a much simpler prescription, we solve the system by taking as input all three slepton masses. For the squark sector, for each generation three physical masses are needed as input to constrain the breaking parametersMQ L for the SU(2) part,Mũ R , Md R for the R-part. In case of mixing, the simplest prescription for the counterterms to the tri-linear couplings A u , A d derives from two conditions on the renormalised mixed two-point functions as is done in [37]. Our plan is to replace these conditions by an onshell input such as a decay of the heavy squark to the lighter one and a Z, to conform with a fully on-shell scheme and study further the scheme dependence. Wave function renormalisation is introduced so that the residue at the pole of all physical particles is 1 and no-mixing is left between the different particles when on shell. This applies for all sectors. Dimensional reduction is used as implemented in the FFL bundle at one-loop through the equivalent constrained dimensional renormalisation [38]. Renormalisation of the strong coupling constant and the gluino is not an issue for the examples we study here.
We have verified our codes and schemes with different calculations on the market for a variety of correction to masses and other observables [28,34,36,37]. The code has also been used for corrections to SM processes and also to one-loop induced processes χ 0 1χ 0 1 → γγ, Zγ, gg [11] relevant for indirect detection.

Non-Linear gauge-fixing
We use a generalised non-linear gauge [11,39] adapted to the minimal supersymmetric model. The gauge fixing writes Unlike the other parts of the model, L GF is written in terms of renormalised fields and parameters. G 3 , G ± are the Goldstone fields. We always work with ξ γ,Z,W = 1 so as to deal with the minimal set of loop tensor integrals. More details will be given elsewhere [35].

The different parts of the cross section
The one-loop amplitudes consist of the virtual corrections A EW +QCD 1loop and the counterterm contributions can contain infrared divergences due to photon and gluon virtual exchange. These are regulated by a small photon or gluon mass. For the QCD corrections we study here, this implementation does not pose a problem. The photon and gluon mass regulator contribution contained in the virtual correction should cancel exactly against the one present in the photon and gluon final state radiation. The photonic (gluonic) contribution is in fact split into a soft part, where the photon (gluon) energy is less than some small cut-off k c , A sof t γ,g (E γ,g < k c ) and a hard part with A hard γ,g (E γ,g > k c ). The former requires a photon/gluon mass regulator. We use the usual universal factorised form with a simple rescaling for the case of the gluon correction. We take α s = α s (M 2 Z ) = 0.118.

Checks on the calculation
i) For each process and for each set of parameters, we first check the ultraviolet finiteness of the results. This test applies to the whole set of the virtual one-loop diagrams. The ultraviolet finiteness test is performed by varying the ultraviolet parameter C U V = 1/ε. We vary C U V by seven orders of magnitude with no change in the result. We content ourselves with double precision.
ii) The test on the infrared finiteness is performed by including both the loop and the soft bremsstrahlung contributions and checking that there is no dependence on the fictitious photon mass λ γ or gluon mass λ g .
iii) Gauge parameter independence of the results is essential. It is performed through the set of the seven gauge fixing parameters defined in Eq. (1). The use of the seven parameters is not redundant as often these parameters check complementary sets of diagrams. It is important to note that in order to successfully achieve this test one should not include any width in the propagators. In fact our tree-level results do not include any width. Because of the parameters and the energies we consider, no width is required to regularise the cross sections.
iv) For the bremsstrahlung part we use VEGAS adaptive Monte Carlo integration package provided in the FFL bundle and verify the result of the cross section against CompHep [31]. We choose k c small enough and check the stability and independence of the result with respect to k c .

Boltzmann equation, the small v expansion
Having the collection of cross sections and the masses of the annihilating (and co-annihilating) DM particles we could have passed them to micrOMEGAs for a very precise determination of the relic density based on a careful treatment of the Boltzmann equation. However, to weigh the impact of the corrections on the relic density it is worth to gain insight through an approximation in going from the cross sections to the relic density, especially that we have found these approximations to be, after all, rather excellent for the cases we study, including co-annihilations. Moreover corrections to the cross sections could be incorporated in the case of non-thermal production. All cross sections σ ij where i, j label the annihilating and co-annihilating DM particles i, j, are expanded in terms of the relative velocity v ij , which for neutralino annihilation is v = 2β = 2 1 − 4m 2 χ 0 1 /s. Away from poles and thresholds, it is a very good approximation to write keeping only the s-wave, a ij , and p-wave, b ij coefficients. With T being the temperature, with g 1 = 2 the neutralino spin degree of freedom (sdof), the co-annihilating particle of sdof g i and mass m i contributes an effective relative weight of The total number of sdof isg ef f = ig i,ef f . A good approximation for the relic density is obtained by carrying a simple one dimensional integration a ij , b ij that are needed to compute σ ij in Eq. (4) are given in cm 3 .s −1 .
x F represents the freeze-out temperature. g * (x F ) is the effective degrees of freedom at freeze-out. g * is tabulated in micrOMEGAs and we read it from there. For the examples we will study g * (x F ) ∼ 9.29. In the freeze-out approximation, x F can be solved iteratively from where the neutralino mass is expressed in GeV. The numerical solutions of the density equation and hence the freeze-out suggest [6,40] c = 1.5 is a very good choice in most, but not all, cases. Though we have verified that Eq. (5) converges quickly and agrees well with the result of micrOMEGAs, the results we give use x F as extracted from micrOMEGAs. The loop corrected cross sections should also impact on the value of x F which is not exactly the same as the value extracted from the tree-level cross sections. However, the shift is marginal, though ultimately in a full computation this should be taken into account. Our results therefore use the same value of x F at both tree and one-loop level. On the other hand to derive the relic density we rely, in this letter, on Eq. (4). For The case ofχ 0 1χ 0 1 annihilation, the latter simplifies to The weight of a channel (see the percentages we will refer to later) corresponds to its relative contribution to J.

Choosing points in the MSSM parameter space
Current limits on the relic density, from WMAP and SDSS [3] give the 2σ range 0.092 < Ωh 2 < 0.124.
In this first exploratory study we thought it is best to consider different scenarii without worrying too much about the absolute value of the derived relic density in order to grasp the origin of the large corrections, if any. Our choice of scenarii was motivated by the physics issues, although our choice is biased towards the popular scenarii that emerge in mSUGRA within thermal production. Nonetheless our choice covers annihilation into fermions, gauge bosons and co-annihilation. This said, apart from the annihilation into gauge bosons, the derived relic density is either within this range or not overly outside. For the gauge bosons the motivation was to take a model that singles out the W W and ZZ final states channels. Moreover since the impact of the radiative corrections can be large there is not much sense in picking up a model based on its agreement with the current value of the relic density on the basis of a tree-level calculation. This said we have used micrOMEGAs as a guide, being careful about translations of effective couplings and input parameters. micrOMEGAs was also quite useful in justifying the approximations we use for deriving the relic density from the cross sections. We should also add that in this letter we do not apply the radiative corrections to all the subprocesses that can contribute to the relic density but only to those channels and subprocesses that contribute more than 5% to the relic density. When calculating the correction to the relic density we include these channels at tree-level.

Annihilation of a bino LSP into charged leptons
The first example we take corresponds to the so-called bulk region, with a neutralino LSP which is mostly bino. The latter will couple with the particles of largest hypercharge, the R-sleptons. Therefore annihilation is into charged leptons. Because of the Majorana nature of the LSP, there is no s-wave contribution in the case of massless fermions. In -17% -17% -15%   Table 1. The e + e − , µ + µ − channels contribute for 31% each. The difference between the three channels is accounted for by the contribution of the s-wave of the τ final states and very little from the fact thatτ 1 is slightly lighter than the lightest smuon and selectron.
Let us first comment on the p-wave contribution which gives the bulk of the contribution to the relic density. The total correction is about 18% in this case. It is tempting to parameterize the corrections. In fact, had we used the value of the gauge coupling not at low energy but at the scale of theχ 0 1 mass of order M Z the bulk of the correction would be absorbed. Indeed, (137.04/127.77) 2 − 1 ∼ 15%. In the few other examples we have looked at concerning the annihilation into leptons, we arrive at the same order of correction, see for example the corrections to the p-wave in the case of co-annihilation and even where there is some higgsino component as in the case of annihilation into bb. The other common trend is that the correction does not show any dependence on the tan β scheme we choose when there is a large bino component. This is not the case for the s-wave. In particular the MH scheme differs from the DR and A τ τ . The tan β dependence comes essentially from the Yukawa contribution, see Eq. (4.8) of ref [1]. The latter is also sensitive to the higgsino component of the neutralino that is also affected by the tan β change. The effect is more obvious in the case of scenario iii), see Table 4. Note however that even in the case of massless fermions, there is (though small) a contribution to the s-wave due to hard photon radiation. Hard photon radiation in association with a light charged fermion pair is not subject, for the s-wave amplitude, to the known helicity suppression when no photon is emitted [12]. Taking both the s and p-wave contribution leads to a correction on the relic density of about 17%. As discussed a few lines above, using α ef f (M Z ) reduces the correction to the level of a few percent.

Neutralinoτ co-annihilation
In this scenario the LSP is still the lightest neutralino and we take it to be essentially bino though with a small higgsino component with a compositionχ 0 1 = 0.986b − 0.049w + 0.144h 1 − 0.070h 2 and mass mχ0 1 = 162.34 GeV. We consider a scenario where the NLSP is the lightestτ 1 with mass mτ 1 = 168.42 GeV coming mainly from itsτ R component. The lightest smuon and selectron are given masses so that they are thermodynamically irrelevant, they have a mass of about 195 GeV. Apart from M A = 1TeV, for the squark and Higgs sector, as well as tan β , the parameters are the same as in the example in the bulk region. We want therefore to concentrate on co-annihilation involving onlyτ 1 . With a mass difference between the LSP and NLSP of 6.08 GeV which is about only 4% of the mass of the LSP,χ 0 1τ 1 annihilation is quite efficient with the two channels χ 0 1τ ± 1 → τ ± γ, τ ± Z accounting for as much as half of the total contribution, see Table 2. τ ± 1τ ± 1 takes up a quarter of the total. Neutralino annihilation makes up for about 15%. The rest is made up byτ − 1τ + 1 → γγ, γZ. It is interesting to note thatτ ± 1τ ± 1 ≫τ −τ + . Our approximation based on Eq. (4) gives Ωh 2 = 0.128, with x F = 26.5. -9% -9% -9% Table 2: Neutralinoτ 1 co-annihilation: meaning of the different cells is as in Table 1. Percentages in the first column represent the weight of the corresponding channels.
•χ 0 1χ 0 1 → f f Compared to the bino-case in the bulk region where this channel accounted for the totality of the relic density, here it only makes up for 6%. Nonetheless the effect of radiative corrections on this channel are very similar to what we found in the scenario of section 3. One may be misled by interpreting the 200% relative correction as a large correction to the relic density. This is a relatively large correction to the s-wave contribution, but in absolute terms this correction is totally negligible compared to the correction brought about by the p-wave contribution at the cross section level as well as after taking the thermal average, notwithstanding that the whole channel in the co-annihilation region has little weight which further dilutes the effect of such large relative correction. As pointed out in the previous section this 200% relative correction is due to the smallness of the s-waveχ 0 1χ 0 1 → f f which is offset by the hard photon emission that now allows for an s-wave contribution [12]. As discussed in the first reference in [12], the relative importance of hard radiation reduces fast as the mass of the intermediate slepton increases. This explains why the relative effect is more prominent in the co-annihilation region.
Before getting into the details about each of the main contributing channels that involveτ ± 1 co-annihilation, let us point at some of their common features. The bulk of the contribution comes now from the s-wave, especially after taking into account thermal averaging. Another common feature is that the tan β scheme dependence is hardly noticeable. Moreover, the corrections to the sand p-wave are, within a margin of 4%, the same.
the correction to the p-wave is within only a 1% from the correction to the s-wave. Here the electroweak correction amounts to about 9%. This is about half the correction we find in all other channels. The reason is the following, the effective coupling for the emitted photon should still be taken at q 2 = 0, and therefore effectively since there is only one neutralino coupling, this should be taken at the scale of mχ0 1 . Proper use of the effective couplings here absorbs about 7% of the correction, leaving therefore a 2% correction. For the τ Z final state the use of the effective coupling would leave out a 6% correction to the s-wave after absorbing the correction due to the effective coupling.
•τ ± 1τ ± 1 → τ ± τ ± The radiative correction toτ ± 1τ ± 1 → τ ± τ ± reveals a quite interesting feature as can be seen from Fig. 1 which shows the dependence of the cross section as a function of the relative velocity (or rather the square of it). It is from this velocity dependence that we usually extract the values of the coefficients of the sand p-wave contribution, at both the tree-level and one-loop, that we need to calculate the relic density. The figure extends to v 2 values as large as 0.5 while we could have contented ourselves with a maximal value of v 2 = 0.3, considering the typical value that one obtains for the freeze-out temperature v 2 ∼ 6/x F ∼ 0.2. Even so, the figure shows that in this case a fit to the tree-level cross section in the form a + bv 2 works quite well. For the one-loop correction a polynomial fit does not do for low enough velocities. There is a large negative correction for v → 0. This correction is in fact very easy to understand. It is the perturbative one-loop manifestation of the non relativistic Coulomb-Sommerfeld effect [41]. With the tree-level cross section denoted as σ 0 and σ 0 v = a 0 + b 0 v 2 , the one-loop perturbative cross section for the samesign stau annihilationτ ± 1τ ± 1 , σ 1−loop Coul. , is such that We thus expect for the one-loop cross section, σ 1  One may ask about how to deal with the 1/v singularity. In fact when calculating the relic density, the one-loop 1/v singularity at the level of the cross section is tamed after thermal average, ∝ ∞ 0 (dv v 2 e −v 2 /4x ) (σv), see also Eq. 2. At the end its contribution to the relic density compared to a 1 is approximately In words, non zero temperature of the problem provides a cut-off. One can also ask whether the one-loop result from the Coulomb-Sommerfeld effect is sufficient. As seen, the QED correction is of O(πα/v) ∼ 0.17 with v typical of the freeze-out temperature. Therefore in our case a one-loop treatment seems to be sufficient especially that theτ ± 1τ ± 1 → τ ± τ ± is not the most dominant channel. This said, these non-relativistic QED 1/v threshold corrections can be resummed to all orders. This resummation as originally performed by Sommerfeld [41] has been known for quite a long time in quantum mechanics, see [42] for a textbook treatment, and amounts to solving the Schrödinger equation in the Coulomb potential. With X nr = −2πα/v for the same-sign stau annihilation, the s-wave factor resums to One might question the validity of Eq. 10 in our case whereτ 1 is not stable. Finite decay width can of course act as a cut-off for the 1/v corrections, see the case of W pair production [43] or slepton pair production at threshold [44]. In our case width effects are of no importance since the characteristic time of the Coulomb interaction, 1/mτ 1 v 2 typical of velocities at freeze-out is much smaller than the decay time, 1/Γτ 1 , since in our example Γτ 1 /mτ 1 ≃ 2 10 −5 and v 2 ∼ 0.2. For smaller δm = mτ 1 − mχ0 1 , the width effect is even more negligible whereas for larger δm, theτ 1τ1 channel would be thermodynamically irrelevant. Therefore in our particular case the resummation can be taken from the old Sommerfeld result. Nonetheless, especially after thermal averaging, in our case this type of QED correction is well approximated by the one-loop approximation 2 .
Taking now all the effects and contributions in our specific example we find an overall correction of −8.6% to the relic density with a corrected value of Ωh 2 = 0.117. As we can see this value would not have been approached with a naive overall rescaling of the effective couplings. Nonetheless, apart from some 4% correction, most of the effect seems to be explained in terms of a proper usage of effective couplings and the Coulomb effect. In fact in the total contribution, the Coulomb effect is diluted and changes the results for the relic density by about only 1.5%.

Annihilation of a mixed gaugino-higgsino LSP into vector bosons
Having studied annihilation into fermions, annihilation into the weak vector bosons is quite interesting. In the context of mSUGRA this occurs for example in the so-called focus point region. In order not to mix issues, we do not consider in this letter a scenario where the LSP neutralino is either dominantly higgsino or wino, therefore avoiding that χ 0 1χ + 1 annihilation is of relevance in this case. We seek a scenario with a neutralino where the largest component is still bino but where one has a substantial higgsino and wino component. In our example one hasχ 0 1 = 0.819b − 0.231w − 0.470h 1 − 0.232h 2 with mχ0 1 = 102.89, mχ+ 1 = 125.13, mχ+ 2 = 215.27 GeV. All other masses (outside the charginoneutralino sector) are taken to be very heavy at 1TeV. This is also to avoid contamination from annihilation into fermions. It would however be worth to study the impact of the sfermions on the radiative corrections. In this case we have not made the extra effort of searching for a set with a relic density within the WMAP range. Here annihilation into W W and ZZ accounts for 80%, see Table 3 with a few other channels below 2% each. These involveχ 0 1χ + 1 into light quarks or W Z which we take into account only at tree-level. First of all, we see that the corrections which affect annihilation into ZZ and W W are about the same (within 5% in the 3 tan β schemes). Moreover the correction to the 2 In a situation, which is not the case here, where theτ ±τ ± andτ ±τ ∓ would contribute equally at treelevel, the Coulomb-Sommerfeld correction would cancel after adding the two-channels. The correction iñ τ ±τ ∓ would be attractive and given by changing the sign of X nr in Eq. 10.The non-perturbative effects of the Coulomb-Sommerfeld-like corrections that might occur when coloured states are involved [9,45] need more care because of the strong QCD coupling, bound states effects might be relevant. The effect of the latter is even more important with models with TeV and multi TeV dark matter almost degenerate with a charged component [16,46]. The non perturbative effects on indirect detection in these models is even more dramatic [15,46].
0.787 -30% -6% +39% Ωh 2 0.053 0.068 0.054 0.039 δΩh 2 Ωh 2 +28% +2% -26% s-wave and p-wave are of the same order, see also Fig. 2 for the v dependence of the ZZ cross section and the extraction of the sand p-wave coefficients. However this is not the most important conclusion. The most important lesson is that there is a very large tan β scheme dependence. In some other investigation concerning Higgs decays, we had noticed, as was also pointed in [34] with a similar scheme, that the MH scheme can lead to very large corrections. However in many instances, like what we saw in the case of the bino annihilation or co-annihilation, the δ tan β is screened. Unfortunately in this mixed neutralino scenario, the δ tan β dependence can be potentially enhanced by 1/(µ 2 −M 2 2 ) from the renormalisation of δµ for example. This needs further investigation. In this model, already at tree-level, we had noticed that the cross sections were very sensitive to small changes in the underlying parameters. Apart from the tan β scheme dependence, the corrections in the DR scheme look modest especially for the dominant swave contribution. However one should not forget that these small corrections are within the use of α in the Thomson limit. Switching to a scale of order M Z , the corrections are large of order 20%. Therefore our conclusion is that in such a scenario there are genuine large corrections in all three schemes we have considered. This also confirms the study of the chargino neutralino system at one-loop in [47] which showed that though the corrections to the masses are modest, there can be a large (of order 30%) change in the gaugino-higgsino component and hence a large impact on cross sections.

Annihilation of neutralinos to bb and a not too heavy pseudo-scalar Higgs
We expose this last example for illustrative purposes. Indeed the one-loop perturbative treatment of the Higgs coupling to a bottom pair using the bottom pole mass, here we have taken all along m b = 4.7 GeV, is far from describing the bulk of the radiative QCD corrections which as we know need to be resummed both for the running effective m b purely from QCD and from the so-called ∆m b effects. The latter being more important for high tan β . These effects are already taken in micrOMEGAs [6] for example. The purpose here therefore is to see whether there are other possible effects, though smaller, that are captured in a complete one-loop calculation. Ideally one would like to subtract the known universal one-loop QCD corrections from the full one-loop QCD corrections that can occur for example from box diagrams, for v = 0, outside the Higgs resonance. By the way because of the Majorana nature of the LSP, at the smallest velocities the most important Higgs resonance is the pseudo-scalar Higgs, A. In any case for a precise calculation of the relic density, non resonant contributions should be taken into account as thermal average is to applied and would bring some smearing. Therefore, here we concentrate onχ 0 1χ 0 1 → bb where the scalar resonance is not negligible. Again like we argued for the previous example we have to rely here also on some higgsino component. The composition of the LSP isχ 0 The mass of the gluino is also lowered to be 400 GeV.
At tree-level with m b = 4.7 GeV, the dominant modes are annihilation into W W followed by bb (about 10% smaller). τ τ and ZZ are about the same level but a factor 10 smaller than the W W channel. We show, see Table 4, the W W and ZZ channel in order to make a comparison with the previous case of a mixed bino with a substantial higgsino component. Here the tan β scheme dependence has considerably reduced especially between the A τ τ 0.033 +3% +9% +52% b 0.631 +19% +18% +12%χ 0 1χ 0 1 → bb A τ τ DR MH δa/a EW -1% +3% +31% δa/a QCD -26% -26% -26% δb/b EW -1% +3% +29% δb/b QCD -30% -30% -30% Table 4: Mixed case 2: As in Table 3 for the array on the left. The array on the right gives the relative corrections to the bb channel for the QCD and EW corrections. and DR scheme. Note also that in the case of τ τ channel there is also a discrepancy with the MH scheme and the other two schemes compared to the almost pure bino case. This again is due to the larger contribution of the higgsino tan β dependent part, naturally in the A-exchange not present in the pure bino case but also in theτ exchange.
Let us look at what we obtain for the bb channel. The electroweak corrections do show some tan β scheme dependence. Compared to the electroweak corrections in the DR scheme of tan β , the QCD corrections are larger by an order of magnitude, they amount to about −30% in both the sand p-wave. If one assumes that most of these corrections arise from the A → bb vertex, then we know that there are large logs resulting from the anomalous dimension of the pseudoscalar (and for that matter the scalar). Its one loop part can be found in [48] and amounts to − 4α s π ln 2mχ0 There is also an important SUSY QCD correction termed ∆m b [49], see section B.4.2 in the second paper of [6]. In this scenario, adding these two corrections amounts to about −35%. Subtracting these from the correction we calculate for the s-wave, leaves us with about +9% QCD corrections. The QCD corrections from the anomalous dimension and the ∆m b effect were extracted for the known effect to A → bb. However, since at v = 0, for the s-wave cross section, the neutralino system constructs a pseudoscalar because of its Majorana nature, the same corrections should affect the a coefficient even if the contribution from the pseudo-scalar Higgs is negligible.

Conclusions
We have performed the first electroweak corrections to some important processes relevant for the relic density of neutralinos in supersymmetry. This has become possible thanks to an automated code for the calculation of loop corrections in the MSSM that will allow to perform with the same tools and the same framework (scheme dependence,..) analyses at one-loop at the collider and for dark matter. Our findings suggest that in the case of a dominantly bino neutralino, a large part of the correction can be accounted for through an effective electromagnetic coupling at the scale of the neutralino mass. Even so, complete one loop corrections would be needed to match the foreseen precision of PLANCK. The corrections to the relic density are not sensitive much to the way tan β is renormalised. In the case of co-annihilation of a bino and stau, the conclusion is similar but one has to be wary of possible Coulomb-Sommerfeld corrections. For a neutralino LSP which is strongly mixed, the corrections are large and the tan β scheme dependence not negligible at all. More investigation of such scenarii should be conducted. Some QCD (and SUSY QCD) corrections affecting final states quarks in the case of neutralino annihilation need that one goes beyond one-loop. Some of these corrections have been identified and already implemented in a code such as micrOMEGAs. Apart from these corrections, there remain some additional one-loop corrections that should be taken into account. Before generalising these conclusions, more work is needed. However, the tools exist. The next step is to interface our code for the loop calculations with a dedicated relic density calculator, avoiding double counting of some of the one-loop corrections implemented as effective operators in the relic density calculator. Work in this direction has already begun based on micrOMEGAs as the relic density calculator.