SUSY-QCD effects on neutralino dark matter annihilation beyond scalar or gaugino mass unification

# SUSY-QCD effects on neutralino dark matter annihilation beyond scalar or gaugino mass unification

## Abstract

We describe in detail our calculation of the full supersymmetric (SUSY) QCD corrections to neutralino annihilation into heavy quarks and extend our numerical analysis of the resulting dark matter relic density to scenarios without scalar or gaugino mass unification. In these scenarios, the final state is often composed of top quarks and the annihilation proceeds through -boson or scalar top-quark exchanges. The impact of the corrections is again shown to be sizable, so that they must be taken into account systematically in global analyses of the supersymmetry parameter space.

###### pacs:
12.38.Bx,12.60.Jv,95.30.Cq,95.35.+d

eprint LPSC 09-082

## I Introduction

The search for physics beyond the Standard Model (SM) is no longer restricted to colliders only. In fact, the most compelling evidence for new physics comes today from cosmological observations such as the mission of the Wilkinson Microwave Anisotropy Probe (WMAP), which have determined the matter and energy decomposition of our universe with unprecedented precision. These observations indicate the existence of Cold Dark Matter (CDM) in the universe, which cannot be accounted for by the SM and likely consists of Weakly Interacting Massive Particles (WIMPs) with non-relativistic velocities. A combination of the five-year measurement of the cosmic microwave background by the WMAP mission with supernova and baryonic acoustic oscillation data yields the narrow -interval for the relic density of dark matter WMAP ()

 0.1097<ΩCDMh2<0.1165, (1)

where denotes the present Hubble expansion rate in units of 100 km s Mpc .

The measured relic density of dark matter can be used to constrain extensions of the Standard Model, which provide a viable WIMP candidate. In the Minimal Supersymmetric Standard Model (MSSM) with -parity conservation, this could be the Lightest Supersymmetric Particle (LSP) , if it is neutral and a color singlet. One can then calculate its relic density, compare it with the experimental limits in Eq. (1), and identify the favored regions of the MSSM parameter space. The relic density

 Ω~χh2 = n0m~χρc (2)

is proportional to the present number density and the mass of the LSP. is the critical density of our universe, and is the gravitational constant. The present number density is obtained by solving the Boltzmann equation describing the time evolution of the number density

 dn~χdt = −3Hn~χ−⟨σannv⟩(n2~χ−n2eq). (3)

The first term on the right-hand side corresponds to a dilution due to the expansion of the universe, and the second term corresponds to a decrease due to annihilations and co-annihilations of the relic particle into SM particles Gondolo1991 (). Here, denotes the time-dependent Hubble expansion parameter, and the density of the relic particle in thermal equilibrium. Details of the dark matter interactions enter the Boltzmann equation through the thermally averaged cross section . The cross section takes into account the thermal velocity distribution of the relic particle and is calculated for a given temperature by

 ⟨σannv⟩ = 4m4~χTK22(m~χT)∫dpcmp3cm√sK1(√sT)σann(s), (4)

where and are the modified Bessel functions of the first and second kind, respectively. The center-of-momentum energy is related to the particle mass and the relative momentum of the annihilating pair through Gondolo1991 ().

In order to keep up with current and future experimental improvements, one has to understand and reduce the different uncertainties involved in the analysis, both for the prediction of the dark matter relic density and the extraction of new mass parameters from cosmological data. These uncertainties include, e.g., a modification of the Hubble expansion rate due to quintessence or an effective dark energy density Arbey2009 (), differences in the new physical particle masses obtained with different spectrum codes Belanger2005 (), or a lack of precision in the annihilation cross-section of dark matter particles Baro2007 (). In this paper, we focus on the impact of next-to-leading order QCD and SUSY-QCD corrections on the latter, but other possible uncertainties will also be briefly discussed.

In many scenarios of the MSSM, the lightest neutralino is the LSP and therefore a suitable dark matter candidate. The thermally averaged cross section is then obtained by computing all relevant neutralino annihilation cross sections into SM particles. Most prominent are the processes with two-particle final states such as a fermion-antifermion pair or a combination of gauge () and Higgs bosons () Jungman (); DreesNojiri (). In this paper, we focus on the annihilation into a massive quark-antiquark pair, since the leading order cross section with a fermion-antifermion final state is proportional to the mass of the fermion, which disfavors the light quarks. Moreover, the annihilation into heavy quarks is important in the regions of parameter space allowed by Eq. (1). We now present the full details of our calculation and investigate scenarios with dominant top-quark final state contributions, extending our analysis of Refs. Letter1 (); Letter2 () beyond minimal supergravity (mSUGRA) models. In mSUGRA models, one is constrained by having only five universal high-scale parameters [, , , , and sgn()], and in regions of parameter space, where quark final states are important, the cross section is dominated by an exchange of Higgs bosons. Here, we relax the unification of either the scalar masses or the gaugino masses, one at a time. This allows for scenarios different from mSUGRA, where the annihilation cross section is not necessarily dominated by Higgs-boson exchanges. Apart from a Higgs-boson dominated scenario, we thus analyze scenarios, where -boson exchanges in the -channel or squark exchanges in the - and -channels play an important role. The full QCD and SUSY-QCD corrections in these scenarios turn out to be significant, and we have therefore included them into the public code micrOMEGAs micromegas ().

