Colored Dark Matter

We explore the possibility that Dark Matter is the lightest hadron made of two stable color octet Dirac fermions ${\cal Q}$. The cosmological DM abundance is reproduced for $M_{\cal Q}\approx 12.5$ TeV, compatibly with direct searches (the Rayleigh cross section, suppressed by $1/M_{\cal Q}^6$, is close to present bounds), indirect searches (enhanced by ${\cal Q}{\cal Q}+\bar{\cal Q}\bar{\cal Q}\to {\cal Q}\bar{\cal Q}+{\cal Q}\bar{\cal Q}$ recombination), and with collider searches (where ${\cal Q}$ manifests as tracks, pair produced via QCD). Hybrid hadrons, made of $\cal Q$ and of SM quarks and gluons, have large QCD cross sections, and do not reach underground detectors. Their cosmological abundance is $10^5$ times smaller than DM, such that their unusual signals seem compatible with bounds. Those in the Earth and stars sank to their centers; the Earth crust and meteorites later accumulate a secondary abundance, although their present abundance depends on nuclear and geological properties that we cannot compute from first principles.

1 Introduction Many models of particle Dark Matter (DM) have been proposed; one common feature is that DM is a new neutral and uncolored particle. We challenge this view: can DM be instead colored or charged, and be dominantly present today in the form of neutral bound states kept together by ordinary electromagnetic or strong interactions analogously to hydrogen or neutrons? The answer is no for electric binding: two charged particles with mass M m e form a negligible amount of neutral bound states, when their thermal relic abundance matches the DM cosmological abundance.
On the other hand, colored particles necessarily form hadronic bound states. We add to the Standard Model (SM) a new stable heavy colored particle Q, for simplicity neutral. Q could be a heavy quark in the 3 ⊕3 representation of SU(3) c , or a 'Dirac gluino' in the 8 ⊕ 8 representation, such that Q annihilates withQ, but not with itself. We dub this neutral quark as quorn. Perturbative annihilations and recombination between Q andQ leave a thermal relic density of order Ω Q h 2 ∼ 0.1 M Q /7 TeV. After the quantum chromo-dynamics (QCD) phase transition at temperature T < ∼ Λ QCD ≈ 0.27 GeV colored particles bind in hadrons. Subsequent annihilations among hadrons reduce their relic abundance, increasing the value of M Q such that DM has the observed cosmological abundance, Ω DM h 2 ∼ 0.1 for M Q ≈ 10 TeV.
The quorn-onlyum hadrons made of Q only (QQ if Q ∼ 8, and QQQ if Q ∼ 3) are acceptable DM candidates, as they have a small Bohr-like radius a ∼ 1/α 3 M Q . This scenario is believed to be excluded because it predicts other hybrid hadrons where Q binds with SM quarks q or gluons g. Such hybrids, Qqq, QQq, Qq (if Q ∼ 3) and Qg, Qqq (if Q ∼ 8), have size of order 1/Λ QCD and thereby cross sections of order σ QCD ∼ 1/Λ 2 QCD , can be charged, and are subject to strong bounds. Their cosmological abundance must be orders of magnitude smaller than the DM abundance Ω DM ≈ 0.1, while naively one might expect that cosmological evolution results into Ω hybrid Ω DM , given that quarks and gluons are much more abundant than quorns Q.
We will show that cosmological evolution gives Ω hybrid ∼ 10 −4 Ω DM , such that this scenario is allowed. This is not surprising, taking into account that quorn-onlyum has a binding energy E B ∼ α 2 3 M Q ∼ 200 GeV much larger than hybrids, E B ∼ Λ QCD . Quorn-onlyum thereby is the ground state, reached by the universe if it has enough time to thermalise. This depends on two main factors: i) quorns are much rarer than quarks and gluons: n Q ∼ 10 −14 n q,g when the DM abundance is reproduced; ii) QCD interactions are much faster than the Hubble rate H ∼ T 2 /M Pl : a loose bound state with a σ QCD cross section recombines N ∼ n q,g σ QCD /H ∼ M Pl /Λ QCD ∼ 10 19 times in a Hubble time at temperature T ∼ Λ QCD .
Since 10 19 is much bigger than 10 14 , chromodark-synthesis cosmologically results into quornonlyum plus traces of hybrids. This is analogous to Big Bang Nucleo-synthesis, that leads to the formation of deeply bounded Helium plus traces of deuterium and tritium.
The paper is organised as follows. In section 2 we define the model, and summarize the main features of its QCD interactions. In section 3 we discuss how cosmology leads to dominant formation of Q-onlyum hadrons. In section 4 we show that the abundance of hybrids is small enough to be compatible with bounds. In section 5 we show that Q-onlyum DM is compatible with bounds. A summary of our results is given in the conclusions in section 6.

The model
We consider the following extension of the SM: 1 The only new ingredient is Q: a Dirac fermion with quantum numbers (8, 1) 0 under SU(3) c ⊗ SU(2) L ⊗ U(1) Y i.e. a neutral color octet. The only free parameter is its mass M Q . Like in Minimal Dark Matter models [4] Q is automatically stable, as no renormalizable interaction with SM particles allows its decay, which can first arise due to dimension-6 effective operators such as QDDU and QLDQ where Q (L) is the SM quark (lepton) doublet, and U (D) is the right-handed SM up-type (down-type) quark. The decay rate is cosmologically negligible if such operators are suppressed by the Planck scale.
After confinement Q forms bound states. For M Q Λ QCD /α 3 states made by Q-only are Coulombian. The QQ bound states are unstable: Q andQ annihilate into gluons and quarks. No such annihilation arises in QQ bound states as we assumed that Q carries an unbroken U(1) dark baryon number that enforces the Dirac structure such that QQ is stable. The DM candidate is the quorn-onlyum QQ ground state, neutral, color-less and with spin-0. 2 As we will see, if QQ is a thermal relic, the observed cosmological DM abundance is reproduced for M Q ∼ 12.5 TeV. This mass is large enough that Q does not form QCD condensates. The QQ potential in the color-singlet channel is V (r) = −3α 3 /r, so the binding energy is E B = 9α 2 3 M Q /4n 2 ≈ 200 GeV/n 2 , which is bigger than Λ QCD up to n ∼ 20. We adopt the value Λ QCD ≈ 0.27 GeV.
The quantum numbers of the hybrid hadrons, Qg and Qqq , are not exotic. We expect that the isospin singlet Qg is lighter than Qqq (isospin 3 ⊕ 1) by an amount of order Λ QCD , 1 Within the SM, QCD could give rise to Dark Matter as 'strangelets' made of many uds quarks [2] or as 'sexaquark' uuddss [3]. However there is no experimental nor lattice evidence that such objects exist. We thereby extend the SM. 2 Other assignments of quantum numbers of Q are possible. A scalar would give similar physics. A fermionic Q ∼ (3 ⊕3, 1) 0 under SU(3) c ⊗ SU(2) L ⊗ U(1) Y would give the QQQ baryon as a viable DM candidate. As the gauge quantum numbers of a neutral color triplet are exotic, the QQq, Qqq and Qq hadrons containing light quarks would have fractional charges. Fractionally charged hadrons are subject to stronger experimental bounds [1]. A Q ∼ (3, 2, 1/6) = (Q u , Q d ), with the same quantum numbers of SM left-handed quarks Q, would give as lightest state the neutral DM candidate Q u Q d Q d . This is excluded by direct detection mediated at tree level by a Z, being a weak doublet with hypercharge Y = 0. Allowing for an additional confining group, a Q ∼ 8 can be build out of Q ∼ 3 obtaining double composite Dark Matter.
which accounts for the relative motion of q andq , where q, q = {u, d}. A lattice computation is needed to safely establish who is lighter. Assuming that Qqq is heavier, then its neutral component Qqq decays to Qg with a lifetime of order 1/Λ QCD . The slightly heavier components Qud and Qdū with electric charges ±1 have a lifetime of order v 4 /Λ 5 QCD . The above DM model has possible extra motivations. The fermion Q appears as a 'Dirac gluino' in some N = 2 supersymmetric models [5], where sfermions can mediate its decay, if R-parity is broken. Alternatively, the heavy quarks Q could be identified with those introduced in KSVZ axion models [6]. In such a case our U(1) symmetry gets related to the Peccei-Quinn symmetry. Corrections to the Higgs mass squared proportional to M 2 Q arise at 3 loops and are comparable to its measured value for M Q ≈ 10 TeV [7].

