Gluequark Dark Matter
Abstract
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 freezeout is followed by a nonperturbative reannihilation 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.
backgroundcolor=gray!10, skipabove = 10pt, skipbelow = 20pt, innertopmargin=20pt, innerbottommargin=20pt, linewidth = 0pt, innerleftmargin = 20pt, innerrightmargin = 20pt, roundcorner=30pt
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 cutoff, 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 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. 1503.08749 (); 1707.05380 () of theories where the dark sector comprises new quarks transforming as real or vectorlike representations under both the SM and dark gauge groups, and where the Higgs field is elementary. Such framework, also known as VectorLike Confinement 0906.0577 (), provides a safe nontrivial 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 Barducci:2018yer ().
Previous studies of accidental DM focused on baryons or mesons of the dark dynamics as DM candidates, see 1604.04627 () for a review. A systematic analysis of models with baryonic DM was performed in Refs. 1503.08749 (); 1707.05380 (). In the regime where dark quark masses are below the dark color dynamical scale, , it was found that a correct relic abundance can be obtained for DM masses of order 100TeV. Lighter DM masses, down to , are instead allowed if . Mesonic DM candidates can be even lighter and can arise for example as pseudo NambuGoldstone 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 UVsuppressed dimension6 operators thus ensuring cosmologically stable gluequarks for sufficiently large cutoff 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, , 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. Boddy:2014yra (); Boddy:2014qxa () in the context of supersymmetric gauge theories. There they arise as the partners of glueballs after confinement and were consequently called glueballinos. Ref. Boddy:2014yra () 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 nongravitational probes. In this case the dark and visible sectors stay in thermal equilibrium until relatively low temperatures, of the order of GeV, and the thermal history of the Universe is rather different than that described in Refs. Boddy:2014yra (); Boddy:2014qxa (). 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. Gudnason:2006ug (); Kouvaris:2007iq (). 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}^{1}1Reference Kouvaris:2007iq () for example considered gluequark DM in the context of the socalled Minimal Walking Technicolor model, but its estimate of the thermal relic abundance focuses on the perturbative freezeout and does not include any of the nonperturbative effects described in this work..
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.
2 The models
We consider the scenario in which the SM is extended by a new confining gauge group (dark colour), and by a multiplet of Weyl fermions (dark quarks) transforming in the adjoint representation of and as a (possibly reducible) representation under the SM group :
(2.1) 
In particular, we consider models where the dark quarks have quantum numbers under but are singlets of color. New fermions colored under 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. DeLuca:2018mzn ()). We assume to be a real or vectorlike representation, so that the cancellation of 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 and the value of the dark gauge coupling at . A violating term can also be included but does not play an important role in what follows. The renormalizable lagrangian thus reads
(2.2) 
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 , dark parity, is an accidental invariance of the renormalizable Lagrangian. It is what remains at the quantum level of the anomalous 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. Cirelli:2005uq () 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 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 cutoff 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 and by the dimension of the SM representation , i.e. by the number of Weyl flavors . Models with too large or imply too low Landau poles for , and are thus excluded from our analysis. The list of minimal blocks that satisfy our requirements is reported in Table 1 for , and dark color groups.
Quantum numbers  Accidental  Classical  