This paper is organized as follows: In Sec. II and appendices A and B, we give all necessary details of the calculation related to the virtual loop corrections, the renormalization procedure, and the calculation of real gluon emission, in particular the subtraction of the induced soft and collinear singularities. We then continue in Sec. III with a discussion of the MSSM models beyond scalar or gaugino mass unification. In Sec. IV, we analyze the impact of the radiative corrections on the relic density in these models. Finally, our results are summarized in V.

## Ii Calculation details

The annihilation of neutralinos into quarks proceeds at tree-level through an exchange of a -boson and Higgs-bosons in the -channel and through the exchanges of scalar quarks in the - and -channels (see Fig. 1). By definition, a WIMP can only have electroweak interactions, and in order to reach a sufficient annihilation rate to explain the dark matter relic density, the cross section has to be enhanced, e.g., by a resonance. One often finds that there is one contribution which dominates the whole cross section. This allows us, by choosing different scenarios, to study contributions from various channels. Moreover, it allows us also to isolate effects of radiative corrections coming from different sources.

We have computed the full QCD and SUSY-QCD corrections to neutralino annihilation into quarks. The next-to-leading order cross section contains virtual contributions stemming from loop diagrams and real contributions, which are due to the radiation of an additional gluon. Symbolically, the next-to-leading order (NLO) cross section can be written as

 σNLO=∫2→2dσV+∫2→3dσR, (5)

where denote the virtual and real emission parts integrated over the two- and three-particle phase space, respectively. The cross section , which includes the tree-level and the one-loop virtual corrections, is given by

 σVann(s)=∫116π212κ(s,m2~χ,m2~χ)κ(s,m2q,m2q)2s[|Mtree|2+2R(M1−loopM∗tree)]dΩ, (6)

where . For our notation and conventions, we refer the reader to App. A. There are two types of contributions to the one-loop amplitude , those coming from the loop diagrams and those coming from the counterterms. The loop diagrams, depicted in Fig. 2, contain ultraviolet (UV) and infrared (IR) divergent loop integrals. We regulate both types of divergences dimensionally () and evaluate the loop integrals in the dimensional reduction scheme (). The full analytic results for the loop diagrams are given in App. B.

The UV divergences are compensated by counterterms (see Fig. 3) related to quarks and their scalar super-partners, the squarks. All counterterm vertices with quark or squark legs contain wave-function renormalization factors , and , that follow from replacing the (s)quark fields by

 (qLqR)→(1+12δZL001+12δZR)(qLqR),~qi→(δij+12δZij)~qj. (7)

Although in principle only the wave-function renormalization constants of the external quarks have to be included, we also include the squark renormalization constants, as they allow us to perform simpler UV-convergence checks. In the full calculation, the squark wave-function renormalization constants cancel out. The wave-function renormalization constants are determined by requiring the residues of the propagators to remain at unity even at one-loop order. This condition gives

 δZL = R[−ΠL(m2q)−m2q(˙ΠL(m2q)+˙ΠR(m2q))+12mq(ΠSL(m2q)−ΠSR(m2q)) (8) −mq(˙ΠSL(m2q)+˙ΠSR(m2q))], δZR = δZL(L↔R), (9) δZii = −R[˙Π~qii(m2~qi)],δZij = 2R[Π~qij(m2~qj)]m2~qi−m2~qjfori≠j, (10)

where and stand for the vector and the scalar parts of the two-point Green’s function as defined in Ref. Kovarik2005 () and . After performing the wave-function renormalization, the remaining divergences are canceled by renormalizing the coupling constants. In our case, the coupling counterterms that receive contributions proportional to the strong coupling constant are the quark Yukawa couplings through the masses of the quarks, the squarks masses and the squark mixing angle . A very important contribution, in particular in scenarios with a dominant Higgs-boson exchange, comes from renormalizing the Yukawa couplings of the quarks Letter2 (). In our calculation, we use the Yukawa couplings for both the top and the bottom quarks. As the quark masses that serve as inputs are defined in different schemes, we take two distinct approaches for top and bottom quarks. For top quarks, the input is the on-shell mass measured at the Tevatron CDFD02008 (), and the -mass of the top quark (and hence the Yukawa coupling) is obtained by subtracting the finite on-shell counterterm

 δmq = 12R[mq(ΠL(mq)+ΠR(mq))+ΠSL(mq)+ΠSR(mq)]. (11)

On the other hand, the input mass for bottom quarks is extracted in the renormalization scheme from the Standard Model analysis of sum rules MBmass (). In order to obtain the appropriate bottom Yukawa coupling in the renormalization scheme within the MSSM, we first use the Standard Model next-to-next-to-leading order (NNLO) renormalization group evolution to obtain the mass of the bottom quark at the scale Baer2002 (). Still in the SM, we then convert to Baer2002 (), and finally we apply the threshold corrections including also contributions from SUSY particles in the loop. For the last step, we take into account the fact that the sbottom-gluino and stop-chargino one-loop contributions are considerably enhanced for large or large and can be resummed to all orders in perturbation theory Carena2000 (); Spira2003 (). Denoting the resummable part by and the finite one-loop remainder by , the bottom quark mass is then given by

 m¯¯¯¯¯¯¯DR,MSSMb(Q) = m¯¯¯¯¯¯¯DR,SMb(Q)1+Δb−Δmb. (12)