Confinement
QCD confinement happens in cosmology through a smooth crossover. In Cornell parametrisation [8] the QCD potential between two quarks in the F undamental representation at temperature T in the singlet configuration is approximated as V qq (r) ≈ −α F eff /r + σ F r. In the perturbative limit one has α F eff = C F α 3 where C F = (N 2 c − 1)/2N c = 4/3 is the quadratic Casimir and α 3 is renormalized around 1/r. At r ∼ 1/Λ QCD lattice simulations find α F eff = 0.4 and σ F ≈ (0.45 GeV) 2 [9]. The potential between two adjoints is similarly approximated by a Coulombian term plus a flux tube: Perturbation theory implies V QQ /C A ≈ V qq /C F [10] where C A = N c = 3. Thereby α eff ≈ 3α 3 and σ(0) ≈ 9σ F (0)/4 ≈ (0.67 GeV) 2 , as verified on the lattice [11]. At finite temperature the Coulombian force gets screened by the Debye mass and the string appears only below the critical temperature T c ≈ 170 MeV as σ(T ) ≈ σ(0) 1 − T 2 /T 2 c [9].

Eigenvalues in a linear plus Coulombian potential
We will need the binding energies of a non-relativistic QQ hadron. We thereby consider the Hamiltonian H = p 2 /2µ + V (r) in 3 dimensions that describes its motion around the center of mass, with reduced mass 2µ M Q . The potential is given by eq. (2). As usual, wave-functions are decomposed in partial waves as ψ(r, θ, φ) = ñ, ,m Rñ (r)Y m (θ, φ) whereñ is the principal quantum number. For each = 0, 1, 2, . . . we define asñ = 1 the state with lowest energy, so thatñ = 1, 2, 3, . . .. The radial wave function Rñ (r) hasñ − 1 nodes. Unlike in the hydrogen atom there are no free states: angular momentum is not restricted to <ñ. In order to match with the Coloumbian limit in its usual notation we define n ≡ñ + such that, at given , only n ≥ + 1 is allowed. The reduced wave function uñ (r) = rRñ (r) obeys the Schroedinger equation in one dimension in the effective potential V eff = V + ( + 1) 2 /2µr 2 . Rescaling arguments imply that Figure 1: Binding energies Eñ in GeV for a QQ in the singlet configuration. States with Eñ < −0.2 GeV (in green) are well approximated by the Coulombian limit. Increasing M Q leads to a larger number of Coulombian states and to a deeper ground state. QQ states are cosmologically mostly produced in the region with larger of the band E ∼ Λ QCD .
energy eigenvalues have the form From [12] we extract the approximation valid at leading order in ε 1 The first term is Coulombian. The second term accounts for the linear potential, and becomes relevant at large n, . In particular, assuming n 1, Coulombian states with negative binding energy exist up to < α 3/4 eff M 1/2 Q /(8σ) 1/4 . The ground state has binding energy E B = −E 10 ∼ 200 GeV for M Q ∼ 10 TeV.
In the opposite limit where the linear force dominates and the Coulomb-like force can be neglected, all energy levels are positive and states with higher have higher energy [12] Eñ ≈ such that thermalisation lowers . The dependence on σ, µ and the ground state energy can also be computed variationally, assuming a trial wave-function ψ(r) = e −r/rc /r 3/2 c , such that the typical size is r c ∼ (µσ) −1/3 . Fig. 1 shows the binding energies for relevant values of the parameters.
We next discuss a bound state B Q made of a heavy Q and a gluon. It cannot be described by non-relativistic quantum mechanics. Nevertheless, its binding energy can roughly be obtained by eq. (5) taking a small reduced mass µ ∼ √ σ. One then expects that such states are in their ground states at T < ∼ Λ QCD , and that their mass is M B Q = M Q + O(Λ QCD ).

Decay rates of excited bound states
Energy losses due to quantum decay of a QQ state with n, 1 into deeper states can be approximated with classical Larmor radiation. This holds in dipole approximation, where a state can only decay to = ± 1.
To see this, we consider a hydrogen-like system with V = −α/r and reduced mass µ. Assuming a circular orbit as in [13] one gets the emitted power having inserted the acceleration a = α/µr 2 and converted the orbital radius into n 2 times the Bohr radius as r = r n = n 2 /αµ. Similarly, the binding energy is E = −α/2r = −α 2 µ/2n 2 . At quantum level, a circular orbit corresponds to a state with maximal = circ = n. In dipole approximation such a state decays only to n = = n − 1, emitting a soft photon with energy ∆E Larmor = |E n − E n−1 | α 2 µ/n 3 , such that the decay rate is This matches the quantum decay rate.
Let us now consider a generic state. Classically, a generic elliptic orbit is parameterized by its energy E and by its angular momentum ≤ circ , where circ = α 2 µ/2E is the value corresponding to a circular orbit. The Larmor radiation power, averaged over the orbit, is Due to the larger acceleration at the point of minimal distance, the radiated energy for circ is much larger than in the circular case: this is why eē colliders are built circular. This classical result for non-circular orbits agrees with the quantum results for n, 1, summarized in appendix A for the hydrogen atom, which can be approximated as In the quantum computation the enhancement at small < n appears after summing over the available final states with small n ≥ − 1 which allows for energy jumps |E n − E n | larger than in the circular case.
In the opposite limit where the linear part of the potential dominates over the Coulombian part, energy losses of highly excited states are again well approximated by classical Larmor radiation, which does not depend on the shape of the orbit, given that the force does not depend on the radius: W Larmor = 8α eff σ 2 /3M 2 Q is negligibly small. This is confirmed by numerical quantum computations.

Cross section for formation of a loose QQ bound state
We here estimate the cross section σ tot (B Q + B Q → B QQ + X) for formation of a loose bound state containing two heavy quarks Q, starting from two bound states B Q containing one Q.
Assuming that B Q = Qg can be approximated as a Q and a gluon kept together by a flux tube with length ∼ 1/Λ QCD , the following geometrical picture emerges. The cross section is σ tot ≈ π 2 ℘ at energies E ∼ M Q v 2 < ∼ Λ QCD such that there is not enough energy for breaking the QCD flux tubes, and the recombination probability of two flux tubes is ℘ ∼ 1, like in string models. Independently from the above geometric picture, the size of the bound state is of order 1/Λ QCD , and thereby one expects a cross section σ QCD = c/Λ 2 QCD , with c ≈ π in the geometric picture. In the following we will consider c = {1, π, 4π}. For example the measured pp cross section corresponds to c ≈ 10.
While this expectation is solid at energies of order Λ QCD , at lower temperatures the cross section might be drastically suppressed if the residual van der Waals-like force has a repulsive component, which prevents the particles to come close enough. We will ignore this possibility, which would result into a higher abundance of hybrid relics.
More in general, processes that only require a small energy exchange E can have large cross sections of order 1/E 2 . 3

