Cosmological Abundance of Colored Relics

The relic cosmological abundance of stable or long-lived neutral colored particles gets reduced by about a few orders of magnitude by annihilations that occur after QCD confinement. We compute the abundance and the cosmological bounds on relic gluinos. The same post-confinement effect strongly enhances co-annihilations with a lighter Dark Matter particle, provided that their mass difference is below a few GeV. Charged colored particles (such as stops) can instead form baryons, which can be (quasi)stable in some models.


Introduction
Extensions of the Standard Model (SM) sometimes predict (quasi)stable colored particles. We show that, due to non perturbative QCD effects, their relic abundance is significantly lower than previously expected, correspondingly reducing the phenomenological constraints.
Weak-scale supersymmetry has been considered the most motivated extension of the SM, as it allows to control quadratically divergent corrections to the Higgs mass keeping them naturally small; improves the prediction for the gauge couplings in SU(5) unification; provides Dark Matter (DM) candidates. The plausibility of the naturalness goal is now endangered by the lack of any new physics in LEP [1] and LHC data [2]. Furthermore the Higgs mass is larger than what predicted by the MSSM with weak-scale sparticles.
Split SuperSymmetry [3,4] (where the new supersymmetric fermions are much lighter than the new supersymmetric scalars) abandoned the naturalness goal, retaining the two other good features, allowing to fit the Higgs mass [5,6], and relaxing the possible supersymmetric flavour problem caused by weak-scale sfermions. If sfermions are very heavy the light gauginos can become long-lived, giving peculiar signatures at colliders and potential cosmological problems. These were explored in [7], where the relic gluino abundance (before late gluino decay in neutralino and colored SM particles) was computed including perturbative gluino annihilations at T ∼ M 3 and arguing that one can neglect non-perturbative effects arising after confinement at T ∼ Λ QCD . Such effects reduce the relic gluino abundance by a few orders of magnitude [8], thereby weakening cosmological bounds.
The relevance of confinement effects has been estimated in [9] in the case of colored charged particles. Unlike in the case of the neutral gluino, QCD bound states of charged particles can be formed or broken by emitting or absorbing photons. We will consider the case of (quasi)stable stopt.
In section 2 we compute the thermal relic abundance of (quasi)stable gluinos and in section 3 we reconsider the cosmological bounds and discuss the associated phenomenology. Conclusions are given in section 4.

Relic gluinos
We consider a Majorana fermion in the adjoint of SU (3). In supersymmetric models this is known as gluino and denoted asg. The gluino can be stable if it is the lightest supersymmetric particle. Otherwise it can decay via squark exchange into a quark, an antiquark and a neutralino or chargino, or radiatively to a gluon and a neutralino, with quarks and squarks in the loop. The resulting lifetime is long if sfermions have a much heavier mass m SUSY [10,11]: where N is an order-one function [11]. A stable or long lived gluino is probed and constrained by cosmology. 1. at tree level in the perturbative expansion;