In the squark sector, we pick five independent quantities, , , , , and , in order to respect the SU(2) symmetry. We renormalize the masses of three squarks in the on-shell scheme, which leads to the counterterm

 δm2~qi=R[Π~qii(m2~qi)]. (13)

The remaining mass of the heavier scalar bottom quark is treated as dependent,

 δm2~b2=1sin2θ~b[δm2~t1cos2θ~t+δm2~t2sin2θ~t−δm2~b1cos2θ~b+(m2~t2−m2~t1)sin2θ~tδθ~t−(m2~b2−m2~b1)sin2θ~bδθ~b]. (14)

The squark mixing angle is renormalized in the scheme and so the corresponding counterterms of the squark mixing matrices contain only the divergent parts. The counterterms can be determined as

 δR~qij=2∑k=114(δZdivik−δZdivki)R~qkj, (15)

using only the divergent parts of the wave-function renormalization constants. This is equivalent to fixing the mixing angle as

 δθ~q = 14(δZdiv12−δZdiv21)=12(m2~q1−m2~q2)R[Π~qdiv12(m2~q2)+Π~qdiv21(m2~q1)]. (16)

After renormalization, the virtual cross section still contains IR divergences, that come from an exchange of a gluon in the loops. In order to compensate them, one has to include the real cross section coming from diagrams with an additional gluon in the final state (see Fig. 4). The symbolic formula of Eq. (5) cannot be directly applied, since the divergences appear as and poles in the number of dimensions in the one-loop amplitude and as a divergent matrix element in certain regions of the parameter space. A convenient way to combine these two cross sections and to cancel the divergences is the dipole subtraction method Catani2002 (). Using this method, we can write the total next-to-leading order cross section as

 σNLO=∫2→2[dσV+∫1dσA]ϵ=0+∫2→3[(dσR)ϵ=0−(dσA)ϵ=0], (17)

where we introduced an auxiliary cross section . The auxiliary cross section does not contribute to the total cross section and serves only as a tool to cancel the IR divergences. It has the same divergence structure as the real cross section and at the same time its structure allows for a partial integration of the gluon phase space, so that it can also be added to the virtual cross section, canceling the divergences in each part. The matrix element leading to the auxiliary cross section for the real part is constructed from two dipole contributions and as

 |M2→3aux|2=D31,2(k1,k2,k3)+D32,1(k1,k2,k3), (18)

where the are the four-momenta of the final state quarks and of the gluon. The dipole contributions and are related by a simple interchange with given by

 D31,2(k1,k2,k3) = CF8παss|Mtree|2× (19) 11−x2⎧⎨⎩2(1−2μ2q)2−x1−x2− ⎷1−4μ2qx22−4μ2qx2−2μ2q1−2μ2q[2+x1−1x2−2μ2q+2μ2q1−x2]⎫⎬⎭,

where , , and Catani2002 (). The leading order matrix element appearing in Eq. (19) is calculated using different kinematics with redefined 4-momenta

 kμ1→˜kμ31=12qμ−√1−4μ2q√x22−4μ2q(kμ2−12x2qμ),kμ2→˜kμ2=12qμ+√1−4μ2q√x22−4μ2q(kμ2−12x2qμ). (20)

The auxiliary matrix element that cancels the infrared divergences of the virtual matrix element is

 |M2→2aux|2 = 2CFαs2π(4π)ϵΓ(1−ϵ)|Mtree|2× (21) [(μ2s12)ϵ(Vq(s12,mq,mq;ϵ)−π23)+1CFΓq(mq;ϵ)+32lnμ2s12+5−π26],

where . The function is composed of a singular part,

 V(S)(s12,mq,mq;ϵ)=1+β22β[1ϵln1−β1+β−12ln21−β1+β−π26+ln1−β1+βln21+β2] (22)

with , and a non-singular part given by

 V(NS)q(s12,mq,mq) = 32ln1+β22+1+β22β[2ln1−β1+βln2(1+β2)(1+β)2+2Li2(1−β1+β)2−2Li2(2β1+β)−π26] (23) +ln(1−12√1−β2)−2ln(1−√1−β2)−1−β21+β2ln√1−β22−√1−β2 −√1−β22−√1−β2+21−β2−√1−β21+β2+π22.

The function is defined as

 Γq(mq;ϵ)=CF[1ϵ+12lnm2qQ2−2], (24)

where is the renormalization scale. Apart from convergence checks of our calculation, we also performed checks of the finite part of the one-loop diagrams against the results obtained by automatic computer packages FeynArts and FormCalc FeynArts () and parts of the calculation also against existing results, e.g. in Ref. PhDs (). The UV and IR divergent parts of our loop integrals were checked and agree with the results given in Ref. Ellis2007 ().

## Iii SUSY models beyond minimal supergravity