Cross section for formation of a un-breakable QQ bound state
We can finally compute the quantity of interest for us: the thermally averaged cross section σ fall (T ) for collisions between two Qg states which produce an unbreakable QQ hadron. This happens when the loose bound state discussed in the previous section radiates more energy than ∼ T in the time ∆t before the next collision, such that it becomes un-breakable and later falls down to its deep ground state.
In view of the previous discussion, we proceed as follows. A large total cross section σ QCD ∼ π/Λ 2 QCD needs a large impact parameter b ∼ 1/Λ QCD , and thereby the QQ state is produced with large angular momentum ∼ M Q vb.
The issue is whether a bound state with large gets broken or radiates enough energy becoming un-breakable [13]. As discussed in section 2.3, abelian energy losses are well approximated by classical Larmor radiation, and it is crucial to take into account that non-circular 3 The authors of [14] propose a quantum mechanical model where processes analogous to σ(B Q + B Q → B QQ + X) are computed in terms of cross sections suppressed by 1/M Q . This large suppression seems to derive from their arbitrary assumption that the cross section should be dominated by an s-channel resonance.  orbits radiate much more than circular orbits. The QQ potential is given by eq. (2), with a large α eff ≈ 3α 3 (μ) renormalized atμ ∼ 1/r ∼ Λ QCD . 4 The cross section for falling into an un-breakable QQ bound state is computed as follows. We simulate classical collisions, averaging over the velocity distribution at temperature T and over the impact parameter b. We numerically solve the classical equation of motion for the QQ system, starting from an initial relative distance b and an orthogonal relative velocity v. From the solution x(t) we compute the radiated energy ∆E by integrating the radiated power W Larmor ∼ 2α eff¨ x 2 /3 for a time ∆t. We impose ∆E > ∼ T where ∆t is the average time between two collisions at temperature T . We estimate it as ∆t ∼ 1/n π v π σ QCD where n π is the pion number density and σ QCD = c/Λ 2 QCD such that ∆t Λ 2 QCD /T 3 at T m π , while the pion density is Boltzmann suppressed at lower T .
The resulting σ fall (T ) is plotted in fig. 2, computed varying the uncertain QCD parameters as α eff , c = {1, π, 4π}. We see that even for α eff ∼ 1 the fall cross section σ fall (T ) equals to the total cross section σ QCD at temperatures below (0.1 − 0.3)Λ QCD , and it is mildly smaller at T ∼ Λ QCD . If instead α eff ∼ 4π one would have σ fall = σ QCD even at T ∼ Λ QCD . The value α eff ∼ 4π can account for non-perturbative QCD effects: it is not unreasonable to think that the bound state can quickly radiate the maximal binding energy E B ∼ 200 GeV by emitting in one shot a hundred of gluons with energy E ∼ 2 GeV each.
A rough analytical estimate for σ fall (T ) can be obtained as follows. As discussed above, states that radiate fast enough arise only in the Coulombian part of the potential. In view of eq. (8), their energy loss rate is W Larmor ∼ α 7 eff M 2 Q / 8 , which can be big enough only for where the order one numerical value was added by roughly fitting to fig. 2, for the values of the total QCD cross section there assumed. The fall cross section is only suppressed by a small power of M Q , explaining why we find a large σ fall ∼ σ tot for M Q ∼ 12.5 TeV. In the analytic estimate we neglected the fact that m π ∼ Λ QCD : this is taken into account by the relatively large ad hoc numerical factor added to eq. (10) such that it provides a better agreement with the numerical result in fig. 2 for M Q ∼ 12.5 TeV.

Cosmological relic densities
We can now compute how strong QCD interactions lead to an abundance of the Q-onlyum DM candidate QQ much larger than the severely constrained hybrid bound states Qg. We describe what happens during the cosmological evolution, from the usual decoupling of free Q at T ∼ M Q /25 (section 3.1), to recoupling (section 3.2) at T > ∼ Λ QCD , to T ∼ Λ QCD (section 3.3), to redecoupling at T < ∼ Λ QCD (section 3.4), to nucleosynthesis at T ∼ 0.1 MeV (section 3.5).

Q decoupling at T ∼ M Q /25
As usual, at T > ∼ M Q the free Q annihilate into SM particles much faster than the Hubble rate, remaining in thermal equilibrium until they decouple at T = T dec ≈ M Q /25, leaving the usual relic abundance, determined by their annihilation cross-section in this decoupling phase. The non-relativistic s-wave cross section reads where the strong coupling is renormalized around M Q , while it is renormalized around α 3 M Q in the Sommerfeld factors S n corresponding to the various color channels: We define Y Q ≡ (n Q + nQ)/s, where s is the entropy density, and assume no dark baryon asymmetry, n Q = nQ. Bound state formation gives an order one correction to the relic abundance, as discussed in [16] that considered Majorana gluinos. The bound states made by our 'Dirac gluinos' can be divided into stable QQ orQQ states that carry two units of dark baryon number, and unstable QQ states, where Q andQ annihilate. The latter come into spin-0 and spin-1 combinations, while the stable states have only the spin allowed by Fermi statistics: in particular the singlet ground state has spin 0. Among the unstable bound states the most relevant for the relic abundance at T Λ QCD are the ones that decay faster and have larger binding energy. These are listed in table 1. The corresponding effective rates are plotted in fig. 3. We only estimated the annihilation widths of those states that exist only as QQ; they are suppressed by O(α 2 3 ) making these states negligible (the formation cross section does not depend on spin) unless numerical factors compensate for the suppression.
These rates determine a network of Boltzmann equations for the abundance of free Q and for the abundances Y I = n I /s of the various bound states I as function of z = M Q /T . In the notations of [16] such equations are Here γ I is the thermal-equilibrium space-time density of formations of bound state I, related to the thermal average Γ Ibreak of the breaking rate Γ Ibreak as described in [16]. Furthermore Γ Iann is the decay rate of bound state I due to annihilations between its Q andQ constituents: it vanishes for the QQ andQQ states. Finally, Γ I→J = −Γ J→I is the decay rate from state I to state J. Taking into account that the annihilation and decay rates are much larger than the Hubble rate, ref. [16] used thermal equilibrium conditions to substitute the network of Boltzmann equations with a single equation for the total DM density, in terms of an effective annihilation rate γ eff ann . This strategy needs to be extended including the QQ andQQ states. Their annihilation rates Γ ann vanish, so we can now only reduce the network of Boltzmann equations to two equations: one for Y Q (density of free Q) and one for Y QQ = I∈QQ Y I (total density of stable bound states, that satisfies Y QQ /Y eq QQ = Y I /Y eq I for all stable states I). The equations are where γ eff ann includes the effects of QQ bound states and is given by the same expression as in [16]. The total fall rate that accounts for the cumulative effect of all QQ andQQ bound made of color S n states is given by the sum of the formation rates of all such states, γ fall = I∈QQ γ I , which equals n eq QQ Γ break ≡ I∈QQ Γ Ibreak n eq I . Notice that Y Q + 2Y QQ remains constant when a QQ bound state is formed.
We now derive an approximated analytic solution by computing the deviation from equilibrium of the stable bound states. First, we appreciate that at temperatures at which the quorn annihilation goes out of equilibrium the second of the above equations is still in equilibrium and thus the effect of stable bound states can be ignored in the solution for the first equation. The asymptotic solution in this phase is where z dec ≈ 25 and 1/λ = H/s| T =M Q . Expanding in small 1/λ one finds Y 1 QQ (z) and determines the temperature at which Y 0 This gives the asymptotic solution for Λ QCD T M Q : Using the specific rates for the main perturbative bound states listed in table 1 we obtain the values of Y Q and of Y QQ at temperatures T Λ QCD . The result is shown in fig. 4b, where they are denoted as 'perturbative'. We see that such effect can be neglected. At confinement, non-perturbative QCD effects force all free Q to bind with SM quarks and gluons to form strongly interacting hadrons, as discussed in the following.

