Effects of QCD bound states on dark matter relic abundance

We study scenarios where there exists an exotic massive particle charged under QCD in the early Universe. We calculate the formation and dissociation rates of bound states formed by pairs of these particles, and apply the results in dark matter (DM) coannihilation scenarios, including also the Sommerfeld effect. We find that on top of the Sommerfeld enhancement, bound-state effects can further significantly increase the largest possible DM masses which can give the observed DM relic abundance, by $\sim 30 - 100\%$ with respect to values obtained by considering the Sommerfeld effect only, for the color triplet or octet exotic particles we consider. In particular, it indicates that the Bino DM mass in the right-handed stop-Bino coannihilation scenario in the Minimal Supersymmetric extension of the Standard Model (MSSM) can reach $\sim 2.5$ TeV, even though the potential between the stop and antistop prior to forming a bound state is repulsive. We also apply the bound-state effects in the calculations of relic abundance of long-lived or metastable massive colored particles, and discuss the implications on the BBN constraints and the abundance of a super-weakly interacting DM. The corrections for the bound-state effect when the exotic massive colored particles also carry electric charges, and the collider bounds are also discussed.


Introduction
In scenarios of physics beyond the Standard Model (SM), the early Universe may have been inhabited by exotic particles charged under QCD. Due to their strong interactions with SM particles, they were initially in thermal equilibrium and later froze out.
If the colored particle is stable, stringent constraints for a strongly interacting dark matter (DM) particle apply (see [1] and references therein). On the other hand, if metastable, it may first freeze out and then decay. This can happen in, e.g., the R-parity conserving Minimal Supersymmetric extension of the Standard Model (MSSM) where the next-to-lightest supersymmetric particle (NLSP) is colored, and the lightest supersymmetric particle (LSP), which is a DM candidate, is extremely weakly interacting (the superWIMP), such as the gravitino or axino. Typically, observational and experimental bounds are applicable only to a metastable colored particle with a lifetime 0.1 sec. In particular, important constraints can be derived from Big-Bang nucleosynthesis (BBN) (see e.g. [2]). Moreover, the decays of the colored particles into superWIMPs can contribute non-thermally to the relic density of DM.
The main purpose of this paper is to study the effects of exotic massive colored particles on DM relic abundance, assuming that they share the same discrete symmetry stabilizing DM (e.g., R-parity in supersymmetric models and KK-parity in UED models). We calculate the DM relic abundance in scenarios where the colored particle coannihilates with the WIMP. We also discuss implications of a metastable colored particle on BBN and the DM relic abundance in the superWIMP scenario.
The major update of our work compared to previous calculations is the inclusion of the effects of QCD bound states of the exotic colored particles, X 1 and X 2 . The formation of a QCD bound state, η, occurs via the process X 1 X 2 → ηg, in which the gluon, g, takes away the binding energy. η can then decay to SM particles via the annihilation of its constituents, X 1 and X 2 . The net effect is to remove X 1 and X 2 from the thermal bath. As will be shown in the following, these processes increase significantly the effective annihilation cross section of the X's or DM. We note that due to color charge conservation, the QCD potential between X 1 and X 2 before they form a bound state can be different from the one after they become the constitutes of η, with a gluon emitted. Moreover, in contrary to the Sommerfeld effect [32] which can enhance the effective annihilation cross section when the potential between the incoming X 1 and X 2 is attractive, bound-state effects can enhance the effective annihilation cross section no matter whether the potential between the incoming X 1 and X 2 is attractive, repulsive or zero. In addition, we calculate the corrections to the bound-state effects when the colored particle is electrically charged.
Furthermore, we study the compatibility of these scenarios with collider constraints, especially those derived from the Large Hadron Collider (LHC) experiments. Studying LHC constraints is particularly timely as the undergoing LHC experiments are capable of excluding masses of colored particles up to the TeV scale. It is worth noting that while the DM relic abundance and BBN constraints often imply upper limits on the masses of colored particles or DM, collider results impose lower limits on the masses of colored particles. Studying both BBN/DM relic density and collider constraints therefore probe the experimentally allowed parameter "window" of the long-lived colored particles or the coannihilating DM scenarios.
The rest of the paper is organized as follows. In Section 2, we discuss the formation, dissociation and annihilation decay of bound states formed by a pair of massive colored particles in the early Universe, focusing in particular on the cases of complex scalars (S3) or Dirac fermions (F3) in the color SU(3) fundamental representation, and real scalars (S8) or Majorana fermions (F8) in the adjoint representation. In Section 3, we calculate the DM relic abundance in massive colored particle -WIMP coannihilation scenarios and the contributions from decays of metastable massive colored particles on the relic density of a superWIMP DM, by including the Sommerfeld and bound-state effects in the Boltzmann equation. We also consider the constraints from BBN on the long-lived massive colored particle scenarios, and discuss the electric charge corrections for the bound-state effect. The thermally-averaged s-and p-wave annihilation cross sections needed in the calculations are listed in Appendix A. In Section 4, we study the collider limits on the exotic massive colored particles. Finally, we summarize our conclusions in Section 5.