In our previous publications, we studied the impact of the SUSY-QCD corrections to neutralino annihilation within minimal supergravity (mSUGRA) models, defined at the grand unification scale in terms of a universal scalar mass , a universal gaugino mass , a common trilinear coupling , the ratio of the Higgs doublet vacuum expectation values, and the sign of the higgsino mass parameter Letter1 (); Letter2 (). In mSUGRA models, the lightest neutralino is often a -ino and annihilates preferably through resonant Higgs-boson exchanges. For example, in the focus-point region at high scalar masses the cross section is dominated by heavy -even Higgs bosons decaying into top quarks. Regions with important bottom quark final states include those with small gaugino masses , where light -even Higgs-boson exchanges dominate, and the -funnel region at high values of , where the bottom Yukawa coupling is enhanced and the annihilation proceeds through pseudo-scalar Higgs-boson exchanges. At large , the bottom-quark contribution can also become sizable in the focus-point region.

In this work, we study in detail numerically the annihilation of neutralinos into top quark-antiquark pairs through the exchange of -bosons in the -channel and of top squarks in the - and -channels (see Fig. 1). These channels can be enhanced in models, where some of the unification conditions have been relaxed, which is well motivated theoretically Anderson:1996bg (); Anderson:1999uia (); Olechowski:1994gm (). To be concrete, we focus on two sets of models based on the SO(10) Grand Unified Theory (GUT): models with non-universal Higgs masses (NUHM), and models without gaugino mass unification. SO(10) theories are particularly promising, as they involve complete 16-dimensional matter multiplets with a right-handed neutrino and can be embedded in string theories involving larger groups like E or SO(32) Ratz2007 (); Nilles2008 ().

### iii.1 SUSY models with non-universal Higgs masses

In SO(10) SUSY GUTs, the matter superfields of one generation belonging to a 16-dimensional representation are completely mass degenerate, if the SUSY-breaking masses are acquired above the SO(10)-breaking scale. Flavor-blind mechanisms can furthermore lead to universal masses of all matter scalars. However, if the Higgs doublets and belong to different, e.g. 10-dimensional, representations, the corresponding SUSY-breaking masses need not be the same, i.e. the Higgs masses need not be universal Baer:2005bu (); Ellis2002 (); Ellis2002b (); Ellis2002c (). In the scalar part of the general MSSM Lagrangian,

the trilinear scalar couplings and the SUSY-breaking scalar masses are still unified to and at the GUT scale, but the SUSY-breaking Higgs mass parameters and are in general different. Models with non-universal Higgs masses (NUHM) are therefore defined by the parameters , , , , , , and . Note that in our NUHM models the gaugino masses remain unified to and electroweak symmetry-breaking (EWSB) is still achieved radiatively, albeit through modified renormalization group equations (RGEs). This leads to a more constrained parameter space in the plane as compared to mSUGRA models and to the fact that the minimum conditions of the tree-level Higgs potential

 sinβ = −2bm2Hu+m2Hd+2μ2, (26) m2Z2 = m2Hd−m2Hutan2βtan2β−1−μ2, and (27) m2A = m2Hu+m2Hd+2μ2 (28)

allow to replace the parameter by , but not to determine the superpotential parameter and the pseudo-scalar Higgs boson mass as functions of and , as was the case in mSUGRA. It is possible to replace the free parameters and by the low-scale parameters and . Although we don’t make use of this fact in our analysis, it can be useful to consider and as free parameters of the model, since they strongly influence the annihilation of neutralinos into quarks and especially the relative weight of each channel. Having as a free parameter means that one can find scenarios, where and the Higgs-boson exchange dominates the cross section even for smaller values of . In mSUGRA, such a scenario was only allowed for values of . We study this case by choosing the benchmark point I given in Tab. 1. Furthermore, the higgsino parameter influences the higgsino component of the neutralino. By making it larger, one can enhance the contributions from the -boson exchange. This can be achieved, as discussed in Ref. Baer:2005bu (), by starting with large positive values of and at the GUT scale. By virtue of the RGE evolution, is driven to small negative values, while remains positive, which results in a small value of as given by Eq. (27). This in turn gives rise to a larger higgsino fraction of the neutralino and enhances the coupling to -bosons. Such a scenario corresponds to our point II in Tab. 1. Finally, one can make use of the NUHM RGEs to reduce the values of and . This can be obtained by choosing a large and negative difference . On top of that, by choosing to be large and negative, one induces a larger splitting in the third-generation sfermion masses, leading to a scenario (point III in Tab. 1) with a small top squark mass and enhanced - and -channel exchanges of top squarks.

### iii.2 SUSY models without gaugino mass unification

In the gaugino sector, the situation is somewhat similar to the one in Sec. III.1. The breaking of the SO(10) symmetry can proceed via a step involving its SU(5) subgroup. The SU(5) later breaks into the Standard Model gauge groups SU(3)SU(2)U(1) Ellis1985 () and also determines the properties of the SUSY-breaking mechanism. As pointed out in Refs. Anderson:1996bg (); Anderson:1999uia (), the breaking of supersymmetry itself is induced by an -term, while the gaugino masses are generated through a chiral superfield , whose auxiliary component acquires a vacuum expectation value. Assuming that the gauginos belong to the adjoint representation 24 of SU(5), the fields and can in principle belong to any of the irreducible representations appearing in the symmetric product

 (24⊗24)sym = 1⊕24⊕75⊕200, (29)

or any linear combination thereof. The relations between the gaugino masses at the unification scale are given by the embedding coefficients of the Standard Model groups in SU(5). Note that these possibilities are all compatible with gauge coupling unification, but only the case where the SUSY-breaking field is taken to be a pure singlet (1) leads to the gaugino mass universality featured by the mSUGRA model.