Symmetry  
1  All  All  All  6  
3  6  
4  5  
6  6 
Each block is characterized by its accidental symmetry (that can be larger than the dark parity) and by the dimensionality of the lowestlying operator which violates it. The latter has the form
(2.3) 
where is a SM composite operator matching the quantum numbers of the dark quark . 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 stronglycoupled IR fixed point at energies , the dimension of is simply given by , as reported in the sixth column of Table 1. Among the minimal blocks, the model has 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 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 and on the value of the dark coupling at the scale . One can identify two regimes. In the first, is perturbative and this necessarily implies a confinement scale smaller than the quark mass, ; we will call this the ‘heavy quark’ regime. In this case, depending on the value of , there are three possible behaviours. For 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. ^{2}^{2}2Notice that the value of , in the case of adjoint fermions, does not depend on the gauge group. For , where is the nonperturbative 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 the coupling grows strong quickly in the infrared and confinement is triggered without delay. Only for this latter range of values of the confinement scale can be larger than the quark mass, ; 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 gauge theories with fermions in the adjoint representation was extensively studied through lattice simulations, see for example Catterall:2011zf (); Bergner:2016hip (); Bergner:2017gzw (); Athenodorou:2014eua (); Bergner:2015adz (); DeGrand:2015zxa (); Ali:2016zke (); Bergner:2017bky (); Nogradi:2016qek (); DeGrand:2013uha () and references therein. There seems to be sufficient evidence for an infrared conformal phase of theories with colors and massless Weyl flavors, while results with are more uncertain though still compatible with a conformal regime. Theories with are supersymmetric and have been shown to be in the confining phase. The case with 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 , while the existence of a weaklycoupled infrared fixed point can be established for 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 () can be obtained from lattice simulations with adjoint static sources, see for example Refs. Foster:1998wu (); Bali:2003jq ().
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 theories show that the state is the lightest with mass , see for example Morningstar:1999rf (). Similar values are found for with different number of colors. The gluequark is expected to be the lightest state made of quarks, with a mass . 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 parityviolating operator has dimension 6 the main decay channel for the lightest gluequark is (where indicates a glueball and ). In the model of Table 1 with three dark flavors transforming as an EW triplet, the dim6 operator
induces the decay of the gluequark with inverse lifetime
(2.4) 
Similar results apply for the and 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 dimension8 operators of the form generated at the scale . For all the minimal models in Table 1 this is the lowestdimensional operator which induces glueball decay. The partial width into photons is estimated to be 0903.0883 (); 0911.5616 ()
(2.5) 
When phase space allows, the decay channels , and open up producing one orderofmagnitude smaller lifetime. Relatively longlived 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 1loop radiative effects at the scale generate the dimension6 effective operator , inducing a much shorter lifetime. If their mass is high enough, , glueballs predominantly decay to two Higgs bosons with a decay width
(2.6) 
where , and , are respectively the Yukawa couplings and masses of the dark quarks circulating in the loop. Lighter glueballs can decay through the mixing with the Higgs boson; as for the Higgs, the dominant channel for GeV is that into bottom quarks, with a corresponding partial width ^{3}^{3}3The scaling is approximately correct for , though eq. (2.7) is a good numerical estimate for as well.
(2.7) 
Light quark regime:
If dark quarks are lighter than , the physical spectrum is radically different and one expects spontaneous breaking of the global symmetry down to . The lightest states are thus the (pseudo) NambuGoldstone bosons , while the DM candidate is the gluequark, accidentally stable and with a mass of the order of the confinement scale . As discussed in section 4, and similarly to the baryonic DM theories of Ref. 1503.08749 (), reproducing the correct DM relic density in this regime fixes . The NGBs with SM quantum numbers get a mass from 1loop electroweak corrections, which is predicted to be for the value of of interest. Besides such a radiative correction, the quark mass term breaks explicitly the global symmetry and gives an additional contribution. Including both effects, the NGB mass squared is given by
(2.8) 
where is the weak isospin of the NGB and are coefficients.
For fermions in the adjoint representations, only models with light quarks can be in the regime , 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 model and the model. The model has a global symmetry breaking which leads to five NGBs transforming as an electroweak quintuplet. In the model one has and nine NGBs transforming as of .
The limit of very small quark masses, , 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 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 higherdimensional operators, possibly resulting into longlived particles. An example of this kind is given by the model, where NGBs made of or constituents have number 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 ()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 , as required to reproduce the DM relic abundance 0906.0577 (); 1707.05380 (); Barducci:2018yer ().
Besides the NGBs, the physical spectrum comprises additional bound states with mass of order . 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}^{4}4The existence of stable baryons in theories with adjoint fermions was investigated in Refs. Bolognesi:2007ut (); Auzzi:2008hu (), where stable skyrmion solutions were identified and conjectured to correspond to composite states with mass of , interpolated from the vacuum by nonlocal operators. We will not include these hypothetical states in our analysis. In the light quark regime they are expected to annihilate with a geometric cross section and contribute a fraction of DM relic density comparable to that of the gluequarks.. 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 dimension6 operator , the most important decay channels of the gluequark are and . The twobody decay dominates at large and gives a lifetime of order
(2.9) 
where is the decay constant of the gluequark ^{5}^{5}5This has been defined by , and scales as in the large limit..
To summarize our discussion on models, Table 1 reports the minimal blocks which have a potentially viable DM candidate and a sufficiently high cutoff, above , as required 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 Boddy:2014yra (); Boddy:2014qxa (), and it will not be considered further in this work. We find that in all the other minimal blocks of Table 1 the 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 model of Table 1, i.e. the minimal block with an triplet. For a discussion of the model see Appendix C.
3 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. 1503.08749 (). Dark color confines when dark quarks are relativistic and in thermal equilibrium. After confinement the gluequarks annihilate into NGBs with a nonperturbative cross section , while glueballs are heavy and unstable. At temperatures 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 (), free dark quarks annihilate into dark gluons and undergo perturbative freezeout at (see section 3.1). At confinement (), the vast majority of the remaining dark quarks hadronizes into gluequarks, while the plasma of dark gluons is converted into a thermal bath of nonrelativistic 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 Scherrer:1984fd (); Kamionkowski:1990ni (); McDonald:1989jd (), and recently reconsidered in hepph/0005123 (); 1506.07532 (), decays of nonrelativistic particles with a large and nonthermal 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 is modified into . During this early epoch of matter domination, the Universe expands faster than in the radiationdominated 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 process ^{6}^{6}6An analogous mechanism was first discussed in Ref. hepth/0405159 (); hepph/0611322 (); 0712.2681 () and more recently by Ref. DeLuca:2018mzn () in the context of DM charged under SM strong interactions.
(3.1) 
where is an excited bound state of dark quarks and 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 with size 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 antihydrogen scattering Morgan:1970yz (). 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 . Once formed, the can either decay () or be dissociated back into two gluequarks by interactions with the SM and glueball baths (). 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, , in which the is produced with a large angular momentum, . Bound states with take more time to deexcite 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).
3.1 Thermal freezeout
Thermal freezeout 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 ) or of gluequarks (for ) is reduced until it becomes so low that chemical equilibrium is no longer attained and freezeout takes place. The number density at freezeout is approximately given by
(3.2) 
where 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 freezeout temperature is of order . A general expression for the annihilation cross section is reported in Appendix A, see eq. (A.2). For the model with analysed in the next section, the annihilation cross section into dark gluons and SM fields is
(3.3) 
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: , , refer respectively to the 1, 8 and 27 color channels and are given by 1402.6287 (); 1702.01141 ()
(3.4) 
In the light quark regime gluequarks annihilate into NGBs with a cross section that is expected to scale naively as
(3.5) 
in analogy with nucleonnucleon scattering in QCD Zenoni:1999st (). NambuGoldstone bosons are unstable and later decay into SM particles.
3.2 Dilution
As well known, the number density of DM particles today is related to the number density at freezeout by
(3.6) 
This relation is usually rewritten in terms of temperatures assuming that between freezeout and today the standard scaling 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 grows faster than . In this case the relation between and is given by:
(3.7) 
where and defines the temperature interval during which the nonstandard 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 latetime decays of dark glueballs can give rise to a nonstandard scaling of the form with . The corresponding suppression factor thus reads:
(3.8) 
After dark color confinement, the energy density of the Universe can be divided into a relativistic component, , containing all the SM relativistic particles, and a nonrelativistic one, , containing all the darksector longlived 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 nonrelativistic species, and this can lead to an early epoch of matter domination. The evolution of is governed by
(3.9) 
where is the glueball decay rate and the Hubble parameter is given by the Friedmann equation:
(3.10) 
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 are given by
(3.11) 
where and 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 is related to the relativistic energy density by:
(3.12) 
The evolution during the early matterdominated epoch, if the latter exists, can be described by solving analytically eq. (3.9) at leading order in for cosmic times 1506.07532 ():
(3.13a)  
(3.13b) 
Here and denote the initial conditions at some time much after the beginning of the matterdominated epoch. The relativistic energy density is given by the sum of (first term in eq.(3.13b)), diluted as , and the energy injected by glueball decays (second term in eq.(3.13b)), diluted as . Initially the first term dominates and the standard scaling is obtained; as long as the glueball lifetime is long enough, the second term will start to dominate at some temperature , implying a nonstandard scaling (see Fig. 2). The value of can be found by equating the first and second terms of eq.(3.13b) and by using eqs.(3.10),(3.11):
(3.14) 
The nonstandard scaling ends when almost all the glueballs are decayed, i.e. around , where is the time at dark confinement. Using eqs. (3.10) and (3.11), one can translate this condition in terms of a temperature finding:
(3.15) 
From eq.(3.8) it follows that latetime decays of glueballs dilute the naive relic density by a factor
(3.16) 
where 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 is very well approximated by:
(3.17) 
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.
3.3 Reannihilation
At 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 hadronize into gluequarks. These can then collide and recombine in excited 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 , one expects naively a recombination cross section of order . 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.
Once formed, states with mass will promptly decay back to two gluequarks. Lighter states, on the other hand, can either deexcite and thus decay into SM particles through the emission of a SM vector boson or a glueball (), or be dissociated by interactions with the glueball and SM plasmas (), see Fig. 3.
If deexcitation occurs faster than dissociation, a second era of efficient DM annihilation can take place, reducing the gluequark number density. While reannihilation processes can be active over a long cosmological time interval, it is the last stage during which the reannihilation cross section gets its largest value that is most important to determine the final gluequark density. This last stage happens relatively quickly and can be characterized by a reannihilation temperature . The exact value of depends on the rate of dissociation and is difficult to estimate. The largest uncertainties arise from the calculation of the deexcitation 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 reannihilation 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 states occurs faster than deexcitation, as long as the glueball bath originating from dark gluons confinement is present; therefore, the reannihilation temperature is approximately equal to the one at which glueballs decay (). Besides this most probable scenario, in the following we will also consider the other extreme possibility where reannihilation occurs right after confinement (). The comparison between these two opposite scenarios will account for the theoretical uncertainties intrinsic to the determination of the nonperturbative dynamics characterizing our DM candidate.
In both benchmark scenarios considered above the last stage of the reannihilation 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 deexcitation or glueball decay rates. This reads
(3.18) 
where , and is the entropy density of the Universe. The equilibrium term can be neglected since . Assuming a reannihilation cross section which is constant and velocity independent ^{7}^{7}7As explained in Appendix B, the last stage of reannihilation can be effectively described by a constant cross section; the latter turns out to be also velocity independent in the relevant region of the parameter space of our theories., eq. (3.18) can be easily integrated analytically; one obtains (for )
(3.19) 
Latetime annihilation significantly affects the gluequark relic density when the second term in the above equation dominates, i.e. roughly when
(3.20) 
in agreement with a naive expectation. When condition (3.20) is met, any dependence from the previous stages of cosmological evolution, encoded in , is washed out and the asymptotic value of the relic density is set only by reannihilation. For temperatures sufficiently smaller than (but higher than a possible subsequent period of dilution, in the case ), eq.(3.19) can be recast in terms of the gluequark relic density as follows:
(3.21) 
4 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 and . For each point in the plane 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 (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 nonperturbative nature of the processes involved. As anticipated in section 3.3, one of the largest uncertainties comes from the identification of the reannihilation temperature . We will consider the two previously discussed benchmarks: , the most plausible according to our estimates, and . We reconstruct for each of them the different possible cosmological evolutions obtained by varying and . Our estimate of the DM abundance for both benchmarks is reported in Fig. 4, where we show the isocurve reproducing the experimentally observed density.
Let us consider first the case . There are three possible thermal histories that can be realized (they are correspondingly indicated in the left plot of Fig. 4):

For very large the Universe undergoes a first perturbative freezeout at , then dark confinement occurs at followed by an epoch of dilution between and ^{8}^{8}8Here we are implicitly assuming that the reheating temperature at the end of inflation is larger than , so that the number density of dark quarks after the perturbative freezeout is thermal.. Glueballs decay at , and the number density of gluequarks is too small, as a consequence of the dilution, to ignite a phase of nonperturbative reannihilation. The DM density is therefore given by
(4.1) where the number density at freezeout 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 one obtains
(4.2) 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.

For smaller values of (but still with ), the dilution between and is not enough to prevent reannihilation (i.e. condition (3.20) is met). The latter thus occurs at , washing out any dependence of from the previous stages of cosmological evolution; the corresponding DM relic density is
(4.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 reannihilation cross section by using the semiclassical model described in Appendix B.1; this gives
(4.4) The parameters and are smaller than 1 and encode the suppression from energy and angular momentum conservation respectively for the recombination processes and . While eq.(4.4) is the result of a rather sophisticated analysis of the reannihilation dynamics and represents our best estimate for , it is subject to large theoretical uncertainties, as discussed in Appendix B.1. We thus also consider the extreme situation where the reannihilation cross section is always large and saturated by its geometric value
(4.5) Varying between the values in eqs.(4.4) and (4.5) will quantify the uncertainty on . By using eq.(3.21) and setting , eq. (4.3) takes the form:
(4.6) This formula describes the intermediate part of the isocurve in the left plot of Fig. 4. Initially (i.e. for ) the reannihilation is dominated by the process and ; in this case the last factor in (4.6) can be well approximated with 1 (the electroweak contribution to is small) and the estimated uncertainty on the gluequark relic density is negligible. For larger reannihilation into plus a glueball becomes kinematically forbidden in our semiclassical model, and quickly drops to zero (see Appendix B.1). In this region and varying between and spans the gray region. The extension of the latter quantifies the uncertainty of our estimate of the relic density.

When , the perturbative freezeout does not take place. If is bigger than , then the Universe undergoes a first epoch of annihilation of dark quarks for , followed after confinement by the annihilation of gluequarks, until thermal freezeout of these latter occurs at . If , 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 freezeout at . Afterwards 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 . Setting , one obtains
(4.7) For the nonperturbative 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 . This implies , so that the DM relic density turns out to be independent of . If instead the reannihilation cross section is estimated by eq.(4.5), then by continuity with the previous cosmological evolution one must take , which also corresponds to a relic density independent of . Varying 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, , the annihilation of gluequarks proceeds through the direct annihilation of their constituents (the theory at is nonperturbative) with a cross section , where is an order 1 coefficient. We vary 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 . The observed relic density in this regime is reproduced for TeV, similarly to the light quark regime in baryonic DM models 1503.08749 ().
Let us turn to the case . As for one can identify three possible thermal histories (correspondingly indicated in the right panel of Fig. 4):

For the Universe goes first through a perturbative freezeout of dark quarks at , then reannihilation occurs right after confinement for . Finally, dilution takes place between and the temperature of the glueball decay . The DM relic density is given by the expression in eq.(4.3) times the dilution factor . Numerically one has
(4.8) In this case, our semiclassical model estimates throughout the parameter space of interest. By varying between and we thus obtain the upper portion of the gray region in the right plot of Fig. 4.

For smaller (but still with ), the glueballs are too short lived to ignite the dilution, and the DM relic density is given by eq.(4.3). Setting one obtains
(4.9)
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 .
5 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. 1503.08749 (); 1707.05380 (). 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.
5.1 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 global symmetry breaking in the light quark regime. In the case of the model, the five NGBs form a multiplet with weak isospin 2, and one expects . The phenomenology of a quintuplet of NGBs was studied recently in Ref. Barducci:2018yer (). These states are pair produced at hadron collider in DrellYan processes through their electroweak interactions, and decay to pairs of electroweak gauge bosons through the anomalous coupling
(5.1) 
A promising discovery channel studied by Ref. Barducci:2018yer () is ; the doubly charged states decay into samesign pairs and are somewhat more challenging to see experimentally. The LHC has an exclusion reach up to TeV masses, while a 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 higherdimensional operators, and are rather elusive at colliders. In models without Yukawa couplings, where interactions with the SM occur through dimension8 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 TeV collider is of order fb for and GeV). In models with Yukawa couplings the glueballs mix with the Higgs boson and production via gluonfusion 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 highintensity experiments like SHiP can only probe light glueballs in a region of parameter space that is already excluded by EW precision tests 1707.05380 ().
Mesons can give interesting signatures in both light and heavy quark regimes. Bound states made of a pair of dark quarks, , can be singly produced through their EW interactions. While the production of spin0 mesons is suppressed since they couple to pairs of EW gauge bosons, spin1 resonances mix with the SM gauge bosons of equal quantum numbers and can be produced via DrellYan processes. In the narrow width approximation the cross section for resonant production can be conveniently written in terms of the decay partial widths as
(5.2) 
where is the dimension of the representation, the spin, the parton producing the resonance and are the dimensionless parton luminosities.
In the heavy quark regime the bound state is perturbative and its decay width can be computed by modelling its potential with a Coulomb plus a linear term. For the decay width of the lowestlying swave bound states scales as
(5.3) 
where originates from the nonrelativistic Coulombian wavefunction. When , 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 1707.05380 (). For example, in the model the decay width of the swave spin1 resonance (isospin 1 in light of the Majorana nature of ) into a lefthanded fermion doublet is
(5.4) 
where refers to the radial quantum number. The tiny energy splitting between levels is irrelevant at colliders and the total rate is dominated by the lowestlying Coulombian ones. The branching ratio into pairs of leptons is about and the strongest bounds currently arise from searches of spin1 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 TeV depending on the ratio (or equivalently on the value of ).
In the light quark regime the lightest spin1 state is the meson with mass . The widths scale as Contino:2006nn ()