Computing the relic gluino abundance
2. taking into account Sommerfeld corrections the s-wave annihilation cross-section is (eq. (2.24) of [12], where the Sommerfeld S factors are defined) 3. taking into account also a related effect: formation of bound states [13].
These effects reduce by about 1 order of magnitude the gluino abundance, controlled by the Boltzmann equation Figure 1: Predicted gluino abundance. Relic stable gluinos exceed the DM density if Mg > ∼ PeV. The bands show the non-perturbative analytic result for σ QCD = 1/Λ 2 QCD (blue) and σ QCD = 4π/Λ 2 QCD (red). The thin (thick) lines assume that only singlet bound states (octet bound states too) can form with QCD size; similarly, the small (large) dots show our numerical computation for some values of the gluino mass.
where z = Mg/T , Yg = ng/s, s is the entropy density at temperature T ; H(T ) is the Hubble constant.
If τg < M Pl /Λ 2 QCD ∼ µsec gluinos decay before the QCD phase transition leaving no cosmological effects. Otherwise gluinos recouple as the temperature approaches the QCD scale, and their relic abundance is determined by a re-decoupling at temperatures mildly below the QCD phase transition. At this point gluinos have formedgg and/orgqq hadrons which scatter with large cross sections σ QCD = c/Λ 2 QCD where c ∼ 1, making about M Pl /Λ QCD ∼ 10 19 scatterings in a Hubble time. For comparison, the proton-proton elastic scattering cross section at low energy is known to be σ el ≈ 100 mb, corresponding to c ≈ 23.
Although gluinos are much rarer than gluons and quarks, occasionally, two gluino hadrons meet forming agg bound state. Classically such state has angular momentum ≈ µv rel b where b ≈ 1/Λ QCD is the impact parameter; µ Mg/2 is the reduced mass; v rel ∼ (T /Mg) 1/2 is the relative velocity. Thereby ∼ (MgT ) 1/2 /Λ QCD , is large for Mg Λ QCD > ∼ T . The quantummechanical total QCD cross section for forminggg bound states is large because many partial waves contribute. This can be parameterized defining the maximal angular momentum as max ≡ c/2πMgv rel /Λ QCD obtaining (see e.g. [14]) where the phase shifts average to sin 2 δ 1/2. This expectation is consistent with numerical results in toy calculable models [15].
The cross section relevant for reducing the gluino abundance is not σ QCD , but the smaller cross section σ ann for forminggg states which annihilate into SM particles before being broken. Assuming that agg with angular momentum and energy ∼ T annihilates before being broken with probability ℘ (T ), one has 1 A large cross section needs large , but ℘ can be small at large . We compute ℘ as the probability that thegg bound state radiates an energy big enough to become unbreakable (bigger than ≈ T ) before the next collision, after a time ∆t ∼ 1/n π v π σ QCD . In such a case it becomes unbreakable and keeps radiating untilgg annihilate.
The key quantity to be computed is thereby the power radiated by the relevant bound states which have n, 1. In the abelian case, this is well approximated by its classical limit: Larmor radiation. Having assumed neutral constituents, we can neglect photon radiation. Similarly, gravitational radiation has cosmologically negligible rates Γ grav ∼ E 3 B /M 2 Pl . The dominant radiation mechanism is gluon radiation, which differs from abelian radiation because gluons are charged under QCD. This makes a difference when (as in our case) particles are accelerated because of the strong force itself. While a photon can be soft and its emission leaves the bound state roughly unchanged, an emitted gluon has its own QCD potential energy, and its emission changes the QCD potential among gluinos by an order one amount (in particular, a singlet bound state becomes octet). As the classical limit of gluon emission is not known, we apply the quantum formulae. 1 This intuitive picture can be formally justified writing a network of Boltzmann equations, one for each bound state I with different and n. Such equations contain the formation rates γ I , the thermally averaged breaking rates Γ break I , the annihilation rates Γ ann I , the decay rates among the states Γ IJ . This is unpractical, given that hundreds of states play a relevant role. To get some understanding, we consider a toy system where only one state 1 can be produced, and only one state 3 can annihilate. The state 1 can decay to state 2, which can decay to state 3. Then, assuming that the rates are faster than the Hubble rate, one can reduce the network of Boltzmann equations [13] to the single Boltzmann equation eq. (3) for the total gluino density, controlled by an effective annihilation rate equal to ℘γ 1 where where the last term takes into account that 2 can upscatter to 1. We see that ℘ does not depend on Γ ann 3 and has the expected physical meaning. In view of QCD uncertainties we cannot compute all order unity factors, such that it is appropriate to employ the simpler intuitive picture. We need to compute the power radiated by highly excited bound states, with sizes of order 1/Λ QCD . Smaller bound states can be approximated by the Coulomb-like non-relativistic limit of the QCD potential, and can have various color configurations, in particular singlets and octets. At large distances, they appear as color singlets because they are surrounded by a soft gluon cloud at distance of order 1/Λ QCD , which acts as a spectator when computing their inner behaviour. In the opposite limit, states larger than 1/Λ QCD can only be color-singlet hadrons.
For our purpose what is needed are QCD-size bound states which are the most challenging, as confinement effects are starting to be relevant. We will estimate their effect into two opposite limits: 8) assuming that color octet bound states are relevant, such that radiation is dominated by single-gluon emission (pion emission after hadronization) into singlet states. This is computed in section 2.2.
1) assuming that only color singlets exists, such that radiation is dominated by color-singlet double-gluon emission (pion emission after hadronization) among singlets. This is computed in section 2.3.
While the two cases are analytically very different (e.g. different powers of the strong coupling), QCD is relatively strongly coupled so that the numerical final results in the two limiting cases will be similar. Before starting the computations, we summarize generic results for QCD bound states.