Bound state formalism
In this section, we discuss the formation and dissociation of the QCD bound state, X 1 X 2 ↔ ηg, in the early Universe. DM bound-state formation due to some new abelian gauge force was considered in [33], and a systematic study of its effects on DM relic abundance was performed in [34]. Using the SM non-abelian color SU(3) group, gluino bound-state effects on neutralino DM coannihilation in the MSSM have been investigated in [23]. A field-theoretic framework for the computations of bound-state effects was established in [35]. A formalism for bound states based on non-relativistic effective theories was given in [36,37]. Bound-state formation due to a Yukawa potential was studied in [38]. See e.g. [39][40][41][42] for other phenomenological implications of DM bound-state formation.
We are concerned with massive colored particles with masses m X 1 , m X 2 Λ QCD . In the early Universe before the quark-hadron phase transition era, the long-range interaction between two massive colored particles can be described by a perturbative one-gluon exchange Coulomb-like potential, which has the form: in which ζ is determined by the quadratic Casimir coefficients of the color representations of the individual colored particles, X 1 and X 2 (C X 1 and C X 2 , respectively), as well as of the one by taking X 1 and X 2 together in a specific color state (C X 1 X 2 ): where α s > 0 is the QCD coupling strength. A positive, negative or zero value of ζ gives an attractive, repulsive or zero potential, respectively. The colored particles we consider in this paper include a complex scalar and a Dirac fermion in the color SU(3) fundamental representation, a real scalar and a Majorana fermion in the adjoint representation. The X 1 X 2 combinations are S3S3, F 3F 3, S8S8 and F 8F 8, abbreviated in the following as S3, F3, S8 and F8, respectively, and hence m X 1 = m X 2 ≡ m X . Examples of S3 and F8 are a squark-antisquark pair and a gluino-gluino pair, respectively, in the MSSM. A KK quarkantiquark pair in models of UED is a realization of F3. One can also build models for the S8 case [43,44]. The product of a color triplet and an anti-triplet is decomposed as and the product of two color octets is decomposed as where the subscripts "S" and "A" indicate symmetric and anti-symmetric color states, respectively. Therefore, the relevant quadratic Casimir coefficients of the color representations for our calculations are C 1 = 0, C 3 = 4/3, C 8 = 3, C 10 = 6 and C 27 = 8. 1 In principle, a bound state can form as long as the potential for it is attractive. In this paper, we focus on the color-singlet bound state, since it is expected to be the deepest bound (i.e., the ground 1 We note that at a temperature T of the Universe, the screening effect from the quarks and gluons in the thermal plasma induces a thermal mass m th ∼ √ α s T to the gluon, modifying the QCD Coulomb potential to a Yukawa one. However, as emphasized in [38], the Coulomb potential is a good approximation as long as the momentum transfer between the two incoming particles, ∼ m X T m X , is larger than m th . We can see that this condition is well satisfied at the usual freeze-out temperature T ∼ m X /20, and further better satisfied for the bound-state effect calculation since the effect of which is important at even lower temperature T ∼ α 2 s m X , as will be shown in the next section. state) and the most copiously produced one, in analogy to atomic physics 2 . Therefore, the coefficient, ζ, in Eq. (1) for the bound states in the S3 and F3 cases is 1/2 × (4/3 + 4/3 − 0)α s = (4/3)α s , and is 1/2 × (3 + 3 − 0)α s = 3α s for the S8 and F8 cases. In all the four cases, we consider that the color-singlet bound states have total orbital angular momentum L = 0 and spin S = 0 3 . The normalized spatial wave function of such a bound state is where a is the Bohr radius, given as a = (ζµ) −1 , where µ ≡ m X 1 m X 2 /(m X 1 + m X 2 ) = m X /2 is the reduced mass. The binding energy of the bound state is 2.1 Bound-state formation, dissociation, and annihilation decay We start with a general description of the formation, dissociation and annihilation decay processes of bound states without specifying the color representation of the particles. First of all, due to color charge conservation, the emission (absorption) of a gluon during boundstate formation (dissociation) makes the color representation of the bound state η not necessarily be the same as the free pair X 1 X 2 . Therefore, the coefficient in the Coulomb potential for the free pair, denoted as ζ , is not necessarily equal to the one for the bound state. In particular, ζ can be negative, so that the potential is repulsive for the free pair. As will be shown, at high temperature in the early Universe, the massive colored particles can have enough kinetic energy to overcome a repulsive potential to form a bound state.
We follow [23] to calculate the bound-state formation and dissociation cross sections, where the method is adapted from the calculations of the photoelectric effect for an atom [45]. The essence of the calculation is to evaluate the transition matrix element between the bound state and the free pair state: 2 In [41], the formation of bound states at excited energy levels is discussed for bound-state effects in the late Universe, and it is found that the total bound-state formation cross section is dominated by levels with principle quantum numbers n < ζ/v rel , where v rel is the relative velocity of the incoming particles. Compared to the ground state, the contribution from the excited states enhances the total bound-state formation cross section by a logarithmic factor ∼ log(ζ/v rel ), which is significant for v rel ∼ 10 −3 in the galactic halo. However, in the early Universe at temperatures relevant for the dark matter relic abundance calculation, v rel is of order 10 −1 , so that ζ/v rel ∼ 1. Moreover, compared to the ground state, the excited states are easier to be dissociated by gluons in the thermal bath, while the dissociation is not a concern for bound states in the late Universe. Therefore, the contribution from the excited states is not significant for the relic abundance calculation. Nevertheless, in the next section we will also show results with a factor of 2 enhancement of the bound-state effect from the considerations of the excited states contribution as well as other uncertainties in our calculations. 3 The total wave function of the bound state is a product of the spatial, spin and color parts of the wave functions. For the S8 (F8) case, because of the nature of identical particles, the total wave function should be symmetric (antisymmetric). L = 0 gives symmetric (symmetric) spatial wave function, and S = 0 gives symmetric (anti-symmetric) spin wave function. Together with the symmetric color wave function of the color-1 S state, indeed the requirement for the total wave function is satisfied.
where φ f is the wave function of the free pair and φ i ≡ φ η (r). The gluon has a momentum k and a polarization c , where "c" is the color index.
For the free pair, the normalized spatial part of the wave function is (see Section 136 of [46]) where | p| is the relative momentum of the free pair, expressed in terms of the reduced mass and their relative velocity as | p| = µv rel . P L ( p· r | p|r ) is the Legendre polynomial and δ L (a real number) is the phase shift. Note that the form of the radial function R pL (r) for an attractive potential between the free pair is different from the one for a repulsive potential (see Section 36 of [46] for details). Consider the spatial part of the wave functions only, the differential dissociation cross section is given as where ω ≡ | k| is the energy of the gluon. The explicit α s factor in the above equation comes from the coupling between the emitted gluon and the massive colored particle. In the Lagrangian of the quantum field theory, this coupling is from the covariant derivative of the kinetic term of the massive colored particle, and it takes the form of ig s T c , where g s = √ 4πα s is the strong coupling, and T c are the generator matrices for the color representation in which the massive colored particle lies. We will specify T c in the next subsection when we consider the color part of the wave functions for the four cases of our interest.
We use the dipole approximation 4 , e i k· r 1, to calculate the transition matrix element Eq. (8). This means that only the L = 1 term in φ f has a non-zero contribution. Also, considering that it is the absolute square of the transition matrix element that appears in Eq. (10), we can drop the phase factors and rewrite φ f as Defining dimensionless quantities ν ≡ |ζ |/v rel (12) and we can write down the integrated dissociation cross section, averaged over the incoming gluon spin polarizations. The result depends on whether the free pair feels an attractive (denoted by the subscript "a") or a repulsive (denoted by the subscript "r") Coulomb potential: σ 0 dis,r = 2 9 π 2 3 α s a 2 E B ω In the case that the free pair feels no potential (denoted by the subscript "f ree", and see Section 33 of [46] for the radial function), we find One can check that in the ζ → 0 limit, Eq. (14) and Eq. (15) both become Eq. (16). The superscript "0" in the above three equations indicates that we have considered the spatial part of the wave functions only, while the full wave functions are products of spatial, color and spin wave functions. Also, if the particles are identical, one needs to symmetrize or anti-symmetrize the wave functions. The full dissociation cross section, σ dis , after taking into account color, spin and the symmetry factors, is related to the bound-state formation cross section, σ bsf , via the Milne relation: where g g,X 1 ,X 2 ,η are the degrees of freedom of gluon, X 1 , X 2 and η, respectively. Note that if X 1 and X 2 are identical, the left-hand side of Eq. (17) has to be multiplied by 1/2 to avoid double counting the number of bound-state formation reactions. Bound state can be destroyed not only by the dissociation process, but also by decays. Moreover, the decays can happen in two ways: the constituent particles inside the bound state can annihilate between themselves (annihilation decay) or an individual constituent particle can decay by itself. The effects of these two kinds of decays on the relic density of the metastable colored particles or DM are different. Since we assume that the constituent particles in the bound state have the same discrete symmetry as the DM particle, the annihilation decay to SM particles removes, for example, two R-odd numbers in supersymmetry, while the individual constitute particle decay does not change the R-odd number. For the colored particle coannihilating in the WIMP DM scenario and the metastable colored particle in the superWIMP DM scenario, the individual constituent particle decay rate is suppressed either by the small mass difference or by the very small coupling between the massive colored particle and the DM particle 5 , while the annihilation decay rate is not suppressed and is proportional to the large mass of the colored particle. Therefore, we will hereafter neglect the individual constitute particle decay rate compared to the annihilation decay rate.

