Gluino Coannihilation Revisited

Some variants of the MSSM feature a strip in parameter space where the lightest neutralino is identified as the lightest supersymmetric particle (LSP), the gluino is the next-to-lightest supersymmetric particle (NLSP) and is nearly degenerate with the LSP, and the relic cold dark matter density is brought into the range allowed by astrophysics and cosmology by coannihilation with the gluino NLSP. We calculate the relic density along this gluino coannihilation strip in the MSSM, including the effects of gluino-gluino bound states and initial-state Sommerfeld enhancement, and taking into account the decoupling of the gluino and LSP densities that occurs for large values of the squark mass. We find that bound-state effects can increase the maximum LSP mass for which the relic cold dark matter density lies within the range favoured by astrophysics and cosmology by as much as ~ 50% if the squark to gluino mass ratio is 1.1, and that the LSP may weigh up to ~ 8 TeV for a wide range of the squark to gluino mass ratio \lesssim 100.


Introduction
In the absence of any signal for supersymmetry during Run 1 of the LHC [1], it is natural to ask how and where supersymmetry may be hiding. Perhaps it is hiding in plain sight with a compressed spectrum [2] that the conventional missing-energy searches at the LHC have been unable to resolve? Or perhaps R parity is violated, in which case supersymmetry may be 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 detected during Run 1 of the LHC?
There are two issues with this last possibility. One is the accentuation of the problem of the naturalness (or fine-tuning) of the electroweak scale that low-scale supersymmetry was postulated to mitigate, and the other is the cosmological cold dark matter density. The cold dark matter may well not consist only, or even predominantly, of the lightest supersymmetric particle (LSP). However, even if the cold dark matter density is considered only as an upper limit on the relic LSP density, it imposes an upper bound on the LSP mass that depends on the specific LSP candidate under consideration.
If R parity is conserved and the LSP is present in the Universe today as a relic from the Big 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 plausible candidates are the lightest neutralino χ and the gravitino [3]. Here we consider the neutralino case, 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 which it 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 rates are maximised, which happens when the LSP is (nearly) degenerate with other particles. If there is only one such coannihilating sparticle species, the coannihilation rate will in general be maximised for a coloured sparticle. There have been analyses in the literature of the cases where the coannihilating particle is a squark, specifically the lighter stop squark [5][6][7][8][9][10][11], and also the case of the gluino [10,[12][13][14][15][16][17]. In general, one would expect that the heaviest LSP will be allowed when it coannihilates with the particle with the largest colour charge, namely the gluino.
We study here the question how heavy the neutralino LSP χ could be, if it is nearly degenerate with, and coannihilates with, the gluinog. This is of relevance to assessing, for example, what centre-of-mass energy would be needed for a proton-proton collider to be 'guaranteed' to detect R-conserving supersymmetry. There can of course be no cast-iron guarantee, even within the MSSM. For example, even in the gluino coannihilation case studied here the neutralino LSP mass limit depends on the squark masses, and the LSP mass limit could be substantially modified if the squarks were degenerate with the neutralino LSP and 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-gluino coannihilation [10,[12][13][14][15][16][17], and the main new elements here are in our discussions of the effects of gluino-gluino bound states and of the issue whether coannihilations can be maintained in the presence of a large squark to gluino mass ratio. Here, we will restrict our attention to the coannihilation processes and leave their application to more complete models (with for example, radiative electroweak symmetry breaking) for future work [18]. As we discuss in detail, bound states can remove from the primordial plasma gluino pairs that may subsequently annihilate into Standard Model particles, before they can decay into the LSP as is usually assumed in discussions of coannihilation. We present numerical estimates of the bound-state production rate, and find that, for fixed sparticle masses, the relic dark matter density is substantially reduced compared with the cases where bound-state formation is neglected. Conversely, the cosmological relic density may lie within the cosmological range for substantially larger LSP masses than would have been estimated in the absence of bound-state effects: this effect is ∼ 50% for mq/mg = 1.1, falling to ∼ 20% for mq/mg ∼ 10 to 50. Another effect we discuss is that, if mq/mg > ∼ 100, the densities of neutralinos and gluinos decouple and coannihilation effects freeze-out early, leaving a significantly higher relic density, thereby reducing the possible LSP mass. There is also a reduction in the possible LSP mass for small mq/mg → 1, due to cancellations between s-, t-and u-channel diagrams that tend to reduce annihilation rates.
Taking these effects into account, we find a maximum LSP mass ∼ 8 TeV if it is the Bino, which may be attained for 10 < ∼ mq/mg < ∼ 100. If the LSP is the neutral Wino, the upper limit is reduced to ∼ 7 TeV, and for a neutral Higgsino the upper limit it becomes ∼ 6 TeV. The layout of this paper is as follows. In section 2, we review the Sommerfeld enhancement in the relevant gluino-gluino annihilation processes. In Section 3, we discuss the formation of gluino-gluino bound states, considering also dissociation processes in the early Universe. In Section 4, we consider the rates for conversion between gluinos and neutralinos. In Section 5, we present and discuss the coupled Boltzmann equations for neutralinos χ, gluinos and gluino-gluino bound states. Section 6 contains some numerical results for the gluino coannihilation strip and a discussion of its endpoint. Section 7 summarises our conclusions and discusses their significance for future colliders. Finally, Appendices present technical aspects of the computation of the 2 → 2 cross sections needed for solving the Boltzmann equations.

Sommerfeld Enhancement
Before discussing the formation and effects of gluino-gluino bound states, we first discuss briefly Sommerfeld effects in gluino-gluino annihilation, which may enhance annihilation rates at low velocities, and are particularly relevant in the case of the strongly-interacting gluino. As a general rule, initial-state interactions modify s-wave cross-sections by factors [19,20] where α is the coefficient of a Coulomb-like potential whose sign convention is such that the attractive case has α < 0, and β is the velocity of one of the annihilating particles in the centre-of-mass frame of the collision. In the cases of strongly-interacting particles, the Coulomb-like potential has the form [21] where α s is the strong coupling strength, C f is the quadratic Casimir coefficient of a specific final-state colour representation, and C i and C i are the quadratic Casimir coefficients of the annihilating coloured particles. In the case of octet annihilating particles such as gluinos, The relevant final states are in singlet, octet, or 27 s representations, for which 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 [22] for calculating the relic dark matter density. In the coannihilation region of interest, this code uses a non-relativistic expansion for annihilation cross-sections: where ... denotes an average over the thermal distributions of the annihilating particles, the coefficient 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 α < 0 in (1), modifying the leading term in (3) so that it acquires a singularity ∝ √ x. In this paper we have included these enhancements in thegg → gg andgg → qq cross sections.
The procedure for obtaining a thermally averaged cross section is given in Appendix A. The expressions for the matrix elements for the coannihilation processes are given in detail in Appendix B.

Gluino-Gluino Bound-State Formation
Gluino-neutralino coannihilations may increase the effective annihilation cross section and thereby lower the final neutralino relic abundance. The Sommerfeld enhancement discussed above further increases the cross section in specific channels and again lowers the abundance of neutralinos allowing for larger masses at the tip of the coannihilation strip defined by ∆m = 0 where ∆m is the gluino-neutralino mass difference [16]. Gluino-gluino bound states can further serve to remove gluinos from the thermal bath and thereby lower the relic density by a factor that is non-negligible relative to the Sommerfeld enhancement, and much larger than 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 isg +g ↔R + g. These processes become important when the plasma temperature falls low enough for typical thermal energies to become comparable to the binding energy of theR state, namely T E B ≡ 2mg − mR. In principle, one may form colour-octet states as well as singlets, but the latter are expected to be more deeply bound with larger wave functions at the origin. Here we focus on the production of the lightest colour-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 considering gluinos weighing several TeV, we expect the leading order of QCD perturbation theory to be a useful approximation, and assume the Coulomb potential V (r) = −3α s /r for the 1 s state, with binding energy E B (3α s /2) 2 mg. The normalised spatial part of the wave function for this 1 s bound state is where a ≡ 2/(3α s mg) is the Bohr radius. The 1 s bound state decays predominantly to a pair of gluons and the leading order decay rate is

Dissociation
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 following Section 56 of [23], where the photoelectric effect for an atom is calculated. The central part of the calculation is the evaluation of the transition amplitude given in Eq. (56.2) of [23]: where φ f is the wave function of the freeg agb pair and φ i ≡ φ bs (r), and c and k are the polarisation and momentum vectors of the gluon, respectively.
We use the dipole approximation, e i k· r ≈ 1, which is justified because the bound-state wave function φ bs (r) is exponentially suppressed for r > a, and because the gluon momentum | k| = ω, where its energy ω satisfies the conservation condition where | p| is the momentum of one of the annihilating gluinos. (Note that | p| is the same as the relative momentum, (mg/2)v rel , and the second term on the LHS of (7) can be neglected.) We find ωa E B a = (3αs) 2 1 for v rel = 0 and α s = 0.1, and more generally ωa < 1 when v rel < 0.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 agb ) = C(R)C(g c ), where the 1 s ground state with L = 0 and S = 0 has J P C = 0 −+ . The Cparity of the colour anti-symmetric 8 A state is the same as that of the gluon [24], 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 parity is conserved in this case, because P (φ f ) = 1 and the gluon has P = −1.) The normalised spatial part of the wave function for the free pairg agb is Only the L = 1 term survives, due to the selection rule from the dipole approximation. Since we wish to calculate |M f i | 2 , we may discard the phase shift factor e −iδ L (δ L is real) and the factor i L . Therefore we write so that where M f i is calculated following Section 56 of [23].
Since φ f is the wave function for an 8 A state, the Coulomb potential is V f (r) = − 3 2 α s /r, whereas φ i is a wave function for the Coulomb potential V (r) = −3α s /r, the result is different from Eq. (56.12) of [23], namely where ξ = 3 2 α s /v rel and κ = 2. This equation is averaged over the gluon polarizations and would reduce to Eq. (56.12) if κ = 1.
The total wave functions for the freeg agb pair and the bound stateR are products of the spin, colour and spatial parts of the wave functions. In view of the Majorana nature of the gluinos, the total wave functions should be anti-symmetric. Concerning the spin part of the wave function, since both the bound state and the free gluino pair are in an S = 0 state, the spin wave functions are both and the spin parts of the wave functions do not introduce any extra factor in σ dis . As for the colour part of the wave function, according to [25,26] (the latter because f abc f abd = 3δ cd ).
The (−i ∇· c (mg/2) )e i k· r factor in the transition amplitude (6) is calculated from the gluinogluino-gluon interaction Lagrangian which can be compared with the corresponding QED interaction Lagrangian We simply replace the electric charge factor Q f in the transition amplitude (6) by the colour factor f abc , since the factor 1/2 in (13) is compensated by a factor of 2 due to the Majorana nature of the gluino. Putting the above colour factors together, we obtain Note that all color indices are summed over.
Concerning the spatial part of the wave function, we need to take into account the fact that both the initial and final states contain two identical particles. In the case of the bound state, they are in the symmetric L = 0 state, and the wave function needs to be symmetrised as in Eq. (2.14) of [25]: On the other hand, the final free pair is in the antisymmetric L = 1 state, and the wave function needs to be antisymmetrised: The coefficients in these two equations introduce an extra factor of | √ 2 √ 2| 2 = 4 into the modulus-squared of the spatial wave function factors. Finally, recall that we have averaged over the polarisations of the gluon, but we also need to average over its colour. This gives a factor of 1/8. Therefore, the final dissociation cross section is where the final factor of 1/2 is to avoid double counting of gluinos in the final-state phasespace integration.