Q recoupling at T > ∼ Λ QCD
DM annihilations recouple below the decoupling temperature T dec if the thermally averaged DM annihilation cross section σ ann (T ) grows at low temperatures faster than 1/T . In such a case DM recouples, and its abundance n DM is further reduced. A tree-level cross section T leads to order one effects, but not to recoupling (unless enhanced by some resonance). Formation of bound states with small quantum number n ∼ 1 give other similar effects. In the previous section we included such order one corrections, adapting the results of [16]. 5 At this stage Q can form relatively deep bound states with heavy quarks, which eventually decay.
The QCD coupling grows non-perturbative at T > ∼ Λ QCD giving a more dramatic recoupling effect: bound states with binding energy E Bn ∼ (α 3 /n) 2 M Q can be formed through a large cross section σ ann ∼ 1/E 2 Bn , having omitted powers of the strong coupling. The increase of the cross section as n → ∞ is tamed by a competing effect: only bound states with E Bn > ∼ T are actually formed at temperature T (as better discussed in appendix B), leading to a re-coupling cross section that grows as σ ann ∼ 1/T 2 for T > ∼ Λ QCD .

Chromodark-synthesis at T ∼ Λ QCD
This effect culminates after confinement. Cosmological effects of confinement begin when the Coulombian force α eff /r 2 becomes weaker than the string tension σ(T ) at the typical distance r ∼ 1/T . Given that gluons and quarks are much more abundant than Q, the free Q form Qg and Qqq bound states, which have a binding energy of order Λ QCD and scatter among themselves and with other hadrons with cross sections of typical QCD size, σ QCD = c/Λ 2 QCD with c ∼ 1. In this stage H ∼ Λ 2 QCD /M Pl ∼ 10 −20 Λ QCD , such that a Qg hadron experiences 10 20 QCD scatterings. Given that the relative abundance of Q is Y Q ∼ 10 −14 , two Qg will meet, forming either deep QQ hadrons (which remain as DM) or QQ hadrons (which annihilate into SM particles). The abundance of Q-only hadrons gets dramatically suppressed, until they decouple.
While most DM particles form in this phase, a precise description is not needed to compute the final abundances, which are dominantly determined by what happens during the final redecoupling, where the dominant SM degrees of freedom are semi-relativistic pions, while the baryon abundance is negligible, in view of the Boltzmann factor e −mp/T and of the small asymmetry.

Q redecoupling at T < ∼ Λ QCD
We need a precise description of the final redecoupling which occurs at temperatures of tens of MeV. One might think that the simplified Boltzmann equations for the density of free Q and of QQ bound states, eq. (14), can be replaced with corresponding equations for the total density of B Q bound states (Qg and Qqq ) and for the total density of B QQ bound states. A slightly different strategy is needed. Indeed, the simplification that allowed to reduce the network of Boltzmann equations (one for each bound state) to two is valid under the following conditions: all B Q bound states are in thermal equilibrium among them; all B QQ bound states are in thermal equilibrium among them. Bound states are subject to QCD interactions, with large σ QCD cross sections, such that the corresponding interaction rates are much faster than the Hubble rate. However, as discussed in section 2.3, non-perturbative QCD interactions now lead to the formation of a large variety of bound states, with large n and quantum numbers which suppress the decay rates among them. Some decay rates can be slower than the Hubble rate. This issue was solved in section 2.5 where we computed an effective cross section for the formation of all unbreakable QQ bound states, that later fall to the QQ ground state. The same cross section, almost as large as the QCD cross section, holds for the formation of unbreakable QQ, that later annihilate: The equality of the classical non-perturbative total cross section for forming QQ bound states with the total cross section for forming QQ bound states, is compatible with the perturbative quantum cross sections computed in section 3.1. Indeed, because of Fermi anti-symmetrisation in the QQ case cross sections are twice bigger, while the number of QQ states is twice bigger (after restricting to colour-singlet bound states and averaging odd with even ). One extra process can take place: annihilations between QQ andQQ in their ground states. In section 5.2 we will compute its cross section, finding that it can be neglected in our present cosmological context. Together with eq. (18) this implies a simple result: half of the Q and Q present before redecoupling annihilate, and half end up in our DM candidates, the QQ and QQ ground states. Boltzmann equations are only needed to compute how small is the residual fraction of Q in loose hybrid hadrons, which are phenomenologically relevant in view of their large detection cross sections.
We thereby group bound states in two categories. We define Y QQ as the density of all unbreakable QQ bound states, produced with cross section σ fall . We define Y Q as the density of Q in loose bound states: the Q in bound states containing a single Q (Qg, Qqq ), and those in loose QQ and QQ bound states at relative distances ∼ 1/Λ QCD , that get broken by QCD scatterings.
The relevant Boltzmann equation are: valid for T < ∼ Λ QCD i.e. z > ∼ z QCD ≡ M Q /Λ QCD . In the non-relativistic limit the space-time density of interactions is determined by the cross sections as 2γ (n eq B Q ) 2 σv rel . The asymptotic solutions to the this system of equations are with the last term roughly equals Y Q (z QCD )/4. Fig. 4 shows our final result: the DM abundance and the hybrid abundance as function of the only free parameter, M Q . The left panel shows the mass abundances Ω = ρ/ρ cr ; the right panel shows the number abundances Y = n/s. The hybrid abundances are plotted as bands, given that they are affected by QCD uncertainties; smaller values are obtained for larger c = σ QCD Λ 2 QCD and for larger α eff . Varying them between 1 and 4π, the hybrid abundance changes by a factor 100. The DM abundance, less affected by QCD uncertainties, is plotted as a blue curve. The right panel shows that the DM QQ abundance is mostly made at non-perturbative level; the perturbative bound states computed in section 3.1 only play a significant role in enhancing QQ annihilations.
The observed DM abundance is reproduced for M Q ≈ (12.5 ± 1) TeV (21) and the hybrid mass abundance is about 10 4 smaller that the DM abundance (between 10 3 and 10 5 within our assumed range of QCD parameters). For such mass, fig. 5 shows the cosmological evolution of the abundances. It also shows how large uncertainties at T ∼ Λ QCD before redecoupling have a negligible impact on the final abundances, which is dominantly determined by redecoupling. An analytic argument that shows that Ω hybrid Ω DM is unavoidable and that gives the dependence of the final abundances on M Q , M Pl , Λ QCD (eq. (57)) is confined to appendix B because it follows a logic different from the one used in the more accurate numerical computation presented here.