The soft SUSY-breaking Lagrangian contains mass terms for the -ino, -ino, and gluino,

 Lsoft ⊂ −12(M1~B~B+M2~W~W+M3~g~g+h.c.). (30)

As argued above, the values of , , and at the unification scale can be considered as independent parameters. Here, we adopt a commonly used parametrization and introduce the dimensionless parameters

 x1 = M1M2andx3 = M3M2, (31)

which will be used together with the -ino mass parameter to describe the gaugino sector. The case reproduces the mSUGRA model with the five parameters , , , , and at the unification scale.

Models with non-universal gaugino masses have been shown to favor annihilation processes where neutralino annihilation into quarks is mediated by a -boson exchange and squark exchanges Orloff2002 (); Martin2007 (); Baer2007 (). It has also been shown that the key parameter in this context is the gluino mass , since it influences practically all sectors of the low-energy mass spectrum through the renormalization group evolution. A decrease in induces a decrease of the mass squared , which through electroweak symmetry breaking conditions induces a decrease in the higgsino mass parameter as well. This in turn increases the higgsino fraction of the neutralino and lowers the pseudo-scalar Higgs mass . Moreover, having independent of the other gaugino mass parameters, it is straightforward to obtain lighter squarks, in particular the scalar tops. This effect can still be enhanced by decreasing the scalar mass parameter and by adjusting the trilinear coupling , which influences the squark mass splitting. With scalar tops being light (even becoming eventually the next-to-lightest SUSY particle), the squark exchange dominates the cross section for a -ino-like neutralino, where a low value of suppresses the Higgs-boson exchange (our point IV in Tab. 2). As the higgsino fraction of the lightest neutralino increases, the squark exchange is enhanced due to the large Yukawa couplings. In the case of a large higgsino fraction, also the -boson exchange becomes important and can take over if the lightest stop is not too close in mass to the lightest neutralino. We analyze the consequences of such a scenario by choosing the point V in Tab. 2.

## Iv Numerical Results and Discussion

Starting from the high-scale parameters, we use the public computer program SPheno 2.2.3 Spheno () for the numerical evaluation of the renormalization group running in order to obtain the SUSY-breaking parameters at the electroweak scale. The relic density of the neutralino is then evaluated numerically using the public code micrOMEGAs 2.1 micromegas (), where we have included our calculation of the neutralino annihilation cross section into third-generation quarks as discussed in Sec. II. For the Standard Model input parameters, such as masses and couplings, we refer the reader to Ref. PDG2008 (), except for the value of the top quark pole mass,  GeV, which has been taken from Ref. CDFD02008 ().

We have chosen five typical parameter points shown in Tabs. 1 and 2, that have a dominant neutralino annihilation into top or bottom quarks through Higgs-boson, -boson or squark exchanges and whose relic density lies reasonably close to the WMAP range of Eq. (1). In all scenarios we only consider relatively low values of . In Sec. IV.1 we analyze the first three parameter points (Tab. 1) where we make use of the possible non-universality of the Higgs-boson masses to construct scenarios with cross sections dominated by Higgs-boson, -boson, and squark exchanges, respectively. In Sec. IV.2, we investigate the point IV, which corresponds to one of the scenarios without gaugino mass unification discussed in Ref. Orloff2002 (), and the point V, which is motivated by one of the “compressed SUSY” scenarios proposed in Ref. Martin2007 (). Our chosen parameter points also satisfy electroweak precision and low-energy constraints such as the measurements of the -parameter, the anomalous magnetic moment of the muon, and the branching ratio of the decay . Due to the large experimental error on the measurement of PDG2008 (), however, only regions featuring very high SUSY masses are excluded at the -level, so that this constraint does not affect our analysis. We have taken into account the anomalous magnetic moment of the muon by taking the higgsino mass parameter . Negative values are disfavored, since they increase the gap between the recent experimental value and the Standard Model prediction Moroi1995 (). The most stringent constraint is given by the inclusive branching ratio of the decay . Recent experimental measurements from BaBar, Belle, and CLEO lead to the combined value BSGamma ()

 BR(b→sγ) = (3.52±0.25)⋅10−4. (32)

The theoretical prediction of the SUSY contribution is particularly sensitive to the masses of the chargino, the charged Higgs boson, and the lightest scalar quark. We have verified that our parameter points lead to values that lie within of the above limit using the public codes FeynHiggs 2.6.5 FeynHiggs () and SusyBSG 1.3 SusyBSG (). Note that also the direct mass limits from collider searches are fulfilled PDG2008 ().

Although the focus of the present paper is the improvement of the annihilation cross section through SUSY-QCD corrections, one should keep in mind that there are also other sources of uncertainties in the calculation of the relic density. From the particle physics side, it is well known that differences in the low-energy mass spectrum may occur when using different spectrum generators. This, in consequence, may induce a sizable difference in the prediction of the dark matter relic density or other observables Belanger2005 (). Note that, in this context, also the numerical value of the Standard Model parameters, in particular the top-quark mass, can influence the favored regions of parameter space, in particular in the case of dominant annihilation into top-quark final states. Corrections to the annihilation cross section are also induced by the electroweak interaction Baro2007 (), but these are generally smaller than the strong corrections and beyond the scope of this work. Concerning cosmology, it has been shown that modifications of the standard cosmological model may also affect the prediction of the dark matter relic density. Including, e.g., an energy content such as quintessence or an effective dark density due to extra dimensions modifies the expansion rate or the entropy density of the early universe (see e.g. Arbey2009 ()) and thus enters into the calculation of the relic density.