Results for S3, F3, S8 and F8
Here we present the full bound-state formation and dissociation cross sections for the cases of S3, F3, S8 and F8, as well as the annihilation decay rates.

S3 and F3
Since we consider that the bound state is a color-singlet state, the emission (absorption) of a gluon in the bound-state formation (dissociation) process dictates that for both S3 and F3 the free pair state must be in a color-octet state (see Eq. (3)), due to color charge conservation. The normalized color wave function is δ kj / √ 3 for the bound state, and λ b ij / √ 2 for the free pair, where λ b ij are the Gell-Mann matrices, and the color indices i, j, k = 1 − 3, b = 1 − 8. The generator T c takes the form λ c ki /2 . Therefore, the color part of the wave functions contributes to σ dis as For S3, there is no spin wave function to worry about. While for F3, without considering the bound state, a pair of heavy colored fermion and anti-fermion can have 3/4 chance in spin-triplet configurations with S = 1 and 1/4 chance in a spin-singlet configuration with S = 0. Since we only consider a bound state with S = 0, then by neglecting the spin-orbit interaction we will only consider a free pair that is also in S = 0 state. Therefore, we consider that both the bound state and the free pair have the same spin wave function, given as so that the spin part of the wave functions does not introduce a factor for σ dis . However, in the next section we will see that when including the bound-state formation and dissociation cross sections in the Boltzmann equation, we need to introduce an additional factor of 1/4 for F3 compared to S3 to take into account the fact that we have only considered the S = 0 possibility in the former.
Putting the factor of 1/8 from the incoming gluon color averaging, the factor of 4/3 from the color part of the wave functions, and by noticing that the free pair has a repulsive potential with (2)), we get the full dissociation cross section for S3 and F3, in which the quantities inside σ 0 dis,r are given in Eqs. (6), (7), (12) and (13) with ζ = (4/3)α s . From Eq. (17), the bound-state formation cross sections are and where the degrees of freedom are written explicitly. For the bound-state annihilation decay, we consider the dominant decay mode only, which is the two-gluon final state (see e.g. [47]), and the results are and where ζ = (4/3)α s . In the above two equations, the α s factor explicitly written is evaluated at the scale of 2m X , while the α s inside ζ is evaluated at the scale of the inverse Bohr radius, a −1 .

S8 and F8
Due to the nature of identical particles, the total wave functions need to be symmetric for S8 whereas anti-symmetric for F8.
The F8 case was studied in detail in [23], and the result for the gluon dissociation of a color-1 S bound state with (S = 0, L = 0) into a free pair in an 8 A state with (S = 0, L = 1) is where the factor 3 comes from the color part of the wave functions together with the generator T c = −if cde in the coupling between the gluon and the massive colored particle, where f cde are the SU(3) structure constants. The factor 4 comes from symmetrization of the spatial part of the bound-state wave function (L = 0) and anti-symmetrization of the spatial part of the free pair wave function (L = 1). 1/8 comes from the color averaging of the incoming gluon. The factor 1/2 is introduced to avoid double counting of the two identical massive colored particles in the outgoing free pair phase-space integration. The spin part of the wave functions do not introduce any extra factor. The quantities inside σ 0 dis,a are given in Eqs. (6), (7), (12) and (13) with ζ = 3α s and The S8 case is exactly the same as the F8 case, namely, a transition from the color-1 S bound state with (S = 0, L = 0) into the 8 A free pair state with (S = 0, L = 1). The only difference is that while for the F8 case the S = 0 state means that the spin wave function is anti-symmetric (i.e., a spin-singlet configuration), for the S8 case there is no spin to worry about, so that for the latter the total wave functions for both bound state and free pair state are symmetric, as they should be. Therefore, we have The corresponding bound-state formation cross sections are and where the factor of 2 in the front of the right-hand side of the two equations is actually the factor of 1/2 which would have been in the left-hand side of Eq. (17) to avoid double counting the number of bound-state formation reactions. In Eqs. (27) and (28), the degrees of freedom are written explicitly. Again, we consider the dominant two-gluon annihilation decay channel only, and the annihilation decay rates are (see e.g. [47]), and where ζ = 3α s . The scales of evaluating the α s explicitly written and the one inside ζ in the above two equations are understood similarly as in the S3 and F3 cases.

Thermal averaging
To study the bound-state effects on the relic abundance of the massive colored particles or the DM, we need the thermally-averaged bound-state dissociation and annihilation decay rates, as well as the formation cross section times the relative velocity of the free pair, since it is these quantities that appear in the Boltzmann equation which determines the evolution of the density of the massive colored particles or the DM with temperature, T . By defining two dimensionless variables, z ≡ E B /T and u ≡ 1 2 µv 2 rel /T , we can rewrite σ dis and σ bsf as functions of z and u, together with factors not changing with T . In particular, the relevant quantities are expressed as The thermally-averaged bound-state dissociation rate is The thermally-averaged bound-state formation cross section times relative velocity is where f (v rel ) is the Maxwell-Boltzmann distribution function for v rel , given as The factor 1 e ω/T −1 in Eq. (35) accounts for the stimulated emission due to the gluons in the thermal bath. Using Eqs. (33) and (7), σv bsf can be rewritten as The thermally-averaged bound-state annihilation decay rate is where m η is the mass of the bound state, given as m η = m X 1 + m X 2 − E B . In the above formula, we have assumed the Maxwell-Boltzmann approximation for the bound-state equilibrium distribution, and K 1,2 (m η /T ) are the modified Bessel functions of the second kind. At m η T , Γ η ≈ Γ η .

