The Thermal History of Composite Dark Matter

We study the thermodynamic history of composite Dark Matter models. We start with classifying the models by means of the symmetries partially protecting the composite Dark Matter decays and constrain their lifetimes. For each model, we determine the impact of number-changing and number-conserving operators on its thermal history. We also develop the analytic formalism to calculate the asymptotic abundance of stable relics. We show how the relative strength between number-changing and number-conserving interactions together with the dark plasma lifetime affect the thermal fate of the various composite models. Additionally, we discover that the final dark relic density of composite particles can be diluted due to an entropy increase stemming from dark plasma decay. Finally, we confront the models with experimental bounds. We find that indirect detection experiments are most promising in testing this large class of models.

Here we focus on the thermal history and ultimate fate of some of the prime candidates of composite dark matter. A generic property of composite models of dark matter is that, at low energies, they feature a complex multicomponent structure. When taken into account this feature plays a crucial role in thermal history influencing the thermal bath.
The paper is structured as follows. In section II we identify the composite dark matter candidates arising in gauge theories. We classify them according to the symmetries protecting their decays. We further discuss higher dimensional operators inducing composite dark matter decay that can emerge in a full ultraviolet theory including gravity. We use these operators to estimate the composite dark matter candidates lifetime. Section III is devoted to the thermal history of the composite dark sector. For each model we elucidate the impact of both number-changing and number-conserving operators. In section IV we construct solutions to the Boltzmann equations allowing us to determine the asymptotic abundance of stable relics. This preparatory phase is then exploited in section V to determine the relic abundance in each model. Here we provide analytic expressions for the asymptotic relic abundance stemming from each scenario. We observe that the relative strength between number-changing and conserving interactions determines the thermal fate of the various composite models and further depends on the dark plasma lifetime. We find that the final composite dark relic density can be diluted due to an increase in entropy stemming from dark plasma decay. We constrain the models via the experimentally observed relic density. The entropy ratio between the standard model and the dark sector is another parameter affecting the ultimate relic abundance which we fix by assuming the two plasma to be in thermal equilibrium in the early universe. Section VI confronts the model predictions with experiments such as H.E.S.S. and Fermi-LAT. Finally, we offer our conclusions in section VII. In the appendix A we provide the effective Lagrangian for the dark pion dynamics, in appendix B and C we offer a glossary of multi-fluid thermodynamics and in appendix D we summarize the Boltzmann equation for number-changing processes.  Table I: Field content of the lightest fermions in the model.