Nucleodark-synthesis
Redecoupling is completed at temperatures T ∼ 10 MeV. Later nucleons bind into light nuclei at the Big Bang Nucleosynthesis (BBN) temperature T BBN ∼ 0.1 MeV. Various authors tried to compute what happens to SIMPs during BBN, and how SIMPs affect ordinary BBN [18][19][20][21] 6 . Our predicted amount of Strongly Interacting Massive Particles, Y SIMP ∼ 10 −18 , has negligible effects on ordinary BBN, which constrains Y SIMP < ∼ 10 −12 . Such studies however disagree on what happens to SIMPs during BBN. Do SIMPs bind with (some) nuclei? Does a significant fraction of SIMPs remain free?
We present our understanding, but we cannot provide a safe answer. Indeed, nuclear forces are not understood from first principles, not even for ordinary p and n [22]. Long-range nuclear properties are determined by couplings to pions, known thanks to chiral perturbation theory [23]. Heavier QCD states contribute to short-range nuclear forces: however QCD is here only used as inspiration to write phenomenological nuclear potentials to be fitted to p, n data, see e.g. [24].
In our scenario there are two types of SIMPs with distinct properties. The Qg hybrids are a isospin singlet and thereby do not couple to pions. The Qqq hybrids form an isospin triplet (with charges 0, ±1) coupled to pions.
Presumably Qqq are heavier and decay promptly to Qg. Then, the Qg singlet states, which do not feel the pion force, are expected to behave similarly to the Λ baryon, which does not bind to protons to form heavy deuterons [25]. Maybe such SIMPs do not bind with any nuclei, or maybe they find a way to form bound states with big enough nuclei. An attractive force can be provided by exchange of an isospon-singlet scalar meson, such as the σ (mass M ∼ 0.6 GeV) or glueballs (mass M ∼ 1.5 GeV) provided that their effective Yukawa couplings y SIMP and y N to the SIMP and to nucleons are large enough and have the same sign. In spherical well and Born approximation and for M Q M , the hybrid can form a bound state in a nucleus with atomic number A if [26] y SIMP y N > 12π If SIMPs bind to light nuclei, after BBN they dominantly end up in Helium or free, with a relatively large amount in Beryllium, according to [19,20]. The Qqq states, which feel the pion force, have an interaction potential of approximately 2 fm. If they are the lighter stable bound states, during BBN they get incorporated into nuclei with an efficiency close to 100% [21]. In the Milky Way, SIMPs in charged nuclei can loose a significant fraction of their energy by interactions with ambient matter.
No SIMP searches have yet been performed in galactic clouds, which would probe the SIMP primordial abundance. After BBN, SM matter forms stars and planets: primordial SIMPs sink to their center before that these objects possibly solidify. Stars (rather than BBN) later produce the observed elements heavier than He. In the next section we estimate the present geological abundance of SIMPs.

Signals of relic hybrid hadrons
In our model Q-onlyum DM is accompanied by hybrid hadrons, containing heavy colored Q bound together with SM quarks or gluons. In this section we discuss their signals. While SIMP DM has been excluded long ago, in our model SIMPs have a sub-dominant abundance, f SIMP ≡ ρ SIMP /ρ DM below 10 −3 , possibly a few orders of magnitude smaller. Such small value of ρ SIMP makes indirect SIMP detection signals negligible (f 2 SIMP σ QCD < ∼ 10 −24 cm 3 /sec) despite that SIMPs interact with matter nucleons and with themselves through large cross sections of order σ QCD ∼ 1/Λ 2 QCD . See also [28]. In some models SIMPs can have electric charge (fractional in exotic models).
As discussed in section 4.1, galactic SIMPs are stopped by the upper atmosphere of the Earth and slowly sink. Thereby SIMPs are not visible in direct detection experiments performed underground. Their later behaviour depends on whether SIMPs bind with nuclei: if yes they indirectly feel atomic forces; otherwise they sink even within solid bodies, such as the present Earth. In section 4.2 we summarize bounds on the SIMP abundance, to be compared with their present abundance, estimated in sections 4.3 and 4.4.

Direct detection of hybrid hadrons
Despite their reduced abundance, SIMPs would be excluded by a dozen of orders of magnitude, if they reach the underground direct detection detectors with enough energy to trigger events. This is not the case. The energy loss of a neutral SIMP in matter is [29] where n A is the number density of nuclei with atomic number A and mass m A ≈ Am p ; 2m A /M Q is the fractional energy loss per collision and σ A ≈ σ p A 2 (m A /m p ) 2 is the SIMP cross section on a nucleus [30], written in terms of the SIMP scattering cross section on protons, σ p ≈ π/Λ 2 QCD ≈ 1.6 10 −26 cm 2 . The cross section σ A is coherently enhanced at the energies of interest for us, The densities n A in the Earth crust can be written as n A = f A ρ/m A where ρ is the total mass density and f A is the mass fraction of material A, A f A = 1. The energy loss following from eq. (23) is Thereby SIMPs with M Q ≈ 10 TeV thermalize in the Earth atmosphere, which has a column depth of 10 4 kg/m 2 and A 4 1/4 ≈ 16.6, before reaching the crust with A 4 1/4 ≈ 31 and density ρ ≈ 3 g/cm 3 . SIMPs do not reach direct detection experiments, situated about a km underground. Some direct detection searches have been performed by balloon experiments at high altitudes. The authors of [31] claim that it is questionable whether such experiments exclude a SIMP with density ρ SIMP = ρ DM . Our predicted abundance ρ SIMP ∼ 10 −4 ρ DM is allowed.
After thermalisation, SIMPs diffuse with thermal velocity v thermal ≈ 6T /M Q ≈ 40 m/s at temperature T ≈ 300 K. In the Earth gravitational field g = 9.8 m/s 2 , SIMPs not bound to nuclei sink with a small drift velocity that can be estimated as follows. Each collision randomises the SIMP velocity because v drift v thermal . Thereby the drift velocity is the velocity v drift ≈ gτ /2 acquired during the time τ ≈ d/v thermal between two scatterings, where d = 1/( A n A σ A ) ∼ 0.1 mm in the Earth crust. Thereby the sinking velocity is The mean free path in Earth of a SIMP in a charged state is thereby L ± ∼ M Q β 4 /K ∼ 2 10 −5 cm (β/0.001) 4 . Again, SIMPs do not reach underground detectors. The main difference is that SIMPs bound in nuclei sink in the ocean and in the primordial Earth, but not in the solid crust, where electric atomic forces keep their positions fixed on geological time-scales.

Searches for accumulated hybrid hadrons
Experimental searches for accumulated SIMPs consist in taking a sample of matter, and searching if some atom has an anomalous mass or charge, see [35] for a recent review. The results, detailed below, imply relative abundances smaller than O(1/N A ) (inverse of the Avogadro number) in the selected samples.
The searches often involve a first phase of sample enrichment in hybrids (for example centrifuge treatment of a sample of water, or use of radioactive materials), followed by a second  Table 2: Experimental bounds on the density of Strongly Interacting Massive Particles with non-exotic electric charges, compared to the expected abundance of our hybrid, roughly estimated assuming that it binds in nuclei (otherwise they sink), and assuming f SIMP ≈ 10 −5 .
phase of hybrid detection, with the most successful being the mass spectroscopy and Rutherford backscattering [34]. Limits on the SIMP fraction in the sample depend on the SIMP mass: in the range GeV to TeV best bounds are derived from mass spectroscopy of enriched sea water samples [36]. Here the hypothetical particle is a positively charged SIMP, which could form heavy water replacing a proton. The bounds on the relative abundance are of order N SIMP + /N N < 10 −27 where N N is the number of nuclei.
For heavier SIMPs, mass spectroscopy seems to provide weaker limits. Stringent limit stems from studies of material from meteorites. In [34] the Rutherford backscattering technique was used to set a limit on the SIMP-to-nucleon number density in the tested meteorites that covers the range 100 GeV < M SIMP < 10 7 GeV. This technique does not depend on the SIMP charge and thus also applies to neutral SIMPs. For M SIMP ∼ 10 TeV the limit is [34] N SIMP N n < ∼ 3 10 −14 10 TeV M SIMP (meteorites) (27) where N n is the number of nucleons.
These bounds should be compared with the predicted SIMP abundance in the selected samples. If the tested samples were representative of the average cosmological composition, our model would predict having used the cosmological density of baryonic matter, Ω b h 2 ≈ 0.022, and of DM, Ω DM h 2 ≈ 0.12. The predicted abundance in the selected samples is much lower than the cosmological average and depends on their geological history.