In the left panels of Figs. 5, 7 and 10, we show the cross section for the annihilation of a neutralino pair into a bottom and/or top quark-antiquark pair as a function of the center-of-momentum energy calculated at the tree-level with Yukawa couplings (solid lines) as well as the leading contributions from the individual annihilation channels (dashed lines) and their interferences (dotted lines). Note that some of the latter have been multiplied by in order to fit into the logarithmic plot. We also show, in arbitrary units, the thermal velocity distribution function, involved in the calculation of the thermal average evaluated at the freeze-out temperature (shaded areas).

In the right panels of Figs. 5, 7 and 10, we show again the total annihilation cross section into bottom and/or top quark-antiquark pairs at different levels of precision, i.e. at leading order with couplings (dash-dotted lines), in the approximation included in the micrOMEGAs code (dashed lines), which uses effective couplings to absorb the leading loop effects, and with our full one-loop QCD and SUSY-QCD corrections (solid lines). The shaded area indicates again the thermal velocity distribution of the neutralinos at the freeze-out temperature in arbitrary units.

In order to generalize our results, the remaining Figs. 6, 8, 9, 11, and 12 show scans in various two parameter planes around our parameter points in Tabs. 1 and 2 and display the contours allowed by WMAP calculated at tree-level (green), with the approximation implemented in micrOMEGAs (red), and with our full SUSY QCD corrections (blue). We also show dependence of the relic density on various physical parameters, e.g. the mass of the lightest neutralino.

### iv.1 Relic density in models with non-universal Higgs masses

We begin our detailed numerical discussion by analyzing a scenario, where the annihilation cross section into quarks is dominated by an exchange of a Higgs boson. This scenario corresponds to our parameter point I from Tab. 1. Here the heavy CP-even and the CP-odd Higgs boson resonances coincide with both Higgs boson masses at GeV. Due to the dominance of the Higgs-boson exchange, the neutralino annihilates predominantly into bottom quarks. This is a consequence of the fact that, although the bottom quarks have a much smaller mass compared to the top quarks, their couplings to the Higgs bosons are enhanced by or . As opposed to mSUGRA, the neutralino in this scenario is rather heavy, which allows for an effective annihilation into top quarks as well. In the bottom panels of Fig. 6, a small discontinuity indicates the place, where the top quark final state starts to contribute. These circumstances result in a mixture of top and bottom quark final states, which cannot be reached in mSUGRA for such a low value of .

A few general features can be observed in Figs. 5 and 6. First, on-resonance annihilation of neutralinos would reduce the relic density too much, so that regions allowed by Eq. (1) sit on each side of the resonance peak. In Fig. 5 one can clearly see that the resonance is not aligned with the kinematic region where the biggest contribution to the relic density comes from. The top panels in Fig. 6 display the bands where quark final states dominate the annihilation. In addition, we show a line that corresponds to the position of the Higgs resonance peak. The WMAP allowed regions are situated on each side of the resonance where the width of the these regions reflects the steepness of the peak. The SUSY-QCD corrections to the relic density calculated here are a combination of corrections to the processes with either bottom or top quarks final states. The corrections to the annihilation into bottom quarks are essentially confined to the Higgs-quark-antiquark vertex. The bulk of the correction comes from the gluon exchange and from the SUSY corrections which become large at large and have to be resummed. As we take a moderate even the non-resummable corrections play a role here. These corrections amount to the difference between our full SUSY-QCD result and the effective coupling approximation implemented in micrOMEGAs in this scenario.

The situation with the top quarks is more complicated. As can be seen from Fig. 5, the squark exchanges in the - and -channels and their interference with the Higgs-boson exchange become now also sizable, which renders also the corrections to squark-quark-neutralino vertex and to the squark propagator important. Overall for this scenario, the effect of the full SUSY-QCD corrections amounts to a 20% shift in the prediction of the relic density with respect to the calculation implemented in micrOMEGAs, as can be seen in the bottom panels of Fig. 6. In consequence, the cosmologically favored regions of parameter space are shifted away from the position of the Higgs pole in order to compensate the larger annihilation cross-section. These shifts amount to up to 5 GeV for the common scalar and gaugino masses and a few GeV for the masses of the Higgs doublets at the unification scale (see Fig. 6). This is about two times larger than the experimental precision, and therefore distinct bands appear.