Relic abundance calculation
We need to set up a Boltzmann equation (or a coupled set of Boltzmann equations) to solve for the relic abundance of the massive colored particles in superWIMP scenario or of the DM in WIMP scenario. The necessary ingredients for the bound state in the Boltzmann equation, i.e., the thermallyaveraged bound-state formation cross section, dissociation and annihilation decay rates, are given in the previous section. There is another important ingredient -the Sommerfeld effect -needs to be considered for massive annihilating particles feeling a long-range force. The Sommerfeld effect in the calculations of thermal relic abundance has been studied extensively in the literature (see e.g. [48][49][50][51][52][53][54][55]). Let us first briefly describe the Sommerfeld corrections to the annihilation cross sections of the massive colored particles of our interest.

Sommerfeld effect
We consider the Sommerfeld effect resulting from the long-range Coulomb potential between the two annihilating massive colored particles. Under the influence of a Coulomb potential of the form V (r) = −α/r, the Sommerfeld corrected s-wave annihilation cross section can be written as where An attractive potential (α > 0) results in an enhancement (S > 1), while a repulsive potential (α < 0) results in a suppression (S < 1). In Eq. (39), the perturbative s-wave cross section, a, does not depend on temperature. Therefore, the thermally-averaged Sommerfeld corrected s-wave cross section is a S(α/v rel ) , where with f (v rel ) given in Eq. (36). For the massive colored particles of our interest, we expect that the dominate annihilation channels are S3S3, F 3F 3, S8S8 and F 8F 8 annihilation into a pair of gluons, gg, and into quark-antiquark pairs, qq. We consider the Sommerfeld corrected s-wave cross sections and the tree-level p-wave cross sections in our calculation. The relevant expressions are collected in Appendix A. To consider the Sommerfeld effect, we need to decompose an s-wave cross section into partial cross sections contributed from each of the two-body states in different color representations, as given in Eqs. (3) and (4), because different representations correspond to different Coulomb-like potentials. We follow the decompositions given in [54]. The thermally-averaged Sommerfeld factors are The s-wave cross sections vanish for S3S3 → qq and S8S8 → qq.

Boltzmann equation
We now have all the ingredients to write down and solve the Boltzmann equation. The general formulae for N species of exotic particles in calculating the thermal relic densities are given in [25].
In that formalism, the possibility of when the rates for interconverting some of the species are not sufficiently large is taken into account, and in such case a coupled set of Boltzmann equations, rather than a single Boltzmann equation, is needed. Also, a simple method for implementing the boundstate effects in the Boltzmann equations is given, and a detailed example for the gluino-gluino bound state is shown in [23].
In this paper, we consider two exotic species -the massive colored particle and the DM, and by considering that they share the same discrete symmetry stabilizing DM, we assume that the decay of one massive colored particle produces one DM together with some SM particles. In the superWIMP scenario, because the DM has long been out of the thermal bath when the freeze-out of the massive colored particle happens, we only need to solve the Boltzmann equation for the latter. In the WIMP scenario, we are interested in the coannihilation between the two exotic species. As highlighted in [23,24], coannihilation is effective only if the interconversion rate between the two species is sufficiently large compared to the Hubble expansion rate, otherwise the two species would freeze out separately. Without committing to specific particle theory models, in this work we assume the interconversion rate is sufficiently large so that the two species freeze out together and we can use a single Boltzmann equation to calculate the DM density. We emphasize that this condition needs to be checked when one consider coannihilations in specific DM models.
Also, as we are mainly interested in QCD interactions and further model specifications would have been needed, in our calculation we neglect the (co)annihilation cross sections for DM -DM and DMmassive colored particle, in comparing to the annihilation cross sections between the massive colored particles. In the WIMP DM scenario, indeed this is usually a good approximation. For example, for the right-handed stop-Bino coannihilation in the MSSM, the stop -antistop annihilation to gg channel dominates the effective annihilation cross section in the coannihilation region for an O(TeV) Bino [13], if the ratio of the lighter stop -heavier stop -Higgs coupling to heavier stop mass is not very large [11]; for the gluino-neutralino coannihilation, the gluino pair annihilations to gg and qq dominate over the neutralino-gluino and neutralino pair (co)annihilation cross sections.
With these assumptions and approximations, in the WIMP DM scenario the evolution of the total yield, that is, the total number density of exotic particles (i.e., the DM χ and the massive colored particles X) over the entropy density,Ỹ ≡ Y X + Y χ , where Y X ≡ n X /s and Y χ ≡ n χ /s, can be described by a single Boltzmann equation as follows: where and G N is the gravitational constant, m χ is the DM mass, g * s and g * are the numbers of effectively massless degrees of freedom associated with the entropy density and the energy density, respectively. Y eq is the equilibrium value of the total yield, given asỸ eq = Y eq X + Y eq χ = (n eq X + n eq χ )/s, where are the thermal equilibrium number densities [56][57][58]. g χ is the DM degrees of freedom. g X is the degrees of freedom of massive colored particles, being 6, 12, 8 and 16 for the cases of S3, F3, S8 and F8, respectively. We note that for the S3 and F3 cases n X and n eq X take into account both the massive colored particles and anti-particles, assuming that there is no asymmetry between the number densities of them. The condition of a sufficiently large interconversion rate between χ and X comparing to the Hubble expansion rate guarantees that n X /n χ ≈ n eq X /n eq χ , to a very good approximation 6 . By further defining the thermally-averaged effective annihilation cross section, σ ef f v , can be written as where σv XX→SM is the thermally-averaged total annihilation cross section of X's into SM particles, with the Sommerfeld and bound-state effects taken into account, given as For the S3 and F3 cases, since n X includes both particles and anti-particles, we need to time a factor of 1/2 to the σv XX→SM term in the right-hand side of Eq. (47), if we are using the usual spin and color averaged cross sections, as explained in the Appendix of [56]. Furthermore, for the F3 case, we need to time another factor of 1/4 to the σv bsf and Γ dis in Eq. (48), since we only consider the bound-state formation and dissociation for the S = 0 state, which occurs with 1/4 of the possibilities of the total numbers of the spin configurations for two spin-1/2 fermions.
The physical interpretation of the above formulae is as follows. First of all, since the mass of the bound state, m η = 2m X − E B , is much larger than m X (note that E B ∼ O(10 −2 ) m X for α s ∼ 0.1), the equilibrium number density of the bound state is negligible compared to the one for X, for m X T . Also, since the bound-state annihilation decay rate, Γ η , is many orders of magnitude larger than the Hubble expansion rate for a massive colored particle with a mass of O(1 − 100) TeV , which is of our interest in this paper, the yield of the bound state keeps decreasing and is always negligible compared to the one for X. These explain why we do not have to worry about the yield of the bound state in the Boltzmann equation (43).
The second term on the right-hand side of Eq. (48) accounts for the bound-state effect: once formed, the bound-state can be destroyed either by dissociation back to individual X's or by annihilation decay to SM particles (as explained in Sec. 2.1, we neglect the individual X decay rate compared to the annihilation decay rate); only the latter process reduces the total number of X's. At high temperatures, a large fraction of the gluons in the thermal bath have enough energy to break up the bound states before they can annihilation decay, so that the bound-state effect is negligible. With the temperature falling, this energetic fraction becomes smaller, and eventually when the temperature becomes comparable and then is below the bound-state binding energy, the annihilation decay process dominates over the dissociation process, and the bound-state effect helps to reduce the total number of X's. Since we assume that X and the DM particle share the same discrete symmetry which stabilizes DM, the net effect is to reduce the DM relic abundance. In other words, the formation and the subsequent annihilation decay of bound states enhance as a new channel of the annihilation of X's. We note that to have the bound-state effect reduce the total number of X's efficiently, the bound-state formation rate should be larger or comparable to the Hubble expansion rate, otherwise bound states essentially cannot form in the first place. Indeed eventually the boundstate effect will cease to be effective as the building block of the bound state, n X , is decreasing with the expanding of the Universe, as long as the bound state formation cross section does not increase too fast with the decrease of temperature so that it compensates the decreasing of n X .
Following the usual procedure to obtain the DM relic abundance [56][57][58] by integrating Eq. (43) from a small value of x (whenỸ =Ỹ eq ) to the current temperature of the Universe (when essentially x → ∞), we getỸ 0 , where the subscript "0" indicates today, and the DM relic abundance is The above formulae can describe the evolution of the yield of the massive colored particles in the superWIMP scenario by the following modifications: in Eqs. (43), (44) and (49), changeỸ to Y X , Y eq to Y eq X , m χ to m X , and σ ef f v to σv XX→SM , which is given in Eq. (48). After solving for the relic abundance of X would have today, if it had not decayed, the amount of DM produced from the decays of X's is obtained as where m SW is the mass of the superWIMP DM, and the superscript "non-th" indicates that there can be also thermally produced amount, Ω th SW h 2 , which contributes to the total DM relic abundance together with Ω non-th SW h 2 .