Abundance of hybrid hadrons in the Earth
Testing a sample of sea water does not lead to bounds, because the atoms that contain heavy hybrid hadrons sink to the bottom. Similarly, the Earth once was liquid, so that the primordial heavy hybrids sank to the core of the Earth. 7 Objects made of normal matter accumulate SIMPs due to collisions with SIMP relics in the interstellar medium. Heavy hybrids accumulated in the Earth crust, if captured by nuclei, presumably stopped sinking after that the crust solidified. In order to set bounds, we thereby consider the SIMPs captured by the Earth in the time ∆t ∼ 4 Gyr passed since it is geologically quasi-stable. We ignore convective geological motion. The Earth is big enough to stop all SIMPs, so that the total mass of accumulated SIMPs is having inserted the escape velocity from the Galaxy v ∼ 10 −3 and assumed that the SIMP galactic density follows the DM matter halo density ρ DM ≈ 0.3 GeV/ cm 3 as n SIMP = f SIMP ρ DM /M SIMP . The rate of QQ annihilations of stopped SIMPs is negligible, because suppressed by e −M Q r where r is the macroscopic distance between Q andQ. 8 The number of SIMPs accumulated in the Earth is If SIMPs are not captured by nuclei and sink as in eq. (25), their present density in the crust is negligibly small, N SIMP /N n ∼ 10 −23 . If SIMPs get captured in nuclei, a significant fraction of such SIMPs could be in the crust, with a local number density higher by some orders of magnitude. In fig. 4 we plot the bound from Earth searches assuming that all SIMPs stop in the atmosphere and sink slowly through earth until captured by a nucleus, which might happen in the upper 10 km. The capture cross section with nuclei is discussed below.

Abundance of hybrid hadrons in meteorites
Meteorites result from accumulation of interstellar dust and contain heavy elements. The tested meteorites consist mainly of carbon and/or iron. These elements have not been produced by Big-Bang-Nucleosynthesis, which produced H and He (Z ≤ 2), nor by cosmic ray fission, which produced Li, Be, B (Z ≤ 5). Heavier elements have been synthesized from nuclear burning in stars and have later been dispersed away through various explosive processes: corecollapse supernovae, accretion supernovae, merging neutron stars and r-process nucleosynthesis.
Primordial SIMPs would have sunk to the center of stars, and would have presumably remained trapped there, undergoing QQ annihilations. Thereby, the SIMP relative abundance in meteorites made of heavy elements is expected to be significantly smaller than the average relative cosmological abundance.
In order to set bounds we compute the amount of SIMPs accumulated in meteorites. Meteorites are the oldest objects in the solar system and are so small that heavy hybrids do not sink in them. While the Earth is large enough that it captures all SIMPs intercepted by its surface, we consider meteorites small enough that the opposite limit applies: SIMPs are captured by all nuclei within the volume of the meteorite. Thus we need to estimate the probability ℘ that a nucleus captured a SIMP in a time ∆t: This value is roughly two orders of magnitude above the meteorite bound in eq. (27).
However, the capture cross sections of SIMP by nuclei are very uncertain. Taking into account that they are not coherently enhanced, the maximal value is the area of the nucleus, σ capture ∼ A 2/3 /Λ 2 QCD [37]. The measured capture cross sections of neutrons by nuclei are smaller: in most cases σ capture ∼ 0.01/Λ 2 QCD at MeV energies. Assuming this capture cross section we obtain the possible meteorite bound plotted in fig. 4 and summarized in table 2. Our SIMPs have MeV energies, but the longdistance attractive force mediated by pions (present for neutrons, where it is the only effect understood from first principles) is absent for Qg SIMPs, which are isospin singlets. Their capture cross section could be much smaller, and possibly our SIMPs do not form bound states with nuclei, such that meteorite bounds are not applicable.

Neutrinos from SIMP annihilations in the Sun
DM accumulates in the center of the Sun and annihilates to neutrinos, giving a detectable signal in IceCube [38]. Given that equilibrium holds between DM capture and DM annihilation, the neutrino rate depends on the cross section for DM direct detection. The IceCube bounds are weaker than those from direct detection experiments, and satisfied in our model [38]. Annihilations of SIMPs accumulated in the center of the sun provide an extra neutrino signal. The capture rate does not depend on the SIMP cross section, given that it is so large that all SIMPs that hit the Sun get captured, such that where R sun ≈ 7 10 8 m is the solar radius. Around the relevant mass, IceCube provides the bound Γ ann < ∼ 7 10 20 sec −1 on DM annihilating to bb [38]. Our Q dominantly annihilates to gluons and quarks, providing a slightly smaller neutrino flux [39]. We thereby conclude that the IceCube bound is satisfied even assuming a SIMP annihilation rate in equilibrium with the capture rate, Γ ann ≈ Γ capt /2.

Dark matter signals
In our model DM is a QQ hadron. In this section we discuss the DM signals: direct detection (section 5.1), indirect detection (section 5.2) and collider (section 5.3).

Direct detection of DM
Direct detection of DM is a low energy process, conveniently described through effective operators. Composite DM gives operators which can be unusual with respect to those characteristic of elementary DM with tree-level-mediated interactions to matter. For example, a fermionic bound state can have a magnetic dipole moment, which is strongly constrained. In our case DM is a non-relativistic scalar bound state QQ made of two colored neutral fermions Q. Its dominant interaction with low-energy gluons is analogous to the Rayleigh scattering of photons from neutral hydrogen. Describing our QQ bound state as a relativistic field B with canonical dimension one, the effective Lagrangian is The first expression employs the conventional basis of operators where (G a µν ) 2 = 2( B a2 − E a2 ) and O g µν ≡ G aρ µ G a νρ − 1 4 η µν G a ρσ G aρσ . In the second expression we rewrote them in terms of the chromo-electric E a i = G a 0i and chromo-magnetic B a components, such that c E is 4π times the chromo-electric polarizability of the bound state, c E ∼ 4πa 3 where a = 2/(3α 3 M Q ) is its Bohr-like radius. Furthermore c B c E is suppressed by the velocity v ∼ α 3 of the Q in the bound state. Neglecting the chromo-magnetic interaction, the coefficients renormalized at the high scale (that we approximate with M Z ) are The low energy effective coupling of DM to nucleons is where f g = 0.064 and g(2, M Z ) = 0.464. The spin-independent direct detection cross-section is This is close to the Xenon1T bound [41], σ SI < ∼ 3 10 −44 cm 2 × M DM /20 TeV, that holds at M DM 100 GeV up to the standard assumptions about the DM galactic halo.
Thereby we perform a dedicated computation of the c E coefficient, which is possible in perturbative QCD. Adapting the techniques developed for the hydrogen atom and for bottomonium [42], the effective Lagrangian of eq. (34) also describes the shift in the QQ ground-state energy induced by external chromo-electric and chromo-magnetic fields: The external field E a adds a chromo-dipole interaction to the non-relativistic Hamiltonian of the QQ bound state, as well as the associated non-abelian effects. Perturbation theory at second order then gives a shift in the ground state energy E 10 , which allows to reconstruct c E as where |B is the QQ ground state, N c = 3 and C is the Casimir coefficient, defined by Cδ ij = (T a T a ) ij and equal to 3 for our assumed octet representation. Summing over all allowed intermediate states with free Hamiltonian H 8 in the octet channel we find (see appendix C) where the first (second) contribution arises from intermediate bound (free) states. The nonabelian nature of QCD manifests in the fact that the allowed intermediate states are p-wave color octets: they are less bound (relatively to the ground state) than in the hydrogen atom case, such that our c E coefficient is significantly smaller than what would be suggested by a naive rescaling of the abelian result. Eq. (41) is the coefficient used as a reference value in the cross section of eq. (38). Higher order QCD interactions and relativistic effects conservatively amount up to a 50% uncertainty. As plotted in fig. 6a our predicted DM mass M DM ≈ 25 TeV is higher than the DM mass excluded by direct detection, M DM > ∼ 14 TeV.