II. COMPOSITE DARK MATTER CANDIDATES
A particle description of dark matter (DM) requires a lifetime above τ > 10 26 −10 28 s [20,21] which is orders of magnitude the age of the universe (10 17 s). A near exact symmetry is a time-honored way to ensure long lifetimes. A famous example is the baryon symmetry that protects the proton from decaying. For this reason, we investigate here models of composite DM stemming from the dynamics of Yang-Mills gauge theories featuring fermionic matter. In these models, we identify potentially long-lived composite relics and systematically investigate their cosmological consequences.
Historically, models of composite DM, such as the ones inspired by Technicolor [7,[22][23][24], linked the baryon asymmetry to a potential DM asymmetry in order to achieve the experimentally observed dark relic density. Another class of composite DM models was introduced to generate the desired relic density via number-changing operators in the darkpion sector [25]. Here we consider composite DM models in which the primary source for the observed relic density is due to a thermal production mechanism.

A. The Composite Framework
Assuming dark asymptotic freedom we expect at low energies the theories to generate a On general grounds, part of the dark building blocks will carry standard model (SM) quantum numbers that can bring into thermal equilibrium the dark and the SM sectors.
Since the fermions are vector-like, an explicit common mass term is allowed, so that L ⊃ M Q QQ. Depending on whether M Q Λ D (strongly coupled regime) or M Q Λ D (weakly coupled regime) the composite states will behave differently. In the latter case the low energy theory is a pure Yang-Mills with the confinement scale that can be estimated to be: • The near Nambu-Goldstone bosons of the broken global symmetry (pions).
• The dark gluon bound states (glueballs), with the lightest state being a scalar singlet.
Another class of interesting composite states is constituted by composite baryons. The reason being that, even though they are not LCPs, the dark baryon number is protected.
As a teaser, we notice that in both limits the LCPs will feature number-changing operators generating intriguing thermodynamics.

B. Composite Particles and their Lifetimes
Here we summarise the composite particle spectra and discuss their lifetimes due to the possible breaking of the accidental global symmetries via higher order operators. The spectrum and operators are summarised for the reader's convenience in Table II.
As mentioned above, among the composite particles we have dark baryons Q N D ∼ Λ (for this case we consider N D odd), dark mesonsQQ ∼ Λ 2 D π and last but not the least dark glueballs G µν G µν ∼ Λ 3 D S. Their decays can occur via higher dimensional operators such as the ones envisioned in Table II and generated by new dynamics or gravitational interactions.
We also introduced the interpolating low energy fields with the help of dimensional analysis.
Due to the dark flavor structure of the models, the mesons can carry further accidental symmetries as well. Examples are: The species number, which would make π ± stable in the SM were not for the presence of weak interactions that violate this symmetry; Generalised G-parity [26], which would stabilize π 0 but it is broken by the chiral anomaly that induces  Here Λ B / , Λ G / and Λ S / are the scales at which baryon number, G-parity and species number are violated, respectively.
its decay into photons. Additionally, new physics can lead to induced higher dimensional operators that further break these symmetries. These, however, are suppressed by the scale of the new physics. Using dimensional analysis for the dark sector and constraining ourselves to dimension five operators we naively estimate the decay rates to scale as Γ ∝ m 3 π /Λ 2 for the dark mesons. Intriguingly even for Λ ≈ M Pl one discovers that the experimentally required lower bound on the dark lifetime induces an upper bound [20,21] of m π 10 MeV for this type of DM.
In general, within this framework, several relics can be long-lived on cosmic scales. Their abundances have to satisfy the overclosure bound i Ω i ≤ Ω DM . Particles with lifetimes shorter than the current on DM (≈ 10 28 s) must decay before the production of light elements in the early universe i.e. τ < 1 s as constrained by successful Big Bang Nucleosynthesis (BBN) [27].
We find four interesting scenarios: • If M Q Λ D the glueball is the LCP and we take it to be the DM candidate. The dark baryon mass is independent on the glueball mass m B ≈ N D M Q . The dark baryons can be made heavy and decay before BBN.
• Within the above mass hierarchy one can find a parameter space of the theory where dark baryons are the DM candidate. Here the glueballs will have to decay before BBN via, for example, SM residual interactions. Nevertheless being still the LCPs they now provide a thermal bath of massive particles with a finite lifetime.
• If there are light fermions in the spectrum with M Q Λ D , the LCPs are dark pions, which generically will not be long-lived on cosmological time scales. Thus the pions have to decay before BBN and provide a plasma with finite lifetime as well, taking up the glueball role. Here the dark baryons are strongly coupled states that we take them to be the DM candidate with a mass of the order of N D Λ D .
• A small fraction of the above parameter space can exhibit long-lived pions when their masses are below the GeV scale. This is indeed the case of the simplest miracle [25] where the baryons, however, are not present in the spectrum of the theory because of the choice of the underlying dynamics. If dark baryons are present one will have to consider the presence of additional stable relic states.
As sketched in Fig. 1, we observe that the ingredient which is common to all scenarios, is the existence of an LCP plasma made of massive particles (glueballs/pions). This plasma will generically lose thermal contact with the SM and undergo a distinct thermal history once the LCPs become non-relativistic [28,29]. Due to the finite lifetime of the LCPs their decay will partially restore the thermal contact with the SM at late times. To investigate the resulting intriguing phenomenology, we will briefly review the unusual thermodynamics of a non-relativisitc dark plasma.

III. DARK THERMODYNAMICS
In DM theories with a nontrivial spectrum with a mass gap, the freeze-out will show strong deviations from the standard scenario, as discussed in [28,29]. Here we will summarise the conditions which lead to non-standard freeze-out and that have to be satisfied before the DM annihilation reactions decouple.
• First of all dark sector particles must not be in kinetic equilibrium with SM particles at the freeze-out temperature, meaning that at that temperature scale the reactions Figure 1: Schematic overview of the thermal history, of a non-abelian dark sector. At high temperatures, we expect the systems to be in thermal contact. After this connection is lost, the dark force confines and forces the dark quarks in dark-color neutral bound states. Baryons (B) and glueballs (S) are non-relativistic at confinement, pions (π) on the other hand become non-relativistic substantially later. Relativistic particles are denoted by red ovals, non-relativistic particles are represented by purple ovals. The number-changing interactions can heat up the LCP plasma, but eventually, they decouple at T c and a phase begins, where dark particle number in the plasma is conserved. Dark baryon freeze-out can happen in the number-changing and number-conserving era of the dark plasma and the freeze-out process will be different. Finally, the decay of the dark plasma to SM final states partially re-introduces a thermal link of the systems at T decay . The decay of the LCPs can furthermore lead to entropy injection into the SM bath and dilute the dark baryon relics.
χSM → χSM have to be weak. This condition is generically satisfied if the lightest composite particles are SM singlets. We will also discuss a few exceptions to this scenario.
• Furthermore, the lightest particle of the dark sector (here the LCPs) is stable on cosmological scales i.e. Γ LCP H(T ).
• The freeze-out temperature is below the LCP mass so that the LCP is non-relativisitc.
• If the number-changing processes are strong enough to maintain chemical equilibrium at freeze-out i.e. σv 2 3→2 n 2 LCP H(T ) the resulting hotter dark plasma leads to a non-standard freeze-out.
• If number-changing processes are inefficient at the dark freeze-out, a chemical potential develops in the dark sector. This is a new phenomenon that must be taken into account for the proper determination of the DM relic abundance.
Thermodynamics will be extremely helpful to describe the evolution of this systems, for a summary of multi-fluid thermodynamics, see appendix B. In particular entropy conservation plays an important role since the dark sector entropy is dominated by the lightest dark degrees of freedom comprising the dark plasma. The ratio of the entropy densities is fixed ζ = s SM /s D , and generically set by weak interactions that keep the SM and the dark sector in kinetic equilibrium at high temperatures. After those interactions decouple the entropies are separately conserved, as we assume both sectors to evolve adiabatically. While this assumption is justified for the SM, for the dark sector one has to be more careful. For example within the scenario in which M Q Λ the pure dark Yang-Mills theory at low energies undergoes a first order phase transition as a function of the temperature when N D ≥ 3. In this case the entropy ratio established at high temperature only provides an upper bound ζ < s SM /s D | T high .
We can now differentiate two profoundly different regimes depending on the relative strength of the number-changing process.
In the following, we will refer to the generic LCP as π although the formalism applies to dark LCP glueballs as well. Of course, the difference in the detailed dynamics will be taken into account when concrete models are considered.

A. Non Conserved Dark Particle Number
Generically at large temperatures, the number-changing processes will maintain chemical equilibrium in the dark sector and enforce a vanishing chemical potential. The crucial feature of this regime is that while the temperature in the SM sector scales as T SM ∝ 1/a the dark sector temperature scales as T D ∝ 1/ log (a), which leads to an exponential temperature difference between the two sectors. This is in strong contrast to the case where the dark sector is dominated by relativistic degrees of freedom (say a dark photon) and the temperatures are linearly dependent. We thus refer to the above possibility as the hot dark plasma phase.
The relevant thermodynamic features follow from entropy conservation. Since the LCPs are a non relativistic fluid with zero pressure P = 0, we have s D = ρ D /T D ≈ m π n π /T D . At the same time the SM entropy density is s SM = 2π 2 /45 g SM T 3 SM . Now, from s SM = ζ s D we get the relation This shows that the dark sector gets exponentially hotter than the SM at scales where the relevant degrees of freedom are the LCPs counted by g π . Clearly a dramatic feature of the hot dark plasma era.
An important quantity we will need is the Hubble rate as a function of the dark sector temperature. Solving the Friedmann equation we obtain with z π = m π /T D .
The Hubble rate is an essential ingredient for the DM freeze-out computation, we are ultimately interested in. We now discuss what happens at temperatures at which the numberchanging processes cannot maintain the chemical equilibrium in the plasma.

B. Conserved Dark Particle Number
The number-changing operators are bound to switch off when the reaction rate drops below the Hubble rate. This happens at a critical temperature T c . Below this temperature, the dark particle number is conserved along with its entropy.
Keeping into account the conservation of the comoving entropy for a non-relativistic pressureless gas we have the following conservation law: Here µ is the chemical potential for the LCP. Solving the first order differential equation we obtain: where we used as a initial condition µ(T c ) = 0.
Furthermore, using with H 2 ∼ 1/a 3 , we verify that the dark temperature quickly drops as a function of the scale factor: T D ∝ 1/a 2 . This scaling holds independently whether SM or dark sector dominates.
The critical temperature, at which the number-changing processes decouple is found by solving H(T D ) = σ 3→2 v 2 n 2 π , which is equivalent to where we defined z π c = m π /T c .
Given that the dark sector entropy is dominated by the LCP, imposing entropy conservation ζ = s SM /s π and taking into account the chemical potential one derives the new relation between the SM temperature and the dark temperature of the two fluids: The interesting feature is that now the two temperatures are no longer exponentially related.
Furthermore, we have the following expression for the Hubble rate as a function of dark temperature, expressed in terms of z π , in this regime Finally, to determine when either of the two fluids dominates the energy budget of the universe we first note that the SM is always relativistic and thus s SM = (ρ SM + p SM )/T SM ≈ 4/3 ρ SM T SM . The temperature T e D (and corresponding SM temperature T e SM ) for which the two energy densities agree is derived by imposing ρ SM (T e SM ) = ρ D (T e D ) with the definition ζ = s SM /s D we obtain the equation 4/3ζ −1 = T e SM /T e D . This can be expressed in terms of the dark sector inverse temperature z e π = m π /T e D using eq. (2) or eq. (7): which defines the temperature at which the dark sector starts to dominate the energy budget of the universe. In the case in which the SM dominates the energy budget this is given by the solution to: while if the dark sector dominates: z e π ≈ 0.64 For all z π > z e π the dark sector dominates the energy budget of the universe. This will be crucially important, when possible entropy dilution of frozen-out relics will be considered [30].
The results of this section are summarized in appendix C.

IV. DARK FREEZE-OUT
We are now ready to discuss the dark freeze-out following reference [29]. One can envision the following possibilities: 1. If the SM energy density dominates the universe at freeze-out the DM relic density is enhanced by a factor proportional to T D /T SM 1. Depending on whether the dark sector particle number is conserved or not the actual functional dependence of T D on T SM changes.
2. If the dark sector dominates the universe energy density at freeze-out the corresponding DM relic density is enhanced by (T D /T SM ) 3/2 1. As for the point above the actual functional dependence of T D on T SM depends on whether the dark sector particle number is conserved at freeze-out.
3. If the LCPs dominate the universe energy budget and decay after freeze-out the entropy injection dilutes DM and suppresses it by a factor D ≈ T e SM /T RH , where T e SM is the SM temperature at which LCPs start dominating and T RH ≈ 0.55 g * SM Γ π M pl is the generated reheating temperature.
The dark baryons annihilate in a rearrangement process, as discussed in Ref. [31]. Since a baryon in an SU (N c ) theory is a Q Nc object it carries baryon number one or equivalently fermion number N c . If a baryon and antibaryon collide they can rearrange into a state with N c − 1 quark-anti-quark pairs by meson emission. The remaining multiquark state has baryon and fermion number zero and thus consequently annihilates via sequential meson emission. We will discuss this process in detail in Section II A. The cross section of this rearrangement process is central for the dark baryon freeze-out.
In references [28,29] the sudden freeze-out approximation was adopted in order to determine the final DM relic abundance. Here we will demonstrate that the sudden freeze-out approximation is not accurate because it does not take into account the late time annihilation processes which count here.
We start with the full Boltzmann equation for the DM baryons which reads: Here the space-time interaction density is 2γ ≈ σ ann v n eq2 B and further assumes that the baryon-antibaryon annihilation takes place via a re-arrangement of the underlying degrees of freedom into N c pions [31].
It is convenient to rewrite the left-hand side in terms of a dimensionless quantity as with Y B = n B /s D ≈ n B /s π , z = m B /T D = z π m B /m π = z π /r and the scale factor evolves as a ∝ t m . The Jacobian factor in (13) contains most of the thermodynamic information. We now move to systematically solve the Boltzmann equations for the different cases.

A. Boundary Layer Solution of the Boltzmann Equations
A direct application of the general formula (D1), derived in appendix D, gives the Boltzmann equation for dark Baryons and LCP abundances: Moreover, we parametrised the rearrangement process of baryons into LCPs with the parameter w, with w = N c for pions and w = 2 for glueballs.
These interactions freeze-out respectively at z f , z c , and we will assume z f << z c . This approximation allows us to study the system in terms of decoupled Boltzmann equation: the baryon freeze-out when µ π = 0 and the (relativistic) LCPs do not deviate from the equilibrium distribution. Afterwards, the 3 → 2 process among the LCPs can freeze-out.
We now discuss boundary-layer solutions to the Boltzmann equation following [32]. Given a general differential equation, a boundary-layer type solution is a solution that satisfieṡ Y (z) Y (z) everywhere apart from a finite number of narrow regions in which the opposite is true. These narrow regions are called boundary layers, and the size is typically controlled by a small parameter in the differential equation. A boundary layer solution is found by matching solutions in the regions outside the layers with solutions found inside the layers.

Baryon Freeze-Out
The general template for the Baryon freeze-out equation iṡ where we kept a general z-dependent function. In the case in which the LCP is the most abundant specie and drives the expansion of the universe we have the ordinary result f (z) = Az −2 , with A constant for s-wave process only. We will find a solution with a boundary layer at z = z f to be determined. For 1/λ 1 and z < z f we have The position of the boundary layer can be estimated by finding the value of z for which the LO approximation in 1λ breaks down: On the other side of the Boundary layer, z > z f the inhomogeneous term can be dropped and the equation readṡ These two solutions have to be matched inside the boundary layer, where can change variable where we used the consistent balance κ = 1/(λf (z f )) and dropped higher orders. The solution is easily found and matched with the z > z f patch: as well as with the z < z f patch: Finally, the asymptotic relic density is calculated as:

Number-changing process freeze-out
We now apply the boundary layer method to the freeze-out of particles interacting with number changing operators among themselves. First let us assume that the pions are relativistic. WeẎ The general structure of the equation iṡ Taking in the outer region z < z f (with z 1) an ansatz Y = Y 0 + 1 λ Y 1 we get Y n = 0. This trivially implies that the freeze out cannot happen while the pions are relativistic, consistently with the Γ 32 /H estimate.
Freeze-out happens when pions non-relativistic, the general equation is theṅ In the outer region z < z f (z 1), we proceed as in the previous section with ansatz .
The freeze-out temperature is determined as In the post-freeze-out region the equation reduces to: Inside the boundary layer, choosing κ = 1/(λf (z f )) the solution reads: Matching in the pre-freeze-out region instead, we have leading to the following asymptotic relic abundance: Equipped with the general analytic or semi-analytic solution for the Boltzmann equations, we will study the asymptotic solutions for baryons and pions/glueballs in concrete regimes, assuming different hierarchies between z f and z c .

B. Application to Freeze-out in the Number-changing Era
In the following sections we will discuss the baryon freeze-out and consider stable LCPs in section V A 4.
First, we assume that the dark freeze-out temperature T f (due to dark baryons annihilation processes) is higher than the temperature T c at which the dark number-changing operators freeze-out. Here the number-changing processes keep the chemical potential at zero while heating up the dark plasma. Furthermore, we assume that the plasma degrees of freedom do not decay in this era i.e. Γ π ≈ H(T 0 ) and T f > T c > T 0 .
We use (3) to evaluate the Jacobian factor and the fact that the LCPs are in thermal equilibrium with vanishing chemical potential. These assumptions yield The large constant in front of the derivative ensures that asymptotic methods will give good results and is given by The temperature dependent function is given by with r = m π /m B being the ratio of LCP and dark baryon masses. The explicit form of the equilibrium baryon number density, normalised to dark entropy is: We determine the freeze-out temperature expanding the solution of the Boltzmann equation in powers of 1/λ. The freeze-out is defined as the temperature at which the next to leading order in 1/λ is of the same size as the leading order contribution. Thus defining Y B ≈ we can insert this ansatz in eq. (33) and solving Y B, eq = Y 1 B /λ, we obtain for the case that the SM dominates the energy budget and for the case when the dark sector dominates we have Following the boundary layer method suggested in [32] we find the asymptotic solution to In the case of σ ann. v ≈ const. the integral can be found analytically (assuming z c where E n (x) is the n−th exponential integral. We see at this point again that the post freeze-out regime dominates over the standard scenario where

C. Application to Freeze-out in the Number-conserving Era
Next, we consider the case T c > T f such that freeze-out occurs when number-changing processes are inactive and a chemical potential leads to conservation of dark particle number.
We still assume that the plasma degrees of freedom do not decay in this era i.e. T c > T f > T 0 .
There will still be two possibilities, depending on whether the SM or the dark sector dominates the universe at DM freeze-out. We use (8) to evaluate the Jacobian factor together with the chemical potential µ = m π (1 − T D /T c ), which implies for n π ≈ n eq π e µ/T D . The Boltzmann equation (12) reads now The equilibrium baryon number, normalised to dark entropy, reads The temperature dependent function reads Following the procedure outlined in the previous subsection the freeze-out temperature is in the regime where the SM dominates the energy budget. The case in which the dark sector dominates the energy density gives The boundary layer method leads to the following asymptotic solution for (39) In the case of σ ann. v ≈ const. the integral can be solved analytically (assuming We find that also in this scenario, the late time annihilations play a significant role and corrections to the sudden freezeout approximation are substantial. These solutions are applicable if freeze-out happens on timescales on which the LCPs are stable.

D. Application to Freeze-out with Decaying LCPs
Here we consider the decay of the LCPs and consider two scenarios, with different effects on the asymptotic dark baryon abundance.

LCP Decay before Baryon Decoupling
Here we assume that the lifetime of the universe becomes larger than the LCP lifetime, while dark baryons are still in chemical equilibrium with the LCP plasma. In this case the set of Boltzmann equations governing the system reads where the space time interaction densities are given by 2γ ≈ σv n 2 B eq and γ ann ≈ Γ ann n π eq and 3γ 32 ≈ σ 3→2 v 2 n 3 π eq . Here we allow for annihilation rate of the LCPs into SM particles. The LCP plasma begins to decay away once the Hubble rate drops below the annihilation rate H(T Decay ) < Γ ann , which defines the decay temperature T Decay . The dimensionless number densities are now conveniently defined with respect to the SM entropy density Y = n/s SM during the radiation dominated era.
As discussed in [33], once the annihilation rate is faster than the Hubble rate, the left hand side of (45) can be neglected at leading order and the equation becomes algebraic. Now, the effective equation for the DM number density can be obtained by inserting the solution of (45) into (44). We obtain where the effective interaction density is for N c = 3 with the branching ration defined as BR(z) = σv n 2 B eq 1 + N c Γ ann n 2 B eq n 4 π eq σv σ 3→2 v 2 Γ ann n π eq + N c σv n 2 In fact baryons in large color groups can be constructed by means of the two index antisymmetric representation [34], in which case this formulae apply to arbitrary N c .

This final Boltzmann equation can be integrated and yields
This effective description clearly shows how the decay of the LCP generates a thermal link between the SM and the dark sector. The crucial factor is the interaction proportional to Γ ann (Y π − Y π eq ) in the Boltzmann equation, which establishes the thermal link, provided that Γ ann H(T Decay ). Below the temperature T Decay , the LCPs can not be regarded as stable any longer, as the Hubble time starts exceeding their life-times.

LCP Decay after Baryon Decoupling
If the dark baryons have already decoupled from the dark plasma at the time of LCP decay, the only effect which can be relevant is entropy injection and DM dilution [35,36]. This is and can rewrite: The reheating temperature is determined by the condition that the LCP lifetime equals the Hubble time τ H = Γ −1 ann . Which for a matter dominated universe is We make use of the fact that all of the dark sector energy density is converted to SM energy density in the reheating process and rewrite the Friedmann equation to find the reheating energy density The reheating temperature is therefore given by The dark baryon relic abundance is then given by Y B = Y B (∞)/∆. In the case that the LCP energy is a subdominant fraction of the total energy density in the universe, no dilution takes place.
The questions regarding whether the dark baryon freeze-out takes place during the numberchanging or preserving era, as well as what is the actual LCP lifetime can be answered in a concrete model. Note that the number-changing interactions are generically present among LCPs (glueballs or pions) at low energies. We will study the most economical and appealing models of this type in the following sections.

V. WEAK AND STRONG COUPLED REGIMES OF DARK BOUND STATES
Here we employ the formalism developed in the previous sections to investigate the asymptotic abundance of stable relics in different regimes of non-abelian gauge theories. First, we focus on the strong coupling regime, where the gauge coupling at typical momenta inside the bound states is non-perturbative.

A. The Strongly Coupled Bound-state Limit
Here we assume the vector-like fermion masses to be small compared to the dark confining scale. In this limit, the dark bound states are relativistic and interact non-perturbatively among each other. This situation mirrors ordinary QCD.

DM Self Annihilation
Because of the QCD resemblance we can use it as an analog computer to determine DM self annihilation by properly rescaling the measured QCD proton-proton annihilation cross section. The annihilation is exothermic since σv rel = constant and the s-wave contribution dominates. Thus the dark baryon annihilation cross section is [37] σv rel. = 4π

Dark Pion Interactions
In the current setup the lightest dark particles are goldstone bosons described by the chiral Lagrangian, see for details appendix A. Because the coset space of the spontaneously broken global symmetry has a non-trivial fifth homotopy group we can write, for SU (N f ) R ⊗ SU (N f ) R broken to SU (N f ) V the following term induced by the WZW [38,39] action: where λ a are the broken generators. The resulting cross section for the 3 → 2 prodcesses is For a breaking pattern of the form SU (N f ) ⊗ SU (N f ) → SU (N f ) we have (see [25]) N π = N 2 f − 1 and t 2 = 4/3 (N 2 f − 1)(N 2 f − 4). Within this framework of chiral perturbation theory the dark baryon abundance and the effects of the dark plasma can be computed.
In Fig. 2   Depending on the parameter space the resulting dark baryon density will come from either number preserving or number-changing type of interactions outlined above. Therefore, the dynamics of the system is either governed by the Boltzmann equations from Section IV B or Section IV C. The pion lifetime is varied as a model independent free parameter and for comparison the two regimes are contrasted. In the case that the pions are stable on cosmic time-scales during the confinement and freeze-out process (τ π > τ DC H ) the relic abundance depends on the initial entropy ratio ζ. If the pion decay is significant during the freeze-out (τ π < τ DC H ) the established connection to the SM plasma erases the ζ dependence. The Boltzmann equation used for the relic abundance is in this case described in Section IV D.

Pion Decay Effects
To discuss the pion lifetime effects more concretely, we specify a model, by choosing quantum numbers for the constituent quarks. The simplest viable model, as discussed in [40] contains The dark Goldstones enjoy a G-parity symmetry (charge-conjugation followed by a weak isospin flip) that would stabilize π 0 . However, the symmetry can be broken by a dimension five operator induced by physics at the scale Λ G / . This physics may induce, for example, the operator −κ (α EM /Λ G / ) π 0 F µνF µν enabling the decay with The coefficient κ is O(1) and Λ D < Λ G / < M P l , which defines the minimal and maximal lifetime of the pions respectively.
It is obvious, that in order for a pion to be stable on cosmological scales its mass has to be below the GeV scale. Thus, we have a severely limited parameter space for a long-lived (a) No entropy dilution due to pion decay. Short lived or stable pions are assumed.
Entropy dilution in the case that pions are unstable and their lifetime is τ π ≈ 1 s. pion as DM candidate. In the generic case, as mentioned before, the baryon will be the DM candidate and thus the pions have to decay before Big Bang Nucleosynthesis, in order to leave the abundance of light elements unchanged i.e. τ π < 1s. In the case that the dark pions have such long lifetimes, the entropy dilution, as discussed in Section IV D 2, takes place. The resulting relic abundance is shown in Fig. 3b. Note that in this case, the dark baryon mass can be significantly higher than expected from unitarity considerations [41,42].

Stable Pion Abundance
As discussed in the introduction, it is possible to construct models with global symmetries, in which the pions are long-lived DM candidates. This implies that the pion cannot be heavier than a few GeV. In this section, we predict the relic abundance of the long-lived pions as a function of the model parameters.
While for the dark baryons, the freeze-out and production depend to a lesser degree on the couplings of the constituent quarks to the SM particles, the situation is different in the pion case. This is due to the fact that viable models feature light pions and therefore, depending on the concrete realization of the model, the pion-SM elastic scattering might be strong enough to maintain kinetic equilibrium between the sectors, leading to the equality of T D = T SM . Whether temperature equality is achieved changes dramatically the resulting freeze-out phenomenology.

Kinetic Equilibrium
In Section III, it is shown that a dark pion sector, which is not in kinetic equilibrium with the SM bath cools slower than normal non-relativistic matter. Thus, structure formation excludes a secluded, self-heating pion sector as the sole component of the DM fluid in the universe [43].
Therefore, we first consider a scenario in which the temperature equality between the dark pion and the SM sector is maintained throughout the freeze-out process, as in [44]. Here the general Boltzmann equation Eq. (26) greatly simplifies and the temperature dependent function is just f (z) = z −n π . The WZW interaction cross section includes the z −2 π scaling of the pion scattering cross section for n = 7.
The asymptotic number density with n = 7 is thus given by where λ = 2g 2 π π 5/2 135 √ 5g is the large parameter enabling the boundary layer method to work. Here we have defined σ 0 32 = σ 32 v 2 rel z 2 . In the case of 3 → 2 freeze-out and equal temperatures in the SM and dark sectors, the critical temperature at which the number-changing interactions decouple can be found an-alytically to be where W (x) is the ProductLog function.
Since the relic density is given by Ω π ≈ s 0 SM m π Y ∞ π /ρ c ≈ 2.6×10 3 (N 2 f − 1)m π σ 0 32 GeV 3/2 −1 this scenario alone does not point to a fixed target cross section, as in the case of the well known thermal WIMP [45]. Fixing the relic density to be the maximally allowed relic abundance of DM we constrain the constraint (N 2 f − 1) 2 m 2 π σ 0 32 > 10 8 GeV −3 . Given the scaling of σ 0 32 ≈ 2.3 × 10 6 N 2 c r 10 m −5 π N −1 f , we find that the relic density is Ω π ≈ O(1) r −5 N −1 c (m π /(N f GeV)) 3/2 and m π ≤ 0.3 r 10/3 N  There is an other important point which is crucial and very general to take into account here. The pion self scattering cross section is given by [46] Note that this cross section only weakly depends on the number of fermion flavors. Observations of galaxy clusters provide upper bounds on the DM self scattering cross section demanding that σ 2→2 /m π < cm 2 /g ≈ 4.5 10 3 GeV −3 . Given the upper bound, we derived on the pion mass from the relic density constraint, we get a lower bound on this quantity

Kinetically Decoupled Pions
Even though Ref. [43] shows that self-heating secluded pions cannot be the DM of the universe, Ref. [48] discusses the possibility that a sub-dominant DM fraction of self-heating light particles is acceptable and might alleviate the observed tension in the matter power spectrum. To this end, the fraction has to be of order 1% of the total DM relic density. Fig. 4 shows the numeric result for the secluded dark pion relic density as the green shaded region. Since in analogy to the baryon freeze-out the Hubble rate has a complicated dark temperature dependence there is no such simple analytical expression for the final relic density. It is interesting to observe, that in this case the mass of the dark pions is expected to be between 10 and 100 keV, with a mild dependence on the ratio r = m π /Λ D ≈ m π /m B , with the dark baryon mass m B .
As shown above, in order to obtain a symmetric relic abundance of dark baryons in accordance with the observed DM relic density, those baryons have to have mass of order 100 TeV, which implies an extremely small value of r < 10 −10 , which makes the pion self scattering cross section tiny. In the case that dark baryons are produced with an asymmetry and a smaller mass, the ratio r can be much larger. However, since the dark pions are a subdominant DM fraction, the DM self-interaction constraint does not apply to them.

B. The Weakly Coupled Bound-state Limit
In this subsection, we consider the interaction regime, in which M Q Λ D , thus bound states are perturbative objects similar to heavy quarkonia. In this regime, the annihilation takes place in two ways.  The rearrangement cross section is given by σ ann ≈ πR 2 B / E kin /E B and σ ann v rel ≈ π/( √ N c C N α D m 2 Q ), where C N is the quadratic Casimir of the fundamental representation of the considered group [40]. The resulting cross section is thus a temperature dependent quantity, given by The group factor can be computed in a given model, for example for an SU (N c ) group it is κ = (N 4 c − 3N 3 c + 2)/16 N 3 c , assuming the quarks are SM singlets [40]. The value of the rearrangement cross section quoted above can be estimated by numerical simulations of the classical rearrangement process, where the color charge of the dark quarks is taken into account. This is done by augmenting the Maxwell equations, by a classical evolution equation for the color vector Q a (t) with the adjoint index a, which is constructed from the bilinear Q a = c † i λ a ij c j with the color vectors of the fundamental fields c i and the corresponding group generators λ a ij [49]. The parallel transport equation dictates the evolution of this color vector v · DQ a = 0 [50,51], where v µ is the four velocity of the quark. This equation can be written and approximated in the non-relativistic limit aṡ where f abc is the structure coefficient of the gauge group. At leading order in the gauge coupling the equation for the potential in the Lorentz gauge is just A µ a = g D J µ a , where the sources are J µ a = i dτ Q a (τ )v µ i δ(x − x i (τ )). We can use the Lienart-Wiechert ansatz and get at leading non-relativistic order With the Lorentz force equation m Q¨ x i = g D Q a i ∂A a 0 for the positions of the color charges we obtain a closed system of differential equations. As an example we show a simulated classical meson, anti-meson rearrangement event in Fig. 5.

Glueball Interactions
The glueball spectrum is well known from lattice studies [52] and contains the 0 ++ = S state as lightest particle. The interactions among the glueballs can be described by a low energy effective Lagrangian, obtained in the large N c expansion [53][54][55] with a i of order unity. From this Lagrangian the cross section for the 3 → 2 processes can be computed

Glueball Decay
Glueballs are not protected from decay by any accidental symmetry. Therefore, generically they will have a limited lifetime. If the theory contains quarks, which carry electroweak quantum numbers, effective operators will be induced. If no Yukawa couplings with the standard model fermions are present, the leading operator will be of dimension eight [40] with the coefficient induced by the quark loop For the simple modle, we have introduced in the previous section, the only quarks in the theory are in the Q = (3, 1, 3, 0) of SU (3) D × SU (3) c × SU (2) L × U (1) Y and thus T D = 1/2, T 2 = 2, d 2 = 3 and Y = 0. The decay with of the glueballs is therefore The constant f G can be extracted from lattice data as in [56] f G := 0| Tr G µν G µν |0 ++ ≈ 3m 3 GB /4πα D . In the case of heavy glueballs with m GB m W , the kinematical bracket separating the γγ, Zγ, ZZ and W W channels respectively, simplifies to 3 α 2 2 .

Phase Transitions
In this section, we will show how a first order phase transition can lead to a dilution of frozen out relic particles [57]. Numerical simulations of pure gluon dynamics indicate the presence of a first order confinement phase transition for N c larger than two. We can use this analogy for the dark scenario with heavy vector-like quarks. Here we are interested in the discontinuity of entropy density, which would affect the abundance of frozen out relics.
The entropy of the Yang-mills sector just above the critical confinement temperature is given by s D = 4π 2 /90 g g T 3 D − ∂B 0 /∂T D , where g g is the number of gluons. The non-perturbative contribution is encoded in the factor B 0 . This factor can be modeled with the MIT bag model [58] and ascribes constant vacuum energy density to the interior of hadrons. However, just above the phase transition, the plasma can be considered as so dense, that B 0 is a universal global vacuum energy and thus does not affect the entropy function. The entropy density before the phase transition is thus s + ≈ s g ≈ 4π 2 /90 g g T 3 D . After the confinement, the lightest degree of freedom is the 0 ++ glueball with a mass of m GB ≈ 7Λ D . Therefore, at T D < Λ D the glueballs are non-relativistic and have an exponentially suppressed contribution to the entropy density If the dark sector dominates the energy budget of the universe at the phase transition i.e. if Λ D < T E a dilution effect will take place. Since we argue that the non-perturbative entropy generating contributions can be neglected the change in entropy densities is compensated by an expansion of the universe i.e. s + a 3 i = s − a 3 f . The increase in the volume leads to a dilution of any relic already frozen out at the phase transition with the dilution factor given by The confinement temperature can be approximated by T conf. ≈ For the case of N c = 3 the phase transition leads to a dilution factor D ≈ 400. A further interesting aspect of this scenario is that during this first order phase transition dark baryons will be accumulated up in clusters at the collapse points of the false vacuum.
In case of an asymmetry in the dark sector those DM clusters will persist until the present day.

Dark Baryon Abundance
Equipped with our formalism the relic density of weakly coupled dark baryons can now be computed. We begin by exploring the parameter space keeping the SM quantum numbers general and varying the glueball lifetime as a free parameter and then focusing on the predictions of a concrete simple model. Fig. 6a shows the relic abundance of dark baryons, in the case that the dark glueballs are stable at the baryon freeze-out. The final abundance depends on the initial entropy ratio ζ, which is bounded from above by the number of SM degrees of freedom, assuming that the dark and SM sectors were in thermal contact at some higher temperature. Fig. 6b shows the relic abundance under the assumption that dark glueballs are not stable at the time of the dark baryon freeze-out and establish a partial thermal link between the systems. Furthermore, we show that given an entropy ratio ζ > 30 at the first order phase transition in the case of low mass dark glueballs additional dilution takes place and leads to larger dark baryon masses.
No effect of the phase transition, since long lived or stable glueballs are assumed.
Effective dark baryon dilution at the phase transition, which is first order in this case.  and strong entropy dilution takes place. This is indicated by the purple contour which goes up to dark baryon masses, which are much larger than the expected unitarity bound [41].
The maximal entropy diluted relic density is shown by the red contour. It takes place if the glueball life-time is assumed to be τ GB ≈ 1 s, so as large as it can get, given the BBN constraint. For comparison, the dark baryon relic abundance result with short-lived glueballs is shown in green. Note that the provided upper bounds on the DM mass from the requirement Ωh 2 < 0.12 also apply if a DM asymmetry is present. Since in that case, the symmetric component has to annihilate away to at least that value in order not to overcolse Figure  the universe.

Stable Glueball Abundance
When dark glueballs are long lived on cosmic time scales, they can be DM candidates as well. The effective Lagrangian for the glueball interactions is not based on approximate global symmetry but it is a result of a 1/N c expansion and works best at large number of colors [53]. The glueball self-scattering cross section can be computed from Eq. (67) and is directly related to the 3 → 2 reaction cross section The associated relic abundance can be computed in two distinct regimes.

Dark Glueballs in Kinetic Equilibrium
The same logic concerning kinetic equilibrium used earlier for dark pions applies also to dark glueballs, thus the stable glueball is only a good DM candidate if kinetic equilibrium is maintained [59].
The structure of the Boltzmann equation is given by Eq. (26) and thus identical to the dark pion freeze-out case. The only difference is that the 3 → 2 annihilation cross section for glueballs is temperature independent. Therefore, the function in Eq. (26) is f (z) = z −5 and hence the asymptotic number density is given by (n = 5): Given the strict relations between the glueball mass and the number-changing 3 → 2 cross section Eq. (68), we have a concrete prediction Ω GB ≈ 8.1 (m GB /GeV) 3/2 N 3 c . This implies an upper bound on the glueball mass m GB < 0.1/N 2 c GeV translating into σ 2→2 > 5 × 10 4 GeV −2 or equivalently σ 2→2 /m GB > 5 ×10 5 N 2 c GeV −3 . This self scattering cross section violates the observational bound of σ 2→2 /m GB < 4.5 × 10 3 GeV −3 by more than two orders of magnitude. Therefore, stable dark glueballs in kinetic equilibrium with the SM plasma can not be the dominant DM component.

Decoupled Dark Glueballs and Multicomponent Dark Matter
The other DM avenue for dark glueballs is that they are not in kinetic equilibrium with the SM plasma and have their own temperature evolution, which leads to a prolonged phase with non-vanishing pressure. The only phenomenologically viable scenario in this case is a multicomponent dark sector [60,61]. The glueballs comprise in this case only a subdominant DM fraction with a much larger free-streaming length than the dominating DM candidates.
This phenomenon could resolve some large scale structure problems if the glueballs make up for about 1% of the total DM relic abundance [48].
The explicit numerical solution of the Boltzmann equations leads to a dark glueball mass between 100 keV and 10 MeV, depending on the initial entropy ratio ζ. A few comments are in order This self-heating glueball component can help to resolve observational tension in the matter power spectrum [48]. The model considered has N c = 3. This scenario is special from an experimental point of view since glueballs are produced during annihilation of dark baryons. Since the glueballs are long-lived and feebly interacting with SM matter those decay channels are experimentally invisible. Therefore, the detectability in this scenario depends strongly on the specific SM quantum numbers of the dark quarks that will lead to subdominant annihilations of dark baryons into SM final states.
Another possibility is that the dark glueballs, despite being long-lived, slowly decay leading to detectable signals in indirect detection experiments. We will now come to the discussion of experimental signatures of composite DM models.

VI. EXPERIMENTAL CONSTRAINTS
The composite DM models we discussed here can be very broadly divided in two classes: • Type I: Theories with composite DM particles which carry residual SM interactions. This is the case, for example, of one weakly charged dark quark in which the lightest baryon is by construction a multiplet of the weak force [37] but still neutral under QCD. Another possibility are models featuring dark quarks charged under QCD, that lead to residual chromo-electric interactions [63,64]. Another time-honored example is a dark technibaryon [65] which has still residual SM interactions due to its SM charged constituents.
• Satellite [68] observations of the CMB [69] play an important role to discover these models [70]. In the near future the upcoming CTA experiment will probe the target cross section regions very efficiently already with the data of a single dwarf galaxy [71]. Furthermore, the details of the spectral shape are less relevant for an order of magnitude limit, because the spectrum is sufficiently softened by cascade decays [72].
In Fig. 9, we present the limits for the minimal SM charge assignment of the dark quarks.
Here, the dark quarks transform according to the adjoint representation of SU (2) L . As discussed above this yields the lightest dark baryon (the DM candidate) to transform as the adjoint representation of SU (2) L . The dark glueballs or dark pions generically decay to photons ore heavy SU (2) L gauge bosons, if kinematically allowed. Note that the limits on the dark baryon annihilation cross section will be in the same ballpark, even if the DM candidate baryon is a SM singlet given that the LCPs decay into photons.
The above scenario does not apply when considering long-lived glueballs yielding invisible dark baryon annihilation processes. Depending on the lifetime one could hope to experimentally test this scenario via the dark glueball decays into SM particles. We postpone a detailed study of this scenario to a later endeavour.

VII. CONCLUSIONS
We investigated the thermodynamic history of confining dark sectors that can feature global accidental symmetries, leading to long-lived particles such as the proton for the Standard We adopted our formalism to investigate two regimes of confining dark sectors, the strongly coupled bound-state regime, and the weakly coupled bound-state regime. We find several viable Dark Matter scenarios: • Stable dark baryons, strongly coupled with masses around the 100 TeV scale.
• Stable dark baryons, weakly coupled with masses between a few and 100 TeV.
• Stable dark pions around the GeV scale and considerable self scattering cross sections.
The effective theory is written in term of a unitary matrix U (x) = exp iΠ(x)/f π where Π = a T a π a ( Tr[T a T b ] = 1 2 δ ab ). Considering the constraint U † U = 1, the lowest order Lagrangian is: where M parametrizes the small breaking of chiral symmetry due to quark masses. For this description to be valid an hierarchy of scales is assumed The Lagrangian A3 has been truncated to the lowest order in momentum expansion, but all higher order terms are present. In particular, the lowest order one providing numberchanging processes is a 5-point interaction. For our specific symmetry breaking pattern this term is present and it's the so-called Wess-Zumino-Witten term: In our treatment this is the main source of number-changing interaction in the piondominated confined dark sector.
We denoted by n(t) the particle number density and by U µ = (1, 0, 0, 0) the fluid fourvelocity in the comoving frame. H(t) =ȧ/a is the Hubble parameter. In general a relation between pressure P and energy density ρ of the fluid is required, for example: Pluggin the FRW ansatz in the Einstein equations one finds: We will consider exclusively the case in which k = 0. These equations uniquely determine the universe scale factor a(t) with his own matter content.

General Multi-fluid Thermodynamics
In this section we provide a compendium of formulas for equilibrium multi-fuid systems.
Considering first a single component fluid of particles, this is in kinetic equilibrium if kinetic energy exchange through particle collision is efficient. In this case, the phase-space distribution is a Bose-Einstein (BE) or Fermi-Dirac (FD): f (p) = 1 e β(E−µ) ± 1 ∼ e −β(E−µ) when β(E − µ) 1 .
Relativistic particles have the dispersion relation E = p 2 + m 2 and β = 1/T . The pressure, energy and particle density are obtained by with g being the number of discrete degrees of freedom (spin, flavour, ecc..). The entropy density of the system is obtained from dS = β(dU + P dV − µdN ) =⇒ s = β(ρ + P − µn) .
Useful expressions are obtained in the relativistic βm ∼ βµ 1 limit : For simplicity we will include the BE/FD factors in the definition of g.
In the non-relativistic regime βm ∼ βµ >> 1 we have instead: Although we include the chemical potential in our treatment due to the presence of numberchanging interaction, we will always assume to be far from the degenerate phase. In general as a consequence of continuity equations we have µ = µ(T ).
Suppose now that the ensemble contains many species of particles, and inter-specie interactions are absent. Then every specie thermalize separately, each one with its own temperature and chemical potential (T i , µ i ).
In our treatment it's convenient to work with dimension-less quantities with respect to a given scale M . In terms of z = M/T, z µ i = µ i /T . The entropy density in relativistic and non-relativistic regimes are, respectively: Where we introduced the mass ratios r i = m i /M . If inter-specie interactions are switched on, these can provide extra channels for thermalization (elastic processes) or specie-changing interactions (inelastic processes). The balancing of the latter reactions with respect to the time-reversed processes is governed by the chemical potentials. When the sum of chemical potentials of reacting particles are equal on both sides of the reaction the system reaches chemical equilibrium, and the fraction of each specie is then determined and not evolving in time.
A multi-fluid system that is both in kinetic and chemical equilibrium is said to be in thermal equilibrium. In this case, every specie follows the equilibrium distribution with a common When all the species are in thermal equilibrium the total entropy density is In square bracket one can read the effective entropy degrees of freedom g * S (z). For all practical purposes, the entropy of the non-relativistic sector can be neglected compared to the relativistic one, and away from mass thresholds g * S (z) ∼ r i z<<1 g i ∼ const.
Studying departure from thermal equilibrium it is useful to introduce the particle number per comoving volume Y i = n i /s, which evaluates to After this point, the temperature of the two sectors evolves differently depending on the degrees of freedom of each sector. We will describe the different phases the system evolves into as the universe expands.
This is self-consistent as long as Λ D m H as expected from a new physics scale. The total energy density is found analogously as in the previous case and the scale factor evolves as The Γ(z) factor contains the microscopic information about the interaction, and it is connected to the matrix element in the zero-temperature QFT in the followgin way: where f eq i ∼ e −β(E i −µ i ) are the particle occupation number for a thermal distribution which is sufficiently diluted. The sum/integral is performed on the phase space of both initial and final states.

a. Calculation of Jacobian factors
In this section we compute the Jacobian factor that appears in the Boltzmann equation expressed in term of dark temperature: In the regime in which the two sector temperature are linearly related and the dark sector is dominated by relativistic pions, we have so that J(z D ) = H(z D )z D . When the pion are non-relativistic and number-changing processes are active, we have instead Finally, once these processes are decoupled we get