Gluequark dark matter

We introduce the gluequark Dark Matter candidate, an accidentally stable bound state made of adjoint fermions and gluons from a new confining gauge force. Such scenario displays an unusual cosmological history where perturbative freeze-out is followed by a non-perturbative re-annihilation period with possible entropy injection. When the gluequark has electroweak quantum numbers, the critical density is obtained for masses as large as PeV. Independently of its mass, the size of the gluequark is determined by the confinement scale of the theory, leading at low energies to annihilation rates and elastic cross sections which are large for particle physics standards and potentially observable in indirect detection experiments.

1 Introduction The striking success of the Standard Model (SM) in reproducing laboratory tests of fundamental interactions and its failure to explain some of the structural features of our Universe motivate the study of extensions based on its same principles of simplicity and elegance.From a modern point of view the SM is understood as an effective field theory with a very high ultraviolet cut-off, which appears renormalizable at energies currently probed in experiments.This feature notoriously gives rise to the SM hierarchy problem, but is also at the very origin of the attractive properties of the SM.In particular, global symmetries arise accidentally in the infrared and explain in the most economical way baryon and lepton number conservation, flavour and electroweak (EW) precision tests.
We consider these remarkable properties as paradigmatic, providing a compelling guidance to build possible extensions of the SM, even at the price of sacrificing the naturalness of the electroweak scale (as hinted anyway by experiments).In particular, the cosmological stability of Dark Matter (DM) can be elegantly explained in terms of accidental symmetries, in analogy with the stability of the proton following from baryon number conservation.This has to be contrasted with SM extensions where global symmetries are imposed ad hoc, like for example the case of R-parity in supersymmetry.A simple way to generate accidental symmetries is to extend the gauge theory structure of the SM by postulating a new confining dark color group.In this paper we will continue the exploration initiated in Refs.[1,2] of theories where the dark sector comprises new quarks transforming as real or vector-like representations under both the SM and dark gauge groups, and where the Higgs field is elementary.Such framework, also known as Vector-Like Confinement [3], provides a safe non-trivial extension of the SM, since it gives small and calculable deviations to EW observables that are in full agreement with current data and potentially observable at future colliders [4].
Previous studies of accidental DM focused on baryons or mesons of the dark dynamics as DM candidates, see [5] for a review.A systematic analysis of models with baryonic DM was performed in Refs.[1,2].In the regime where dark quark masses are below the dark color dynamical scale, M Q ă Λ DC , it was found that a correct relic abundance can be obtained for DM masses of order 100 TeV.Lighter DM masses, down to Op10 TeVq, are instead allowed if M Q ą Λ DC .Mesonic DM candidates can be even lighter and can arise for example as pseudo Nambu-Goldstone bosons (NGBs) of a spontaneously broken global symmetry of the dark dynamics.
In this work we explore scenarios with a new kind of accidental composite DM candidate, the gluequark, which has properties different from those of dark baryons and mesons in several respects.Gluequarks are bound states made of one dark quark and a cloud of dark gluons in theories where the new fermions transform in the adjoint representation of dark color.They are accidentally stable due to dark parity, an anomaly free subgroup of dark fermion number, which is exact at the level of the renormalizable Lagrangian.Depending on the SM quantum numbers of the new fermions, violation of dark parity can arise from UV-suppressed dimension-6 operators thus ensuring cosmologically stable gluequarks for sufficiently large cut-off scales.Contrary to baryons and mesons, the physical size of the gluequark is determined by the confinement scale independently of its mass.In the regime of heavy quark masses, M Q ą Λ DC , this implies a physical size larger than its Compton wavelength, see Fig. 1.The annihilation cross section for such a large and heavy bound state can be geometric, much larger than the perturbative unitarity bound of elementary particles.This in turn modifies the thermal relic abundance and can lead to significant effects in indirect detection experiments.Also, the resulting cosmological history is nonstandard and different from that of theories with baryon or meson DM candidates.
Bound states made of one dark fermion (the dark gluino) and dark gluons were considered in Refs.[6,7] in the context of supersymmetric gauge theories.There they arise as the partners of glueballs after confinement and were consequently called glueballinos.Ref. [6] showed that the observed DM abundance can be reproduced by a mixture of glueballs and glueballinos provided that the dark and SM sectors are decoupled very early on in their thermal history.In such scenario the two sectors interact only gravitationally, the dark gluino being neutral under the SM gauge group.Notice that the stability of glueballs in this case does not follow from an accidental symmetry but is a consequence of the feeble interaction between the SM and dark sectors.In this paper we will focus on the possibility that dark fermions are charged under the SM gauge group, so that the lightest states of the dark sector may be accessible through non-gravitational probes.In this case the dark and visible sectors stay in thermal equilibrium until relatively low temperatures, of the order of 1 GeV, and the thermal history of the Universe is rather different than that described in Refs.[6,7].In particular, we will argue that in our scenario dark glueballs cannot account for a sizeable fraction of DM because of BBN and CMB constraints.
Composite DM candidates from theories with adjoint fermions were also considered in the context of Technicolor models, see for example Refs.[8,9].Those constructions differ from ours in that technicolor quarks are assumed to transform as complex representations under the SM, but they can share common features with some of the models described in this paper 1 .
The paper is organized as follows.Section 2 provides a classification of models with adjoint fermions that can lead to a realistic DM candidate.We outline the cosmological history of the gluequark in section 3 and present our estimate for the thermal relic abundance in section 4. Section 5 discusses a variety of bounds stemming from cosmological and astrophysical data, DM searches at colliders, direct and indirect detection experiments.We summarize and give our outlook in section 6.A discussion of the relevant cross sections can be found in the Appendices. 1 Reference [9] for example considered gluequark DM in the context of the so-called Minimal Walking Technicolor model, but its estimate of the thermal relic abundance focuses on the perturbative freeze-out and does not include any of the non-perturbative effects described in this work.