Bound-state effects
To get an idea of the size of bound-state effects, we plot in Fig. 1 with orange curves the ratio of thermally-averaged Sommerfeld corrected XX → gg, qq annihilation cross section to the one without Sommerfeld correction, i.e., σv rel (XX → gg, qq) Sommerfeld σv rel (XX → gg, qq) w/o Sommerfeld , the thermally-averaged bound-state formation cross section with (solid black curve) and without (dotted black curve) considering bound state dissociation and annihilation decay rates, also normalized to the tree-level XX → gg, qq annihilation cross section, i.e., σv bsf respectively, and with purple curves for the same ones but multiplied by a factor of 2 in the boundstate formation cross section (solid and dotted) and dissociation rate (solid), that is 7 , The left and right panels are for the S3 and F8 cases, respectively 8 . As can be observed from both panels, at early times when E B /T 1, the solid black curves are much smaller than 1, although the dotted black curves have already bigger than 1. This is because for this high temperature range, a large fraction of the gluons in the thermal bath are energetic enough to break up the newly formed bound states before the latter can annihilation decay, that is, Γ dis Γ η . Also, a factor of 2 increase of the bound-state formation cross section and also in the dissociation rate (since they are related via the Milne relation) has little effect. With the decrease of temperature, Γ dis becomes smaller and eventually negligible compared to Γ η when E B /T 1, so that the solid and dotted black (and purple) curves merge. By comparing the orange and solid black (and purple) curves, one can see that when E B /T 1, the size of the full-fledged bound-state effect is of the same order and even larger than that of the Sommerfeld enhancement.
There is a qualitative difference between the S3 and F8 cases at E B /T 1: while in the F8 case the solid black and purple curves keep growing with the decrease of temperature, in the S3 case they are decreasing after achieving maximum values around E B /T ∼ 2. This is due to the fact that 7 To account for the errors involving the evaluations of α s 's [59], the corrections from a more accurate QCD potential and the thermal mass of the gluon [54], as well as the possibility that excited bound states may also contribute to the bound-state effect [33], we plot purple curves in this and following figures with an uncertainty of a factor of 2. 8 We use a common value, α s = 0.1, for all the α s 's appearing in the formulae for the curves in Fig. 1 and the upper left, upper right and lower left panels in Fig. 6. In this way, all the ratios do not depend on m X . However, we note that in other parts of this paper and other plots, α s 's are evaluated differently: the α s 's appearing in the Sommerfeld factors in Eq. (42) are evaluated at βm X which is the typical scale of the momentum transfer of the soft-gluon exchanges which are responsible for the Sommerfeld effect [48], and we take β = 0.3, which is roughly the average of the thermal velocities of the X's at the freeze-out temperature; the α s 's in the bound-state formation cross sections and dissociation rates, as well as in the ζ part of the annihilation decay rates (given in Eqs. (23), (24), (29) and (30)), are evaluated at the bound-state inverse Bohr radius scale; the α s 's in the tree-level XX → gg, qq annihilation cross sections and the ones appearing explicitly in the bound-state annihilation decay rates, are evaluated at the scale of 2m X .  Figure 1: The thermally-averaged Sommerfeld corrected XX → gg, qq annihilation cross section (orange curve), bound-state formation cross section with (solid black curve) and without (dotted black curve) considering bound state dissociation and annihilation decay rates, as functions of E B /T , for the S3 (left panel) and F8 (right panel) cases. The purple curves have the same meaning as the corresponding black curves, but multiplied by a factor of 2 in the bound-state formation cross section (solid and dotted) and dissociation rate (solid). All curves are normalized to the tree-level thermallyaveraged XX → gg, qq annihilation cross section without Sommerfeld correction. The thin vertical black lines correspond to, from left to right, m X /T = 20, 30, 100. an S3 incoming pair feels a repulsive potential prior to forming a bound state, while the potential is attractive for an F8 incoming pair. At lower temperature, there is a smaller fraction of the incoming S3 pairs which have enough kinetic energy to overcome the repulsive potential to form bound states. While for the F8 case, a low temperature (and hence small velocities) favors the formation of bound states, similar to that of the Sommerfeld enhancement.
In order to see the bound-state effect on the evolution of the yield and on the result of final DM relic abundance, we plot in Fig. 2 the change ofỸ with E B /T in the WIMP coannihilation scenarios for the S3 (left panel) and F8 (right panel) cases, by assuming a WIMP DM with the number of degrees of freedom g χ = 2. For S3 (F8), we take m χ = 2 TeV (8 TeV), and the mass difference between DM and the colored particle 5 GeV (15 GeV). The solid red, orange, black and purple curves are results calculated from Eq. (43) without the Sommerfeld and bound-state effects, with the Sommerfeld effect but without bound-state effect, with both the Sommerfeld and boundstate effect, and with Sommerfeld effect and a factor of 2 enlargement of the bound-state effect (i.e., σv bsf → 2 σv bsf and Γ dis → 2 Γ dis ), respectively. In particular, the limiting values of the solid black curves at E B /T → ∞ give roughly the correct DM relic abundance consistent with the observational value, as we will show more in the next subsection by parameter scans. The equilibrium yield,Ỹ eq , is shown by a green dotted line. The thin vertical black lines correspond to, from left to right, m X /T = 20, 30, 100. As can be seen from the plots, the thermal freeze-out, defined when (Ỹ /Ỹ eq − 1) becomes order unity, happens between m X /T ∼ 20 and ∼ 30 9 . After freeze-out,Ỹ keeps decreasing to its limiting value at E B /T → ∞, which can be more than one order of magnitude smaller than its value at freeze-out, especially when considering the Sommerfeld and bound-state effects. At E B /T 1, the black and purple curves are very close to the orange curve, indicating that comparing to the Sommerfeld effect, bound-state effect is not important for that temperature range. However, at E B /T ∼ 1, the black and purple curves depart from the orange curve, and eventually result in significantly lower values ofỸ at E B /T → ∞ (the limiting values for the purple and black curves are ∼ 40% − 60% of the ones for the orange curve, for the two parameter sets shown). We can also see that the purple curves almost overlap with the black curves until E B /T ∼ 1, indicating that enlarge bound-state formation and dissociation by a factor of 2 has little effect on the yield at high temperature. These features are expected from the relative sizes of the Sommerfeld and bound-state enhancement factors illustrated in Fig. 1 10 .
To understand the evolution ofỸ more transparently, we also plot in Fig. 2 approximate solutions ofỸ with dashed red, orange, black and purple curves by using the fact thatỸ eq becomes negligible compared toỸ shortly after freeze-out, so that Eq. (43) can take the approximate form which has the solutioñ We plot the dashed curves by using Eq. (43) until x 1 = 30, at whichỸ is already much larger thañ Y eq , then we input the value ofỸ (x 1 = 30) in Eq. (55) and use it to plotỸ at larger values of x.
One can see that this approximation is very accurate, as the dashed curves completely overlap with the corresponding solid curves. It is easier to see from Eq. (55) that with a large σ ef f v for at least a range of not-too-large x (in other words, when n X has not been too diluted away by Hubble expansion), the contribution from the integration part can be large compared to 1/Ỹ (x 1 ). This explains why for the S3 case the bound-state effect on the result of the final relic abundance is still significant, even though the bound-state formation cross section decreases to zero at x → ∞.