The bound states
We compute the energy levels of thegg bound states assuming the non-relativistic QCD potential where λ = (C R + C R − C Q )/2 for the potential among representation R and R in the Q configuration with C 1 = 0, C 3 = 4/3, C 8 = 3 being the Casimirs. So λ = 3 (3/2) for the potential among octets in the singlet (octet) configuration. Lattice simulations indicate α 3lattice ≈ 0.3 and σ ≈ (0.4 GeV) 2 . The one-loop correction to the perturbative term means that the QCD potential is roughly given by the tree level potential with the strong coupling renormalised at the RGE scaleμ ≈ 1/r. At finite temperature σ(T ) ≈ σ(0) 1 − T 2 /T 2 QCD with T QCD ≈ 170 MeV [17]. The product of two color octets decomposes as such that there are three attractive channels and the gluino bound states exist in the following configurations The energy eigenvalues in a potential V = −α eff /r + σ eff r are [18] E n ≈ µα 2 eff 2 − 1 n 2 t + 12tnεx −µα 2 eff /2n 2 Coulomb limit 3(xσ eff ) 2/3 /2µ 1/3 string limit (10) where µ ≈ Mg/2 is the reduced mass, = {0, 1, . . .} is angular momentum, n ≥ 1 + , x = 1.79(n − ) + − 0.42, ε = σ eff /4α 3 eff µ 2 is a dimension-less number and t is the positive solution to t = 1 − 4n 3 εxt 3 . In the limit where the Coulomb force dominates one has t 1 and ε 0; bound states have size n 2 a 0 where a 0 = 1/µα eff is the Bohr radius. The linear force dominates when n 2 a 0 α eff /σ ∼ 1/Λ QCD . Fig. 2 shows the energy levels with nearly zero energy for Mg = 3 TeV.

The breaking rate
The probabilities ℘ that a given state radiates enough energy before being broken by a collision can be computed in two different ways.
Based on classical intuition, one can simply compare its energy loss rate with the breaking rate. While this simplification holds in the abelian case, we have to deal with a non-abelian dynamics, where gluon emission changes singlet to octet states, and vice versa. This is relevant, as singlet and octet decay rates are significantly different (especially for some singlet states which only decay through higher-order effects, as discussed below). It's not clear what is the classical limit of this system in the limit of large quantum numbers n, .
We then perform a quantum computation, determining the ℘ by simulating transitions among the many different states. This is feasible up to masses Mg ∼ 10 TeV, because it involves a growing number of states at larger Mg.
We then need the breaking rate of the individual bound states. Thermal equilibrium between direct and inverse process (also known as Milne relation) does not allow to infer the breaking rates from the total creation rate, because the latter is cumulative over all bound states. We assume that the breaking rate is given by the thermal average of the pion scattering cross section, assumed to be equal to 1/Λ 2 QCD , and perform the thermal average σ break v rel over the distribution of pions with energies large enough to break the bound states. The number density of pions with enough energy to break a bound state with binding energy E B is such that Γ break ≈ σ break v rel n eq π (E π > E B ).

Color octet states and single gluon emission
We here assume that two collidingg can form agg system with all 64 possible color configurations of eq. (8), and with relative weights determined by combinatorics rather than by energetics. Then the effective annihilation cross section is determined summing over attractive channels as σ ann ∝ 1 64 We fix the proportionality factor to ≈ 4 such that the total cross section is σ QCD = c/Λ 2 QCD , where c ∼ 1 parameterizes our ignorance of the overall QCD cross section. The annihilation cross section is dominated by σ 8 A ann because the state 8 A radiates much more than 1 or 8 S . Indeed, because of selection rules, single-gluon emission allows the following decays with ∆ = ±1: Taking hadronization into account two pions are emitted, such that the binding energy of the final state E B must be larger than E B + 2m π , otherwise the decay is kinematically blocked. If the energy gap is somehow bigger than Λ QCD , inclusive decay rates can be reliably computed treating the gluon as a parton.
Since the 1 state is more attractive than 8 S,A , the above conditions are easily satisfied for the 8 A → 1 decay, while 1 → 8 A decays are kinematically blocked at larger and allowed at small enough (elliptic enough classical orbit), but suppressed with respect to the abelian result.
In our numerical results we sum over all possible final states using wave-functions computed in WKB approximation using the Langer transformation. We also provide a simple approximated analytic result obtained assuming Coulombian wave-functions (which is valid for deep final states, but not for the QCD-size initial states) 2 The decay rate must be compared with the thermal breaking rate, which is given by pion scatterings such as (gg) + π → (gg) + (gg) + π. Since we considered bound states made of neutral gluinos, they are not broken by photon scatterings to leading order. The result is very simple: the 8 A decay rate is so fast that its actual value is irrelevant: all 8 A allowed states have ℘ = 1 at the relevant temperatures T < ∼ Λ QCD . On the other hand, 8 S and 1 states contribute negligibly. Then, the annihilation rate is controlled by a much simpler condition: 8 A bound states with binding energy E B ∼ T only exist up to some maximal ≤ max8 , which can be easily computed. For Mg = 3 TeV fig. 2 shows that max8 ≈ 25. For generic Mg T , max8 is well approximated by imposing the vanishing of E n in eq. (10), finding having approximated t ≈ 1 in the last expression. Using eq. (10), the deepest available singlet state has energy gap ∆E = 9 4 √ 3α 3 σ ≈ 0.9 GeV (see also fig. 2) and can only decay via higher order processes.
The effective annihilation cross section is At low (high) temperatures one has cr max ∝ v rel ( cr max8 ∝ v 0 rel ) such that the thermal average for 1 is σ ann v rel 2σ QCD T /πMg ( σ ann v rel 3πα 3 3 /16MgT σ). Taking the minimum of these two limits (which are equal at T = T cr = π 3α 3 3 /σ/8σ QCD with σ QCD = c/Λ 2 QCD ), we obtain an approximation valid at a generic intermediate T : The Boltzmann equation of eq. (3) is approximatively solved by 2 In the same approximation, the smaller energy radiated into 8 S is given by a Larmor-like formula, given that the initial and final state are equally attractive. Figure 3: The effective annihilation cross section of gluinogg bound states, assuming that they form color-octet 8 A states (left) or only color-singlet states (right). The solid curves are the numerical computation the dashed lines are the maximal geometrical cross sections given by the analytic approximation.
where the dz integral is dominated by T ∼ T QCD : for T cr T QCD the abundance simplifies to The final relic abundance does not have a strong dependence on σ QCD , as it is only relevant at relatively low temperatures. The DM critical density is exceeded if Mg > ∼ PeV. Fig. 3a shows the full numerical result for σ ann v rel , which agrees with the analytic maximal value (apart from some smoothing at T ∼ T cr ) up to about 50 MeV: thereby the numerical abundance is better reproduced lowering T QCD down to 50 MeV in eq. (19). This is done in the analytic estimate plotted in fig. 1.