Formation
We come finally to the bound-state formation cross section, σ bsf , which is related to σ dis through the Milne relation: where the 1 2 on the LHS of the above equation is introduced to avoid double-counting the number of bound-state formation reactions, and the factor 1 e ω/T −1 comes from the enhancement of bound-state formation due to the stimulated gluon emission in the thermal background (similar to the stimulated recombination in e − p ↔ Hγ). Using and (7), we find where For comparison, the Sommerfeld enhanced s-wave cross section forg agb → g c g d is given in Eqs. (2.13) and (2.25) of [16]: where ξ = 3 2 α s /v rel . Therefore, in the v rel → 0 limit we find Therefore, we see that the inclusion ofgg bound states is a non-negligible component in determining the final neutralino relic density.

Conversion Rates
For coannihilation to be effective, the coannihilating species (in this case neutralinos and gluinos) must be in thermal contact. That is, the rates for interconverting the LSP and NLSP must be faster than the Hubble rate. In both the familiar cases of stop and stau coannihilation, connectivity of the two species can be taken for granted, as the conversion rates are mediated by light Standard Model particles and are always fast. This implies that the ratio of densities (n N LSP /n LSP ) is approximately equal to the equilibrium ratio and allows for a simplification in the Boltzmann equations. However, the interconversion of neutralinos and 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 the ratio of squark masses to the gluino mass to be less than approximately 100 as we show below. 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 ↔gq, and the gluino decays and the inverse decaysg ↔ χqq. When the neutralino is a Wino or Higgsino, the processes involving a chargino, χ + d ↔gu, χ +ū ↔gd andg ↔ χ + dū, as well as the corresponding processes for χ − , are also included. We note that q stands here for all six quark flavors, and the 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 1,2 , are both taken into account. For each relevant process, we first calculate the transition matrix element |T | 2 .
We calculate the gluino decay rates forg → χqq,g → χ + dū and its charge conjugation process. The squared transition matrix elements |T | 2 are identical to the corresponding ones for the coannihilation processes given in Appendix B, except that the expressions in Appendix B should be multiplied by a factor of 2, because the statistical factor for the initial spin averaging is 1 2 × 1 2 = 1 4 for the coannihilation processes, whereas it is 1 2 for the gluino decay processes. We note also that the definitions of the Mandelstam variables should also be changed correspondingly as follows: for the coannihilations, s = (p 1 + p 2 ) 2 , t = (p 1 − p 3 ) 2 and u = (p 1 − p 4 ) 2 , whereas for the gluino decays, s = (p 1 − p 2 ) 2 , while t and u do not change.
The gluino decay rates are then obtained by performing the standard 3-body phase space integration. The inverse-decay processes do not have to be calculated separately, because they are taken into account automatically by the Boltzmann equations given in the next section.
To calculate the conversion rates for χq →gq, χ + d →gu, χ +ū →gd and their chargeconjugated processes, we first calculate the cross sections. Again, the squared transition matrix elements |T | 2 are identical to the corresponding ones for the coannihilation processes given in Appendix B, except that the expressions in Appendix B should be multiplied by a factor of 8 3 , because the factor for the initial color averaging is 1 8 for the coannihilations, whereas it is 1 3 for the conversions. Also, compared to the coannihilation processes, the Mandelstam variables for the conversion processes are re-defined as s = ( where the upper signs in the definition of t and u apply ifq B ord B is brought into the initial state, while the lower signs apply if q A or u A is pulled over to the initial state.
For each of the quark flavors, the thermally-averaged conversion rate is obtained by integrating σ c v q over the Fermi-Dirac distribution of the quark in the initial state, where σ c is the conversion cross section for any of the relevant processes discussed above. In the 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 , where q 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 = | p q |/E q . The factors 3 and 2 in (25) count the quark color and spin degrees of freedom, respectively. Again, the inverse conversion rates are taken into account automatically by the Boltzmann equations.

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 considering three separate density components: neutralinos, gluinos and bound states.
To set up a coupled set of Boltzmann equations, it is convenient to rescale the number densities of neutralinos, gluinos and bound states by the entropy density, These are governed by the following set of coupled Boltzmann equations: where and g * s and g * are the total numbers of effectively massless degrees of freedom associated with the entropy density and the energy density, respectively, σv χχ is the relative velocity times the total cross section for the channels for χχ annihilation into Standard Model particles, and σv χg and σv gg are to be understood similarly, q Γ c and Γ g are the total conversion rate and gluino decay rate discussed in the previous section, and all possible quark and anti-quark channels for the χ are summed over, Γ R is the decay rate of theR, and σv bsf is the bound-state formation cross section times the relative velocity of the two incoming gluinos, taking into account the 1/(e ω/T − 1) enhancement factor as discussed in Section 3.2. Finally, σv gR→gg Yg has the same effect as Γ R /s, namely, it converts the bound states to gluons without altering the density of free gluinos. All the quantities bracketed by . . . are thermally averaged, and the superscript 'eq' denotes equilibrium yields.
Eq. (29) can be written in a more intuitive form: One can check that the LHS of Eq. (31) is of order -10, whereas each of the terms on the RHS of Eq. (31) are of order Hence, to a good approximation, we can set the two terms on the RHS equal to each other and solve for where Therefore, we find that Moreover, we note that σv gR→gg ng is much smaller than Γ R for x 10, due to the fact that ng decreases with the decrease of temperature while Γ R is nearly temperature independent. Since the processgR →gg is related to the bound-state formation process by crossing, σv gR→gg should be related to σv bsf by a coefficient not too much different from order 1.
If at least one of the q Γ c and Γ g is sufficiently larger than H(T ) throughout the period during which (Y χ + Yg) changes substantially, which is the case when the squark mass appearing 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 Yg/Y χ ≈ Y eq g /Y eq χ . In this case, Eq. (34) can be recast in the familiar form suitable for coannihilation calculations, and we where we have included the term x 3g * s dg * s dx which takes into account the evolution of g * s with temperature. As we will see, this approximation is valid so long as mq/mg < ∼ 20.
In Eq. (35), Y = n/s, where n is interpreted as the total number density, and Y eq = n eq /s, where n eq is the total equilibrium number density, n eq ≡ i n eq,i = n eq χ + n eq g .
The effective annihilation cross section is As one can see from Eq. (34), the expression for σv gg is the 'standard' term in the first line of (34) combined with the second line involving the bound states. We re-emphasise that this simplification requires a fast interconversion rate as discussed in the previous Section, so that we can set (Yg/Y χ ) = (Y eq g /Y eq χ ), which is true only when mq/mg < ∼ 20. For larger squark 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 solve for the relic density. All we need to do is re-define the following quantities to include the contributions from each of the χ components, χ i , neutral or charged, as where q is the same as q when χ i is a neutralino, and they are different when χ i is a chargino, and the q and q indicate all the possible quark and anti-quark channels for the conversion rates and gluino decay rates. In Eq. (39), r χ i ≡ n eq χ i /n eq χ = n χ i /n χ , where the latter '=' is guaranteed by the fast conversion and/or decay rates among the different χ i 's. For later discussion, it is useful to define an effective number of degrees of freedom for χ: where ∆ χ i ≡ (m χ i /m χ 1 − 1), and we assume χ 1 is the lightest component (i.e., the LSP) among the χ i 's. We can then write r χ i explicitly as In the limit that all the χ components have the same mass, g χ ef f = 2, 6 and 8 for Bino, Wino and Higgsino, respectively.

Numerical Results
We now present some numerical results obtained using the above formalism. Our results in this section are based on simplified supersymmetric spectra defined at the weak scale. We assume degenerate squark masses, mq and for the most part, our results do not depend on supersymmetric parameters such as µ, A 0 , and tan β 1 . We assume that the neutralino is a pure state of either a Bino, Wino, or Higgsino. Thus our free parameters are simply the neutralino mass, m χ , the gluino mass, mg and the squark masses, mq. In future work, we apply these results to more realistic CMSSM-like models (without gaugino mass universality) and pure gravity mediation models with vector-like multiplets [18].
We begin with the case in which the lightest neutralino χ is the Bino.   that the relic cold dark matter density is higher than previously: Ω χ h 2 = 0.21. This is due to the fact that at a low squark to gluino mass ratio, there is a cancellation among the t and u channel annihilations with the s channel leading to a smaller gluino annihilation cross section and hence a larger relic density. The results also change, even more significantly, for large values of mq/mg, as shown in the right panel of Fig. 2, where mq/mg = 120. In this In order to summarize the effects of both the cancellations in the annihilation cross section at low mq/mg and the decoupling of the gluino coannihilations at high mq/mg, we show in Fig. 3, the relic neutralino density as a function of mq/mg for our nominal value of m χ = 7 TeV, and ∆m ≡ mg − m χ = 0, 40, and 120 GeV (black, red, and blue lines, respectively). We see clearly the rise in Ω χ h 2 at small mq/mg as well as the very rapid rise in Ω χ h 2 at high mq/mg 100. In between there is a plateau with lower Ω χ h 2 , as exemplified by the case mq/mg = 10 shown in Fig. 1. In general, there is a shallow minimum in Ω χ h 2 around mq/mg ∼ 50 whose location depends on ∆m. The horizontal band indicates the 3-σ range for the Planck determination of the cold dark matter density of Ωh 2 = 0.1193 ± 0.0014 [27].

Bino LSP
The panels of Fig. 4  The relic cold dark matter density Ω χ h 2 as a function of mq/mg for m χ = 7 TeV and the choics ∆m ≡ mg −m χ = 0, 40, and 120 GeV (from bottom to top, black, red, and blue lines, respectively). The rise at small mq/mg is due to the cancellations between the s-, t-and u-channel diagrams for gluino pair annihilation, and the rise at large mq/mg is due to the decoupling of the gluino and neutralino densities. The horizontal band indicates the 3-σ range for the Planck determination of the cold dark matter density of Ωh 2 = 0.1193 ± 0.0014 [27].
to larger values of m χ . The effect of including bound-state effects is to suppress further the value of Ω χ h 2 for fixed model parameters, so that the corresponding black Ω χ h 2 bands in Fig. 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 than our calculations, as might arise from higher-order QCD or other effects.
The upper left panel of Fig. 4 is for the case mq/mg = 1.1, where the t and u channels partially cancel the s-channel contributions to the gluino annihilation cross section. Here we see that the black band calculated including both the Sommerfeld enhancement and gluino bound-state effects extends to m χ ∼ 6.2 to 6.4 TeV. In this case, the numerical effects of the Sommerfeld enhancement are similar to those of gluino bound-state formation, and both effects are considerably larger than the current observational uncertainties in the dark matter density represented by the breadths of the bands. The purple band, which includes an allowance of a factor 2 uncertainty in the bound-state effects, as might arise from higherorder QCD, excited bound states, etc., extends to larger m χ ∼ 7.2 to 7.5 TeV. In the case mq/mg = 10 (upper right panel of Fig. 4), the effect of bound-state formation is somewhat smaller than the Sommerfeld effect, and the black (purple) band extends to m χ ∼ 8 (9) TeV.   120 (lower right). These results are calculated without the Sommerfeld enhancement factor and gluino bound-state formation (red bands), with the Sommerfeld enhancement factor but without gluino bound-state formation (orange bands), with both the Sommerfeld enhancement factor and gluino boundstate formation (black bands), and allowing for the possibility that the bound-state formation rate is a factor 2 larger than our calculations (purple bands).
These trends are also seen in the case mq/mg = 50 (lower left panel of Fig. 4), where the black and purple bands also extend to m χ ∼ 8 (9) TeV. On the other hand, the results for mq/mg = 120 (lower right panel of Fig. 4) are quite different. The Sommerfeld effect is much larger than the bound-state effect though the latter is still slightly larger than the widths of the 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 ∼ 5.3 TeV (∼ 5.4 TeV allowing for a factor 2 uncertainty in the bound-state effects). including bound-state effects is to increase the range of m χ compatible with the measured value of Ω χ h 2 by ∼ 50% for mq/mg = 1.1, decreasing to ∼ 20% for mq/mg = 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 2 = 0.1193 ± 0.0042 (green band), as a function of mq/mg: the brown and red contours are for Ω χ h 2 = 0.05 and 0.15, respectively. The band and contours exhibit the inverse of the behaviour of the relic density seen previously in Fig. 3. The neutralino mass at low mq/mg is below the maximum value of m χ , which has a shallow maximum around mq/mg = 10 to 50, and falls sharply when mq/mg 100, reflecting the effect of a breakdown ing − χ conversion. We conclude that, within the framework studied here, m χ 8 TeV (rising to ∼ 9 TeV when allowing for a factor 2 uncertainty in the bound-state formation rate) in the Bino LSP case.

Wino LSP
We now consider the case of a Wino LSP. The left panel of Fig. 7 displays the gluino-Wino coannihilation strips for Ω χ h 2 = 0.1193 ± 0.0042 for mq/mg = 10, using the same colour codings as for the Bino case (red with neither the Sommerfeld enhancement nor gluino bound states, orange including the QCD Sommerfeld enhancement but again no boundstate effects, black with both effects included, and purple with the bound-state formation rate enhanced by a factor 2). We see that in this case the black coannihilation strip extends to m χ ∼ 7 TeV. Note that the curves appear to diverge at low m χ . The reason is that even in The left panel of Fig. 8 is the analogue of Fig. 6 for the case of a Wino LSP, with the green band corresponding to Ω χ h 2 = 0.1193 ± 0.0042 and the brown and red contours to Ω χ h 2 = 0.05 and 0.15. We see that Ω χ h 2 is within the preferred range for m χ ∼ 7 TeV over a broad range 5 mq/mg 100. The percentage increase in the allowed range of m χ due to bound-state effects, as a function of mq/mg, is similar to the Bino case. As shown in Fig. 6 for the Bino case, the fall in the Ω χ h 2 to lower values of m χ is due to the breakdown ofg − χ conversion. The curve hits a plateau for mq/mg 300 which represents the decoupling limit at m χ ∼ 2.3 TeV.

Higgsino LSP
We now consider the case of a Higgsino LSP. The left panel of Fig. 9  Hiiggsino annihilations are sufficient to reduce the relic density below the Planck density.
The right panel of Fig. 9 shows how Ω χ h 2 at the endpoints varies with m χ , with the colours of the lines corresponding again to the colours of the strips in the left panel of Fig. 9. The black line crosses the horizontal green band where Ω χ h 2 = 0.1193 ± 0.0042 for m χ ∼ 6 TeV. As seen in the right panel of Fig. 8, similar values of m χ are found for a range 5 mq/mg 100, with the drops in the Ω χ h 2 contours to lower values of m χ again being due to cross section cancellations at low mq/mg and due to the breakdown ofg − χ conversion at high mq/mg.
As in the case of the Wino, the curves drop to a plateau for mq/mg 300 representing the decoupling limit. In this case, the asymptotic value of m χ is ∼ 1.2 TeV. The percentage increase in the allowed range of m χ due to bound-state effects is again similar to the Bino case.
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 [28], namely that coannihilations may, under some circumstances, increase the relic abundance by coupling 'parasitic' degrees of freedom. In the Bino (Wino) (Higgsino) case, there are 2 (6) (8) electroweak degrees of freedom, linked by coannihilation to the gluinos, that contribute incrementally to the relic abundance. This effect is compensated by the decreases in the maximum values of m χ that we find in the Wino and Higgsino cases.

Summary
We have studied in this paper MSSM scenarios in which the LSP is (almost) degenerate with the gluino, exploring the characteristics and locating the endpoints of the gluino-LSP coannihilation strip in the cases where the LSP is the Bino, the neutral Wino or a neutral Higgsino. Important ingredients in our analysis are the Sommerfeld enhancement of gluino annihilation rates, gluino-gluino bound-state formation and gluino-neutralino conversion. As we show, these can affect significantly the preferred range of the gluino-LSP mass difference along the coannihilation strip, and also the position of the endpoint.
In the Bino LSP case, we find that at the endpoint the LSP mass ∼ 8 TeV, increasing to ∼ 9 TeV if we allow for a factor 2 increase in the bound-state formation rate above our calculations. These values are decreased by ∼ 1 TeV if the LSP is a Wino, and by a further ∼ 1 TeV if it is a neutral Higgsino. The upper limit on the LSP mass of ∼ 8 TeV is weakly sensitive to the squark mass for 10 < ∼ mq/mg < ∼ 50, but is substantially reduced for either smaller or larger values of mq/mg. In all cases, the percentage increase in the allowed range of 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 that it is substantially higher than what is possible along the stop coannihilation strip, reflecting the larger annihilation rates that are possible for the gluino because of its larger colour charge. 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 in upper limit on the LSP mass for small mq/mg is due to cancellations in the annihilation matrix elements, whilst the decrease at large mq/mg is due to the breakdown of gluino-LSP conversion. However, we have not studied the limit mq/mg → 1, where many more coannihilation processes would come into play, as might also be the case in non-minimal supersymmetric models.
Nevertheless, our analysis does show that a large mass reach to at least 8 TeV will be necessary to explore conclusively the possibility of supersymmetric dark matter within the MSSM and a conventional cosmological framework.
The usual partial-wave expansion can be obtained by expanding |T | 2 in powers of p 1 (s)/m 1 . The odd powers vanish upon integration over θ CM , while the zeroth-and second-order terms correspond 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 | 2 at two different values of cos θ CM .
The proper procedure for thermal averaging has been discussed in [30,34] for the case m 1 = m 2 , and in [31,33] for the case m 1 = m 2 , so we do not discuss it in detail here. One finds the coefficients a and b in the expansion (3) of the thermal-averaged cross-sections for the processes of interest: where x ≡ m 1 /T (assuming m 1 < m 2 ) by following the prescription given in [29], using the transition amplitudes listed in Appendix B for each final state. When the conversion rates are large compared to the Hubble rate, these amplitudes can be used to compute the total effective coefficients a eff and b eff by performing the sum over initial states as in (38), and we then integrate the rate equation (35) numerically to obtain the relic density.
Here we list the |T | 2 for each of these processes, separating the contributions from s-, t-and u-channel diagrams. In the following expressions, final spins and colors are summed over, and initial spins are averaged over. A factor c ini is used to average over initial colors. Therefore, |T | 2 takes the form We note that there is also the charge-conjugated process for the chargino,gχ − j →ū A d B , which we do not list separately.

gg → gg
There is an s-channel gluon-exchange diagram, and t-and u-channel gluino-exchange diagrams. We note that, because there are two identical gluons in the final state, an extra factor of 1/2 is needed when performing the momentum integration in (42).
We find 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.
gg → q AqB ,gχ 0 i → q AqB ,gχ + j → u AdB These three processes all have t-and u-channel squark-exchange diagrams, andgg → q AqB 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 have the same structures as for T t ×T t , T u ×T u and T t ×T u . In the case ofgχ + j → u AdB , because the quark CKM matrix is involved in the chargino-quark-squark vertex, the indices A and B can be different even if we restrict to the case of no generation mixing with only left-right mixing in the third generation for the up-type and down-type squarks. Therefore, it is convenient to write a 6 × 6 up-type squark mixing matrix, ZŨ , which relates the interaction eigenstates and mass eigenstates of the up-type squarks as follows: where θt is the stop left-right mixing angle. The mass eigenvalues are correspondingly defined as mŨ 1 = mũ 1 , mŨ 2 = mũ 2 , mŨ 3 = mc 1 , mŨ 4 = mc 2 , mŨ 5 = mt 1 and mŨ 6 = mt 2 . A similar mixing matrix, ZD, is introduced for the down-type squarks, by changing θt to the sbottom left-right mixing angle, θb. The mass eigenvalues mD 1−6 are also defined similarly. Forgg → q AqB we find where m 1 = mg and m 3 = m f A . The indexf =Ũ ,D, the index f = U, D. m U 1,2,3 = m u,c,t , m D 1,2,3 = m d,s,b , and T s ×T u is related to T s ×T t by exchanging t ↔ u.
For all three processes, T t ×T t , T u ×T u and T t ×T u take the following forms: where N is the 4 × 4 neutralino mixing matrix as defined in [35], and g 2 is the Standard Model SU(2) L coupling constant. For up-type quark final states, the indexf =Ũ , and tan β ≡ v 2 /v 1 , and the vacuum expectation values of the two Higgs doublets are defined as, For down-type quark final states, the indexf =D, and where the K matrix is the quark CKM matrix, and U and V are the 2×2 chargino mixing matrices as defined in [35].
Finally, we give the s-wave result (i.e., the coefficient a in Eq. (43)) for thegg → q AqB channel, in the limit of a common squark mass and massless quarks, with no generation or left-right mixing in the squark mixing matrices (the case considered in the main body of the text). In this limit, the contributions from each of the six quark flavor final states are the same, and the result of putting all the six quark flavors together is agg →q AqB limit value = 9πα 2 s m 2 g − m 2 q 2 8m 2 g m 2 g + m 2 where mq is the common squark mass. When mq mg, only the s-channel gluon-exchange diagram contributes, and the above expression is proportional to m −2 g . On the other hand, when mq → mg, 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 mq/mg that are commented upon in the main body of the text.