Coannihilations in the WIMP DM scenario
Let us now study scenarios in which the WIMP DM χ has a mass close to a certain massive colored particle X, such that χ − X coannihilation is important in determining the DM relic abundance. We assume that the DM has degrees of freedom g χ = 2, corresponding to, for example, the Bino in the MSSM. The relic abundance of DM then depends only on the DM mass, m χ , and the mass splitting between DM and the colored particle, m X − m χ . We plot in Fig. 3 in the (m χ , m X − m χ ) planes the contour bands of the DM relic abundance falling within the 3-σ range of the Planck determination of the cold DM density, Ω CDM h 2 = 0.1193 ± 0.0014 [60]. These bands are calculated using Eq. (43) without the Sommerfeld and bound-state effects (red), with the Sommerfeld effect but without bound-state effect (orange), with both the Sommerfeld and bound-state effect (black), and with Sommerfeld effect and a factor of 2 enlargement of the bound-state effect (purple). We can see that for the S3, S8 and F8 cases, on top of the Sommerfeld enhancement, bound-state effects  further push upwards the largest mass splittings which can result in correct DM relic density. Also, the largest possible DM masses achieved at the endpoints of the coannihilation strips when the mass splittings approach zero, increase by ∼ 50%, 100% and 30% with respect to the Sommerfeldenhanced-only values, reaching ∼ 2.5, 11 and 9 TeV for the S3, S8 and F8 cases, respectively 11 . For the F3 case in the upper right panel, however, the Sommerfeld and bound-state effects are much smaller compared to the other three cases. As can be seen from the positions of the red and orange bands, the Sommerfeld effect gives a slightly suppressed rather than an enhanced σ ef f v 12 . Fig. 4 shows in the (m χ , Ω χ h 2 ) planes the locations of the endpoints of the coannihilation strips for different values of Ω χ h 2 , achieved when m X − m χ = 0, for S3 (upper left), F3 (upper right), S8 (lower left) and F8 (lower right) coannihilating with a WIMP DM which has g χ = 2. The color conventions are the same as in Fig. 3. The horizontal green band shows the 3-σ range determined by Planck, 0.1151 < Ω χ h 2 < 0.1235. We can see that for the S3, S8 and F8 cases, for a given value of m χ the Sommerfeld effect greatly reduces the calculated Ω χ h 2 compared to the one without the inclusion of Sommerfeld factors. Also, the calculated Ω χ h 2 is further significantly reduced after including the bound-state effect, in particular for a DM mass of TeV scale or larger. On the other hand, for the F3 case, again we see that the Sommerfeld and bound-state effects are small, and the Sommerfeld effect is opposite compared to the other three cases. 11 The numerical differences for the F8 case in Figs. 3 and 4 compared to Figs. 4 and 5 in [23] are due to a different use of α s in the bound-state formation and dissociation cross sections, as well as the effect from the squark masses in the tree-level cross section of gluino pair annihilation into quark-antiquark pairs. 12 The red and orange bands and curves in Figs. 3 and 4 are consistent with the red and light green bands and curves in Figs. 1 and 2 in [54] for the S3, S8 and F8 cases. For the F3 case, the red band and curve presented here are also consistent with the ones in [54], but the orange band and curve are different.