Color-singlet states and two gluon emission
Single-gluon emission switches the color of the bound state as 1 ↔ 8 and its angular momentum by ±1: as a consequence kinematics blocks single-gluon decays of various color-singlet bound states, roughly all the ones in fig. 2 which don't have nearby octet states. In particular, decays of singlet states with maximal are blocked, and octet states with maximal can (but need not) decay to singlets with blocked decays.
We thereby take into account two-gluon emission, which allows for 1 → 1 decays with ∆ = {0, ±2}. The rates of 2g transitions are mildly suppressed by O(α 3 3 ) compared with the 1g decay rates. If the energy difference ∆E is much bigger than Λ QCD , gluon hadronization proceeds with unit probability and the 2g decay widths can be computed using 2nd order non-relativistic perturbation theory [19]: Hadronization is possible down to the kinematical limit ∆E ≈ 2m π . However the energy difference between two singlet states with maximal , |∆ | = 2 and nearby n is ∼ σ 3/4 α which, in view of the Mg suppression, can be smaller than 2m π . In such a case the decay can still proceed through off-shell pions, which produce photons and leptons. We estimate these suppressed decays following section 5.6 of [20]. We neglect multi-gluon emission, which allows bigger jumps in .
The 2g rates are included in numerical computations which assume that QCD-scale color octets exist. The result was discussed in the previous sub-section, as 2g decays give a relatively minor correction.
We consider the opposite extreme possibility that octet states with QCD-size do not exist, and that only color singlets exist. We can again obtain an analytic lower bound on the final g abundance by assuming that all singlet levels fall fast. Then the cross section σ ann ≈ σ 1 ann is only limited by max1 = √ 2 max8 such that where now T cr = π 3α 3 3 /σ/4σ QCD . The resulting relic gluino abundance is 2 times lower than in eq. (18), and with the new value of T cr . Fig. 3b shows that this limit only holds at T < ∼ 20 MeV, such that the analytic expression reproduces the numerical value for Yg by reducing T QCD down to ∼ 20 MeV. Figure 4: Cosmological constraints on long-lived gluinos. Left: As a function of the gluino lifetime. Right: As a function of the sfermion mass scale m SUSY , which in Split SuperSymmetry determines the gluino lifetime.