Indirect detection of DM
Two DM particles in the galactic halo can annihilate into gluons and quarks giving rise to indirect detection signals. The energy spectra of the resulting final-state stable particles (p, e, γ, ν) is well approximated by the general results of non-relativistic annhilations computed in [44]. We need to compute the annihilation cross section between the DM = QQ Coloumbian bound state and DM =QQ. It is enhanced and dominated by the recombination process (QQ) + (QQ) → (QQ) + (QQ) (42) followed by later QQ annihilations to SM particles. This is similar to what happens for hydrogen/anti-hydrogen annihilation, which proceeds through recombination (ep) Figure 6: Left: Direct detection signals of QQ dark matter, as computed in section 5.1. We also show the neutrino floor, which will eventually limit future direct searches. Right: Indirect detection signals as computed in section 5.2. We show the current dwarf galaxy constraints by FermiLAT, which have only a mild systematic uncertainty due to the dark matter J-factor, and the future sensitivity of the CTA [43] experiment to photons from dwarf galaxies.

�������� ���������
(eē) + (pp) followed by later eē and pp annihilations, giving rise to a large σ ann , of atomicphysics size, rather than of particle-physics size, σ ann ∼ α 2 /m 2 e,p . The DM recombination cross section can be estimated as follows. At relative velocities v rel comparable to the orbital velocity α 3 (which corresponds to a center-of-mass kinetic energy K = M Q v 2 rel /2 comparable to the binding energy E B ) the cross section is given by the Bohr radius squared, At larger velocities v rel > ∼ α 3 the cross section gets progressively suppressed by (E B /K) 2 , reducing at v rel ∼ 1 to the particle-physics cross section σ ann ∼ α 2 3 /M 2 Q . At smaller v rel α 3 recombination is enhanced by a classical Sommerfeld effect which can be estimated as follows. The interaction between two neutral atoms at distance r a is given by the non-abelian Van der Waals electric attraction, V el ≈ −0.7a 6 /r 7 [45,42,46], having used eq. (41) for the numerical coefficient. impact parameter b max , and thereby the cross section 9 At astrophysically low velocities v rel ∼ 10 −3 < ∼ α 5/2 3 the magnetic dipole interaction V mag ∼ α 3 /r 3 M 2 Q becomes as important as the electric interaction, giving σ ann v rel ∼ α rel . Detailed quantum computations suggest that a reasonable estimate consistent with Wigner's threshold law is [47] σ ann v rel ∼ πa 2 v rel /2 The 1/α 3 enhancement is numerically mild, given that α 3 ∼ 0.1. As a consequence indirect detection signals are below present bounds, as shown in fig. 6b. We plotted bounds on gamma ray emission from dwarfs, given that searches in the galactic center region are subject to large astrophysical uncertainties, and other bounds are weaker.

Collider signals of DM
While DM usually gives missing-energy signals which are hardly detectable at hadron colliders, DM made of colored quorns Q gives very visible signals. Indeed, DM constituents Q are pair 9 A more precise result can be obtained from a classical computation. Focusing on the color singlet channel, we numerically compute the classical motion of a QQ bound state in its ground state (circular orbit with radius a in some plane) which collides with relative velocity v rel and impact parameter b with a similarQQ system. When the two bound states get closer and interact they can produce two QQ bound states, which later annihilate. Confinement takes place at larger distances and plays a negligible role. By averaging over the relative orientations of the two systems and over the impact parameter gives the classical probability for this process, encoded into a velocity dependent cross section. produced at colliders via QCD interactions. After hadronization they form hadrons. Presumably the neutral Qg is stable, and the charged Qqq are long-lived on collider time-scales, giving rise to tracks. This is the dream of LHC experimentalists [48]. Experiments at the LHC pp collider at √ s = 13 TeV set the bound M Q > ∼ 2 TeV [49]. A larger √ s ∼ 85 TeV is needed to discover the quorn with the mass expected from cosmology, M Q ∼ 12.5 TeV. A pp collider with √ s = 100 TeV would be sensitive up to M Q < ∼ 15 TeV [50], as long as the detector can see the signal.
Furthermore, we explore the possibility of detecting collisions of protons in collider beams with ambient QQ DM. The QQ binding energy is E B ∼ 200 GeV. Protons with energies much larger than E B see the QQ system as two free Q and the QCD cross section is suppressed by the energy squared. Protons with energies comparable to E B see the system as a ball with Bohr radius a = 2/3α 3 M Q . The cross-section for the excitation of the ground state through the absorption of a gluon can be estimated as the cross-section for ionization computed in [16,51] where E g is the gluon energy and ζ = α 3 /v rel = 1/(3ap DM ) parametrises the momentum of Q in the final state. Energy conservation implies E g ≈ E B + M DM v 2 rel /4. Fig. 7 shows the proton-DM cross-section obtained convoluting with parton distribution functions. The event rate in a beam containing N p protons is small, QQ dark matter excitation by cosmic rays is negligible on cosmological time-scales.