Implications of metastable particles on BBN and superWIMP DM abundance
It is well-known that the concordance of the standard BBN predications of the primordial lightelement abundances with the values inferred from observational data, provides strong constraints on the abundance, mass, lifetime, and decay spectra of a massive particle decaying during or after BBN (see e.g. [2]). Since our focus in this work is to study the impacts of the bound-state effect of massive colored particles, we want to see how much the bound-state effect can change the BBN constraints on the abundance and mass comparing to in particular the Sommerfeld effect, for a given lifetime and decay spectra which depend on other details of a specific particle theory model. With this in mind, we simply use the parametrization given in Eq. (56) of [59] for the BBN constraints obtained by [2] for a massive metastable particle with a lifetime of ∼ 0.1 − 10 2 sec and assuming that its hadronic decay branching ratio is 1 (which can be a good approximation for a massive colored particle), given as The above constraint comes from the would be overproduction of 4 He, due to new proton ↔ neutron interconversion reactions induced by the hadronic shower from X decays. In Eq. (56), Y X is the sum of the yield for particle and anti-particle in our convention. We show in Fig. 5 Y X as functions of m X for the S3 (left panel) and F8 (right panel) cases, calculated using Eq. (43) with the meanings of the variables understood as mentioned at the end of Sec. 3.2 for massive colored particles. As before, the red, orange, black and purple lines are results without the Sommerfeld and bound-state effects, with the Sommerfeld effect but without bound-state effect, with both the Sommerfeld and bound-state effects, and with the Sommerfeld effect and a factor of 2 enlargement of the bound-state effect, respectively. The blue dashed line is given by Eq. (56), and the parameter region above this line is excluded for an X with a lifetime of ∼ 0.1 − 10 2 sec. We can see that the bound-state effect pushes the allowed regions of m X to larger values compared to the ones with the Sommerfeld effect included only, namely, ∼ 1.1 → 2.1 TeV and ∼ 8 → 11 TeV for S3 and F8, respectively. As will be discussed more in the next section, since the LHC is pushing the exclusion limit of long-lived colored particles to TeV scale, it is useful to update the exclusion limit from BBN as well by including the previously omitted bound-state effect, so that we can be more confident to close or to leave open the mass window a long-lived colored particle can have.
We now consider the contribution to superWIMP DM relic abundance from the out-of-equilibrium decays of the massive colored particles. The gravitino [61,62] and axino [63,64] LSP in R-parity conserving supersymmetric models serve as good examples of superWIMP DM. Using Eq. (50) and including the Sommerfeld and bound-state effects, the contribution can be parameterized approximately as We note that in our calculations of the relic abundance of massive colored particles, we have not included the possible further reduction of the abundance due to annihilations of heavy exotic color-neutral hadrons, which are formed by not-yet-decayed massive colored particles together with quarks and gluons after the quark-hadron phase transition in the early Universe [65,66]. In this sense, the BBN constraints given here are conservative. The contribution to Ω non-th SW h 2 from the massive colored particle decays becomes smaller as well if there is further reduction of Ω X h 2 after the quark-hadron phase transition.

Electric charge corrections for the bound-state effect
So far we have focused on bound-state effects with a gluon being emitted/absorbed in the boundstate formation/dissociation process. If a massive colored particle also carries some electric charge, for example, the squark in supersymmetry, a bound state can form (or be dissociated) by emitting (or absorbing) a photon, i.e., X 1 X 2 ↔ ηγ. Also, the previously calculated bound-state formation/dissociation cross sections associated with gluon emission/absorption are modified due to the change of the potentials between the massive colored particles. To see the impacts of the electric charge on the bound-state effect, we use the S3 case as an example by assigning a charge Q (−Q) to S3 (S3), and consider the processes S3S3 ↔ ηg and S3S3 ↔ ηγ.
For S3S3 ↔ ηg, we still consider the transition between the (color-octet, L = 1, S = 0) free pair state and the (color-singlet, L = 0, S = 0) bound state. By modifying the coefficients of the Coulomb potentials, ζ → ζ + α EM Q 2 and ζ → ζ + α EM Q 2 , where α EM is the electromagnetic fine structure constant, the formulae given in Sec. 2 still apply. Quantities depending on ζ and/or ζ , e.g., a, E B , κ, and the cross sections and rates into which they enter, therefore all change. With an electric charge, for the bound state the previous attractive potential becomes more attractive. On the other hand, for the free pair state the previous repulsive potential becomes less repulsive, and even it can become attractive when |Q| is large enough, i.e., (−1/6)α s + α EM Q 2 is positive when |Q| 3/2. The bound-state annihilation decay rate changes as well due to the change of ζ.
For S3S3 ↔ ηγ, the bound state and the free pair state are in the same color state, so that ζ = ζ . Using dipole approximation 13 , we consider the transition between the (color-singlet, L = 1, S = 0) free pair state and the (color-singlet, L = 0, S = 0) bound state. The calculation is the same as for S3S3 ↔ ηg, except that there is no color factors to worry about and the explicit coupling factor α s in Eq. (10) is changed to α EM Q 2 . The bound-state dissociation and formation cross sections are where the superscript "γ" indicates photon emission/absorption. The quantities a, E B etc. are evaluated taking into account the change of potential due to the electric charge as mentioned above.
In the thermally-averaged dissociation rate given in Eq. (34), g g is changed to g γ = 2. The formula for the thermally-averaged formation cross section times relative velocity given in Eq. (35) stays the same.
We show in Fig. 6 where the superscript "g" indicates that only S3S3 ↔ ηg is considered. We can see that larger |Q| makes the bound-state effect stronger. Also, for a large enough |Q| such that the free pair potential becomes attractive (the blue dotted line and especially the brown dashed line), the behavior of the ratio at large E B /T becomes more like the black line in the right panel of Fig. 1 where the free pair potential is attractive. The upper right panel shows the ratio Again, larger |Q| leads to larger bound-state effect with photon emission/absorption. By comparing with the corresponding lines in the upper left panel, we see that for |Q| < 2, at E B /T 10 the bound-state effect due to gluon emission/absorption is much larger than the one due to photon emission/absorption, while they become comparable for larger |Q|. The lower left panel shows the ratio  in which the numerator is the one entering into the Boltzmann equation as the second term in Eq. (48). The shape of curves at smaller E B /T is controlled by the gluon emission/absorption bound-state effect, while the one from photon emission/absorption becomes important at larger E B /T . We plot in the lower right panel the locations of the endpoints of the coannihilation strips for different values of Ω χ h 2 , after taking into account the total impacts of electric charge on the boundstate effects. As before, the horizontal green band shows the 3-σ range of the Planck determination of the cold DM relic density, 0.1151 < Ω χ h 2 < 0.1235, and we assume a WIMP DM with the degrees of freedom g χ = 2. Comparing to the black line which corresponds to Q = 0, the |Q| = 1/3 and 2/3 cases (corresponding to the charges of squarks) only slightly increase the endpoint values of m χ for a given Ω χ h 2 , and by ∼ O(10) GeV on the Planck band. For larger |Q|, the increase becomes significant, and the endpoint on the Planck band reaches ∼ 3 (4) TeV for |Q| = 2 (3).