The remaining analyzed scenarios were chosen so as to investigate effects not related to Higgs-boson exchanges. The Higgs-boson masses are very heavy in these scenarios ( GeV) and suppress the importance of the Higgs-boson -channel exchange. The parameter point II in Tab. 1 was chosen as described in Sec. III.1 so that the -boson exchange is enhanced. The only viable way to do this is to increase the higgsino fraction of the neutralino, which in turn increases the coupling of neutralinos to the boson. Another possibility would be to sit on the -boson resonance, which would, however, lead to very small neutralino and chargino masses that are already ruled out by the LEP experiments. The parameter point III in Tab. 1 was picked to use the Higgs potential high-scale parameters to drive down the squark masses via the renormalization group evolution. On top of that, by choosing GeV we induced a large mixing of the third generation squarks. The mixing is biggest for scalar tops increasing the contribution of their exchange. A feature common to both scenarios is a distinct destructive interference between the squark and the boson which links these two scenarios (see Fig. 7). This fact makes corrections to both -boson and squark exchange important in each of the two scenarios. In the case of the parameter point II, the corrections are about 20% in the cross section as compared to the cross section in micrOMEGAs. This is not reflected in all regions of the WMAP allowed regions in Fig. 8, since the contribution of the quark-antiquark final state falls to only 60%-40%, as we lower the higgsino parameter . This decreases the masses of the lightest chargino and of the second-lightest neutralino, which, through their -channel exchange, increase the annihilation cross sections into and final states. Nevertheless, the full SUSY QCD corrections shift the contour in the plane by GeV in and in all instances shift the contour by more than the current experimental precision. The effect of the corrections to the cross section is not screened by other final states in the case of the scenario with a dominant squark exchange (point III in Tab. 1). The top quark final state accounts in certain regions for more than 90% of the annihilation cross section, and the mass of the scalar top is not light enough to allow for efficient co-annihilations. The hyperbolic shape of the WMAP allowed regions in the plane is governed by the Higgs-mass parameter combination . The effect of the corrections on the cross section is about 20% and causes a shift of the preferred value of the pseudo-scalar Higgs boson mass by about GeV.

### iv.2 Relic density in models without gaugino mass unification

The relevant masses of our points III and IV featuring very light stops are very similar, as can be seen in Tabs. 1 and 2. The same holds for the annihilation cross sections. For the benchmark point IV, the exchange of a stop in the - or -channel is therefore favored, which is well visible in the top left panel of Fig. 10. The contribution from the -boson exchange is here one order of magnitude lower than the one from squark exchange. Again, an important role is played by the squark -boson interference effect, which is here even stronger than in the case of point III. The differences are the Higgs boson masses, which were more than  GeV for point III, but are  GeV here (see Fig. 10). The mass difference can be traced back to the mechanism by which we lowered the squark masses. For point III, we had to induce a big difference between the high-scale Higgs potential parameters, which drove the Higgs-boson masses higher, whereas for point IV we changed freely the gluino mass parameter influencing the RGEs for the squark masses, which has not such a big effect on the masses of the Higgs bosons. Due to the large contribution from the exchange of the light squark, the annihilation into top quarks accounts for up to 80% of the total annihilation cross section (see Fig. 11). The subdominant final states here are again and , but also pairs of bottom quarks or leptons. Note that also here the mass difference between the lightest neutralino and the lightest stop is not small enough for efficient co-annihilations. One small, but important difference is that the destructive interference is stronger for point IV than for point III. This leads to negative correction effects when including the full SUSY QCD corrections (when compared to the cross section of micrOMEGAs). The reason is that for this particular point the corrections to the squark exchange are effectively taken into account by including Yukawa couplings and the only difference between our full calculation and the micrOMEGAs approximation are the corrections to the interference which push the full result down (as seen in the top right part of Fig. 10).

Let us now turn to the impact of the corrections in this “compressed SUSY” scenario. The corresponding favored regions in the - and the - planes, the projection on the - plane as well as the prediction of the neutralino relic density as a function of the gaugino mass parameter around our parameter point IV are shown in Fig. 11. Again, in the regions where annihilation into a top quark-antiquark pair is kinematically allowed, the effect of the SUSY-QCD corrections is sizable, resulting in an important shift of the favored regions in the parameter space. The preferred region of parameter space is shifted to higher values of the scalar mass parameter and consequently to higher stop masses . As already discussed above, a heavier stop mass compensates the increasing effect of the additional loop diagrams on the annihilation cross-section. In the right bottom panel of Fig. 11 one sees again that the prediction of the relic density is decreased by the order of 15% with respect to the tree-level calculation and by about 10% with respect to the approximation implemented in micrOMEGAs.

