Gluino Coannihilation Revisited
KKCL-PH-TH/2015-12, LCTS/2015-04, CERN-PH-TH/2015-049UMN–TH–3424/15, FTPI–MINN–15/10
Gluino Coannihilation Revisited
John Ellis , , Feng Luo and Keith A. Olive , Theoretical Particle Physics and Cosmology Group, Department of Physics,King’s College London, London WC2R 2LS, United Kingdom Theory Division, CERN, CH-1211 Geneva 23, Switzerland School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA William I. Fine Theoretical Physics Institute, School of Physics and Astronomy,University of Minnesota, Minneapolis, MN 55455, USA
Abstract
Some variants of the MSSM feature a strip in parameter space where the lightest neu-tralino χ is identified as the lightest supersymmetric particle (LSP), the gluino ˜ g is thenext-to-lightest supersymmetric particle (NLSP) and is nearly degenerate with χ , and therelic cold dark matter density is brought into the range allowed by astrophysics and cos-mology by coannihilation with the gluino NLSP. We calculate the relic density along thisgluino coannihilation strip in the MSSM, including the effects of gluino-gluino bound statesand initial-state Sommerfeld enhancement, and taking into account the decoupling of thegluino and LSP densities that occurs for large values of the squark mass m ˜ q . We find thatbound-state effects can increase the maximum m χ for which the relic cold dark matter den-sity lies within the range favoured by astrophysics and cosmology by as much as ∼
50% if m ˜ q /m ˜ g = 1 .
1, and that the LSP may weigh up to ∼ m ˜ q /m ˜ g (cid:46) a r X i v : . [ h e p - ph ] A ug Introduction
In the absence of any signal for supersymmetry during Run 1 of the LHC [1], it is natural toask how and where supersymmetry may be hiding. Perhaps it is hiding in plain sight witha compressed spectrum [2] that the conventional missing-energy searches at the LHC havebeen unable to resolve? Or perhaps R parity is violated, in which case supersymmetry maybe hiding among the jets and leptons produced by Standard Model processes? Or perhaps R parity is conserved, but supersymmetric particles are too heavy to have been detectedduring Run 1 of the LHC?There are two issues with this last possibility. One is the accentuation of the problem ofthe naturalness (or fine-tuning) of the electroweak scale that low-scale supersymmetry waspostulated to mitigate, and the other is the cosmological cold dark matter density. The colddark matter may well not consist only, or even predominantly, of the lightest supersymmetricparticle (LSP). However, even if the cold dark matter density is considered only as an upperlimit on the relic LSP density, it imposes an upper bound on the LSP mass that depends onthe specific LSP candidate under consideration.If R parity is conserved and the LSP is present in the Universe today as a relic from theBig Bang, it is expected to be electromagnetically neutral and have only weak interactions.In the minimal supersymmetric extension of the Standard Model (MSSM), the most plausiblecandidates are the lightest neutralino χ and the gravitino [3]. Here we consider the neutralinocase, and the cosmological upper bound on its mass.The relic LSP density depends not only on the LSP mass, but also on the rates at whichit annihilated with itself and coannihilated with other sparticles in the early Universe [4].Other things being equal, the largest LSP mass is allowed when such coannihilation ratesare maximised, which happens when the LSP is (nearly) degenerate with other particles. Ifthere is only one such coannihilating sparticle species, the coannihilation rate will in generalbe maximised for a coloured sparticle. There have been analyses in the literature of the caseswhere the coannihilating particle is a squark, specifically the lighter stop squark [5–11], andalso the case of the gluino [10, 12–18]. In general, one would expect that the heaviest LSPwill be allowed when it coannihilates with the particle with the largest colour charge, namelythe gluino.We study here the question how heavy the neutralino LSP χ could be, if it is nearlydegenerate with, and coannihilates with, the gluino ˜ g . This is of relevance to assessing,for example, what centre-of-mass energy would be needed for a proton-proton collider tobe ‘guaranteed’ to detect R -conserving supersymmetry. There can of course be no cast-1ron guarantee, even within the MSSM. For example, even in the gluino coannihilation casestudied here the neutralino LSP mass limit depends on the squark masses, and the LSP masslimit could be substantially modified if the squarks were degenerate with the neutralino LSPand the gluino. However, a complete analysis of this case lies beyond the scope of this paper.As already mentioned, there have been several previous analyses of neutralino-gluinocoannihilation [10, 12–18], and the main new elements here are in our discussions of theeffects of gluino-gluino bound states and of the issue whether coannihilations can be main-tained in the presence of a large squark to gluino mass ratio. Here, we will restrict ourattention to the coannihilation processes and leave their application to more complete mod-els (with for example, radiative electroweak symmetry breaking) for future work [19]. Aswe discuss in detail, bound states can remove from the primordial plasma gluino pairs thatmay subsequently annihilate into Standard Model particles, before they can decay into theLSP as is usually assumed in discussions of coannihilation. We present numerical estimatesof the bound-state production rate, and find that, for fixed sparticle masses, the relic darkmatter density is substantially reduced compared with the cases where bound-state forma-tion is neglected. Conversely, the cosmological relic density may lie within the cosmologicalrange for substantially larger LSP masses than would have been estimated in the absence ofbound-state effects: this effect is ∼
50% for m ˜ q /m ˜ g = 1 .
1, falling to ∼
20% for m ˜ q /m ˜ g ∼ m ˜ q /m ˜ g > ∼ reducing the possible LSP mass. There is also a reduction in the possibleLSP mass for small m ˜ q /m ˜ g →
1, due to cancellations between s -, t - and u -channel diagramsthat tend to reduce annihilation rates.Taking these effects into account, we find a maximum LSP mass ∼ < ∼ m ˜ q /m ˜ g < ∼ ∼ ∼ χ , gluinos and gluino-gluino bound states. Section 6 contains some numerical results forthe gluino coannihilation strip and a discussion of its endpoint. Section 7 summarises ourconclusions and discusses their significance for future colliders. Finally, Appendices present2echnical aspects of the computation of the 2 → Before discussing the formation and effects of gluino-gluino bound states, we first discussbriefly Sommerfeld effects in gluino-gluino annihilation, which may enhance annihilationrates at low velocities, and are particularly relevant in the case of the strongly-interactinggluino. As a general rule, initial-state interactions modify s-wave cross-sections by factors[20, 21] F ( s ) ≡ − πs − e πs : s ≡ αβ , (1)where α is the coefficient of a Coulomb-like potential whose sign convention is such thatthe attractive case has α <
0, and β is the velocity of one of the annihilating particles inthe centre-of-mass frame of the collision. In the cases of strongly-interacting particles, theCoulomb-like potential has the form [22] V = α s r [ C f − C i − C (cid:48) i ] , (2)where α s is the strong coupling strength, C f is the quadratic Casimir coefficient of a specificfinal-state colour representation, and C i and C (cid:48) i are the quadratic Casimir coefficients of theannihilating coloured particles. In the case of octet annihilating particles such as gluinos, C i = C (cid:48) i = C = 3. The relevant final states are in singlet, octet, or 27 s representations, forwhich C f = 0, C f = 3, or C f = 8.As discussed in [9], Sommerfeld effects such as these have been implemented in the SSARD code [23] for calculating the relic dark matter density. In the coannihilation region of interest,this code uses a non-relativistic expansion for annihilation cross-sections: (cid:104) σv (cid:105) = a + bx − + . . . , (3)where (cid:104) ... (cid:105) denotes an average over the thermal distributions of the annihilating particles, thecoefficient a and b represent the contributions of the s- and p-wave cross-sections, x ≡ m/T ,and the dots represent terms of higher order in 1 /x . A Sommerfeld enhancement occurs when α < ∝ √ x . Inthis paper we have included these enhancements in the ˜ g ˜ g → gg and ˜ g ˜ g → q ¯ q cross sections.The procedure for obtaining a thermally averaged cross section is given in Appendix A. Theexpressions for the matrix elements for the coannihilation processes are given in detail inAppendix B. 3 Gluino-Gluino Bound-State Formation
Gluino-neutralino coannihilations may increase the effective annihilation cross section andthereby lower the final neutralino relic abundance. The Sommerfeld enhancement discussedabove further increases the cross section in specific channels and again lowers the abundanceof neutralinos allowing for larger masses at the tip of the coannihilation strip defined by∆ m = 0 where ∆ m is the gluino-neutralino mass difference [17]. Gluino-gluino bound statescan further serve to remove gluinos from the thermal bath and thereby lower the relic densityby a factor that is non-negligible relative to the Sommerfeld enhancement, and much largerthan the uncertainty in the cosmological cold dark matter density.The dominant process for the formation and dissociation of gluino-gluino bound states˜ R in the thermal plasma is ˜ g + ˜ g ↔ ˜ R + g . These processes become important when theplasma temperature falls low enough for typical thermal energies to become comparable tothe binding energy of the ˜ R state, namely T (cid:46) E B ≡ m ˜ g − m ˜ R . In principle, one may formcolour-octet states as well as singlets, but the latter are expected to be more deeply boundwith larger wave functions at the origin. Here we focus on the production of the lightestcolour-singlet state, 1 s , with orbital angular momentum L = 0 and spin angular momentum S = 0, which is expected to be the most copiously produced. Since we are consideringgluinos weighing several TeV, we expect the leading order of QCD perturbation theory to bea useful approximation, and assume the Coulomb potential V ( r ) = − α s /r for the 1 s state,with binding energy E B (cid:39) (3 α s / m ˜ g . The normalised spatial part of the wave function forthis 1 s bound state is φ bs ( r ) = ( πa ) − / e − r/a , (4)where a ≡ / (3 α s m ˜ g ) is the Bohr radius. The 1 s bound state decays predominantly to apair of gluons and the leading order decay rate isΓ ˜ R = 2434 α s m ˜ g . (5) In order to calculate bound-state formation and dissociation via the dominant processes˜ g a + ˜ g b ↔ ˜ R + g c , we first calculate the bound-state dissociation cross section, σ dis followingSection 56 of [24], where the photoelectric effect for an atom is calculated. The central partof the calculation is the evaluation of the transition amplitude given in Eq. (56.2) of [24]: M fi = (cid:90) φ ∗ f ( − i (cid:126) ∇ · (cid:126)(cid:15) c ( m ˜ g /
2) ) e i(cid:126)k · (cid:126)r φ i d (cid:126)r , (6)4here φ f is the wave function of the free ˜ g a ˜ g b pair and φ i ≡ φ bs ( r ), and (cid:126)(cid:15) c and (cid:126)k are thepolarisation and momentum vectors of the gluon, respectively.We use the dipole approximation, e i(cid:126)k · (cid:126)r ≈
1, which is justified because the bound-statewave function φ bs ( r ) is exponentially suppressed for r > a , and because the gluon momentum | (cid:126)k | = ω , where its energy ω satisfies the conservation condition ω + ω m ˜ g − E B ) = | (cid:126)p | m ˜ g + E B , (7)where | (cid:126)p | is the momentum of one of the annihilating gluinos. (Note that | (cid:126)p | is the same asthe relative momentum, ( m ˜ g / v rel , and the second term on the LHS of (7) can be neglected.)We find ωa (cid:39) E B a = (3 α s )2 (cid:28) v rel = 0 and α s = 0 .
1, and more generally ωa < v rel < .
6, so that the dipole approximation should be sufficient for our purposes.The dipole approximation imposes a selection rule on φ f , which needs to be in an L = 1 state. Further, charge conjugation ( C -parity) conservation requires that C (˜ g a ˜ g b ) = C ( ˜ R ) C ( g c ), where the 1 s ground state with L = 0 and S = 0 has J P C = 0 − + . The C -parity of the colour anti-symmetric 8 A state is the same as that of the gluon [25], while the C -parity of colour-symmetric 8 s state is opposite of that of the gluon, for all color indices.Therefore, the only possible state for φ f is 8 A , with L = 1 and S = 0. (Note also that parityis conserved in this case, because P ( φ f ) = 1 and the gluon has P = − g a ˜ g b is φ f = 12 | (cid:126)p | ∞ (cid:88) L =0 i L (2 L + 1) e − iδ L R pL ( r ) P L ( (cid:126)p · (cid:126)r | (cid:126)p | r ) . (8)Only the L = 1 term survives, due to the selection rule from the dipole approximation. Sincewe wish to calculate |M fi | , we may discard the phase shift factor e − iδ L ( δ L is real) and thefactor i L . Therefore we write φ f = 32 | (cid:126)p | R p ( r ) P ( (cid:126)p · (cid:126)r | (cid:126)p | r ) . (9)so that dσ dis = α s ( m ˜ g / | (cid:126)p | πω |M fi | d Ω (cid:126)p , (10)where M fi is calculated following Section 56 of [24].Since φ f is the wave function for an 8 A state, the Coulomb potential is V f ( r ) = − α s /r ,whereas φ i is a wave function for the Coulomb potential V ( r ) = − α s /r , the result is differentfrom Eq. (56.12) of [24], namely σ dis = 2 π α s a (cid:18) E B ω (cid:19) (1 + ξ )[1 + ( κξ ) ] e − ξ arctan(1 /κξ ) − e − πξ κ − , (11)5here ξ = α s /v rel and κ = 2. This equation is averaged over the gluon polarizations andwould reduce to Eq. (56.12) if κ = 1.The total wave functions for the free ˜ g a ˜ g b pair and the bound state ˜ R are products ofthe spin, colour and spatial parts of the wave functions. In view of the Majorana nature ofthe gluinos, the total wave functions should be anti-symmetric. Concerning the spin part ofthe wave function, since both the bound state and the free gluino pair are in an S = 0 state,the spin wave functions are both ( ↑↓ − ↓↑ ) / √ , (12)and the spin parts of the wave functions do not introduce any extra factor in σ dis . As for thecolour part of the wave function, according to [26,27] φ i,color = δ de / √
8, and φ f,color = f hae / √ f abc f abd = 3 δ cd ).The ( − i (cid:126) ∇· (cid:126)(cid:15) c ( m ˜ g / ) e i(cid:126)k · (cid:126)r factor in the transition amplitude (6) is calculated from the gluino-gluino-gluon interaction Lagrangian L = i g s g cµ f abc ¯˜ g a γ µ ˜ g b , (13)which can be compared with the corresponding QED interaction Lagrangian L = − eQ f A µ ¯ f γ µ f . (14)We simply replace the electric charge factor Q f in the transition amplitude (6) by the colourfactor f abc , since the factor 1 / | √ f had √ f cda | = | √ √ δ ch | = 3 . (15)Note that all color indices are summed over.Concerning the spatial part of the wave function, we need to take into account the factthat both the initial and final states contain two identical particles. In the case of the boundstate, they are in the symmetric L = 0 state, and the wave function needs to be symmetrisedas in Eq. (2.14) of [26]:1 √ φ bs ( r, θ, φ ) + φ bs ( r, π − θ, φ + π )] = √ × ( πa ) − / e − r/a . (16)On the other hand, the final free pair is in the antisymmetric L = 1 state, and the wavefunction needs to be antisymmetrised:1 √ φ f ( r, θ, φ ) − φ f ( r, π − θ, φ + π )] = √ × | (cid:126)p | R p ( r ) P ( (cid:126)p · (cid:126)r | (cid:126)p | r ) . (17)6he coefficients in these two equations introduce an extra factor of |√ √ | = 4 into themodulus-squared of the spatial wave function factors. Finally, recall that we have averagedover the polarisations of the gluon, but we also need to average over its colour. This gives afactor of 1 /
8. Therefore, the final dissociation cross section is σ dis = 3 × × × × σ dis , (18)where the final factor of 1/2 is to avoid double counting of gluinos in the final-state phase-space integration. We come finally to the bound-state formation cross section, σ bsf , which is related to σ dis through the Milne relation:12 n eq ˜ g n eq ˜ g σ bsf v rel (cid:18) e ω/T − (cid:19) f ( v rel ) dv rel = n eq ˜ R σ dis dn eqg , (19)where the on the LHS of the above equation is introduced to avoid double-counting thenumber of bound-state formation reactions, and the factor e ω/T − comes from the enhance-ment of bound-state formation due to the stimulated gluon emission in the thermal back-ground (similar to the stimulated recombination in e − p ↔ Hγ ). Using dn eqg = g g π (2 π ) ω dωe ω/T − ,f ( v rel ) = (cid:18) m ˜ g / πT (cid:19) / πv rel e − ( m ˜ g / v rel / T ) (20)and (7), we find σ bsf = 2 g ˜ R g g ω g g [( m ˜ g / v rel ] σ dis , (21)where g ˜ R g g g g = 1 × (2 × × = 116 . (22)For comparison, the Sommerfeld enhanced s-wave cross section for ˜ g a ˜ g b → g c g d is given inEqs. (2.13) and (2.25) of [17]: S ann ( σ ann v rel ) = (cid:18)
16 2 π (2 ξ )1 − e − π (2 ξ ) + 13 2 πξ − e − πξ + 12 2 π ( − ξ )1 − e π ( ξ ) (cid:19) (4 πα s ) m g π , (23)7here ξ = α s /v rel . Therefore, in the v rel → σ bsf v rel S ann ( σ ann v rel ) = 323 e ≈ . . (24)Therefore, we see that the inclusion of ˜ g ˜ g bound states is a non-negligible component indetermining the final neutralino relic density. For coannihilation to be effective, the coannihilating species (in this case neutralinos andgluinos) must be in thermal contact. That is, the rates for interconverting the LSP andNLSP must be faster than the Hubble rate. In both the familiar cases of stop and staucoannihilation, connectivity of the two species can be taken for granted, as the conversionrates are mediated by light Standard Model particles and are always fast. This implies thatthe ratio of densities ( n NLSP /n LSP ) is approximately equal to the equilibrium ratio and allowsfor a simplification in the Boltzmann equations. However, the interconversion of neutralinosand gluinos must proceed via squarks, leading to a suppression if the squarks are heavy.The relevance of the coannihilation process relies on fast conversion rates, and requires theratio of squark masses to the gluino mass to be less than approximately 100 as we showbelow. For larger squark masses, the gluino and neutralino abundances evolve separately,and coannihilation effects are essentially shut off independent of the mass difference.The interconversion processes we consider are χq ↔ ˜ gq , χ ¯ q ↔ ˜ g ¯ q , and the gluino decaysand the inverse decays ˜ g ↔ χq ¯ q . When the neutralino is a Wino or Higgsino, the processesinvolving a chargino, χ + d ↔ ˜ gu , χ + ¯ u ↔ ˜ g ¯ d and ˜ g ↔ χ + d ¯ u , as well as the correspondingprocesses for χ − , are also included. We note that q stands here for all six quark flavors, andthe u, d stand for all the three generations of up-type and down-type quarks. Also, when χ is a Higgsino, the two lightest mass-degenerate neutralino components, ˜ H , , are both takeninto account. For each relevant process, we first calculate the transition matrix element |T | .We calculate the gluino decay rates for ˜ g → χq ¯ q , ˜ g → χ + d ¯ u and its charge-conjugatedprocess . The squared transition matrix elements |T | are identical to the corresponding Since we treat quarks as massless, three-body decays will always occur as long as m ˜ g > m χ . Due tologarithmic corrections, the gluino two-body decay into a neutralino and a gluon becomes important whenthe squark-to-gluino mass ratio is very large ( (cid:29) (cid:46) g and χ are so close in mass that all the gluino three-body decays would bekinematically forbidden, if quarks are not treated as massless. × = for the coannihilation processes, whereas it is for the gluinodecay processes. We note also that the definitions of the Mandelstam variables should alsobe changed correspondingly as follows: for the coannihilations, s = ( p + p ) , t = ( p − p ) and u = ( p − p ) , whereas for the gluino decays, s = ( p − p ) , while t and u do notchange.The gluino decay rates are then obtained by performing the standard 3-body phase spaceintegration. The inverse-decay processes do not have to be calculated separately, becausethey are taken into account automatically by the Boltzmann equations given in the nextsection.To calculate the conversion rates for χq → ˜ gq , χ + d → ˜ gu , χ + ¯ u → ˜ g ¯ d and their charge-conjugated processes, we first calculate the cross sections. Again, the squared transitionmatrix elements |T | are identical to the corresponding ones for the coannihilation processesgiven in Appendix B, except that the expressions in Appendix B should be multiplied bya factor of , because the factor for the initial color averaging is for the coannihilations,whereas it is for the conversions. Also, compared to the coannihilation processes, theMandelstam variables for the conversion processes are re-defined as s = ( p − p ) , t =( p ± p ) and u = ( p ∓ p ) , where the upper signs in the definition of t and u apply if ¯ q B or ¯ d B is brought into the initial state, while the lower signs apply if q A or u A is pulled overto the initial state.For each of the quark flavors, the thermally-averaged conversion rate is obtained byintegrating σ c v q over the Fermi-Dirac distribution of the quark in the initial state, (cid:104) Γ c (cid:105) = (cid:90) σ c v q dn q = (cid:90) + ∞ E qmin σ c v q · · π (2 π ) | (cid:126)p q | d | (cid:126)p q | e E q /T + 1 , (25)where σ c is the conversion cross section for any of the relevant processes discussed above. Inthe initial neutralino or chargino rest frame, σ c is a function of the incoming quark energy E q . In this reference frame, t or u is the squared center-of-mass energy, and is given by m χ + m q + 2 m χ E q , where the lower limit of E q is E q min = [( m ˜ g + m q (cid:48) ) − m χ − m q ] / (2 m χ ),where q (cid:48) represents the quark in the final state. Here v q is the velocity of the incoming quark,and it is related to the energy and 3-momentum of the quark by v q = | (cid:126)p q | /E q . The factors3 and 2 in (25) count the quark color and spin degrees of freedom, respectively. Again, theinverse conversion rates are taken into account automatically by the Boltzmann equations.9 Boltzmann Equations
We are now in a position to put all of the components discussed above into a rate equation(or set of equations) in order to solve for the relic density. To do so, we begin by consideringthree separate density components: neutralinos, gluinos and bound states.To set up a coupled set of Boltzmann equations, it is convenient to rescale the numberdensities of neutralinos, gluinos and bound states by the entropy density, Y χ ≡ n χ s , Y ˜ g ≡ n ˜ g s , Y ˜ R ≡ n ˜ R s . (26)These are governed by the following set of coupled Boltzmann equations : dY χ dx = xsH ( m χ ) (cid:104) − (cid:104) σv (cid:105) χχ (cid:0) Y χ Y χ − Y eqχ Y eqχ (cid:1) − (cid:104) σv (cid:105) χ ˜ g (cid:0) Y χ Y ˜ g − Y eqχ Y eq ˜ g (cid:1) − (cid:88) q (cid:104) Γ c (cid:105) s (cid:18) Y χ − Y eqχ Y ˜ g Y eq ˜ g (cid:19) + (cid:104) Γ (cid:105) ˜ g s (cid:18) Y ˜ g − Y eq ˜ g Y χ Y eqχ (cid:19) (cid:105) , (27) dY ˜ g dx = xsH ( m χ ) (cid:104) − (cid:104) σv (cid:105) ˜ g ˜ g (cid:0) Y ˜ g Y ˜ g − Y eq ˜ g Y eq ˜ g (cid:1) − (cid:104) σv (cid:105) χ ˜ g (cid:0) Y χ Y ˜ g − Y eqχ Y eq ˜ g (cid:1) + (cid:88) q (cid:104) Γ c (cid:105) s (cid:18) Y χ − Y eqχ Y ˜ g Y eq ˜ g (cid:19) − (cid:104) Γ (cid:105) ˜ g s (cid:18) Y ˜ g − Y eq ˜ g Y χ Y eqχ (cid:19) −(cid:104) σv (cid:105) bsf (cid:32) Y ˜ g Y ˜ g − Y eq ˜ g Y eq ˜ g Y ˜ R Y eq ˜ R (cid:33) (cid:105) , (28) dY ˜ R dx = xsH ( m χ ) (cid:104) − (cid:104) Γ (cid:105) ˜ R s (cid:16) Y ˜ R − Y eq ˜ R (cid:17) − (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg Y ˜ g (cid:16) Y ˜ R − Y eq ˜ R (cid:17) + 12 (cid:104) σv (cid:105) bsf (cid:32) Y ˜ g Y ˜ g − Y eq ˜ g Y eq ˜ g Y ˜ R Y eq ˜ R (cid:33) (cid:105) , (29)where x ≡ m χ T , s = 2 π g ∗ s T , H ( m χ ) ≡ H ( T ) x = (cid:18) π G N g ∗ (cid:19) m χ , (30)and g ∗ s and g ∗ are the total numbers of effectively massless degrees of freedom associated withthe entropy density and the energy density, respectively, (cid:104) σv (cid:105) χχ is the relative velocity timesthe total cross section for the channels for χχ annihilation into Standard Model particles, and (cid:104) σv (cid:105) χ ˜ g and (cid:104) σv (cid:105) ˜ g ˜ g are to be understood similarly, (cid:80) q (cid:104) Γ c (cid:105) and (cid:104) Γ (cid:105) ˜ g are the total conversionrate and gluino decay rate discussed in the previous section, and all possible quark and A set of coupled Boltzmann equations for the photino and a gluino R -hadron was studied in [29]. In adifferent context, Boltzmann equations involving a bound state can be found, for example, in [30]. χ are summed over, (cid:104) Γ (cid:105) ˜ R is the decay rate of the ˜ R , and (cid:104) σv (cid:105) bsf is the bound-state formation cross section times the relative velocity of the two incominggluinos, taking into account the 1 / ( e ω/T −
1) enhancement factor as discussed in Section 3.2.Finally, (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg Y ˜ g has the same effect as (cid:104) Γ (cid:105) ˜ R /s , namely, it converts the bound states togluons without altering the density of free gluinos. All the quantities bracketed by (cid:104) . . . (cid:105) arethermally averaged, and the superscript ‘eq’ denotes equilibrium yields.Eq. (29) can be written in a more intuitive form: d ln Y ˜ R d ln x = − (cid:104) Γ (cid:105) ˜ R + (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g H ( T ) (cid:32) − Y eq ˜ R Y ˜ R (cid:33) + (cid:104) σv (cid:105) bsf n ˜ g (cid:16) Y ˜ g Y ˜ R (cid:17) H ( T ) (cid:34) − (cid:18) Y eq ˜ g Y ˜ g (cid:19) (cid:32) Y ˜ R Y eq ˜ R (cid:33)(cid:35) . (31)One can check that the LHS of Eq. (31) is of order -10, whereas each of the terms on the RHSof Eq. (31) are of order α s M P /m χ , where M P = G − / N . Hence, to a good approximation, wecan set the two terms on the RHS equal to each other and solve for (cid:16) Y ˜ R Y eq ˜ R (cid:17) : Y ˜ R Y eq ˜ R = (cid:104) Γ (cid:105) ˜ R + (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g + (cid:104) Γ (cid:105) dis (cid:16) Y ˜ g Y eq ˜ g (cid:17) (cid:104) Γ (cid:105) ˜ R + (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g + (cid:104) Γ (cid:105) dis , (32)where (cid:104) Γ (cid:105) dis = 12 (cid:104) σv (cid:105) bsf ( n eq ˜ g ) /n eq ˜ R . (33)Therefore, we find that d ( Y χ + Y ˜ g ) dx = xsH ( m χ ) (cid:104) − (cid:88) i,j = χ, ˜ g (cid:104) σv (cid:105) ij (cid:0) Y i Y j − Y eqi Y eqj (cid:1) −(cid:104) σv (cid:105) bsf (cid:104) Γ (cid:105) ˜ R + (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g (cid:104) Γ (cid:105) ˜ R + (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g + (cid:104) Γ (cid:105) dis (cid:0) Y ˜ g Y ˜ g − Y eq ˜ g Y eq ˜ g (cid:1) (cid:105) . (34)Moreover, we note that (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg n ˜ g is much smaller than (cid:104) Γ (cid:105) ˜ R for x (cid:38)
10, due to thefact that n ˜ g decreases with the decrease of temperature while (cid:104) Γ (cid:105) ˜ R is nearly temperatureindependent. Since the process ˜ g ˜ R → ˜ gg is related to the bound-state formation process bycrossing, (cid:104) σv (cid:105) ˜ g ˜ R → ˜ gg should be related to (cid:104) σv (cid:105) bsf by a coefficient not too much different fromorder 1.If at least one of the (cid:80) q (cid:104) Γ c (cid:105) and (cid:104) Γ (cid:105) ˜ g is sufficiently larger than H ( T ) throughout theperiod during which ( Y χ + Y ˜ g ) changes substantially, which is the case when the squark massappearing in the denominators of the matrix elements for these processes is not too large,Eq. (34) can be solved by using the very good approximation Y ˜ g /Y χ ≈ Y eq ˜ g /Y eqχ . In this case,11q. (34) can be recast in the familiar form suitable for coannihilation calculations, and wecan write dYdx = − xsH ( m χ ) (cid:18) − x g ∗ s dg ∗ s dx (cid:19) (cid:104) σ eff v rel (cid:105) (cid:0) Y − ( Y eq ) (cid:1) , (35)where we have included the term x g ∗ s dg ∗ s dx which takes into account the evolution of g ∗ s withtemperature. As we will see, this approximation is valid so long as m ˜ q /m ˜ g < ∼ Y = n/s , where n is interpreted as the total number density, n ≡ (cid:88) i n i = n χ + n ˜ g . (36)and Y eq = n eq /s , where n eq is the total equilibrium number density, n eq ≡ (cid:88) i n eq ,i = n eqχ + n eq ˜ g . (37)The effective annihilation cross section is (cid:104) σ eff v rel (cid:105) ≡ (cid:88) ij n eq ,i n eq ,j n (cid:104) σ ij v rel (cid:105) . (38)As one can see from Eq. (34), the expression for (cid:104) σv (cid:105) ˜ g ˜ g is the ‘standard’ term in the firstline of (34) combined with the second line involving the bound states. We re-emphasise thatthis simplification requires a fast interconversion rate as discussed in the previous Section,so that we can set ( Y ˜ g /Y χ ) = ( Y eq ˜ g /Y eqχ ), which is true only when m ˜ q /m ˜ g < ∼
20. For largersquark masses, we use the coupled set of Boltzmann equations to solve for the relic density.When the LSP is a Wino or a Higgsino, we can still use all the above equations to solvefor the relic density. All we need to do is re-define the following quantities to include thecontributions from each of the χ components, χ i , neutral or charged, as n χ ≡ (cid:88) χ i n χ i ,n eqχ ≡ (cid:88) χ i n eqχ i , (cid:104) σv (cid:105) χχ ≡ (cid:88) χ i ,χ j (cid:104) σv (cid:105) χ i χ j r χ i r χ j (cid:104) σv (cid:105) χ ˜ g ≡ (cid:88) χ i (cid:104) σv (cid:105) χ i ˜ g r χ i , (cid:104) Γ c (cid:105) ≡ (cid:88) χ i (cid:104) Γ c (cid:105) χ i q → ˜ gq (cid:48) r χ i , (cid:104) Γ (cid:105) ˜ g ≡ (cid:88) χ i (cid:104) Γ (cid:105) ˜ g → χ i q ¯ q (cid:48) , (39)12here q (cid:48) is the same as q when χ i is a neutralino, and they are different when χ i is a chargino,and the q and q (cid:48) indicate all the possible quark and anti-quark channels for the conversionrates and gluino decay rates. In Eq. (39), r χ i ≡ n eqχ i /n eqχ = n χ i /n χ , where the latter ‘=’ isguaranteed by the fast conversion and/or decay rates among the different χ i ’s. For laterdiscussion, it is useful to define an effective number of degrees of freedom for χ : g χ eff ≡ (cid:88) χ i g χ i (1 + ∆ χ i ) / exp( − ∆ χ i m χ /T ) , (40)where ∆ χ i ≡ ( m χ i /m χ − χ is the lightest component (i.e., the LSP)among the χ i ’s. We can then write r χ i explicitly as r χ i = g χ i g χ eff (1 + ∆ χ i ) / exp( − ∆ χ i m χ /T ) . (41)In the limit that all the χ components have the same mass, g χ eff = 2, 6 and 8 for Bino,Wino and Higgsino, respectively. We now present some numerical results obtained using the above formalism. Our results inthis section are based on simplified supersymmetric spectra defined at the weak scale. Weassume degenerate squark masses, m ˜ q and for the most part, our results do not depend onsupersymmetric parameters such as µ , A , and tan β . We assume that the neutralino isa pure state of either a Bino, Wino, or Higgsino. Thus our free parameters are simply theneutralino mass, m χ , the gluino mass, m ˜ g and the squark masses, m ˜ q . In future work, weapply these results to more realistic CMSSM-like models (without gaugino mass universality)and pure gravity mediation models with vector-like multiplets [19].We begin with the case in which the lightest neutralino χ is the Bino. Fig. 1 compares a naive calculation using a single Boltzmann equation for the total relicabundance (left panel) with a treatment of the three coupled Boltzmann equations for thegluino, Bino and gluino-gluino bound state abundances (right panel). These results arefor the representative case m χ = 7 TeV, ∆ m ≡ m ˜ g − m χ = 40 GeV and m ˜ q /m ˜ g = 10.The dashed line in the left panel shows the total relic abundance as would be given by the The exception is the case of Higgsino-gluino coannihilations for which the vertices do depend on tan β . (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) m Χ (cid:144) T n (cid:144) s
20 50 100 500200 10 (cid:45) (cid:45) (cid:45)
10 20 50 100 200 50010 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) m Χ (cid:144) T n (cid:144) s Figure 1:
Left panel: The evolution of the total supersymmetric particle abundance n/s as afunction of m χ /T for the representative case m χ = 7 TeV, ∆ m ≡ m ˜ g − m χ = 40 GeV, and m ˜ q /m ˜ g = 10 using the single Boltzmann equation (34). The dashed line exhibits the naivethermal equilibrium abundance, and the solid line shows the numerical solution of equation(34) for the sum of the neutralino and gluino densities, exhibiting the familiar freeze-outwhen m/T ∼ . Right panel: the separate evolutions of the abundances (solid lines) usingthe full set of equations given by Eqs. (27 - 29) for the gluino (blue lines), Bino (red lines),gluino-gluino bound-state (green lines) and sum of the Bino and gluino densities (black lines)as functions of m χ /T . Details of the evolutions of the abundances are shown in the inset.The dashed lines again exhibit the naive thermal equilibrium abundances. Boltzmann distribution if thermal equilibrium were maintained, and we see a clear departurefor m/T (cid:38)
30, as expected from a freeze-out calculation. In the right panel of Fig. 1, the blueline shows the evolution of the gluino abundance, the red line that of the Bino abundance,the green line that of the gluino-gluino bound states, and the black line that of the sumof the gluino and Bino densities. The inset shows details of the evolutions of the gluon,Bino and bound-state abundances. The dashed lines again show the corresponding naivethermal equilibrium abundances. With the stated choices of m ˜ g , m χ and m ˜ q , the relic darkmatter density using the single Boltzmann equation is Ω χ h = 0 . χ h = 0 . m ˜ q /m ˜ g for the same representative values m χ = 7 TeVand ∆ m ≡ m ˜ g − m χ = 40 GeV. The left panel is for m ˜ q /m ˜ g = 1 . , in which case we findthat the relic cold dark matter density is higher than previously: Ω χ h = 0 .
21. This is dueto the fact that at a low squark to gluino mass ratio, there is a cancellation among the t We do not show results for smaller values of m ˜ q /m ˜ g , since then squark-Bino coannihilations should alsobe taken into account. (cid:45) (cid:45) (cid:45)
10 20 50 100 200 50010 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) m Χ (cid:144) T n (cid:144) s
20 50 100 500200 10 (cid:45) (cid:45) (cid:45) (cid:45)
10 20 50 100 200 50010 (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) m Χ (cid:144) T n (cid:144) s Figure 2:
As in the right panel of Fig. 1, but for the choices m ˜ q /m ˜ g = 1 . (left), (right). and u channel annihilations with the s channel leading to a smaller gluino annihilation crosssection and hence a larger relic density. The results also change, even more significantly, forlarge values of m ˜ q /m ˜ g , as shown in the right panel of Fig. 2, where m ˜ q /m ˜ g = 120. In thiscase, we find a much larger value of Ω χ h = 6 .
0. In both panels, the insets show details of theevolutions of the gluino, Bino and bound-state abundances. At still larger values of m ˜ q /m ˜ g the relic density grows very sharply: for example, for m ˜ q /m ˜ g = 150, we find Ω χ h = 780, andfor m ˜ q /m ˜ g = 200 we find Ω χ h = 1 . × assuming the same gluino-Bino mass differenceof 40 GeV. These large numbers reflect the failure of gluino-neutralino conversion to keeppace with the Hubble expansion for large m ˜ q /m ˜ g .In order to summarize the effects of both the cancellations in the annihilation crosssection at low m ˜ q /m ˜ g and the decoupling of the gluino coannihilations at high m ˜ q /m ˜ g , weshow in Fig. 3, the relic neutralino density as a function of m ˜ q /m ˜ g for our nominal valueof m χ = 7 TeV, and ∆ m ≡ m ˜ g − m χ = 0 ,
40, and 120 GeV (black, red, and blue lines,respectively). We see clearly the rise in Ω χ h at small m ˜ q /m ˜ g as well as the very rapid rise inΩ χ h at high m ˜ q /m ˜ g (cid:38) χ h , as exemplified bythe case m ˜ q /m ˜ g = 10 shown in Fig. 1. In general, there is a shallow minimum in Ω χ h around m ˜ q /m ˜ g ∼
50 whose location depends on ∆ m . The horizontal band indicates the 3- σ rangefor the Planck determination of the cold dark matter density of Ω h = 0 . ± . m χ , ∆ m ) plane where 0 . < Ω χ h < . σ below and above the current central value [31]) for a selection of values of m ˜ q /m ˜ g , ascalculated in various dynamical approximations. The red bands were calculated droppingboth the Sommerfeld enhancement factor and the effect of gluino bound-state formation. Aswas already noted in [9, 17] the Sommerfeld enhancement causes a significant suppression of15
10 1005 50 5000.050.550.11 m q (cid:142) (cid:144) m g (cid:142) (cid:87) Χ h Figure 3:
The relic cold dark matter density Ω χ h as a function of m ˜ q /m ˜ g for m χ = 7 TeVand the choics ∆ m ≡ m ˜ g − m χ = 0 , , and 120 GeV (from bottom to top, black, red, and bluelines, respectively). The rise at small m ˜ q /m ˜ g is due to the cancellations between the s -, t - and u -channel diagrams for gluino pair annihilation, and the rise at large m ˜ q /m ˜ g is due to thedecoupling of the gluino and neutralino densities. The horizontal band indicates the 3- σ rangefor the Planck determination of the cold dark matter density of Ω h = 0 . ± . [31]. Ω χ h for fixed values of the model parameters, so long as the gluino-Bino conversion rates arelarge enough that the gluino coannihilation is effective. Correspondingly, the orange Ω χ h bands, calculated including the Sommerfeld factor, appear at larger values of ∆ m and extendto larger values of m χ . The effect of including bound-state effects is to suppress further thevalue of Ω χ h for fixed model parameters, so that the corresponding black Ω χ h bands inFig. 4 extend to even larger values of ∆ m and m χ . We also show in Fig. 4 (coloured purple)the bands that would be found if the bound-state formation rate were a factor 2 larger thanour calculations, as might arise from higher-order QCD or other effects . We evaluate the α s appearing in the Sommerfeld enhancement factor at a scale βm ˜ g that is typicalof the momentum transfer of the soft-gluon exchanges responsible for the Sommerfeld effect [32], and take β = 0 .
3, which is comparable to the thermal velocities of the gluinos at the freeze-out temperature. Becausethe Sommerfeld enhancement is a precursor to the formation of bound states, for simplicity we take thesame α s in evaluating the bound-state effects. The full QCD potential and the thermal mass of the gluonwere considered in computing the Sommerfeld enhancement in [17], and they should be relevant also to thecalculation of the bound-state effects. We estimate that these effects could result in a shift of our orangebands inward by ∼ ∼
000 4000 6000 8000 10 000050100150200 m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) Figure 4:
The ( m χ , ∆ m ≡ m ˜ g − m χ ) planes for a Bino LSP, exhibiting bands where . < Ω χ h < . (3 σ below and above the current central value), for different valuesof m ˜ q /m ˜ g = 1 . (upper left), (upper right), (lower left) and (lower right). Theseresults are calculated without the Sommerfeld enhancement factor and gluino bound-state for-mation (red bands), with the Sommerfeld enhancement factor but without gluino bound-stateformation (orange bands), with both the Sommerfeld enhancement factor and gluino bound-state formation (black bands), and allowing for the possibility that the bound-state formationrate is a factor 2 larger than our calculations (purple bands). m ˜ q /m ˜ g = 1 .
1, where the t and u channelspartially cancel the s -channel contributions to the gluino annihilation cross section. Herewe see that the black band calculated including both the Sommerfeld enhancement andgluino bound-state effects extends to m χ ∼ . m χ ∼ . m ˜ q /m ˜ g = 10 (upper right panel of Fig. 4), the effect of bound-state formation is somewhatsmaller than the Sommerfeld effect, and the black (purple) band extends to m χ ∼ m ˜ q /m ˜ g = 50 (lower left panel of Fig. 4), where theblack and purple bands also extend to m χ ∼ m ˜ q /m ˜ g = 120 (lower right panel of Fig. 4) are quite different. The Sommerfeld effect is muchlarger than the bound-state effect though the latter is still slightly larger than the widths ofthe coloured bands corresponding to the 3- σ ranges for the cold dark matter density. Also,the allowed range of the LSP mass is greatly reduced, extending only to ∼ . ∼ . m ˜ g = m χ ) of the gluino coannihilation strips as functions of the assumed value of Ω χ h , again fromcalculations with neither the Sommerfeld enhancement nor gluino bound states (red lines),with the Sommerfeld enhancement but without the bound states (orange lines), with botheffects included (black lines), and allowing for a factor 2 uncertainty in the bound-state effects(purple lines). The horizontal green bands again show the 3- σ band 0 . < Ω χ h < . χ h would be relevant if the LSP provides only part of thedark matter density, e.g., if there is also a significant axion component, and parameter choicesyielding higher values of Ω χ h in conventional Big Bang cosmology (as assumed here) couldbe relevant in models with non-standard cosmological evolution . As was also seen in Fig. 4,the smallest value of m ˜ q /m ˜ g = 1 . m χ for any fixedvalue of Ω χ h , as compared to the m ˜ q /m ˜ g = 10 case (upper right). The choice m ˜ q /m ˜ g = 50 obtained by considering only the ground state, but the color charge in the gluino case may change this value.We therefore plot purple bands with an uncertainty of a factor 2 to allow for these uncertainties. It is worth recalling that for m ˜ g ∼ O (50) GeV and the freeze-out temperature in conventional Big Bang cosmology is hundreds of GeV, so the assumption made in thispaper of standard cosmological evolution is rather different from that made in more conventional thermaldark matter scenarios where the freeze-out temperature may be in the GeV range.
000 4000 6000 8000 10 0000.000.020.040.060.080.100.120.14 m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) Figure 5:
The locations of the endpoints of the gluino coannihilation strips for different valuesof Ω χ h , using the same colour conventions as in Fig. 4. As in that Figure, the calculationsassume the different values m ˜ q /m ˜ g = 1 . (upper left), (upper right), (lower left) and (lower right). The horizontal green bands show the 3- σ band . < Ω χ h < . . (lower left) leads to a marginally smaller value of m χ , and the choice m ˜ q /m ˜ g = 120 (lowerright) leads to significantly lower values of m χ for any fixed value of Ω χ h . The effect ofincluding bound-state effects is to increase the range of m χ compatible with the measuredvalue of Ω χ h by ∼
50% for m ˜ q /m ˜ g = 1 .
1, decreasing to ∼
20% for m ˜ q /m ˜ g = 10 to 50.Finally, we show in Fig. 6 the value of m χ at the endpoint of the coannihilation strip when∆ m = 0 and Ω χ h = 0 . ± . m ˜ q /m ˜ g : the brown andred contours are for Ω χ h = 0 .
05 and 0 .
15, respectively. The band and contours exhibit theinverse of the behaviour of the relic density seen previously in Fig. 3. The neutralino massat low m ˜ q /m ˜ g is below the maximum value of m χ , which has a shallow maximum around m ˜ q /m ˜ g = 10 to 50, and falls sharply when m ˜ q /m ˜ g (cid:38) g − χ conversion. We conclude that, within the framework studied here, m χ (cid:46) ∼ m q (cid:142) (cid:144) m Χ m Χ (cid:64) G e V (cid:68) Figure 6:
The value of m χ at the endpoint of the gluino coannihilation strip when ∆ m = m ˜ g − m χ = 0 in the Bino LSP case, as a function of m ˜ q /m ˜ g . The drop at small m ˜ q /m ˜ g is due tothe cancellations between the s -, t - and u -channel diagrams for gluino pair annihilation, andthat at large m ˜ q /m ˜ g is due to the decoupling of the gluino and neutralino densities. The greenband corresponds to the current - σ range of the dark matter density: Ω χ h = 0 . ± . ,and the brown and red contours are for Ω χ h = 0 . and . , respectively. We now consider the case of a Wino LSP. The left panel of Fig. 7 displays the gluino-Winocoannihilation strips for Ω χ h = 0 . ± . m ˜ q /m ˜ g = 10, using the same colourcodings as for the Bino case (red with neither the Sommerfeld enhancement nor gluinobound states, orange including the QCD Sommerfeld enhancement but again no bound-state effects, black with both effects included, and purple with the bound-state formationrate enhanced by a factor 2). We see that in this case the black coannihilation strip extendsto m χ ∼ m χ . The reason is that even inthe absence of coannihilation (large m ˜ g − m χ ), Wino-Wino annihilations are strong enoughto suppress the relic density below the density indicated by Planck and other experiments20
000 4000 6000 8000 10000050100150200 m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) Figure 7:
As in Fig. 4 and 5, but for a Wino LSP. when m χ < ∼ . The right panel of Fig. 7 shows how Ω χ h at the endpoints of the stripsvaries with m χ . As previously, the colours of the lines correspond to the colours of the stripsin the left panel of Fig. 7. We see that the black line crosses the horizontal green band whereΩ χ h = 0 . ± . m χ ∼ χ h = 0 . ± . χ h = 0 .
05 and 0 .
15. We see that Ω χ h is within the preferred range for m χ ∼ (cid:46) m ˜ q /m ˜ g (cid:46) m χ due tobound-state effects, as a function of m ˜ q /m ˜ g , is similar to the Bino case. As shown in Fig. 6for the Bino case, the fall in the Ω χ h to lower values of m χ is due to the breakdown of ˜ g − χ conversion. The curve hits a plateau for m ˜ q /m ˜ g (cid:38)
300 which represents the decoupling limitat m χ ∼ We now consider the case of a Higgsino LSP. The left panel of Fig. 9 displays the gluino-Higgsino coannihilation strips for Ω χ h = 0 . ± . m ˜ q /m ˜ g = 10 using the same The Sommerfeld enhancements in the calculations of the thermal relic abundance of a Wino- or Higgsino-like LSP were discussed in detail in [34]. We have included an estimate of the Sommerfeld enhancementfactor for the Wino-Wino annihilations, and we get a similar result as [16], namely m χ ∼
10 100 10005 50 500200040006000800010 000 m q (cid:142) (cid:144) m Χ m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m Χ m Χ (cid:64) G e V (cid:68) Figure 8:
As in Fig. 6, but for a Wino LSP (left panel) and a Higgsino LSP (right panel). colour codings as for the Bino and Wino cases (red with neither the Sommerfeld enhancementnor gluino bound states, orange including the QCD Sommerfeld enhancement but againno bound-state effects, black with both effects included, and purple with the bound-stateformation rate enhanced by a factor 2). In this case the black strip extends to m χ ∼ m = 0. The Higgsino couplings depend on tan β and in our calculationswe have taken tan β = 10. Our results are very weakly dependent on this choice. Once againwe see a divergence of the contours at low m χ . In this case, when m χ < ∼ . χ h at the endpoints varies with m χ , with the colours ofthe lines corresponding again to the colours of the strips in the left panel of Fig. 9. The blackline crosses the horizontal green band where Ω χ h = 0 . ± . m χ ∼ m χ are found for a range 5 (cid:46) m ˜ q /m ˜ g (cid:46) χ h contours to lower values of m χ again being due to cross sectioncancellations at low m ˜ q /m ˜ g and due to the breakdown of ˜ g − χ conversion at high m ˜ q /m ˜ g .As in the case of the Wino, the curves drop to a plateau for m ˜ q /m ˜ g (cid:38)
300 representing thedecoupling limit. In this case, the asymptotic value of m χ is ∼ . m χ due to bound-state effects is again similar to the Binocase.The decreases in the maximum values of m χ allowed in the Wino and Higgsino cases,compared to the Bino case, are due to the effect noted in [35], namely that coannihilations22
000 4000 6000 8000 10000050100150200 m Χ (cid:64) GeV (cid:68) m g (cid:142) (cid:45) m Χ (cid:64) G e V (cid:68) m q (cid:142) (cid:144) m g (cid:142) (cid:61) m Χ (cid:64) GeV (cid:68) (cid:87) Χ h m q (cid:142) (cid:144) m g (cid:142) (cid:61) Figure 9:
As in Fig. 4 and 5, but for a Higgsino LSP. may, under some circumstances, increase the relic abundance by coupling ‘parasitic’ degreesof freedom. In the Bino (Wino) (Higgsino) case, there are 2 (6) (8) electroweak degrees offreedom, linked by coannihilation to the gluinos, that contribute incrementally to the relicabundance. This effect is compensated by the decreases in the maximum values of m χ thatwe find in the Wino and Higgsino cases. We have studied in this paper MSSM scenarios in which the LSP is (almost) degeneratewith the gluino, exploring the characteristics and locating the endpoints of the gluino-LSPcoannihilation strip in the cases where the LSP is the Bino, the neutral Wino or a neutralHiggsino. Important ingredients in our analysis are the Sommerfeld enhancement of gluinoannihilation rates, gluino-gluino bound-state formation and gluino-neutralino conversion. Aswe show, these can affect significantly the preferred range of the gluino-LSP mass differencealong the coannihilation strip, and also the position of the endpoint.In the Bino LSP case, we find that at the endpoint the LSP mass ∼ ∼ ∼ ∼ ∼ < ∼ m ˜ q /m ˜ g < ∼
50, but is substantially reduced for eithersmaller or larger values of m ˜ q /m ˜ g . In all cases, the percentage increase in the allowed range23f m χ due to bound-state effects may be as large as 50%.We are loath to claim that our upper limit on the LSP is absolute, but we do note thatit is substantially higher than what is possible along the stop coannihilation strip, reflectingthe larger annihilation rates that are possible for the gluino because of its larger colourcharge. However, these annihilation rates also depend on the masses of other sparticles,notably the squarks in the gluino NLSP case studied here. As have shown, the decrease inupper limit on the LSP mass for small m ˜ q /m ˜ g is due to cancellations in the annihilationmatrix elements, whilst the decrease at large m ˜ q /m ˜ g is due to the breakdown of gluino-LSP conversion. However, we have not studied the limit m ˜ q /m ˜ g →
1, where many morecoannihilation processes would come into play, as might also be the case in non-minimalsupersymmetric models.Nevertheless, our analysis does show that a large mass reach to at least 8 TeV will benecessary to explore conclusively the possibility of supersymmetric dark matter within theMSSM and a conventional cosmological framework.
Acknowledgements
The authors would like to thank Brian Batell, Mathias Garny, Matthew McCullough andHua-Sheng Shao for very helpful discussions. The work of J.E. was supported in part by theLondon Centre for Terauniverse Studies (LCTS), using funding from the European ResearchCouncil via the Advanced Investigator Grant 267352 and from the UK STFC via the researchgrant ST/J002798/1. The work of F.L. was also supported by the European Research CouncilAdvanced Investigator Grant 267352. The work of K.A.O. was supported in part by DOEgrant de-sc0011842 at the University of Minnesota.
Appendix A
In this appendix, we first recall our procedure for computing the thermal-averaged velocity-weighted cross-sections (cid:104) σ v rel (cid:105) for the process 1 + 2 → SSARD [23]. More details of our approach can be found in [6, 36].We start with the squared transition matrix elements |T | (summed over final spins andcolors, averaged over initial spins and colors) for the coannihilation processes of interest,which are given here in Appendix B, expressed as functions of the Mandelstam variables s , t , u . We then express |T | in terms of s and the scattering angle θ CM in the centre-of-mass24rame, as described in [36]. Next we define w ( s ) ≡ (cid:90) d p (2 π ) E d p (2 π ) E (2 π ) δ ( p + p − p − p ) |T | = 132 π p ( s ) s / (cid:90) +1 − d cos θ CM |T | . (42)The total annihilation cross section σ ( s ) is given in terms of w ( s ) by σ ( s ) = w ( s ) /s / p ( s ) .The usual partial-wave expansion can be obtained by expanding |T | in powers of p ( s ) /m .The odd powers vanish upon integration over θ CM , while the zeroth- and second-order termscorrespond to the usual s and p waves, respectively. We can therefore evaluate the s - and p -wave contributions to w ( s ) simply by evaluating |T | at two different values of cos θ CM .The proper procedure for thermal averaging has been discussed in [37, 41] for the case m = m , and in [38, 40] for the case m (cid:54) = m , so we do not discuss it in detail here. Onefinds the coefficients a and b in the expansion (3) of the thermal-averaged cross-sections forthe processes of interest: (cid:104) σ v rel (cid:105) = a + b x − + O ( x − ) , (43)where x ≡ m /T (assuming m < m ) by following the prescription given in [36], using thetransition amplitudes listed in Appendix B for each final state. When the conversion ratesare large compared to the Hubble rate, these amplitudes can be used to compute the totaleffective coefficients a eff and b eff by performing the sum over initial states as in (38), and wethen integrate the rate equation (35) numerically to obtain the relic density. Appendix B
The 2 → g ˜ g → gg , ˜ g ˜ g → q A ¯ q B , ˜ gχ i → q A ¯ q B and˜ gχ + j → u A ¯ d B , where the indices A, B = 1 , , i = 1 , ..., j = 1 , |T | for each of these processes, separating the contributions from s -, t - and u -channel diagrams. In the following expressions, final spins and colors are summedover, and initial spins are averaged over. A factor c ini is used to average over initial colors.Therefore, |T | takes the form |T | = c ini ( T s ×T s + T t ×T t + T u ×T u + T s ×T t + T s ×T u + T t ×T u ) . (44) Our w ( s ) is also the same as w ( s ) in [37–39], which is written as W/
25e note that there is also the charge-conjugated process for the chargino, ˜ gχ − j → ¯ u A d B ,which we do not list separately. ˜ g ˜ g → gg There is an s -channel gluon-exchange diagram, and t - and u -channel gluino-exchange dia-grams. We note that, because there are two identical gluons in the final state, an extrafactor of 1 / c ini = 164 , T s ×T s = 1152 π α s s (cid:2) s − ( t − u ) (cid:3) , T t ×T t = − π α s s (cid:0) m g − t (cid:1) (cid:110) m g (cid:2) s ( t + 3 u ) + 2 s ( t + 2 u ) + 2( t + u ) (cid:3) + m g (cid:2) s − s ( t + 2 u ) − t + u ) (cid:3) + 2 m g [ s + 4( t + u )] − m g − tu (cid:2) s + 2 su + 2( t + u ) (cid:3) (cid:111) , T s ×T t = − π α s s (cid:0) m g − t (cid:1) (cid:2) s ( t − u ) (cid:0) m g − t + u (cid:1) + s + s ( u − t ) + ( t − u ) (cid:3) , T t ×T u = − π α s s (cid:0) m g − t (cid:1) (cid:0) m g − u (cid:1) (cid:0) m g − tu (cid:1) (cid:2) − t + u ) m g + 8 m g + ( t − u ) (cid:3) , and T s ×T u and T u ×T u are related to T s ×T t and T t ×T t , respectively, by exchanging t ↔ u in the corresponding expressions. ˜ g ˜ g → q A ¯ q B , ˜ gχ i → q A ¯ q B , ˜ gχ + j → u A ¯ d B These three processes all have t - and u -channel squark-exchange diagrams, and ˜ g ˜ g → q A ¯ q B also has an s -channel gluon-exchange diagram, whereas the other two processes do not (hence T s ×T s = T s ×T t = T s ×T u = 0 for them). Apart from the couplings, these three processes havethe same structures as for T t ×T t , T u ×T u and T t ×T u . In the case of ˜ gχ + j → u A ¯ d B , because thequark CKM matrix is involved in the chargino-quark-squark vertex, the indices A and B canbe different even if we restrict to the case of no generation mixing with only left-right mixingin the third generation for the up-type and down-type squarks. Therefore, it is convenientto write a 6 × Z ˜ U , which relates the interaction eigenstates26nd mass eigenstates of the up-type squarks as follows: ˜ u L ˜ c L ˜ t L ˜ u R ˜ c R ˜ t R = θ ˜ t − sin θ ˜ t θ ˜ t cos θ ˜ t ˜ u ˜ u ˜ c ˜ c ˜ t ˜ t , (45)where θ ˜ t is the stop left-right mixing angle. The mass eigenvalues are correspondingly definedas m ˜ U = m ˜ u , m ˜ U = m ˜ u , m ˜ U = m ˜ c , m ˜ U = m ˜ c , m ˜ U = m ˜ t and m ˜ U = m ˜ t . A similarmixing matrix, Z ˜ D , is introduced for the down-type squarks, by changing θ ˜ t to the sbottomleft-right mixing angle, θ ˜ b . The mass eigenvalues m ˜ D − are also defined similarly.For ˜ g ˜ g → q A ¯ q B we find T s ×T s = 384 π δ AB α s s (cid:2) m (cid:0) m + s − t − u (cid:1) + 2 m ( s − t − u ) + 2 m + 2 m + t + u (cid:3) , T s ×T t = (cid:88) p =1 π δ AB α s s (cid:16) t − m f p (cid:17) (cid:110) (cid:2) m (cid:0) m + s − t (cid:1) + m ( s − t ) + m + m + t (cid:3) (cid:16) | Z ˜ f ( A +3) p | + | Z ˜ fAp | (cid:17) − m m (cid:0) m + 3 m − t − u (cid:1) (cid:104) Z ˜ f ( A +3) p (cid:16) Z ˜ fAp (cid:17) ∗ + Z ˜ fAp (cid:16) Z ˜ f ( A +3) p (cid:17) ∗ (cid:105) (cid:111) , where m = m ˜ g and m = m f A . The index ˜ f = ˜ U , ˜ D , the index f = U, D . m U , , = m u,c,t , m D , , = m d,s,b , and T s ×T u is related to T s ×T t by exchanging t ↔ u .27or all three processes, T t ×T t , T u ×T u and T t ×T u take the following forms: T t ×T t = (cid:88) p,q =1 π c tt α s (cid:0) t − m t p (cid:1) (cid:0) t − m t q (cid:1) (cid:2) m m (cid:0) f tL ( A, p ) f tR ( A, q ) ∗ + f tR ( A, p ) f tL ( A, q ) ∗ (cid:1) + (cid:0) m + m − t (cid:1) (cid:0) f tL ( A, p ) f tL ( A, q ) ∗ + f tR ( A, p ) f tR ( A, q ) ∗ (cid:1)(cid:3) × (cid:2) m m (cid:0) g tL ( B, p ) g tR ( B, q ) ∗ + g tR ( B, p ) g tL ( B, q ) ∗ (cid:1) + (cid:0) m + m − t (cid:1) (cid:0) g tL ( B, p ) g tL ( B, q ) ∗ + g tR ( B, p ) g tR ( B, q ) ∗ (cid:1)(cid:3) , T t ×T u = − (cid:88) p,q =1 π c tu α s (cid:40) (cid:16) u − m u p (cid:17) (cid:0) t − m t q (cid:1) × (cid:26) (cid:16) f tL ( A, q ) ∗ f uL ( A, p ) (cid:110) m g tR ( B, q ) ∗ (cid:2) m (cid:0) m + m − t (cid:1) g uL ( B, p ) + m (cid:0) m + m − t − u (cid:1) g uR ( B, p ) (cid:3) + g tL ( B, q ) ∗ (cid:2)(cid:0) m m + m m − tu (cid:1) g uL ( B, p ) + m m (cid:0) m + m − u (cid:1) g uR ( B, p ) (cid:3) (cid:111) + m f tL ( A, q ) ∗ f uR ( A, p ) (cid:110) m g tL ( B, q ) ∗ (cid:2)(cid:0) m + m − u (cid:1) g uL ( B, p ) + 2 m m g uR ( B, p ) (cid:3) + g tR ( B, q ) ∗ (cid:2) m (cid:0) m + m − t − u (cid:1) g uL ( B, p ) + m (cid:0) m + m − t (cid:1) g uR ( B, p ) (cid:3) (cid:111)(cid:17) + ( L ↔ R ) (cid:27) + ( t ↔ u, m ↔ m ) (cid:41) , where the ( L ↔ R ) in T t × T u applies to all the L and R in the indices, and the ( t ↔ u )applies to both the t and u in the indices and the Mandelstam variables. Again, T u ×T u isrelated to T t ×T t by exchanging m ↔ m and t ↔ u in both the indices and the Mandelstamvariables.The couplings and masses involved in the above expressions are listed below.28or ˜ g ˜ g → q A ¯ q B : c ini = 164 ,c tt = c uu = 2563 ,c tu = − ,m = m = m ˜ g ,m = m = m f A ,m t p = m u p = m ˜ f p ,f tL ( A, p ) = f uL ( A, p ) = Z ˜ f ( A +3) p ,f tR ( A, p ) = f uR ( A, p ) = − Z ˜ fAp ,g tL ( B, p ) = g uL ( B, p ) = (cid:16) Z ˜ fBp (cid:17) ∗ ,g tR ( B, p ) = g uR ( B, p ) = − (cid:16) Z ˜ f ( B +3) p (cid:17) ∗ , where the index ˜ f = ˜ U , ˜ D , and the index f = U, D .29or ˜ gχ i → q A ¯ q B : c ini = 18 ,c tt = c uu = − c tu = 8 πα s ,m = m ˜ g ,m = m χ i ,m t p = m u p = m ˜ f p ,f tL ( A, p ) = Z ˜ f ( A +3) p ,f tR ( A, p ) = − Z ˜ fAp ,g tL ( B, p ) = − ig m f B c tL (cid:16) Z ˜ f ( B +3) p (cid:17) ∗ √ m w − i √ g (cid:16) Z ˜ fBp (cid:17) ∗ (cid:2) ( N i1 ) ∗ tan ( θ w ) (cid:0) Q f B − T f B (cid:1) + T f B ( N i2 ) ∗ (cid:3) ,g tR ( B, p ) = i √ g N i1 Q f B tan ( θ w ) (cid:16) Z ˜ f ( B +3) p (cid:17) ∗ − ig m f B c tR (cid:16) Z ˜ fBp (cid:17) ∗ √ m w ,f uL ( A, p ) = i √ g Q f A ( N i1 ) ∗ tan ( θ w ) Z ˜ f ( A +3) p − ig m f A c uL Z ˜ fAp √ m w ,f uR ( A, p ) = − ig m f A c uR Z ˜ f ( A +3) p √ m w − i √ g Z ˜ fAp (cid:2) N i1 tan ( θ w ) (cid:0) Q f A − T f A (cid:1) + N i2 T f A (cid:3) ,g uL ( B, p ) = (cid:16) Z ˜ fBp (cid:17) ∗ ,g uR ( B, p ) = − (cid:16) Z ˜ f ( B +3) p (cid:17) ∗ , where N is the 4 × g is the StandardModel SU(2) L coupling constant. For up-type quark final states, the index ˜ f = ˜ U , and T f A = T f B = 12 ,Q f A = Q f B = 23 ,c tL = c uL = (cid:0) c tR (cid:1) ∗ = ( c uR ) ∗ = csc( β ) ( N i4 ) ∗ ,m = m U A ,m = m U B ,m f A = m U A ,m f B = m U B . β ≡ v /v , and the vacuum expectation values of the two Higgs doublets are defined as, (cid:104) H (cid:105) ≡ (cid:18) v (cid:19) , (cid:104) H (cid:105) ≡ (cid:18) v (cid:19) . For down-type quark final states, the index ˜ f = ˜ D , and T f A = T f B = − ,Q f A = Q f B = − ,c tL = c uL = (cid:0) c tR (cid:1) ∗ = ( c uR ) ∗ = sec( β ) ( N i3 ) ∗ ,m = m D A ,m = m D B ,m f A = m D A ,m f B = m D B . gχ + j → u A ¯ d B : c ini = 18 ,c tt = c uu = − c tu = 8 πα s ,m = m ˜ g ,m = m χ + j ,m = m U A ,m = m D B ,m t p = m ˜ U p ,m u p = m ˜ D p ,f tL ( A, p ) = Z ˜ U ( A +3) p ,f tR ( A, p ) = − Z ˜ UAp ,g tL ( B, p ) = (cid:88) C =1 ig csc( β ) K CB m U C ( V j2 ) ∗ (cid:16) Z ˜ U ( C +3) p (cid:17) ∗ √ m w − ig K CB ( V j1 ) ∗ (cid:16) Z ˜ UCp (cid:17) ∗ ,g tR ( B, p ) = (cid:88) C =1 ig sec( β ) K CB U j2 m D B (cid:16) Z ˜ UCp (cid:17) ∗ √ m w ,f uL ( A, p ) = (cid:88) C =1 (cid:32) ig K AC csc( β ) m U A ( V j2 ) ∗ Z ˜ DCp √ m w (cid:33) ,f uR ( A, p ) = (cid:88) C =1 (cid:32) ig K AC sec( β ) U j2 m D C Z ˜ D ( C +3) p √ m w − ig K AC U j1 Z ˜ DCp (cid:33) ,g uL ( B, p ) = (cid:16) Z ˜ DBp (cid:17) ∗ ,g uR ( B, p ) = − (cid:16) Z ˜ D ( B +3) p (cid:17) ∗ , where the K matrix is the quark CKM matrix, and U and V are the 2 × a in Eq. (43)) for the ˜ g ˜ g → q A ¯ q B channel, in the limit of a common squark mass and massless quarks, with no generation orleft-right mixing in the squark mixing matrices (the case considered in the main body of thetext). In this limit, the contributions from each of the six quark flavor final states are thesame, and the result of putting all the six quark flavors together is a ˜ g ˜ g → q A ¯ q B limit value = 9 πα s (cid:0) m g − m q (cid:1) m g (cid:0) m g + m q (cid:1) , (46)32here m ˜ q is the common squark mass. When m ˜ q (cid:29) m ˜ g , only the s -channel gluon-exchangediagram contributes, and the above expression is proportional to m − g . On the other hand,when m ˜ q → m ˜ g , the above expression approaches zero. This cancellation of the s -, t - and u -channel contributions results in the feature of the plots at small values of m ˜ q /m ˜ g that arecommented upon in the main body of the text. References [1] G. Aad et al. [ATLAS Collaboration], arXiv:1405.7875 [hep-ex]; S. Chatrchyan et al. [CMS Collaboration], JHEP (2014) 055 [arXiv:1402.4770 [hep-ex]].[2] S. P. Martin, Phys. Rev. D , 115005 (2007) [hep-ph/0703097 [HEP-PH]]; B. Bhat-tacherjee, A. Choudhury, K. Ghosh and S. Poddar, Phys. Rev. D , no. 3, 037702(2014) [arXiv:1308.1526 [hep-ph]].[3] H. Goldberg, Phys. Rev. Lett. (1983) 1419; J. Ellis, J. Hagelin, D. Nanopoulos,K. Olive and M. Srednicki, Nucl. Phys. B (1984) 453.[4] K. Griest and D. Seckel, Phys. Rev. D , 3191 (1991).[5] C. Boehm, A. Djouadi and M. Drees, Phys. Rev. D , 035012 (2000) [arXiv:hep-ph/9911496]; J. Edsjo, M. Schelke, P. Ullio and P. Gondolo, JCAP , 001 (2003)[arXiv:hep-ph/0301106]; I. Gogoladze, S. Raza and Q. Shafi, Phys. Lett. B , 345(2012) [arXiv:1104.3566 [hep-ph]]; M. A. Ajaib, T. Li and Q. Shafi, Phys. Rev. D ,055021 (2012) [arXiv:1111.4467 [hep-ph]].[6] J. R. Ellis, K. A. Olive and Y. Santoso, Astropart. Phys. , 395 (2003) [arXiv:hep-ph/0112113].[7] J. L. Diaz-Cruz, J. R. Ellis, K. A. Olive and Y. Santoso, JHEP , 003 (2007)[arXiv:hep-ph/0701229].[8] J. Harz, B. Herrmann, M. Klasen, K. Kovarik and Q. L. Boulc’h, Phys. Rev. D (2013) 5, 054031 [arXiv:1212.5241]; J. Harz, B. Herrmann, M. Klasen and K. Kovarik,Phys. Rev. D (2015) 3, 034028 [arXiv:1409.2898 [hep-ph]].[9] J. Ellis, K. A. Olive and J. Zheng, Eur. Phys. J. C , 2947 (2014) [arXiv:1404.5571[hep-ph]]. 3310] S. Raza, Q. Shafi and C. S. n, arXiv:1412.7672 [hep-ph].[11] A. Ibarra, A. Pierce, N. R. Shah and S. Vogl, arXiv:1501.03164 [hep-ph].[12] S. Profumo and C. E. Yaguna, Phys. Rev. D , 115009 (2004) [hep-ph/0402208].[13] D. Feldman, Z. Liu and P. Nath, Phys. Rev. D , 015007 (2009) [arXiv:0905.1148[hep-ph]]; N. Chen, D. Feldman, Z. Liu, P. Nath and G. Peim, Phys. Rev. D , 035005(2011) [arXiv:1011.1246 [hep-ph]].[14] I. Gogoladze, R. Khalid and Q. Shafi, Phys. Rev. D , 115004 (2009) [arXiv:0903.5204[hep-ph]]; I. Gogoladze, R. Khalid and Q. Shafi, Phys. Rev. D , 095016 (2009)[arXiv:0908.0731 [hep-ph]]; M. Adeel Ajaib, T. Li, Q. Shafi and K. Wang, JHEP ,028 (2011) [arXiv:1011.5518 [hep-ph]].[15] K. Harigaya, M. Ibe and T. T. Yanagida, JHEP , 016 (2013) [arXiv:1310.0643[hep-ph]]; J. L. Evans and K. A. Olive, Phys. Rev. D , no. 11, 115020 (2014)[arXiv:1408.5102 [hep-ph]].[16] K. Harigaya, K. Kaneta and S. Matsumoto, Phys. Rev. D , no. 11, 115021 (2014)[arXiv:1403.0715 [hep-ph]].[17] A. De Simone, G. F. Giudice and A. Strumia, JHEP , 081 (2014) [arXiv:1402.6287[hep-ph]].[18] M. Low and L. T. Wang, JHEP , 161 (2014) [arXiv:1404.0682 [hep-ph]].[19] J. Ellis, J. Evans, F. Luo, and K. A. Olive, in preparation (2015).[20] A. Sommerfeld, Ann. Phys. , 257 (1931).[21] J. Hisano, S. Matsumoto and M. M. Nojiri, Phys. Rev. Lett. , 031303 (2004) [hep-ph/0307216]; J. Hisano, S. .Matsumoto, M. M. Nojiri and O. Saito, Phys. Rev. D ,063528 (2005) [hep-ph/0412403]; J. L. Feng, M. Kaplinghat and H. -B. Yu, Phys. Rev.D , 083525 (2010) [arXiv:1005.4678 [hep-ph]]; A. Hryczuk, Phys. Lett. B , 271(2011) [arXiv:1102.4295 [hep-ph]].[22] W. Fischler, Nucl. Phys. B , 157 (1977); Y. Schroder, Phys. Lett. B , 321 (1999)[hep-ph/9812205]; A. Strumia, Nucl. Phys. B , 308 (2009) [arXiv:0806.1630 [hep-ph]]. 3423] Information about this code is available from K. A. Olive: it contains important con-tributions from T. Falk, A. Ferstl, G. Ganis, F. Luo, A. Mustafayev, J. McDonald,K. A. Olive, P. Sandick, Y. Santoso, V. Spanos, and M. Srednicki.[24] V. B. Berestetskii, E. M. Lifshitz and L. P. Pitaevskii, Quantum Electrodynamics (Perg-amon Press, 1971).[25] M. R. Kauth, J. H. Kuhn, P. Marquard and M. Steinhauser, Nucl. Phys. B , 28(2012) [arXiv:1108.0361 [hep-ph]].[26] J. T. Goldman and H. Haber, Physica , 181 (1985).[27] K. Hagiwara and H. Yokoya, JHEP , 049 (2009) [arXiv:0909.3204 [hep-ph]].[28] M. Toharia and J. D. Wells, JHEP , 015 (2006) [hep-ph/0503175]; P. Gam-bino, G. F. Giudice and P. Slavich, Nucl. Phys. B , 35 (2005) [hep-ph/0506214];M. A. Ajaib, T. Li and Q. Shafi, Phys. Lett. B , 87 (2011) [arXiv:1107.2573 [hep-ph]].[29] D. J. H. Chung, G. R. Farrar and E. W. Kolb, Phys. Rev. D , 6096 (1997) [astro-ph/9703145].[30] W. Detmold, M. McCullough and A. Pochinsky, Phys. Rev. D , no. 11, 115013(2014) [arXiv:1406.2276 [hep-ph]]; B. von Harling and K. Petraki, JCAP , no. 12,033 (2014) [arXiv:1407.7874 [hep-ph]].[31] P. A. R. Ade et al. [Planck Collaboration], arXiv:1502.01589 [astro-ph.CO].[32] H. Baer, K. m. Cheung and J. F. Gunion, Phys. Rev. D , 075002 (1999) [hep-ph/9806361].[33] J. L. Feng, M. Kaplinghat, H. Tu and H. B. Yu, JCAP , 004 (2009)[arXiv:0905.3039 [hep-ph]].[34] J. Hisano, S. Matsumoto, M. Nagai, O. Saito and M. Senami, Phys. Lett. B , 34(2007) [hep-ph/0610249]; M. Cirelli, A. Strumia and M. Tamburini, Nucl. Phys. B ,152 (2007) [arXiv:0706.4071 [hep-ph]]; A. Hryczuk, R. Iengo and P. Ullio, JHEP ,069 (2011) [arXiv:1010.2172 [hep-ph]]; M. Beneke, C. Hellmann and P. Ruiz-Femenia,JHEP , 162 (2015) [arXiv:1411.6930 [hep-ph]].[35] S. Profumo and A. Provenza, JCAP (2006) 019 [hep-ph/0609290].3536] J. Ellis, T. Falk, K. A. Olive and M. Srednicki, Astropart. Phys.
181 (2000)[arXiv:hep-ph/9905484].[37] M. Srednicki, R. Watkins and K. A. Olive, Nucl. Phys. B , 693 (1988).[38] T. Falk, K. A. Olive and M. Srednicki, Phys. Lett. B , 248 (1994) [arXiv:hep-ph/9409270].[39] J. Ellis, T. Falk and K.A. Olive, Phys. Lett. B , 367 (1998) [arXiv:hep-ph/9810360].[40] J. Edsjo and P. Gondolo, Phys. Rev. D , 1879 (1997) [arXiv:hep-ph/9704361].[41] P. Gondolo and G. Gelmini, Nucl. Phys. B , 145 (1991).[42] J. F. Gunion and H. E. Haber, Nucl. Phys. B , 1 (1986) [Erratum-ibid. B402