Cosmological bounds and signatures
Bounds on quasi-stable relics depend on their lifetime τg; on their mass Mg; on their relic abundance, that for gluinos we computed in terms of Mg, and on their decay modes. As mentioned above, we assume that gluinos decay to neutralinos (assumed to be the Lightest Super-symmetric Particle, LSP) plus either a gluon or a quark and an antiquark. Here we assume that half of gluino energy is carried away by the LSP; if the LSP is not much lighter than the gluino, even less energy goes into SM states and one would obtain weaker bounds. Our final result is plotted in fig. 4, using the thick red dashed line of fig. 1: even using updated experimental bounds (discussed below), our bounds on a (quasi)stable gluino are significantly weaker than those derived in [7]. The reason is that our relic density takes into account non-perturbative gluino annihilations, and is much smaller than the 'perturbative' gluino relic density assumed in [7], see fig. 1. In particular, we find that a (quasi)stable gluino just above present collider bounds is still allowed provided that its lifetime is smaller than about 10 12 s or larger than about 10 22 s.
In the rest of this section we summarize the various bounds on decaying relics plotted in fig. 4, moving from smaller to larger lifetimes.

Big Bang Nucleosynthesis
A gluino that decays during BBN can disturb the successful BBN predictions of light element abundances, which get affected in different ways, depending on the gluino lifetime (for more details see [21,22]): 3 • For 0.1 s τg 10 2 s the mesons and nucleons produced by gluino decays quickly reach kinetic equilibrium with the thermal bath of background photons and e ± and thus do not have enough energy to destroy light nuclei. However, the extra pions, kaons and nucleons present in the thermal bath increase the p ↔ n conversion rate, thus increasing the n/p ratio and as a consequence the primordial 4 He mass fraction Y p .
• For τg 10 2 s the gluino decay products do not thermalize before interacting with nuclei, due to the lower temperature of the plasma at these times. The still energetic nucleons (the mesons decay before they can interact) can thus hadrodissociate 4 He which in turn also increases the D abundance (e.g. via p + 4 He → D + 3 He).
• For τg 10 7 s photodissociation of 4 He, which induces increased 3 He and D abundances, becomes relevant. Photodissociation is not relevant at earlier times because the γ-spectrum is cut off at the threshold energy E γ th ≈ m 2 e /(22 T ) [25] for e + e − pair production from energetic γ's with thermal γ's, so that photons are not energetic enough to break up nuclei.
The resulting constraints have been computed in [21] and updated and improved in [22]. The constraints are given in the (τ X , ξ X ) plane for different main decay modes of X, where X is the unstable relic (the gluino in our case) and ξ X = E vis Y X is its destructive power. Since we assume that half of gluinos' energy is carried away by the LSP we have E vis ≈ Mg/2. The bounds for the various hadronic decay modes are similar since in all cases they induce hadronic showers, and our bounds are based on the plot for the tt mode.
The effects from photodissociation depend only on the total injected energy, so that for τg 10 7 s the bounds do not explicitly depend on Mg to a good approximation. At earlier times, the effects depend on the number of hadrons produced in the hadronization process, which scales with a power of Mg. Thus we fit the bounds, given in [22] for M X = 1 TeV, 10 TeV, 10 2 TeV, 10 3 TeV, to a power-law function of Mg.
The left-handed panel of fig. 4 shows the resulting bounds in green. In the right-hand panel we show the same bounds with the gluino lifetime computed as function of the SUSY breaking scale m SUSY .

Distortion of the CMB blackbody spectrum
Gluinos with lifetimes between ∼ 10 7 s and ∼ 10 13 s (the latter corresponds to recombination) can lead to deviations of the CMB spectrum from a blackbody form. When the Universe is 10 7 s old, photon number changing processes such as double Compton scattering are not efficient any more, so that photons injected into the plasma can induce a chemical potential µ 1.41 δ / [26] in the Bose-Einstein distribution of the CMB radiation, where [27] δ 4 × 10 −3 τg 10 6 s MgYgB γ 10 −9 GeV exp − 6.1 × 10 6 s τg After ∼ 4 × 10 11 Ω b h 2 s [27], elastic Compton scatterings do not maintain thermal equilibrium anymore. An injection of photons 'Comptonizes' the spectrum, i.e. it leads to a mixture of blackbody spectra of different temperatures. This is described by the Compton y-parameter, given by y = δ /4 [26]. The 95% CL limits from the FIRAS instrument on the COBE satellite are |µ| < 9 × 10 −5 and |y| < 1.5 × 10 −5 [28,29]. The resulting constraints on the gluino lifetime are shown in pink in fig. 4. Here we assumed that ∼ 45% (see e.g. [30]) of the energy that is not carried away by the LSP goes into photons. The resulting bounds are less constraining than the BBN bounds. However future bounds from PIXIE [31] will be stronger by 2 to 3 orders of magnitude.