Conclusions
We have shown that Dark Matter can be obtained from a colored neutral quark Q (dubbed quorn), that, after the QCD phase transition, forms deeply bound hadrons made of Q only (dubbed quorn-onlyum), plus traces of hybrid hadrons made of Q together with SM gluons or quarks (dubbed Strongly Interacting Massive Particles or SIMP). We explored the simplest model, where Q is an automatically stable neutral Dirac fermion in the adjoint representation of SU(3) c . Such a state could be a Dirac gluino, or appear in natural axion models (see section 2). Fig. 5 shows the cosmological evolution of the DM and hybrid abundances for the value of the quorn mass, M Q ≈ 12.5 TeV, which reproduces the DM cosmological abundance as discussed in section 3. A first decoupling occurs, as usual, at T ∼ M Q /25. Quorns recouple while the universe cools approaching the QCD phase transition at T ∼ Λ QCD . This opens a phase of chromodark-synthesis: quorns fall into QQ singlet bound states, which have a binding energy E B ∼ 200 GeV. The cross sections grow large, up to σ QCD ∼ 1/Λ 2 QCD , because excited states with large angular momenta are formed. Such states efficiently cool falling to the ground state before being broken, as computed in section 2.5 where we show that quantum states with n, 1 are well approximated by classical physics. It is important to take into account that (non-abelian) Larmor radiation from elliptic orbits is much larger than for circular orbits.
Details of this uncertain phase are not much important for the final result: one half of free quorns annihilate, one half end up in QQ DM; the small residual abundance of Qg hybrids, ρ SIMP /ρ DM between 10 −3 and 10 −6 , is mostly determined at T ∼ 30 MeV, when the states decouple again.
In section 5 we studied DM phenomenology. The quorn-onlyum DM state QQ with mass M DM ≈ 2M Q ≈ 25 TeV has small residual interactions suppressed by powers of 1/M Q . The cross section for direct DM detection is of Rayleigh type, suppressed by 1/M 6 Q . In section 5.1 we performed a non-trivial QCD bound-state computation, finding a cross section just below present bounds. The cross section for indirect DM detection is enhanced by recombination, (QQ) + (QQ) → (QQ) + (QQ), and still compatible with bounds (section 5.2). At colliders quorns manifest as (quasi)stable charged tracks: LHC sets the bound M Q > ∼ 2 TeV.
In section 4 we studied the SIMP hybrid states, which have large cross sections of order 1/Λ 2 QCD and a relic abundance 3 or more orders of magnitude smaller than DM. In view of this, they seem still allowed by the experiments which excluded SIMP DM (ρ SIMP = ρ DM ), although a Manhattan-like project would be needed to predict their properties. Our model contains two kind of SIMPs: the isospin-singlet Qg with no interaction to pions; and the isospin triplet Qqq . Presumably the latter are heavier and decay. We do not know whether Qg can bind with (large enough?) nuclei, and how they would bind during Big Bang Nucleosynthesis, given that there is no first-principle understanding of nuclear potentials. The following statements are safe: our predicted SIMP abundance is so small that they negligibly affect ordinary BBN; SIMPs get stopped by the Earth atmosphere and are not visible in underground detectors; SIMP annihilations negligibly heat the Earth.
The interpretation of searches for rare hybrid heavy nuclei in samples of materials depends on the history of SIMPs and of the selected samples: from the Big Bang, to star burning, through Earth geology. The primordial abundance of SIMPs in the Earth and in stars sank down to their centres, undergoing QQ annihilations. Thereby, in order to set bounds, we consider the smaller secondary abundance of SIMPs. Presumably most primordial SIMPs still are in galactic clouds, and the Earth is big enough to capture all SIMPs encountered along its trajectory. The total energy stored in captured SIMPs likely exceeds the energy of the world fossile fuel reserve by 10 4 . What happens after capture is unclear. If SIMPs do not bind in nuclei, they sink in the Earth ocean and crust with drift velocity v ∼ 0.1 km/yr, such that their ground-level abundance is much below existing bounds. They can be searched for through dedicated enrichment processes and Rutherford backscattering experiments. If instead SIMPs bind within nuclei, electromagnetic interactions keep them in the crust since when the crust become geologically stable. Then, the local SIMP density can be comparable to present bounds, depending on the capture cross section by nuclei, which is highly uncertain. SIMP searches have been also performed in meteorites, where SIMPs cannot sink. Despite this, meteorites are made of heavy elements synthesised by stars: primordial SIMPs sank to the center of stars, and never come back. The secondary abundance of SIMPs in meteorites depends on the SIMP capture cross section by individual nuclei, which is highly uncertain and possibly vanishing. Present bounds are satisfied assuming a SIMP capture cross sections comparable to the one of neutrons with similar MeV energy, σ capture ∼ 0.01/Λ 2 QCD . In conclusion, colored DM seems still allowed, altough close to various bounds. Direct detection seems to provide the strongest and safest probe.
We discussed the apparently nicer model of colored DM: a neutral Dirac fermion Q in the adjoint representation of color. A scalar would give a similar phenomenology, and the DM abundance would be reproduced thermally for a similar M Q ∼ 12.5 TeV. A smaller mass would be obtained for quorns in the fundamental of SU(3) c , although the mass of the quorn-onlyum DM state QQQ would be M DM ≈ 3M Q . In models where Q has an asymmetry, the DM abundance can be obtained for lower M Q .
Finally, we notice that the fall of free Q down to deep multi-Q bound states occurs around the QCD phase transition out of thermal equilibrium. It could thereby contribute to baryogenesis, provided that violation of baryon number can be added at an acceptable model, possibly assuming that B is a gauge symmetry spontaneously broken giving rise to processes that violate ∆B = 1. where 2 F 1 is the Hypergeometric2F1 function. A similar formula can be obtained for R n , n, −1 by the interchange of the indices n and n . The total decay rate and energy loss rate from an initial state (n, ) is obtained by summing over all available lower-energy states with n < n.

B Toy redecoupling
We here show that the chromodark-synthesis mechanism is absolutely unavoidable by discussing a toy model that allows to analytically understand some of its features. We consider formation of one bound state B QQ containing two DM quarks Q from two bound states B Q containing one DM quark: 10 where X denotes any other SM particles, such as pions. We define δ ≡ 2M B Q − M B QQ . In the real situation described in section 3, many bound states with a semi-classical discretuum of binding factors δ can be produced. We simplify the problem by considering just one of them, with δ ∼ Λ QCD such that the QCD cross section for the above process is large, σ QQ ∼ 1/δ 2 . One then reaches thermal equilibrium n B QQ n 2 This means that the B Q dominantly form B QQ at the redecoupling temperature is an entropy factor that describes how much formation of B QQ gets delayed by having a plasma with much more particles X than can break it, than particles B Q that can form it. This is analogous to how e, p bind in hydrogen at T < ∼ δ/ ln(n γ /n p ), and to how p, n bind in deuterium at T < ∼ δ/ ln(n γ /n p ), where δ are the binding energies of hydrogen and deuterium respectively. 11 In the toy model, the residual density of B Q is estimated as its thermal equilibrium value at the redecoupling temperature where the interaction rate Γ QQ ∼ n B Q σ QQ v rel for the process of eq. (54) becomes smaller than the Hubble rate. Imposing Γ QQ ∼ H with H ∼ T 2 /M Pl , v rel ∼ T /M Q and This shows that re-annihilation is dominated by bound states with smaller δ ∼ Λ QCD , rather than by deep states. In the full computation many bound states contribute to the depletion of Y B Q , that gets about 2 orders of magnitude smaller than the toy-model estimate of eq. (57). In turn, the unavoidable toy-value is much smaller than what obtained by including only perturbative QCD annihilations at T ∼ T dec Λ QCD .
C Chromo-polarizability of QQ DM Eq. (40) provides the formula for the polarizability of a QCD bound state. We here evaluate it for our DM, the QQ singlet bound state |B = |1, s, α eff with energy E 10 = −α 2 eff M Q /4, where α eff = 3α 3 . By emitting a gluon it becomes a p-wave octet, with free Hamiltonian H 8 = p 2 /M Q − α 8 /r where α 8 = 3α 3 /2, whose eigenvalues are E 8n = −α 2 8 M Q /4n 2 for bound states and p 2 /M Q for positive energy states. To evaluate the matrix element in eq. (40) we insert the completeness relation for the octet eigenstates where the first (second) term is the contribution from bound (free) states. The factor 1/3 is introduced not to double count the angular momentum states. In coordinate space r|n, , m = R n (r)Y m (θ, φ) for bound states and r| p, , m = R p (r)Y m (θ, φ) for continuum positive energy states, where Y m (θ, φ) are spherical harmonics. The Coulombian wave-functions are where 1 F 1 is the Hypergeometric1F1 function; a i = 2/(α i M Q ) are the Bohr radii in each channel with effective coupling α i = {α eff , α 8 } and L 2 +1 n− −1 are Laguerre polynomials. 11 In the numerical computation such entropy factor was accounted in section 2.5 by imposing a small time allowed to radiate enough energy down to an unbreakable state. To keep the argument simple we here ignore the Boltzmann suppression in the π abundance at T < ∼ m π (in the full numerical computation this is taken into account and increases the σ fall computed in section 2.5, consequently suppressing the hybrid abundances). where the matrix element is | 1, s, α 1 | r|n, p, α 8 | 2 = | ∞ 0 dr r 3 R 10α eff (r)R n1α 8 (r)| 2 .