The models
We consider the scenario in which the SM is extended by a new confining gauge group G DC (dark colour ), and by a multiplet of Weyl fermions Q (dark quarks) transforming in the adjoint representation of G DC and as a (possibly reducible) representation R under the SM group G SM : Q " padj, Rq . (2.1) In particular, we consider models where the dark quarks have quantum numbers under SUp2q L ˆUp1q Y but are singlets of SUp3q c color.New fermions colored under SUp3q c are an interesting possibility but are subject to very strong experimental constraints and their analysis deserves a separate study (see for example the recent discussion in Ref. [10]).
We assume R to be a real or vector-like representation, so that the cancellation of G SM anomalies is automatic and mass terms for the dark quarks are allowed.We performed a classification of the minimal models, i.e. those with the smallest representations and minimal amount of fields, which give a consistent theory of DM.We refer to them as 'minimal blocks'.Each block is characterised by two parameters: the dark quark mass M Q and the value of the dark gauge coupling g DC at M Q .A CP -violating θ term can also be included but does not play an important role in what follows.The renormalizable lagrangian thus reads It is possible to combine more than one minimal block; in this case the number of parameters increases: each module can have a different mass and, depending on the SM quantum numbers, Yukawa couplings with the Higgs boson may be allowed.
As long as all dark quarks are massive, the theory described by (2.2) confines in the infrared forming bound states.The symmetry Q Ñ ´Q, dark parity, is an accidental invariance of the renormalizable Lagrangian.It is what remains at the quantum level of the anomalous Up1q dark fermion number.The physical spectrum is characterized by states that are either even or odd under dark parity.The gluequark, denoted by χ in the following, is the lightest odd state and has the same SM quantum numbers of its constituent dark quark, thus transforming as an electroweak multiplet.Radiative corrections will induce a mass splitting among different components, with the lightest state being accidentally stable at the renormalizable level thanks to its odd dark parity.The mass difference computed in Ref. [11] shows that the lightest component is always the electromagnetically neutral one, which therefore can be a DM candidate provided it has the correct relic density.
We select models with a suitable gluequark DM candidate by requiring them to be free of Landau poles below 10 15 GeV.This is a minimal assumption considering that, as discussed below, astrophysical and cosmological bounds on the gluequark lifetime can be generically satisfied only for a sufficiently large cut-off scale.It is also compatible with Grand Unification of SM gauge forces.The ultraviolet behaviour of each model is dictated by the number of dark colors N DC and by the dimension of the SM representation R, i.e. by the number of Weyl flavors N f .Models with too large N f or N DC imply too low Landau Quantum numbers Table 1: Minimal building blocks for models of gluequark DM.We require that a multiplet contains an electromagnetic neutral component and that the gauge couplings do not have Landau poles below 10 15 GeV, assuming a representative mass of 100 TeV for the dark quarks.
poles for G SM , and are thus excluded from our analysis.The list of minimal blocks that satisfy our requirements is reported in Table 1 for SUpN DC q, SOpN DC q and SppN DC q dark color groups.Each block is characterized by its accidental symmetry (that can be larger than the dark parity) and by the dimensionality of the lowest-lying operator O dec which violates it.The latter has the form where O SM is a SM composite operator matching the SUp2q L ˆUp1q Y quantum numbers of the dark quark Q.The operator (2.3) can in general induce the mixing of the gluequark with SM leptons, providing an example of partial compositeness.As long as the theory is not in the vicinity of a strongly-coupled IR fixed point at energies E " M Q , Λ DC , the dimension of O dec is simply given by rO dec s " 7{2 `rO SM s, as reported in the sixth column of Table 1.Among the minimal blocks, the L L model has rO dec s " 5 classically.In this case the naive suppression of the gluequark decay rate is not enough to guarantee cosmological stability, although a stable DM candidate can still be obtained through additional dynamics, see the discussion in Appendix C. In the remaining minimal blocks the classical dimension of O dec is 6 and the gluequark can be sufficiently long lived.Indirect detection experiments and data from CMB and 21 cm line observations set important constraints on these models which will be discussed in section 5.The behaviour of the theory at energies above the confinement scale depends largely on the number of dark flavors N f and on the value of the dark coupling g DC at the scale M Q .One can identify two regimes.In the first, g DC pM Q q is perturbative and this necessarily implies a confinement scale smaller than the quark mass, Λ DC ă M Q ; we will call this the 'heavy quark' regime.In this case, depending on the value of N f , there are three possible behaviours.For N f ě N AF f " 11{2 the theory is not asymptotically free, hence starting from the UV the coupling gets weaker at lower scales until one reaches the quark mass threshold below which the dynamics becomes strong and confines. 2For N c f ď N f ă N AF f , where N c f is the non-perturbative edge of the conformal window, the theory flows towards an IR fixed point at low energies until the quark mass threshold is passed, below which one has confinement.Finally, if N f ă N c f the coupling grows strong quickly in the infrared and confinement is triggered without delay.Only for this latter range of values of N f the confinement scale can be larger than the quark mass, M Q ă Λ DC ; we will refer to this as the 'light quark' regime.The physical spectrum, the phenomenology and the thermal history are rather different in the two regimes.
The infrared behaviour of SUpN DC q gauge theories with fermions in the adjoint representation was extensively studied through lattice simulations, see for example [12][13][14][15][16][17][18][19][20][21] and references therein.There seems to be sufficient evidence for an infrared conformal phase of theories with N DC " 2 colors and N f " 4, 3 massless Weyl flavors, while results with N f " 2 are more uncertain though still compatible with a conformal regime.Theories with N f " 1 are supersymmetric and have been shown to be in the confining phase.The case with N DC " 3 colors is much less studied and no firm conclusion can be drawn on the conformal window.Notice that, independently of the number of colors, asymptotic freedom is lost for N f ě 6, while the existence of a weakly-coupled infrared fixed point can be established for N f " 5 by means of perturbation theory.Besides determining which phase the massless theory is in, simulations give also information on the spectrum of bound states.In particular, information on the gluequark mass in the limit of heavy quark masses (M Q " Λ DC ) can be obtained from lattice simulations with adjoint static sources, see for example Refs.[22,23].
Heavy quark regime: In the heavy quark regime, the lightest states in the spectrum are glueballs, while those made of quarks are parametrically heavier.The value of the glueball mass is close to the one of pure gauge theories.Lattice results for pure glue SUp3q theories show that the 0 ``state is the lightest with mass m 0 ``" 7Λ DC , see for example [24].Similar values are found for SUpN DC q with different number of colors.The gluequark is expected to be the lightest state made of quarks, with a mass M χ " M Q .Other states made of more dark quarks (collectively denoted as mesons) quickly decay to final states comprising glueballs and gluequarks, depending on their dark parity.
The gluequark lifetime can be accurately estimated by computing the decay of its constituent heavy quark, similarly to spectator calculations for heavy mesons in QCD.In the minimal blocks where the dark parity-violating operator has dimension 6 the main decay channel for the lightest gluequark χ 0 is χ 0 Ñ hν `nΦ (where Φ indicates a glueball and n ě 1).In the V model of Table 1 with three dark flavors transforming as an EW triplet, the dim-6 operator the decay of the gluequark with inverse lifetime Similar results apply for the N and T ' T minimal blocks.Glueballs can decay to SM particles through loops of dark quarks.In particular, since the latter are assumed to have electroweak charges, glueballs can always decay to photons through dimension-8 operators of the form G µν G µν W αβ W αβ generated at the scale M Q .For all the minimal models in Table 1 this is the lowest-dimensional operator which induces glueball decay.The partial width into photons is estimated to be [25,26] ΓpΦ Ñ γγq " 1 s ´1 ˆNDC 3 When phase space allows, the decay channels Zγ, W `W ´and ZZ open up producing one order-of-magnitude smaller lifetime.Relatively long-lived glueballs, as implied by the estimate (2.5), are subject to cosmological and astrophysical constraints as discussed in section 5. Models with Yukawa couplings to the Higgs doublet can be obtained combining minimal blocks.In this case 1-loop radiative effects at the scale M Q generate the dimension-6 effective operator |H| 2 G 2 µν , inducing a much shorter lifetime.If their mass is high enough, M Φ ą 2m h , glueballs predominantly decay to two Higgs bosons with a decay width ΓpΦ Ñ hhq " 10 10 s ´1 ˆNDC 3 where MQ " pM Light quark regime: If dark quarks are lighter than Λ DC , the physical spectrum is radically different and one expects spontaneous breaking of the global SU pN f q symmetry down to SOpN f q.The lightest states are thus the (pseudo) Nambu-Goldstone bosons ϕ, while the DM candidate is the gluequark, accidentally stable and with a mass of the order of the confinement scale Λ DC .As discussed in section 4, and similarly to the baryonic DM theories of Ref. [1], reproducing the correct DM relic density in this regime fixes Λ DC " 50 TeV.The NGBs with SM quantum numbers get a mass from 1-loop electroweak corrections, which is predicted to be Op10 TeVq for the value of Λ DC of interest.Besides such a radiative correction, the quark mass term breaks explicitly the SUpN f q global symmetry and gives an additional contribution.Including both effects, the NGB mass squared is given by where I is the weak isospin of the NGB and c 0,1 are Op1q coefficients.
For fermions in the adjoint representations, only models with N f ă 5 light quarks can be in the regime M Q ă Λ DC , since those with more fermions are either IR conformal or IR free.Therefore, among the minimal blocks of Table 1 only two are compatible with the light quark regime, i.e. the V model and the L ' L model.The V model has a global symmetry breaking SUp3q Ñ SOp3q which leads to five NGBs transforming as an electroweak quintuplet.In the L ' L model one has SUp4q Ñ SOp4q and nine NGBs transforming as 3 ˘, 3 0 of SUp2q EW ˆUp1q Y .
The limit of very small quark masses, M Q Λ DC !m 2 ϕ , is experimentally interesting, since NGBs have predictable masses.In general, the lightest NGBs decay to SM final states through anomalous or Yukawa couplings, as in the case of the V model.In some cases, however, some of the NGBs are accidentally stable due to unbroken symmetries of the renormalizable Lagrangian.An explicit violation of such accidental symmetries is expected to arise from higher-dimensional operators, possibly resulting into long-lived particles.An example of this kind is given by the L ' L model, where NGBs made of LL or L L constituents have U p1q number ˘2 and are stable at the renormalizable level, see Appendix C.
Since we assumed the dark quarks to transform as real or vectorlike representations under the SM gauge group, the fermion condensate responsible for the global symmetry breaking in the dark sector can be aligned along an (SU p2q L ˆU p1q Y )-preserving direction.In practice, we assume that no vacuum misalignment arises from perturbative interactions such as Yukawa couplings.As the strong dynamics preserves the EW gauge symmetry of the SM, it also affects electroweak precision observables through suppressed corrections which are easily compatible with current constraints for sufficiently high values of Λ DC , as required to reproduce the DM relic abundance [2][3][4].
Besides the NGBs, the physical spectrum comprises additional bound states with mass of order Λ DC .These include the gluequarks, which are expected to be the lightest states with odd dark parity, and mesons (i.e.bound states made of more than one dark quark) 4 .Except for the lightest gluequark, which is cosmologically stable, all the other states promptly decay to final states comprising NGBs and gluequarks, depending on their dark parity.In the minimal blocks where dark parity is broken by the dimension-6 operator HG µν σ µν Q, the most important decay channels of the gluequark are χ 0 Ñ hν and χ 0 Ñ hν `ϕ.The two-body decay dominates at large-N DC and gives a lifetime of order where f χ is the decay constant of the gluequark5 .
To summarize our discussion on models, Table 1 reports the minimal blocks which have a potentially viable DM candidate and a sufficiently high cut-off, above 10 15 GeV, as re-quired for SM Grand Unification and to suppress the DM decay rate.In particular, the requirement on the absence of Landau poles restricts the list of possible models to a few candidates.As mentioned before, the case of the singlet was studied already in the literature [6,7], and it will not be considered further in this work.We find that in all the other minimal blocks of Table 1 the SU p3q c ˆSU p2q L ˆU p1q Y gauge couplings unify with much lower precision than in the SM.Making the dark sector quantitatively compatible with SM Grand Unification thus requires extending these minimal blocks by including additional matter fields.Also, it would be interesting to explore the possibility of unifying both the visible and dark gauge couplings.We leave this study to a future work.
In the next sections we will discuss the thermal history of the Universe and try to estimate the DM relic density: section 3 explains the general mechanisms at work and is largely independent of the details of the models; section 4 gives a concrete example, adopting as a benchmark the V model of Table 1, i.e. the minimal block with an SU p2q L triplet.For a discussion of the L ' L model see Appendix C.

Cosmological History
The Universe undergoes different thermal histories in the light and heavy quark regimes.We first give a brief overview of such evolution, followed by a more detailed discussion with quantitative estimates.
In the light quark regime the thermal history is relatively simple and similar to that described for baryonic DM in Ref. [1].Dark color confines when dark quarks are relativistic and in thermal equilibrium.After confinement the gluequarks annihilate into NGBs with a non-perturbative cross section σv rel " π{Λ 2 DC , while glueballs are heavy and unstable.At temperatures T " M χ {25 the annihilation processes freeze out and the gluequark start behaving as ordinary WIMPs.
In the heavy quark regime the thermal history is more complex and characterized by three different stages.Before confinement (T Á Λ DC ), free dark quarks annihilate into dark gluons and undergo perturbative freeze-out at T " M Q {25 (see section 3.1).At confinement (T " Λ DC ), the vast majority of the remaining dark quarks hadronizes into gluequarks, while the plasma of dark gluons is converted into a thermal bath of non-relativistic glueballs.The formation of mesons is suppressed by the low density of dark quarks compared to the ambient dark gluons.Glueballs overclose the Universe if they are cosmologically stable, therefore we consider the region of the parameter space where their lifetime is sufficiently short.As first pointed out in [29][30][31], and recently reconsidered in [32,33], decays of non-relativistic particles with a large and non-thermal energy density -like the glueballs -can modify the standard relation between the scale factor and the temperature during the cosmological evolution.If the glueballs are sufficiently long lived and dominate the energy density of the Universe at some stage of the cosmological evolution, the standard scaling a ∝ T ´1 is modified into a ∝ T ´8{3 .During this early epoch of matter domination, the Universe expands faster than in the radiation-dominated era, leading to an enhanced dilution of the DM relic density (see section 3.2).Finally, interactions among gluequarks can lead to a second stage of DM annihilation through the process6 where QQ ˚is an excited bound state of dark quarks and V stands for a SM vector boson or possibly a Higgs boson in models with Yukawa interactions.The process (3.1) proceeds in two steps.Initially, an excited bound state QQ ˚with size Op1{Λ DC q is formed by a collision of two χ's through a recombination of the constituent heavy quarks.This is similar to what happens for example in hydrogen anti-hydrogen scattering [37].As a consequence of the large size of the gluequark (see the discussion in section 3.3), the corresponding recombination cross section is expected to be large σ rec « π{Λ 2 DC .Once formed, the QQ can either decay (QQ ˚Ñ QQ `V Ñ SM) or be dissociated back into two gluequarks by interactions with the SM and glueball baths (Φ{ V `QQ Ñ χ `χ).A naive estimate shows that the latter process typically dominates.This is because the largest contribution to the total cross section comes from scatterings with large impact parameters, b " 1{Λ DC , in which the QQ ˚is produced with a large angular momentum, " M Q vb.Bound states with " 1 take more time to de-excite to lower states, and dissociation can happen before they reach the ground state.The annihilation of gluequarks through recombination is therefore inefficient as long as the glueball bath is present.Only when the glueballs decay away, a second stage of DM annihilation can take place through the process (3.1).

Thermal freeze-out
Thermal freeze-out is the first (only) phase of the cosmological evolution in the regime with heavy (light) quarks.In this stage the number density of free dark quarks (for M Q ą Λ DC ) or of gluequarks (for M Q ă Λ DC ) is reduced until it becomes so low that chemical equilibrium is no longer attained and freeze-out takes place.The number density at freezeout is approximately given by npT f.o.q » HpT f.o.q xσ ann v rel y , where H is the Hubble parameter, and afterwards it is diluted by the Universe expansion.
In the heavy quark regime, free dark quarks annihilate with a perturbative cross section into dark gluons and into pairs of SM particles (vector bosons, Higgs bosons and fermions).The freeze-out temperature is of order T f.o.« M Q {25.A general expression for the annihilation cross section is reported in Appendix A, see eq. (A.2).For the V model with N DC " 3 analysed in the next section, the annihilation cross section into dark gluons and SM fields is neglecting the mass of final states.The term from annihilation into SM particles separately shows the contribution of vectors and fermions/longitudinal gauge bosons.Terms in the second parenthesis encode the Sommerfeld enhancement from dark gluon exchange: S 3 , S 3{2 , S ´1 refer respectively to the 1, 8 and 27 color channels and are given by [38,39] S n " In the light quark regime gluequarks annihilate into NGBs with a cross section that is expected to scale naively as in analogy with nucleon-nucleon scattering in QCD [40].Nambu-Goldstone bosons are unstable and later decay into SM particles.

Dilution
As well known, the number density of DM particles today is related to the number density at freeze-out by This relation is usually rewritten in terms of temperatures assuming that between freezeout and today the standard scaling a ∝ T ´1 holds.However, the validity of the standard scaling relies upon the assumption that entropy is conserved in the SM sector, i.e. that no energy is injected into the SM plasma.In presence of large entropy injection one can have an epoch during which a grows faster than a ∝ T ´1.In this case the relation between n DM pT 0 q and n DM pT f.o q is given by: ˙3 ˆa pT i q a pT f q where T i and T f defines the temperature interval during which the non-standard scaling holds (see Fig. 2).The last term in parenthesis accounts for the suppression with respect to the naive relict density which would be obtained using the standard scaling.In the following we will show that late-time decays of dark glueballs can give rise to a non-standard scaling of the form a ∝ T ´α with α ą 1.The corresponding suppression factor thus reads: After dark color confinement, the energy density of the Universe can be divided into a relativistic component, ρ R , containing all the SM relativistic particles, and a non-relativistic one, ρ M , containing all the dark-sector long-lived degrees of freedom (i.e.dark glueballs and gluequarks).In particular, the energy density of glueballs at confinement is much larger than the corresponding thermal energy density for a non-relativistic species, and this can lead to an early epoch of matter domination.The evolution of ρ M,R is governed by  3.8) imply a scale factor at freeze-out, a f.o., smaller than the one obtained from a standard cosmology, ãf.o. .This in turn leads to a suppression of the relic density by a factor pa f.o.{ã f.o.q 3 " pT f {T i q 3α´3 .Right: Scaling obtained by solving numerically (3.9) and (3.10) for M Q " 100 TeV and Λ DC " 10 GeV.There exists an epoch during which the apT q scaling is very well approximated by a power law a ∝ T ´α with α " 8{3.
where Γ Φ is the glueball decay rate and the Hubble parameter H is given by the Friedmann equation: Since in the relevant region of the parameter space the dark and SM sectors are in thermal equilibrium at dark confinement, the initial conditions at T " T dc « Λ DC are given by ρ M pT dc q " ξ ρ R pT dc q with ξ " g ˚pT dc q g D pT dc q , (3.11) where g D pT q and g ˚pT q count the number of relativistic degrees of freedom in the dark and SM sector respectively.Furthermore, assuming that the decay products thermalize fast enough, the temperature of the Universe below T dc is related to the relativistic energy density by: The evolution during the early matter-dominated epoch, if the latter exists, can be described by solving analytically eq.(3.9) at leading order in ρ R {ρ M for cosmic times t !1{Γ Φ [33]: Here ρM,R and ā denote the initial conditions at some time t much after the beginning of the matter-dominated epoch.The relativistic energy density is given by the sum of ρR (first term in eq.(3.13b)), diluted as a ´4, and the energy injected by glueball decays (second term in eq.(3.13b)), diluted as " a ´3{2 .Initially the first term dominates and the standard scaling a ∝ T ´1 is obtained; as long as the glueball lifetime is long enough, the second term will start to dominate at some temperature T i , implying a non-standard scaling a ∝ T ´8{3 (see Fig. 2).The value of T i can be found by equating the first and second terms of eq.(3.13b) and by using eqs.(3.10), (3.11): The non-standard scaling ends when almost all the glueballs are decayed, i.e. around pt ´tdc q " Γ ´1 Φ , where t dc is the time at dark confinement.Using eqs.(3.10) and (3.11), one can translate this condition in terms of a temperature finding: From eq.(3.8) it follows that late-time decays of glueballs dilute the naive relic density by a factor where Op1q numerical factors omitted in eq.(3.14) and (3.15) have been included.When the glueballs are sufficiently long lived to give a sizeable dilution, the second term in the numerator inside the parenthesis of eq.(3.16) can be neglected and F is very well approximated by: While the analytic formulas (3.13b)-(3.17)turn out to be quite accurate, in our estimate of the relic density performed in section 4 we will solve eq.(3.9) numerically without making any approximation.

Reannihilation
At T " T dc " Λ DC the theory confines and the dark degrees of freedom reorganize into singlets of dark color.In the heavy quark regime, the number density of gluons is much larger than the one of fermions and the vast majority of free quarks Q hadronize into gluequarks.These can then collide and recombine in excited QQ ˚states by emitting an electroweak gauge boson (or a Higgs boson in theories with Yukawa couplings) or a glueball when kinematically allowed, see eq. (3.1).The process goes through a recombination of the constituent heavy quarks, while the direct annihilation of these latter has a small and perturbative rate.Given that gluequarks have a size of order 1{Λ DC , one expects naively a recombination cross section of order σ rec " 1{Λ 2 DC .This value can in fact be reduced by kinematic constraints and the actual total cross section depends ultimately on the temperature at which the process takes place.A detailed discussion and estimates for the recombination cross section are given in Appendix B.
r n l v z 7 i 6 q 9 e u i m R I 7 Z i f s j H n s k t X Z L W u w J h N s z J 7 Y M 3 t x H p 1 X 5 8 1 5 / x l d c Y q d I / Y H z s c 3 j m G S p A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " c K s h m S 5 + r o 9 B r g D r n l v z 7 i 6 q 9 e u i m R I 7 Z i f s j H n s k t X Z L W u w J h N s z J 7 Y M 3 t x H p 1 X 5 8 1 5 / x l d c Y q d I / Y H z s c 3 j m G S p A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " c K s h m S 5 + r o 9 B r g D      Once formed, QQ ˚states with mass M pQQ ˚q ą 2M Q will promptly decay back to two gluequarks.Lighter states, on the other hand, can either de-excite and thus decay into SM particles through the emission of a SM vector boson or a glueball (QQ ˚Ñ QQ `V {Φ Ñ SM), or be dissociated by interactions with the glueball and SM plasmas (Φ{ V `QQ ˚Ñ χ `χ), see Fig. 3.
If de-excitation occurs faster than dissociation, a second era of efficient DM annihilation can take place, reducing the gluequark number density.While re-annihilation processes can be active over a long cosmological time interval, it is the last stage during which the re-annihilation cross section gets its largest value σ rea that is most important to determine the final gluequark density.This last stage happens relatively quickly and can be characterized by a re-annihilation temperature T R .The exact value of T R depends on the rate of dissociation and is difficult to estimate.The largest uncertainties arise from the calculation of the de-excitation rate, which can vary over several orders of magnitude.We performed a thorough analysis taking into account the many dynamical ingredients which play a role in determining both the re-annihilation cross section and temperature.A detailed account is reported in Appendix B. We find that, under the most reasonable assumptions, dissociation of the most excited QQ ˚states occurs faster than de-excitation, as long as the glueball bath originating from dark gluons confinement is present; therefore, the re-annihilation temperature is approximately equal to the one at which glueballs decay (T R « T D ).Besides this most probable scenario, in the following we will also consider the other extreme possibility where re-annihilation occurs right after confinement (T R " Λ DC ).The comparison between these two opposite scenarios will account for the theoretical uncertainties intrinsic to the determination of the non-perturbative dynamics characterizing our DM candidate.
In both benchmark scenarios considered above the last stage of the re-annihilation epoch occurs while entropy is conserved in the Universe and can thus be described by a set of standard Boltzmann equations given in eq.(B.1).They reduce to a single equation for sufficiently large de-excitation or glueball decay rates.This reads where z " M Q {T , Y χ " n χ {s and s is the entropy density of the Universe.The equilibrium term can be neglected since T R ď Λ DC !M Q .Assuming a re-annihilation cross section which is constant and velocity independent7 , eq. (3.18) can be easily integrated analytically; one obtains (for T ă T R ) Late-time annihilation significantly affects the gluequark relic density when the second term in the above equation dominates, i.e. roughly when in agreement with a naive expectation.When condition (3.20) is met, any dependence from the previous stages of cosmological evolution, encoded in Y χ pT R q, is washed out and the asymptotic value of the relic density is set only by re-annihilation.For temperatures T sufficiently smaller than T R (but higher than a possible subsequent period of dilution, in the case T R " Λ DC ), eq.(3.19) can be recast in terms of the gluequark relic density as follows:

Estimate of the Relic density
The cosmological evolution of gluequarks is determined by the interplay of the mechanisms described in the previous section and depends on the two fundamental parameters M Q and Λ DC .For each point in the plane pΛ DC , M Q q one can thus in principle reconstruct the thermal history of the Universe and compute the DM abundance Ω χ .In this section we will sketch the different possible thermal histories and give an estimate for Ω χ .As a reference model we consider the minimal module with a triplet of SUp2q L (see Table 1).We will assume the theory to be outside its conformal window, so that the regime of light dark quarks is well defined.We will discuss at the end how the picture changes for different SM representations and when the theory is in the conformal window or is not asymptotically free.
We will try to quantify the large uncertainties that arise in the determination of the cosmological evolution and of the relic density as a consequence of the non-perturbative nature of the processes involved.As anticipated in section 3.3, one of the largest uncertainties comes from the identification of the re-annihilation temperature T R .We will consider the On the black solid line Ω χ " Ω DM .We also report (dashed line) the isocurve Ω χ " Ω DM for the case where re-annihilation is not considered.The numbers indicate the three thermal histories described in the text.In the yellow region the glueballs are either stable or have a lifetime bigger than 1 s.In the first case they will over-close the Universe while in the latter they will spoil BBN, both cases are therefore forbidden.The blue region is ruled out by indirect searches, namely modifications of the CMB power spectrum, 21-cm line observables and indirect detection (see section 5.3).
two previously discussed benchmarks: T R " T D , the most plausible according to our estimates, and T R " Λ DC .We reconstruct for each of them the different possible cosmological evolutions obtained by varying M Q and Λ DC .Our estimate of the DM abundance for both benchmarks is reported in Fig. 4, where we show the isocurve Ω χ h 2 " 0.119 reproducing the experimentally observed density.
Let us consider first the case T R " T D .There are three possible thermal histories that can be realized (they are correspondingly indicated in the left plot of Fig. 4): 1.For very large M Q {Λ DC the Universe undergoes a first perturbative freeze-out at T f.o." M Q {25, then dark confinement occurs at T " Λ DC followed by an epoch of dilution between T i and T f " T D " T R8 .Glueballs decay at T À T D , and the number density of gluequarks is too small, as a consequence of the dilution, to ignite a phase of non-perturbative re-annihilation.The DM density is therefore given by where the number density at freeze-out is estimated by solving the Boltzmann equations numerically and is approximately given by eq.(3.2).By using the dilution factor reported in eq.(3.17) and setting T f.o." M Q {25 one obtains which describes well the slope of the upper part of the relic density isocurve in the left panel of Fig. 4. Because of the extreme dilution happening during the early epoch of matter domination, the experimental DM abundance is reproduced in this case for very large DM masses, of order of hundreds of TeV or more, much above the naive unitarity bound.
2. For smaller values of M Q {Λ DC (but still with M Q {Λ DC Á 25), the dilution between T f.o. and T D is not enough to prevent re-annihilation (i.e.condition (3.20) is met).The latter thus occurs at T » T D , washing out any dependence of Ω χ from the previous stages of cosmological evolution; the corresponding DM relic density is 3) The first factor corresponds to the gluequark energy density at the end of the reannihilation (given by eq.(3.21)), and the second one encodes the standard dilution due to the Universe expansion.We evaluate the re-annihilation cross section by using the semiclassical model described in Appendix B.1; this gives The parameters ε Φ and ε V are smaller than 1 and encode the suppression from energy and angular momentum conservation respectively for the recombination processes χχ Ñ QQ ˚`Φ and χχ Ñ QQ ˚`V .While eq.(4.4) is the result of a rather sophisticated analysis of the re-annihilation dynamics and represents our best estimate for σ rea , it is subject to large theoretical uncertainties, as discussed in Appendix B.1.We thus also consider the extreme situation where the re-annihilation cross section is always large and saturated by its geometric value Varying σ rea between the values in eqs.(4.4) and (4.5) will quantify the uncertainty on n χ p T q.By using eq.(3.21) and setting T R " T D " ?M Pl Γ Φ , eq. ( 4.3) takes the form: This formula describes the intermediate part of the isocurve in the left plot of Fig. 4. Initially (i.e. for 150 GeV À Λ DC À 800 GeV) the re-annihilation is dominated by the process χχ Ñ QQ ˚`Φ and ε Φ » 1; in this case the last factor in (4.6) can be well approximated with 1 (the electroweak contribution to σ model rea is small) and the estimated uncertainty on the gluequark relic density is negligible.For larger Λ DC re-annihilation into QQ ˚plus a glueball becomes kinematically forbidden in our semiclassical model, and ε Φ quickly drops to zero (see Appendix B.1).In this region ε V » 1{10 and varying σ rea between σ model rea and σ geo rea spans the gray region.The extension of the latter quantifies the uncertainty of our estimate of the relic density.
3. When M Q {Λ DC ď 25, the perturbative freeze-out does not take place.If M Q is bigger than Λ DC , then the Universe undergoes a first epoch of annihilation of dark quarks for T Á M Q , followed after confinement by the annihilation of gluequarks, until thermal freeze-out of these latter occurs at T » M χ {25.If M Q ă Λ DC , on the other hand, the theory is in its light quark regime and the only epoch of annihilation is that of gluequarks after dark confinement, again ending with a freeze-out at T » M χ {25.
Afterwards n χ is diluted by the Universe expansion without any enhancement from the decay of glueballs (these are too short lived to give an early stage of matter domination).The expression for the DM relic density is formally the same as in eq.(4.1) with For 1 À M Q {Λ DC À 25 the non-perturbative annihilation of gluequarks proceeds through the same recombination processes of eq.(3.1).According to the model of Appendix B.1, only the final state with a vector boson is kinematically allowed, and ε V » α 2 {10.This implies σ ann » pα 2 {10q 4π{Λ 2 DC , so that the DM relic density turns out to be independent of M Q .If instead the re-annihilation cross section is estimated by eq.(4.5), then by continuity with the previous cosmological evolution one must take σ ann » 4π{Λ 2 DC , which also corresponds to a relic density independent of M Q .Varying σ ann between these two values gives the largest vertical portion of the gray region in the left plot of Fig. 4.
As soon as one enters the light quark regime, M Q ă Λ DC , the annihilation of gluequarks proceeds through the direct annihilation of their constituents (the theory at M Q is non-perturbative) with a cross section σ ann " 4πc{Λ 2 DC , where c is an order 1 coefficient.We vary 1{5 ă c ă 1 to quantify the uncertainty in this last nonperturbative process.We thus obtain the narrower vertical portion of the gray region in the left plot of Fig. 4, which extends down to arbitrarily small M Q .The observed relic density in this regime is reproduced for Λ DC » 50 TeV, similarly to the light quark regime in baryonic DM models [1].
Let us turn to the case T R " Λ DC .As for T R " T D one can identify three possible thermal histories (correspondingly indicated in the right panel of Fig. 4): 1.For M Q {Λ DC " 25 the Universe goes first through a perturbative freeze-out of dark quarks at T f.o.» M Q {25, then re-annihilation occurs right after confinement for T » Λ DC .Finally, dilution takes place between T i and the temperature of the glueball decay T D .The DM relic density is given by the expression in eq.( 4.3) times the dilution factor F. Numerically one has In this case, our semiclassical model estimates ε Φ » 1{100 throughout the parameter space of interest.By varying σ rea between σ model rea and σ geo rea we thus obtain the upper portion of the gray region in the right plot of Fig. 4.

2.
For smaller M Q {Λ DC (but still with M Q {Λ DC ą 25), the glueballs are too short lived to ignite the dilution, and the DM relic density is given by eq.(4.3).Setting T R " Λ DC one obtains 3. When M Q {Λ DC ă 25 the cosmological evolution of the Universe is the same as thermal history 3 in the case T R " T D .The DM relic density is given by eq.(4.7), corresponding to the vertical gray regions of the right plot of Fig. 4.
The previous qualitative picture is largely independent of the details of the specific model.However, the quantitative results can change significantly in models with Yukawa couplings, where the glueballs lifetime is much shorter.Models that, in the limit of zero quark masses, are infrared free or in the conformal window are constrained to be in the regime M Q ą Λ DC .

Phenomenology and Experimental Constraints
In this section we outline the main phenomenological signatures for collider physics and cosmology of the models with gluequark DM.In general, the phenomenology has analogies to the one of baryonic DM studied in Refs.[1,2].Given the large gluequark masses needed to reproduce the DM relic density both in the light and heavy quark regimes, searches at collider are not promising, whereas cosmological observations provide interesting bounds.

Collider searches
The dark sector has a rich spectrum of states which, in principle, one would like to study at colliders.
The lightest states in the spectrum, with mass given by eq.(2.8), are the NGBs from the SUpN F q Ñ SOpN F q global symmetry breaking in the light quark regime.In the case of the V model, the five NGBs form a multiplet with weak isospin 2, and one expects m ϕ Á Λ DC {5.The phenomenology of a quintuplet of NGBs was studied recently in Ref. [4].These states are pair produced at hadron collider in Drell-Yan processes through their electroweak interactions, and decay to pairs of electroweak gauge bosons through the anomalous coupling 2N DC ? 3 A promising discovery channel studied by Ref. [4] is pp Ñ ϕ 0 ϕ ˘Ñ 3γW ˘; the doubly charged states decay into same-sign W pairs and are somewhat more challenging to see experimentally.The LHC has an exclusion reach up to TeV masses, while a 100 TeV collider would test the light quark scenario approximately up to 5 TeV.In this regime colliders could start probing the thermal region.The lightest states in the heavy quark regime are the glueballs.They couple to the SM only through higher-dimensional operators, and are rather elusive at colliders.In models without Yukawa couplings, where interactions with the SM occur through dimension-8 operators, the production cross section via vector boson fusion (VBF) or in association with a SM vector boson is too small to observe a signal in current or future colliders (for example, the VBF cross section at a 100 TeV collider is of order σppp Ñ Φ `jjq À 10 ´9 fb for M Q {Λ DC " 10 and M Φ ą 500 GeV).In models with Yukawa couplings the glueballs mix with the Higgs boson and production via gluon-fusion becomes also possible.While this leads to larger cross sections, the corresponding rate is too small to see a signal at the LHC and even high-intensity experiments like SHiP can only probe light glueballs in a region of parameter space that is already excluded by EW precision tests [2].
Mesons can give interesting signatures in both light and heavy quark regimes.Bound states made of a pair of dark quarks, QQ, can be singly produced through their EW interactions.While the production of spin-0 mesons is suppressed since they couple to pairs of EW gauge bosons, spin-1 resonances mix with the SM gauge bosons of equal quantum numbers and can be produced via Drell-Yan processes.In the narrow width approximation the cross section for resonant production can be conveniently written in terms of the decay partial widths as where D QQ is the dimension of the representation, J QQ the spin, P the parton producing the resonance and C PP are the dimension-less parton luminosities.
In the heavy quark regime the QQ bound state is perturbative and its decay width can be computed by modelling its potential with a Coulomb plus a linear term.For α DC M Q ą Λ DC the decay width of the lowest-lying s-wave bound states scales as where α 3 DC originates from the non-relativistic Coulombian wave-function.When α DC M Q ă Λ DC , the effect of the linear term in the potential becomes important and eq.( 5.3) gets modified; since confinement enhances the value of the wave function at the origin, the width becomes larger in this regime.Using the Coulombian approximation thus provides conservative bounds.Explicit formulas for the rates are found in [2].For example, in the V model the decay width of the s-wave spin-1 QQ resonance (isospin 1 in light of the Majorana nature of V ) into a left-handed fermion doublet is where n refers to the radial quantum number.The tiny energy splitting between levels is irrelevant at colliders and the total rate is dominated by the lowest-lying Coulombian ones.
The branching ratio into pairs of leptons is about 7% and the strongest bounds currently arise from searches of spin-1 resonances at the LHC decaying into electrons and muons.We show the limits in the left plot of Fig. 5 and find that the LHC excludes masses up to 2 ´3.5 TeV depending on the ratio M Q {Λ DC (or equivalently on the value of α DC pM Q q). Figure 5: Left: ATLAS bounds on the cross section for the direct production of a spin-1 QQ resonance decaying into muons and electrons [41].Right: Estimated reach on gluequark pair production obtained by recasting the limits of Ref. [42] from disappearing tracks searches at the HL-LHC (red), the HE-LHC (green) and a 100 TeV collider (blue).The solid (dashed) lines assume a 20% (500%) uncertainty on the background estimate. .
In the light quark regime the lightest spin-1 state is the ρ meson with mass M ρ " Λ DC .The widths scale as [43] Γpρ Ñ ϕϕq " where g ˚characterizes the interaction strength among bound states.For moderately large g ˚, as suggested by large-N DC counting g ˚" 4π{ ?N DC , the decay into light NGBs dominates while final states with leptons are suppressed.It thus follows a weaker bound than in the heavy quark regime, as illustrated in Fig. 5 for g ˚" 4.
Gluequarks can also be pair produced at colliders through their EW interactions.In the heavy quark regime the energy threshold is much higher than the confinement scale and quarks are produced in free pairs.Because dark quarks are in the adjoint representation of dark color, when they get separated by a distance of Op1{Λ DC q they hadronize producing color singlets that fly through the detector.On the contrary, dark quarks in the fundamental representation would not be able to escape, leading to quirks/hidden valley phenomenology [2,44,45].The phenomenology of the open production is then identical to the one of an elementary electroweak multiplet except that the cross-section is enhanced by the multiplicity of the dark color adjoint representation, i.e.N 2 DC ´1 for SUpN DC q.Such enhancement factor is not present for gluequark pair production near threshold in the light quark regime.In general, an electroweak triplet can be searched for in monojet and monophoton signals or disappearing tracks, the latter being more constraining.We derived the reach of the high-luminosity LHC (HL-LHC), the high-energy LHC (HE-LHC) and the proposed 100 TeV collider by recasting the results of Ref. [42] for the V model in the heavy quark regime, see the right plot of Fig. 5.We find that the HL-LHC could discover gluequark triplets up to " 600 GeV while a 100 TeV collider could reach " 7 TeV.Such bounds are typically weaker than the ones from the production of QQ spin-1 resonances decaying to leptons.

DM Direct Detection
From the point of view of DM direct detection experiments, where the momentum exchanged is less than 100 KeV, the gluequark behaves as an elementary particle with the same electroweak quantum numbers as the constituent quark.The main difference from elementary candidates with same quantum numbers is that the relic abundance is not controlled by the electroweak interaction, leading to a different thermal region.
For a triplet of SUp2q the spin-independent cross-section is σ SI " 1.0 ˆ10 ´45 cm 2 , which is below the neutrino floor for masses M χ ą 15 TeV.For an SUp2q doublet tree-level Zmediated interactions induce a spin-independent cross section on nucleons σ SI « 10 ´40 cm 2 , which is ruled out for M χ À 10 8 GeV [46].Dark quark masses large enough to make the doublet model viable can be obtained only in the scenario where T R " Λ DC , while the scenario with T R " T D is ruled out (see Appendix C for more details).

DM Annihilation and Decay
After freeze-out, decays or residual annihilations of gluequarks can give rise to indirect detection signals.
In the region of parameter space relevant for our purposes, annihilation processes set constraints on theories in the heavy quark regime, and we thus focus on that case to analyse them.As discussed in section 3.3 (and more extensively in Appendix B), the annihilation can be either direct (χχ Ñ nΦ{SM) or mediated by the formation of a QQ bound state that subsequently decays (χχ Ñ QQ ˚Ñ nΦ{SM).In the former case the annihilation cross section is perturbative (see eq.(A.2)) and, given the relatively high mass of the gluequark, it does not lead to any interesting indirect detection signatures.The latter case, instead, because of the enhanced annihilation cross section, could lead to an interesting phenomenology and it is worth further study.
As discussed in section 3.3 and in Appendix B, for angular momenta " 1 the recombination cross section is of order σ rea » ε 4π{Λ 2 DC .However, given the small velocities relevant for indirect searches (v rel " 10 ´6a TeV{M Q at the CMB epoch and v rel " 10 ´3 in our galaxy), the angular momentum of the colliding particles " M Q v rel {Λ DC is of order unity in a large region of the parameter space.In this case only s-wave processes take place and σ rea v rel rather than σ rea is constant.In this regime the value of the cross section is very uncertain, and we chose to estimate it in terms of two benchmark scenarios (see Appendix B.1): Figure 6: Exclusion limits from experimental searches sensitive to the DM annihilation (red regions) and decay (blue regions).The limits from DM annihilation respectively in the left and right panels have been obtained by adopting the two benchmarks of eq.(B.2), while limits from DM decay were derived by setting Λ U {g U V " MPl " 2.4ˆ10 18 GeV and f χ " 3Λ DC {p4πq in eqs.(2.4),(2.9).
Once formed, the QQ ˚bound states de-excite and in general decay into dark gluons, SM gauge bosons or SM fermions.The branching ratios can be derived in terms of the perturbative annihilation cross section of dark quarks into the corresponding final states, see eq.(A.2).In the region of interest α DC ą α 2 and one finds where G denotes a dark gluon.For the specific case of the V model, the tree-level decay into SM fermions and ZG is forbidden (the χ 0 has vanishing coupling to the Z) and use of eq.(A.2) thus gives where the last factor is from the branching ratio of QQ into W W . Similarly to residual annihilations, decays of the gluequark could give rise to indirect signals.The χ 0 decays mostly to hν plus glueballs in the heavy quark regime, and to hν or hν `ϕ in the light quark regime (see eqs.(2.4),(2.9)).Both glueballs and NGBs in turn decay into SM particles and ultimately the gluequark decay leads to the production of light SM species which can be observed experimentally.Bounds can be avoided, on the other hand, if some mechanism is at work that makes the χ 0 absolutely stable or give it a much longer lifetime than the one estimated in eq.(2.4) and (2.9).
Figure 6 summarizes the constraints in the plane pΛ DC , M Q q that arise from experiments probing DM decay and annihilation.The red exclusion regions from DM annihilation have been derived for the two benchmark values of xσ ann v rel y in eq.(B.2), while the blue ones from DM decay were obtained by setting Λ U V {g U V " MPl " 2.4 ˆ10 18 GeV and f χ " 3Λ DC {p4πq when evaluating τ DM from eqs. (2.4),(2.9).Experimental bounds are given in terms of the DM mass M χ ; in order to translate them into the pΛ DC , M Q q plane we set M χ " M Q in the heavy quark regime and M χ " Λ DC in the light quark regime.In the following we discuss the three classes of experiments that probe the DM annihilation and decay, whose limits appear in Fig. 6.

Cosmic Rays Experiments
Given the large gluequark mass needed to reproduce the DM relic density in the heavy quark regime, the strongest indirect detection bounds on DM annihilation come from the ANTARES neutrino telescope [47], HESS [48] and the multi-messenger analysis made by the Fermi gamma-ray telescope and IceCube [49].The bounds can be roughly summarized as follows 9 : xσ ann v rel y À 10 ´7 GeV ´2 (from ANTARES, HESS) xσ ann v rel y À 10 ´5 GeV ´2 (from Fermi-ICeCube) . (5.9) Indirect searches also place bounds on the lifetime of heavy DM candidates.In the highmass range, Ice-Cube provides stringest bounds [50].For M χ " p10 5 ˜10 7 q GeV, they are roughly given by τ pχ 0 q Á 10 26 ˆMχ 100 TeV ˙3{2 s . (5.10)

CMB power spectrum
The energy released by gluequark annihilations and decays around the epoch of recombination modifies the CMB power spectrum.This, similarly to indirect experiments, constrains the lifetime and the annihilation rate of the DM.The annihilation cross section is constrained to be smaller than [51] xσ ann v rel y À 10 ´8 ˆMχ 100 GeV ˙GeV ´2 , (5.11) while the limits on the DM lifetime are [52] τ `χ0 ˘Á 10 24 s . (5.12) These bounds are slightly less stringent than the ones coming from indirect detection, but have the advantage to be free from astrophysical uncertainties.They are provided for DM masses up to 10 TeV, but are expected to be approximately mass-independent for masses above this value [53].The CMB bounds shown in Fig. 6 have been obtained under this assumption.

21 cm line
While CMB is sensitive to sources of energy injection at the epoch of recombination, the cosmic 21-cm spectrum is sensitive to sources of energy injection during the dark ages.The recent observation of an absorption feature in this spectrum, if confirmed and in agreement with standard cosmology, can be used to put bounds on both the lifetime and the DM annihilation cross section.Conservative limits can be derived by neglecting astrophysical heating sources; the one on annihilation is of order [54,55]: while the one on the DM lifetime is [55][56][57] τ `χ0 ˘Á 10 25 s. (5.14 The latter is independent of astrophysical uncertainties on the distribution of DM.
As for the case of CMB, these bounds are provided up to M χ " 10 TeV and to obtain Fig. 6 we assumed that they are constant at higher masses.Differently from the previous case, however, this assumption is not completely justified and further studies are needed to provide solid bounds in the high-mass range [53].

Glueball lifetime
In the region of the parameter space that we consider, Λ DC Á GeV, glueballs with lifetime larger than 1 s are excluded by a combination of bounds.Cosmologically stable glueballs have a too large relic density and overclose the Universe.Long-lived glueballs, on the other hand, are constrained by BBN observations [58] in the range 1 s ă τ Φ ă 10 12 s and by observations of the diffuse gamma-ray spectrum [59] in the range 10 12 s ă τ Φ ă 10 17 s.
These bounds constrain the high-mass region of the V model as shown in Fig. 4. Notice, however, that they could be potentially relaxed if glueballs decay through dimension-6 operators (generated for example in models with Yukawa couplings).

Summary
In this work we continued the systematic study of gauge theories with fermions in real or vector-like representations, initiated in Ref. [1], where a DM candidate arises as an accidentally stable bound state of the new dynamics.We considered in detail the gluequark DM candidate, a bound state of adjoint fermions with a cloud of gluons, stable due to dark fermion number.What makes this scenario special in the context of accidental DM is that the physical size of DM, that controls the low-energy interactions, is determined by the dynamical scale of the gauge theory independently of its mass.In the heavy quark regime the DM mass and size can be vastly separated leading to an interesting interplay of elementary and composite dynamics.In particular, cross sections much larger than the perturbative unitarity bound of elementary particles can arise, modifying the thermal abundance and producing potentially observable signals in indirect detection experiments.
Gluequarks display a rich and non-standard cosmological history and could be as heavy as PeV if their abundance is set by thermal freeze-out.
The low-energy dynamics and the spectrum of gluequarks are non-perturbative and we were only able to give rough estimates of various effects.In particular, the quantitative estimate of the re-annihilation relevant for the thermal relic abundance and indirect detection of DM is highly uncertain, as it depends on the details of the spectrum and on the rate of non-perturbative transitions.A more firm conclusion would require a precise knowledge of binding energies of gluequarks that could be obtained by dedicated lattice studies.
Our estimates show that the observed DM density can be reproduced by gluequarks both in the light and heavy quark regimes.The mass of the DM is of order 100 TeV or larger, which makes the models difficult to be directly tested at present and future colliders.On the other hand, indirect experiments sensitive to the decay and the annihilation of the DM are a powerful probes of gluequark theories.We found that these experiments can already set important limits, excluding part of the curve which reproduces the observed DM density, depending on the value of the annihilation cross section and if the naive estimate for the gluequark decay rate is assumed (see Fig. 6).This suggests that gluequark theories in the very heavy quark regime require non-generic UV completions to ensure the accidental stability of the DM at the level of dimension-6 operators.For example, the dark parity could be gauged in the UV (see Appendix C), or its violation could be generated only by nonperturbative gravitational effects in a weakly-coupled UV completion.Similar arguments are put forward also in the context of axion models concerning the quality of the Peccei-Quinn symmetry, see [60].Assuming that an appropriate UV completion exists, gluequark models are interesting examples where the DM density can be generated thermally after inflation by very heavy particles.This can be contrasted with other scenarios, such as Wimpzillas [61], where ultra heavy DM candidates are never in thermal equilibrium.
In this work we studied gluequarks as thermal relic candidates and focused on the simplest, minimal theories of Tab. 1. Investigating non-minimal models would be certainly interesting and important under several aspects.For example, SM gauge couplings unify at high energies with less precision in the minimal blocks of Tab. 1 than in the SM.Achieving precision unification thus necessarily requires extending the models to include additional matter with SM quantum numbers.Furthermore, while the thermal relic abundance hints to a large DM mass, this conclusion can be modified in more general gluequark theories where the DM is asymmetric (this requires a larger accidental symmetry than dark parity) or where the DM abundance is determined by the decay of unstable heavier states.These theories would have a smaller mass gap and could be tested at the LHC and at future colliders.We leave an investigation along these directions to a future work.

A Dark Quark Annihilation Cross Section
In this section we report the formulas for the annihilation cross section of dark quarks, which are useful to study the perturbative freeze-out and DM indirect detection.
Dark quarks can annihilate into dark gluons and into SM final states (mainly V V , V h, hh and ψ ψ, where V " W, Z, γ).These latter contribute significantly to the total cross section in the case of perturbative freeze-out whenever M Q {Λ DC " 1 and thus the dark color interaction strength does not exceed much the electroweak one.Final states into SM particles are also expected to be important for direct detection even though they have a smaller rate compared to DM annihilation into glueballs.
The tree-level annihilation cross-section of dark quarks χ i χ j in a representation pA DC , R SM q of the dark color and SUp2q L groups into massless vectors at low energy reads, where the generators are written as T " pg DC T DC b 1q ' p1 b g SM T SM q.Selecting the neutral component in the equation above and averaging over dark color gives the perturbative annihilation cross-section of DM today.Averaging over all initial states as required for the thermal freeze-out one finds [2], where and A stands for the adjoint representation.T pRq and CpRq are respectively the Dynkin index and the quadratic Casimir of the representation R, and g χ " 2p4q for real (complex) representations.
Dark quarks can also annihilate into final states with SM fermions and Higgs bosons through their SM gauge and Yukawa interactions.These channels have been included in eq.(3.3).

B Reannihilation
As discussed in section 3.3, a second stage of annihilation involving gluequarks can occur after confinement.The annihilation can proceed in a single step into glueballs or SM vector and Higgs bosons: • χ `χ Ñ nΦ{nV : in the heavy quark regime this process has a perturbative cross section; indeed the exchanged momentum in the interaction is of OpM Q q with M Q " Λ DC , thus the interaction strength is governed by g DC pM Q q which is perturbative.
Alternatively, it can take place in two steps, a non-perturbative recombination followed by de-excitation and decay into SM particles: • χ `χ Ñ QQ ˚ SM: the recombination is two to one and energy conservation implies M QQ ˚ą 2M χ , therefore the opposite decay process is always allowed.The matrix element for the inverse decay is non-perturbative and the corresponding rate is expected to be much larger than the rate of the de-excitation process QQ ˚Ñ QQ `nΦ{nV .
• χ `χ Ñ QQ ˚`Φ{V SM: the recombination takes place with the emission of one electroweak gauge boson or, if kinematically allowed, one glueball.Bound states with M QQ ˚ă 2M χ will in general be formed which cannot decay back into gluequarks.They can de-excite and decay into SM particles.The corresponding re-annihilation rate is expected to be non-perturbative and potentially large.
Only the last of the three processes described above can ignite an epoch of re-annihilation.The dynamics of re-annihilation is described by a set of coupled Boltzmann equations of the form Hz pY Φ ´Y eq Φ q . (B.1) These expressions are simplified in that the actual system of equations involves the number densities of all possible QQ ˚bound states.Furthermore, we have omitted the effect of the recombination into EW vector bosons and of the corresponding inverse process.Equation (B.1) will be however sufficient for our discussion, and the generalization to the full case is straightforward.
In order to annihilate into SM particles with an unsuppressed rate, an excited QQ ˚bound state needs to reach first a state with low angular momentum.Consequently, re-annihilation is efficient only when the rate of de-excitation is larger than the one of dissociation 10 .
Obtaining a precise estimate of the ratio between the de-excitation and dissociation rates is difficult because: i) the dynamics of these processes is non perturbative and the lifetime depends on the different initial and final QQ ˚states considered; ii) the rate of the dissociation process initiated by EW vector bosons depends on their energy, which follows a thermal distribution and thus varies with the temperature.The result is that the re-annihilation process can be efficient for some of the QQ ˚states and inefficient for others, and it becomes more and more efficient as the temperature decreases.
This can be effectively described as a non-perturbative re-annihilation process happening with a temperature-dependent cross section that saturates to a maximal value when dissociation becomes inefficient for all the bound states.Since the evolution of the relic density takes place on relatively short time scales, the final abundance after this second freeze-out can be approximately characterized by two parameters: the final (maximal) value of the cross section, and the temperature at which this final cross section is reached.These two quantities will be dubbed respectively as the re-annihilation cross section, σ rea , and the re-annihilation temperature, T R .
During the last stage of re-annihilation, for sufficiently large Γ Φ or Γ QQ ˚, the system of equations given in eq.(B.1) simplifies.The abundance of gluequarks can be described by a single equation, see eq. (3.19).

B.1 Estimate of the Re-annihilation Cross Section
In this section we try to estimate σ rea using considerations based on energy and angular momentum conservation and simplified phenomenological models.
First of all, it is useful to determine if (depending on value of the temperature, Λ DC and M Q ) the recombination process takes place in a semiclassical or fully quantum regime.If the De Broglie wavelength of the particle λ " h{p is of order or larger than the typical interaction range R " 1{Λ DC the collision is fully quantum mechanical, otherwise a semiclassical picture can be adopted.The condition for a semiclassical behaviour can be recast as l max " M χ vR " 1, where l max is the maximum angular momentum of the process given the short-range nature of the interaction.We find that the processes occurring in the very early Universe (at T " T R ) are always in the semiclassical regime in the region of parameter space where the DM experimental density can be reproduced.Recombination processes occurring at the CMB or at later times, instead, turn out to be quantum mechanical because of the much lower gluequark velocity.
In the quantum regime, the lowest partial wave is expected to dominate in the low momentum limit k Ñ 0. In the case of exothermic reactions11 , as the one considered here, general arguments of scattering theory suggest that the cross section scales as 1{k for k Ñ 0 if the process can take place in s-wave [62].In the non-relativistic limit we thus expect a cross section σ ∝ 1{v.Since the process is non-perturbative it is not possible to compute this cross section from first principles; furthermore, since two different scales (M Q and Λ DC ) appear in the problem, it is not clear what is the cross section scaling 12 .In light of this we adopt two different benchmark scenarios: In the first one the cross-section is controlled by the size of the gluequark while in the latter is the size of the s-wave ground state which fixes the cross section.
In the semiclassical regime we estimate the re-annihilation cross section using a simple dynamical model.We first discuss the process χ `χ Ñ QQ ˚`Φ and then analyse the recombination into EW vector bosons.Our semiclassical model is defined in terms of the following simplified assumptions: • The gluequarks are modelled as hard spheres with radius of order R " 1{Λ DC , colliding with impact parameter b and thermal velocity v.
• The interaction is short range, therefore the maximum impact parameter for which an interaction occurs is b max " 2R " 2{Λ DC .We define a corresponding geometric total cross section For thermal velocities, b max can be converted into a maximum angular momentum l max " b max {M χ v " 2pM χ {Λ DC q a 3T {M χ for the colliding particles.
• Energy conservation implies that only some bound states can be formed.Among these we identify the states with maximum angular momentum l ˚allowed by energy conservation and by the short range constraint l ˚ď l max .
• Angular momentum conservation implies that only interactions with impact parameter smaller than b ˚« pl ˚`1q{pM χ vq can lead to bound state formation 13 .The short range interaction constraint then forces b ˚ď b max .If no bound state is allowed by energy conservation we take b ˚" 0.
• The recombination cross section is estimated by the geometrical value σ " πb 2 ˚.
The model predicts a re-annihilation cross section into glueballs that can be conveniently expressed in terms of a suppression factor ε Φ as follows: where ε Φ is computed to be Notice that ε Φ is a function of M Q , Λ DC and indirectly of the temperature through the value of l ˚and v.
In order to determine l ˚we use the energy balance equation in the center-of-mass frame: where K χ is the kinetic energy of the colliding gluequarks.The gluequark mass can be written in terms of the quark mass plus a binding energy B χ : Similarly, the mass of the di-quark bound state is decomposed as We set the gluequark binding energy to the value computed in QCD lattice simulations of SUp3q gauge theories: B χ " 3.5 Λ DC [23]  and implies a constraint on the maximal angular momentum l ˚(as well as on the principal number).In general one should also impose the additional condition M QQ ˚ă 2M χ , to ensure that the decay of the QQ ˚state back into gluequarks is kinematically forbidden.In terms of binding energies, this condition reads The average gluequark kinetic energy in eq.(B.10) is of order of the temperature, which in turn is smaller than Λ DC .We set the glueball mass to its value computed on the lattice in SUp3q Yang-Mills theories, M Φ » 7Λ DC , and thus find 2K χ ´MΦ ă 0. As a consequence, the condition (B.10) is always weaker than (B.11). 14The bare quark mass and binding energy are renormalization scheme dependent.Here we quote the result of reference [23] valid in the RS scheme which, according to the authors, smoothly converges to the M S scheme in the perturbative regime.Since we are interested in just an order-of-magnitude determination of the relic density, we neglect the scheme dependence of Bχ in what follows.We notice however that our numerical estimate of the re-annihilation cross section is rather sensitive to the value of Bχ, hence the scheme dependence can have a strong impact.We take such theoretical uncertainty effectively into account by considering two different benchmark values of the re-annihilation temperature, TR " TD and TR " ΛDC, as explained in section 3.3.Since eq.(B.10) depends on the gluequark kinetic energy, which we set to K χ " T in our numerical computation, the value of l ˚will have a dependence on T .For illustration we show in Fig. 7 the value of ε Φ as a function of M Q {Λ DC obtained at T " T D for Λ DC " 1 TeV.Changing Λ DC while keeping the temperature fixed leads to very small variations of ε Φ .For T " Λ DC , on the other hand, ε Φ turns out to be small and of order of a fewˆ10 ´2 in the region of interest (100 GeV ă Λ DC ă 10 TeV and 1 ă M Q {Λ DC ă 100).

���
In the case of the recombination with the emission of a vector boson, χ`χ Ñ QQ ˚`V , we expect the re-annihilation cross section to be suppressed by at least a factor α 2 .Clearly, this process becomes relevant only when the recombination with glueball emission is strongly suppressed or forbidden for kinematic reasons.The transitions χ `χ Ñ QQ ˚`V that are relevant for re-annihilation are those where the QQ ˚is sufficiently light so that it cannot decay back in two χ's.Such bound states satisfy the condition (B.11), which requires the kinetic energy of the emitted vector boson to be larger than the sum of the kinetic energies of the colliding gluequarks (K V ą 2K χ ).The re-annihilation cross section can be written as σ rea,V " ε V α 2 4π Λ 2 DC , (B.12) and the suppression factor ε V can be estimated using our semiclassical model by following a procedure similar to the one described for the case of glueball emission.We find that ε V has a behaviour similar to ε Φ as a function of its variables, and a slightly larger absolute value, see Fig. 7.

B.2 Re-annihilation temperature
The temperature at which the re-annihilation cross section saturates and the relic abundance freezes out is determined by two competing processes: de-excitation and dissociation.The cross section saturates when the former dominates over the latter for all the bound states with M QQ ˚ă 2M χ .We will now try to argue that for T ą T D there are states for which the dissociation rate is larger than the de-excitation one.Therefore, the most reasonable scenario is one in which the re-annihilation cross section saturates at T R À T D .
States with M QQ ˚ą 2M χ ´MΦ can be dissociated by glueballs with vanishing kinetic energy.Therefore, these states are the most easily dissociated since all the glueballs present in the Universe contribute to their breaking rate Γ dis " n Φ xσ dis vy , (B. 13) where σ dis " Λ ´2 DC and n Φ is the number density of glueballs which, at T ą T D , is dominated by the population coming form the confinement of dark gluons: n Φ " T 3 .The de-excitation rate can be estimated using the well known result for spontaneous emission15 where ∆E is the difference of energy levels.A reasonable estimate for this rate can be given for transitions between these states and the ground level.In this case ∆E " Λ DC ὰ2 DC M Q , while the matrix element is a fraction of the Bohr radius, r b " 1{pα DC M Q q.This estimate gives Γ QQ ˚smaller than Γ dis and suggests that for T ą T D re-annihilation cannot proceed through the formation of these states.At T " T D glueballs start to decay.Their number density decays exponentially and the dissociation process becomes soon inefficient.Therefore all the states can contribute to the re-annihilation process and the cross section saturates.
After the decay of the primordial glueballs, dissociation processes involving electroweak gauge bosons can play a role.However their cross section is suppressed by an electroweak factor and, moreover, their energy distribution is thermal.At T À T D one needs vector bosons in the tail of the Bose-Einstein distribution in order to have enough energy to dissociate the bound states.As a result, the rate of this process is exponentially suppressed by a factor expr´p2B χ ´BQQ ˚q{T s and, even if it is efficient at T D , it becomes soon inefficient.
For these reasons, we consider the case in which the reannihilation occurs at T D as the most plausible.This is in agreement with what suggested in Ref. [7].Due to the large uncertainties on the estimates of the rates, however (especially for what concerns Γ QQ ˚, where neither ∆E nor the matrix element can be computed from first principles), we do not exclude the possibility that the dissociation processes are never efficient and re-annihilation takes place directly at Λ DC .

C A model with hypercharge
In this article we focused on the minimal block V of Table 1 as a benchmark for our analysis.However, the L ' L model has many peculiarities and deserves a separate discussion.In particular, in this case the DM candidate has non-vanishing hypercharge and interacts at tree level with the Z boson.

Higher-dimensional operators
This model has a U p1q D accidental symmetry, comprising dark parity as a subgroup, under which the dark quarks L and L have charge ˘1.Differently from the other models, this symmetry is broken by higher-dimensional operators with classical dimension rO dec s " 5 of the form O dec " G µν σ µν L .
In order to have a stable DM candidate and make the model viable, one can gauge the U p1q D in the ultraviolet and break it spontaneously to the dark parity subgroup by means of a scalar field.For instance, if a scalar φ with charge 2 acquires a vacuum expectation value the symmetry is broken according to the pattern: At the scale of the spontaneous breaking only operators that are dark-parity even are generated, hence the gluequark is absolutely stable.

Z-boson mediated direct detection
Below the confinement scale, the spectrum comprises a composite Dirac fermion with SM quantum numbers 2 1{2 , whose EM neutral component is identified with the DM.The non-zero hypercharge induces a tree-level interaction with the Z boson which is strongly constrained by direct searches.The corresponding spin-independent elastic cross section on nuclei N is given by [65]: where Z and N are the number of protons and neutrons in the nucleus N and M N is its mass.This cross section is excluded by direct detection experiments for masses M χ À 10 8 GeV [46].This bound rules out the model in the scenario T R " T D , corresponding to the left panel of Figure 4, but can be satisfied in the scenario T R " Λ DC .
In fact, the constraint from direct detection experiments can be also avoided by introducing a heavy singlet gluequark.In this case the presence of Yukawa couplings induces a splitting among the electromagnetically neutral Majorana fermions.The DM is the lightest among these fermions, so that tree-level elastic scattering mediated by the Z boson cannot exist due to its Majorana nature.Inelastic scatterings are kinematically forbidden if the splitting is large enough; this is easily realized for M N À y 2 ˆ10 5 TeV, where y is the Yukawa coupling.This scenario is analogous to Higgsino DM in supersymmetry, see [1,2] for an extensive discussion.

Q ⇤ 1 DCFigure 1 :
Figure 1: Cartoon of the gluequark DM candidate.A heavy fermion in the adjoint of color gives rise to color singlet state surrounded by a gluon cloud of size 1{Λ DC " 1{M Q .

Figure 2 :
Figure2: Left: Sketch of the non-standard apT q scaling.Values of α ą 1 in eq.(3.8) imply a scale factor at freeze-out, a f.o., smaller than the one obtained from a standard cosmology, ãf.o. .This in turn leads to a suppression of the relic density by a factor pa f.o.{ã f.o.q 3 " 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I 3 4 6 K 6 t r 6 x u b p a 3 y 9 s 7 u 3 n 7 l 4 L B l o 8 Q I

Figure 3 :
Figure 3: Cartoon of the re-annihilation processes occurring after dark confinement.First, free quarks combine into color singlets gluequarks.Next, fast collisions form excited QQ states

ΛFigure 7 :
Figure 7: Suppression factors ε Φ (left panel) and ε V (right panel) at T " T D as a function of the ratio M Q {Λ DC for Λ DC " 1 TeV.In red the results of the numerical computation and in black the interpolation used to compute the relic density.Discontinuities in the numerical results are due to the integer nature of l ˚.
[64]he binding energy of the QQ ˚bound state, B QQ ˚, is instead approximated by the energy levels of a confining model with a Coulomb potential plus a linear term[64]QQ ˚is computed numerically as a function of the principal and orbital quantum numbers of the bound state. Th energy balance of eq.(B.6) can be rewritten as B QQ ˚ď 2B χ `2K χ ´MΦ , (B.10)