CMB anisotropies
The electromagnetic energy ejected into the gas at or after recombination by decaying relics modifies the fraction of free electrons and heats the intergalactic medium. This leads to modifications of the CMB angular power spectrum, measured by Planck. The maximally allowed density of a long-lived relic as a function of its lifetime has been computed assuming decay products with fixed energies in the range from 10 keV up to 10 TeV [32] respectively 1 TeV [33]. The e + , e − , γ from hadronic decays do not have fixed energies, and moreover we do not know the energy spectrum of the decay products of relics with a mass significantly larger than 10 TeV. For very large gluino masses the bounds we show are therefore only indicative. We consider the middle of the band in [33] and obtain bounds by assuming that half of gluinos energy goes into SM states and that 60% (see e.g. [30]) of the latter goes into e + , e − , γ. In fig. 4 we show the resulting constraints for a gluino with a lifetime 10 12 s in yellow.

21-cm line
If confirmed, the observation of an absorption feature in the low energy tail of the CMB spectrum [34] allows us to put an upper bound on the temperature of the intergalactic medium (IGM) at redshift z ≈ 17. Decays of relic particles during the dark ages are constrained, mainly because they inject energy in the IGM heating it, erasing the absorption feature. Bounds on decaying DM particles, with masses up to 10 TeV, have been computed in [35][36][37]. We rescale these bounds to a generic abundance, still assuming that half of gluino energy goes into SM states and that 60% (see e.g. [30]) of the latter goes into e + , e − , γ. The result is shown in fig. 4. Similarly to the case of the CMB bounds in the previous section, the 21 cm bounds for very large gluino masses are only indicative and subject to significant uncertainty.

Constraints from gamma-ray telescopes and neutrino detectors
Decaying gluinos with larger lifetimes are constrained by the measurement of cosmic ray spectra, in particular of photons of neutrinos. We adopt the results of [38] who computed limits on the lifetime of DM decaying to bb, from data from the Fermi gamma ray telescope and the neutrino detector IceCube, up to a DM mass of 10 12 GeV. We rescale the bounds of [38] taking into account that the density of our relics differs from the DM density. Ref. [38] derives bounds assuming a relic that decays to bb. We assume that 50% of the gluino's energy goes to the LSP and the rest goes into hadronic decay channels, which lead to similar spectra as bb. Fig. 4 shows the resulting constraints on a long-lived gluino from Fermi (in blue) and IceCube (in orange). The IceCube limits exceed the bounds from Fermi data for Mg 10 7 GeV.

Searches for super-massive nuclei
Coming finally to stable gluinos, lattice simulations indicate that they would form neutralgg hadrons [39], as well as a minor component of baryonic states such asguud (according to [40] the lightest gluino baryon could beguds). They behave as strongly interacting Dark Matter. This is allowed by direct detection experiments performed in the upper atmosphere and by searches for super-massive nuclei in the Earth and in meteorites if their relic abundance is a few orders of magnitude smaller than the cosmological DM abundance, although the precise bound is subject to considerable uncertainties (see the discussion in [8]). In fig. 4 we indicate the tentative constraints that arise from the search for supermassive nuclei in meteorites by Rutherford backscattering of 238 U, N SIMP /N n | meteorites 2 × 10 −12 [41], assuming a heavy nuclei capture cross section of σ capture = 10 −2 /Λ 2 QCD . Presumably, there is still an open window, from TeV masses above the LHC [9] up to about 10 TeV.