For our point V with a dominant -boson exchange, the correction accounts for about 20% of the micrOMEGAs cross-section. Although micrOMGEAs does not include any correction for the exchange of a -boson, the corresponding curve is approximately 20% over the tree-level prediction, which can be explained by the presence of effective couplings for the sub-leading squark-exchange. The difference between the approximation included in micrOMEGAs and our full one-loop calculation originates from the supplementary corrections, especially for the exchange of the -boson, shown in Figs. 2 and 4 and discussed in Sec. II. We study again the influence of our corrections on the regions of parameter space that are favored with respect to the WMAP limits of Eq. (1). The upper left panel of Fig. 12 shows these regions in the - plane for our scenario V and fixed values of and as given in Tab. 2. In the lower left panel, we show the same contours projected on the corresponding plane of the physical neutralino and stop masses and . In each plot, we also indicate the isocontours corresponding to a contribution from top quark-antiquark final states of 50%, 30%, and 10% to the total neutralino annihilation cross-section. The points lying in the grey shaded areas do not allow for physical solutions of the renormalization group equations. The correction to the relic density is reduced from the 20% it was at the cross section level (as compared with the micrOMEGAs cross section), because the quark-antiquark final state constitutes only about 50% of the annihilation cross section. The remaining contributions include mostly and final states. Nevertheless, as it was also the case for the mSUGRA scenarios analyzed in Ref. Letter2 (), the impact of the one-loop SUSY-QCD corrections is larger than the experimental uncertainty. Therefore distinct bands are observed in wide regions of both the - and in the - planes. Note that for GeV, the annihilation of a neutralino pair into top quarks is kinematically forbidden. The dominating channels for that region are mainly the annihilation into combinations of the gauge bosons and the light Higgs boson. In the - plane, shown in the right top panel of Fig. 12, the WMAP-favored points are confined to rather narrow bands corresponding to the different levels of included corrections. This underlines that is one of the key parameters to which phenomenology is very sensitive. It is interesting that for lower values of their positions are almost independent of , while for higher the dependence becomes stronger. The shift from the tree-level prediction to our full one-loop result is in the direction of higher gluino masses . This is explained by the fact that increasing the gluino mass implies an increase in the squark masses, which in turn implies a decrease of the annihilation cross-section. This effect is then combined with the increase of the cross-section due to the one-loop corrections discussed above, so that the relic density remains in the narrow range which is favored by cosmological data. Finally, in the right bottom panel of Fig. 12, we show the prediction for the neutralino relic density as a function of the lightest stop mass . The neutralino mass has been fixed to the value GeV of our point V. The graph corresponds to a cut through the - plane shown in the lower left panel of Fig. 12. The favored region of Eq. (1) is indicated by a shaded area. Since the SUSY-QCD corrections increase the annihilation cross section by about 50%, the prediction for the neutralino relic density is reduced by about the same amount.

## V Conclusions

The theoretical calculation of the dark matter relic density is an interesting tool to obtain rather stringent constraints on the MSSM parameter space, both at the electroweak and at the grand unification scale. Cosmological precision measurements therefore play an important role in the extraction of SUSY mass parameters from experimental data. In times of increasing experimental sensitivity, in particular for the cosmological parameters such as the relic density of cold dark matter, it is essential to increase the accuracy of the theoretical calculation. Consequently, higher-order corrections to the dark matter annihilation cross section become important, since they enter directly into the calculation of the relic density.

We have presented here the full analytic details of our calculation of the corrections to the annihilation of a neutralino pair into a massive quark-antiquark pair, i.e. of the one-loop and real gluon emission corrections. Annihilation processes into heavy quarks have been shown to be important in large cosmologically allowed regions of SUSY GUT models without scalar or gaugino mass unification, and the effects of our SUSY-QCD corrections have been investigated numerically for these two classes of models. We have identified regions of parameter space, which correspond to current cosmological limits on the dark matter relic density and which feature important annihilations of the neutralino pair into a top quark-antiquark pair through the exchange of a -boson or a scalar quark. Higgs boson resonances were shown to be suppressed in those regions. For five selected parameter points, the effects of the corrections were shown to be sizable, since they enhanced the annihilation cross-section by typically 20 and up to 50% with respect to the tree-level calculation.

As a consequence, the theoretical prediction of the neutralino relic density is also affected by the contributions at the one-loop level. Since the cross section is increased by typically 20 and up to 50%, the relic density is reduced by about the same amount. We have shown that the impact of our corrections is more important than the uncertainty on the observational limits at the confidence level. It is therefore essential to take these corrections into account when analyzing cosmological data with the goal of extracting SUSY mass parameters, which may be shifted by typically 5 GeV and up to 50 GeV, and of determining the favored regions of the SUSY parameter space.

###### Acknowledgements.
The authors would like to thank A. Pukhov for his help in implementing the results into the micrOMEGAs code and W. Porod for useful discussions. The work of K.K. is supported by the ANR projects ANR-06-JCJC-0038-01 and ToolsDMColl BLAN07-2-194882. The work of B.H. is supported by the DFG project PO 1337/1-1.

## Appendix A Notation and couplings

We follow closely the notation of the couplings defined in Ref. GunionHaber () and the conventions used in Ref. Kovarik2005 (); PhDs () and begin by listing all necessary couplings of the -boson. The couplings of fermions to the -boson in the MSSM are identical to those in the SM, and the Lagrangian is

 L = −gcWZ0μ¯fγμ(CfLPL+CfRPR)f, (33)

where is the weak coupling constant, , , () is the (co-)sine of the weak mixing angle , and and are the weak isospin and electric charge of the fermion . The Lagrangian for the interaction of the -boson with two neutralinos is given by

 L = g2cWZ0μ¯~χ0iγμ(O′′LijPL+O′′RijPR)~χ0j, (34)

where

 O′′Lij = −12Zi3Zj3+12Zi4Zj4 = −O′′Rij (35)

depend bilinearly on the neutralino mixing matrix . The Lagrangian of the -boson coupling to two sfermions

 L=−igcWZ0μ~f∗iz~fij\lx@stackrel↔∂μ~fj (36)

is proportional to

 z~fij = CfLR~fi1R~fj1+CfRR~fi2R~fj2, (37)

which depends bilinearly on the sfermion mixing matrix

 R=(RiL,RiR)=(cosθ