Collider constraints
As shown in the previous section, the DM relic abundance in the WIMP DM coannihilation scenarios and the BBN constraints on the long-lived massive particle decays impose upper bounds on the masses of exotic massive colored particles. On the other hand, collider experiments are constraining the masses from below.
In the massive colored particle coannihilating with a WIMP DM scenario, the coannihilation region is characterized by a small mass splitting, which typically is difficult to probe at the LHC using the conventional multiple jets plus large missing energy searches. Massive colored particle pair production accompanied by a hard initial-state radiation, i.e., the monojet search, is utilized to constrain scenarios with such a compressed mass spectrum. For simplicity, we focus only on the m X − m χ → 0 region so that the kinematics of decay products of the colored particles can be ignored 14 . The 13 TeV limits for S3 (stop) and F8 (gluino) are 0.32 TeV [67] and 0.63 TeV [68], respectively. To impose limits on F3 and S8, we perform monojet simulations [69,70] using the MadGraph-Pythia-Checkmate [71][72][73][74][75][76][77][78] pipeline. The obtained mass limits are 0.41 TeV and 0.43 TeV for F3 and S8, respectively. The results are summarized in Table 1, together with the endpoint values for the coannihilation strips including the Sommerfeld and bound-state effects read from Fig. 3.
Let us briefly comment on the discovery reach of a prospective 100 TeV proton-proton collider at an integrated luminosity of 3000 fb −1 , in particular for the S3 case. The bound-state effect increases the mass range of DM significantly, making the coannihilation scenario more difficult to probe. With the inclusion of the bound-state effect so that the right-handed stop-Bino coannihilation strip ends at ∼ 2.5 TeV when the Bino accounts for the total DM density (the inclusion of additional electroweak coannihilation channels [13] and a large lighter stop -heavier stop -Higgs coupling to heavier stop mass ratio [11] can shift the endpoint to even larger values), the ending part of the strip may be not within the discovery reach any more, though may be still within the exclusion reach given sufficiently low systematics [79].  In the long-lived massive colored particle scenario, the produced massive colored particles at a collider form R-hadrons and travel through the detector with velocities significantly less than the speed of light, leaving ionization energy dE/dx characteristically higher than that of charged SM particles. Searches for such events have been performed by both ATLAS and CMS collaborations at the LHC [80,81]. The current 13 TeV limits on long-lived (stable at collider scales) S3 (stop) and F8 (gluino) are 890 GeV and 1580 GeV, respectively [82]. While a dedicated simulation of R-hadron events at the LHC is beyond the scope of this work, in order to make a simple estimate of the bounds on long-lived F3 and S8, we assume that the signal efficiency and hadron formation probability of F3 (S8) equal to that of S3 (F8). The pair production cross section of F3 is estimated up to the next-to-next-to leading order (NNLO) using Hathor [83], while a next-to leading order (NLO) K-factor of 2 is used for the S8 scenario [84]. The obtained mass limits are 1.2 TeV and 1.4 TeV for F3 and S8, respectively. The results are summarized in Table 2, together with the upper bounds from BBN for the masses of long-lived massive colored particles assuming lifetimes between 0.1 and 100 sec, including the Sommerfeld and bound-state effects, as discussed in Sec. 3.5.

Summary
We have studied in this paper the bound-state effects of exotic massive colored particles on DM relic abundance calculations in scenarios where the massive colored particles coannihilating with a WIMP DM. In general the bound-state effect increases the effective annihilation cross section through the formation and then annihilation decays of bound states, draining the number of DM particles in the thermal bath, when the massive colored particles and DM share the same discrete symmetry which stabilizes the latter, and provided the interconversion rate between the two particle species is fast enough compared to the Hubble expansion rate. For a given DM relic abundance, this effect allows a larger DM mass and a larger mass splitting between the massive colored particle and the DM. As examples, we consider the massive colored particles being complex scalars (S3) or Dirac fermions (F3) in the color SU(3) fundamental representation, and real scalars (S8) or Majorana fermions (F8) region. The detailed study is left for future work.
in the adjoint representation. We find that the bound-state effect significantly increases the largest possible DM masses which can give the observed DM relic abundance, reaching ∼ 2.5, 11 and 9 TeV for the S3, S8 and F8 cases, respectively. Comparing to the corresponding ones when considering only the Sommerfeld effect but without the bound-state effect, these values increase by ∼ 50%, 100% and 30%, respectively. The increase for the F3 case is smaller, but still the bound-state effect can more than counterbalance the Sommerfeld effect which is a suppression rather than an enhancement in this case.
We note that while the potentials for the bound states are attractive, due to color charge conservation the potential for an incoming massive colored particle pair can be attractive, zero or repulsive. In the early Universe bound states can form when incoming pairs have sufficiently large relative velocities to overcome the repulsive potential. In particular, for the S3 case, we find that although fading at low temperatures, the large bound-state formation cross section achieved when temperatures are comparable to the bound-state binding energy makes the bound-state effect significant enough, such that to probe the entire stop-Bino coannihilation strip can be quite challenging, if possible, even in a prospective 100 TeV proton-proton collider at an integrated luminosity of 3000 fb −1 .
We have also calculated the corrections for the bound-state effect when the massive colored particles carry electric charges. Using the S3 case as an example, we find that larger electric charge makes the bound-state effect stronger, and the enhancement can make the above mentioned ∼ 2.5 TeV coannihilation endpoint to ∼ 3 (4) TeV for |Q| = 2 (3), for which the incoming pair potential changes from being repulsive to attractive. However, for |Q| < 1, the enhancement is quite small.
As we have briefly discussed, bound-state effects should also be included in calculations of su-perWIMP DM relic density from the decays of metastable massive colored particles, as well as when applying BBN constraints on long-lived massive colored particles. Furthermore, considering that BBN constraints and the DM relic abundance in coannihilation scenarios impose upper bounds on the masses of massive colored particles, we have studied the collider limits on the exotic massive colored particles we consider in this paper.
Before we close, we note that many other QCD bound states are possible. For example, in supersymmetry there can be di-squark and squark-gluino bound states [85]. Also, a squark and an antisquark with different flavors can form a bound state. Studying the effects of these bound states in specific supersymmetric models is left for future works.  The coefficients a and b in σv rel for massive colored particle pair annihilation to gg or qq, up to the common factor πα 2 s /m 2 X , for the S3, F3, S8 and F8 cases.
A Thermally-averaged annihilation cross sections for S3, F3, S8 and F8 We follow the procedure given in [56,86] in calculating the thermally-averaged s-and p-wave annihilation cross sections for a pair of massive colored particles into gg or qq, given as σv rel = a + b T /m X + O ((T /m X ) 2 ). Keeping results at leading order in α s , for qq final state we only need to evaluate the s-channel gluon exchange diagram, while for the gg final state we need to evaluate the s-channel gluon exchange, t-and u-channel X exchange diagrams, and the point interaction diagram when X is a scalar. For the S3 case, the relevant squared amplitudes are covered in the stop -antistop annihilation calculations, given in [5]. The results for the F8 case are covered in the gluino -gluino annihilation calculations, listed in [23]. The results for the F3 case are covered in the heavy quark -antiquark annihilation calculations in the SM. For the relatively less familiar S8 case, the relevant part of the Lagrangian is L = 1 [87], with the color index b summed over through 1 to 8. φ b is a real scalar filed in the color SU(3) adjoint representation, and D µ φ b = ∂ µ φ b + g s f abc G a µ φ c , where G a µ is the gluon field. Up to the common factor πα 2 s /m 2 X , the a and b terms in σv rel we found are listed in Table 3. The results for the qq channel have summed over all 6 types of SM quarks, and we have dropped quark mass dependent terms, as we are considering massive colored particles much heavier than the SM quarks.