Higgs mass
In the right panel of fig. 4 we considered Split SuperSymmetry, such that the gluino lifetime is computed as function of the sfermion mass m SUSY , see eq. (1). This scale is further constrained within the split MSSM by the observed Higgs mass, which is reproduced within the green region (for different values of tan β) in the (M 3 , m SUSY ) plane. We computed M h as in [6], assuming that gauginos and Higgsinos are degenerate at the gluino mass M 3 and that all scalars are degenerate at m SUSY . Allowing the masses to vary and taking into account uncertainties on M t and α 3 slightly expands the region. Within the Higgs-allowed region the gluino decays promptly on cosmological time-scales, evading all cosmological bounds.
No prediction for the Higgs mass arises in extensions of the MSSM. However, roughly the same region is obtained imposing the meta-stability bound on Higgs vacuum decay, which implies that the Higgs quartic λ H cannot be too negative, λ H > ∼ − 0.05. A substantially larger m SUSY , such that the gluino is long-lived, is obtained assuming that Higgsinos are heavy (possibly with masses of order m SUSY : in such a case the RGE for the Higgs quartic are those of the SM (with slightly different values of g 2,3 due to the light gluino and wino), and the Higgs quartic can remain positive up to m SUSY ∼ M Pl within the uncertainty range for the top quark mass.

Collider signals
Next, we discuss some aspects of the phenomenology of long-lived gluinos at hadron colliders, in particular LHC. Long-lived gluinos can be pair produced and after hadronization form long-lived hybrid states with SM quarks and gluons, known as 'R-hadrons'. We conservatively assume that the signal at the LHC is just energy deposit in the calorimeter, rather than charged particles in the tracker. It is difficult to trigger on these event and so an initial state jet is required. The LHC places the limit Mg > 1.55 TeV on a Majorana gluino [43].
The other possibility is the production of agg bound state. Assuming that states with = 0 dominate the rates, they are color 8 A with spin S = 1 and color singlets or 8 S with S = 0 (see eq. 9). The production cross sections are given by gluon and quark fusion respectively where L ij is the luminosity of partons ij. The decay rates are given by [8] with F = 2 for the Majorana gluino and F = 1 for a Dirac particle, and with the channel strength λ 1 = 3 and λ 8 = 3/2. Since the resonances annihilate to two gluons or two quarks, we assume a 100% branching ratio to two jets and apply the LHC di-jet bounds [44] to the sum of the cross sections. In fig.  5 we compare the bounds on the resonances to, slightly stronger, the R-hadron bound.
Concerning future colliders, the expected reach of a 100 TeV hadron collider with 1000 fb −1 is 7 (9)TeV for a Majorana (Dirac) gluino, having used [45] to perform an approximate rescaling. The R-hadron search would then reach 10 TeV and 14.5 TeV respectively. Thus a 100 TeV collider would reach the benchmark mass of a thermally produced Dirac gluino, which recently was found to be a dark matter candidate [8].

Implications for Dark Matter co-annihilations
The thermal relic abundance of a particle is affected by co-annihilations with particles of similar mass. One example is co-annihilations of neutralino DM with heavier colored particles, for example gluinos. Co-annihilations can be enhanced by Sommerfeld corrections [12] and bound-state formation [42,13]. We point out that a much bigger effect is produced by nonperturbative QCD effects after the QCD phase transition, if the mass splitting ∆M between  the co-annihilating species is comparable or smaller than Λ QCD . Such a near-degeneracy is unnatural. This is shown in fig. 7a in the neutralino/gluino co-annihilation case, assuming that squarks mediate fast neutralino/gluino rates. We see that the neutralino mass which reproduces the observed DM density gets much higher at ∆M < ∼ GeV. In the limit ∆M GeV the relic abundance is dominantly set by the new QCD annihilations. As a result, the neutralino mass can reach up to a PeV, heavier than the maximal relic DM mass allowed if DM annihilations are dominated by partial waves with low [14].

Quasi-stable squark
In the previous sections we considered a Majorana gluino. A real scalar in the octet of SU(3) c would behave similarly to the Majorana gluino. On the other hand, a (quasi)stable particle in the fundamental 3 of color SU(3) c can behave in a qualitatively different way. Since the 3 is a complex representation, the particle must be a complex scalar or a Dirac fermion, which can carry a conserved charge.
For definiteness, we consider the possibility of a (quasi)stable squark, and more specifically a stopt, as RGE effects tend to maket lighter than other squarks. A stable stop arises ift is the lightest SUSY particle (LSP) and R-parity is conserved. A quasi-stable stop arises if R-parity is almost conserved, or if the stop decays slowly into the LSP: this can happen e.g. when the LSP is a gravitino. Collider bounds on stops [48] tend to ignore the possibility that the lighter stopt is the (quasi)stable LSP, because it is perceived to be already excluded by cosmology.
In cosmology, perturbative QCDtt * → gg annihilations dominate overtt → tt annihilations  and leave a roughly equal amount of relict andt * . Perturbative QCD annihilations are enhanced by Sommerfeld and bound-state effects, computed in [13]. The relict abundance after perturbative annihilations is plotted in fig. 6 and approximated by For Mt < PeV this is smaller than the baryon asymmetry n b /s ∼ 10 −10 , that we neglect given that its effect is model dependent. Indeed, we do not know how the baryon asymmetry is generated: it might be generated at the weak scale such that it would not affect heavier stops. Even if a baryon asymmetry is present at stop decoupling,tt ↔t * t scatterings could easily concentrate the baryon asymmetry to lighter baryons fast enough that the asymmetry is irrelevant for stops. If instead the baryon asymmetry enhances the relic stop abundance, bounds would become stronger.
After the QCD phase transition, stops form hadrons. In view of the large QCD cross sections, the stop hadrons with dominant abundance are deeply-bounded states which contain stops only. They arett * and the charged baryonsttt. Both fall to the ground state and decay through annihilations of the constituents. In particular, a bound state containing two or more stops decays, in its ground state, with a life-time Γtt ∼ α 3 3 M 3 t σttv rel where the cross section for tt → tt can be roughly estimated as σttv rel ∼ i={1,2,3} α 2  Figure 7: Non-perturbative QCD annihilations that take place at T < ∼ Λ QCD significantly increase the DM neutralino mass such that the observed DM abundance is reproduced trough co-annihilations with gluinos (left) or stops (right), if their mass difference with neutralinos is smaller than a few GeV. In the case of stops (right panel), the big effect is only estimated and only present if stop baryons decay to SM particles before decaying to neutralinos; otherwise confinement only gives a O(1) effect.
We expect a roughly equal number oftt * annihilations for each producedttt given that QCD group algebra implies that bothtt * andtt feel an attractive Coulombian QCD force, such that they can form deep, unbreakable, Coulombian bound states. Assuming that at binds with probability ℘ to at and with probability 1 − ℘ to at * and thereby that a deeptt binds with probability 1 − ℘ tot and with probability ℘ to at * , the average number oftt * per produced baryon is Ntt * Nttt + Nt * t * t * = This equals 3 assuming no baryon asymmetry r ≡ Nt/N * t and ℘ = 1/2, namely neglecting thattt * is more attractive thantt. Extra hadrons and mesons that contain quarks have a much smaller abundance, that is not relevant here. If the charge 2 statesttt decay fast on cosmological scales, final abundances and bounds are similar to the gluino case. If (quasi)stable, they are instead subject to strong cosmological constraints. In particular during BBNt * t * t * can bind to 4 He reducing its charge and thereby the Coulomb suppression of nuclear reactions, opening up a new channel for 6 Li production, (t * t * t * 4 He) + D → 6 Li +t * t * t * , which can strongly alter Lithium abundances (see [49] for a brief review). Charge −1 states with lifetime 10 5 are subject to the BBN bound Y 2.5 × 10 −17 [23]. A study of analogous constraints on relics with charge −2 is beyond the scope of this paper.
Next, we study the scenario where a quasi-stable stop co-annihilates with a slightly lighter DM neutralino. Post-confinement effects are relevant if ∆M < ∼ GeV. Roughly half of the stops formtt * mesons, and the others formttt baryons. The impact on the DM abundance is very different, depending on which process dominatesttt decays. If it is dominated by stop annihilations into SM particles, post-confinement effects strongly suppress the DM abundance, similarly to the gluino/neutralino co-annihilation scenario. A much smaller order one effect is obtained if instead stops decay to DM neutralinos and SM particles with rate Γt > ∼ Γtt. The region where the DM abundance is reproduced is estimated in fig. 7b in the two extreme possibilities, having assumed σ QCD = 1/Λ 2 QCD .

Conclusions
We have reconsidered the relic abundance of neutral colored relics, finding that hadron collisions at temperatures below the QCD scale reduce it by a few orders of magnitude. In particular we considered a quasi-stable gluino: fig. 1 shows its relic abundance, and fig. 4 the cosmological constraints, taking into account the new effect and new data. Co-annihilations between gluinos and neutralino DM are similarly strongly affected by confinement, provided that their mass difference is smaller than a few GeV, as shown in fig. 7a.
In section 3.4 we considered charged colored relics, considering in particular the case of a quasi-stable stop. In this case, confinement gives a big contribution to co-annihilations with neutralinos only ifttt baryons decay into SM particles viatt → tt before that stop decays to neutralinos.
We consider emission of a single-vector V a in dipole approximation, such that the angular momenta of the initial and final states differ by ∆ = ±1. We denote with α the non-abelian gauge coupling, with M a the vector mass, and with M the common mass of the two particles which form the bound state.

A.1 Cross sections for bound state formation
Production of a bound state with angular momentum proceeds from initial states with angular momentum ± 1: (σ n bsf v rel ) a = (σ n bsf v rel ) −1→ a + (σ n bsf v rel ) +1→ a . The cross sections are where R n ,ij is the bound state wave-function in the two-particle space |i ⊗ |j , and R p ,ij is the wave-function of the initial free state with relative momentum p and angular momentum . The other cross section is that holds separately for each initial channel J and final channel J , using the notations of [13].

A.2 Bound state decays
The decay widths of a bound state trough single-vector emission are obtained from the previous expressions substituting the free-particle final state wave function R p with the wave-function of the desired final bound states. Assuming again degenerate (or massless) vectors and a bound state in a representation R with dimension d R , we find Γ(n, → n , − 1) = 1 d