Glueball dark matter, precisely

We delve deeper into the potential composition of dark matter as stable scalar glueballs from a confining dark $SU(N)$ gauge theory, focusing on $N=\{3,4,5\}$. To predict the relic abundance of glueballs for the various gauge groups and scenarios of thermalization of the dark gluon gas, we employ a thermal effective theory that accounts for the strong-coupling dynamics in agreement with lattice simulations. We compare our methodology with previous works and discuss the possible sources of discrepancy. The results are encouraging and show that glueballs can account for the totality of dark matter in many unconstrained scenarios with a phase transition scale $20$ MeV$\lesssim\Lambda\lesssim10^{10}$ GeV, thus opening the possibility of exciting future studies.


I. INTRODUCTION
Dark Yang-Mills sectors, which undergo confinement to form stable composite states known as glueballs, may potentially explain the nature of Dark Matter (DM) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] (see also Ref. [16], for a more general class of DM).This type of self-interacting DM has been shown to provide a consistent explanation for the structure of the Universe at small scales and may also help address issues such as the missing satellite problem [17] and the cusp-core problem in the DM distribution at galactic scales [18,19].
The first-order confining phase transition at a critical temperature T c , present in these models [20][21][22][23][24], makes it highly non-trivial to follow the formation of glueball DM.This challenging and interdisciplinary study requires a detailed knowledge of thermal field theory in a nonperturbative domain, and a productive exchange of results with lattice Quantum Chromodynamics (QCD).
In light of previous studies [1], the calculation of glueball relic density is extended to a generic SU (N ) gauge group for N = {3, 4, 5} by using a low-energy effective model for the gluon-glueball dynamics [25].Additionally, different cosmological scenarios are considered to determine the temperature of the dark gluon gas, an important parameter in determining DM glueball abundance.More specifically, the possibility that dark gluons are thermally produced in the primordial plasma or result from a heavy particle (perhaps the inflaton) decay is explored.In a very model-independent way, we determine the glueball models capable of explaining the existence of DM.
In Sec.II we present the effective Lagrangian used to describe the thermal evolution of the dark gluon-glueball plasma.In Sec.III we discuss how this picture can be merged in a cosmological setting to predict the relic glueball DM abundance.In Sec.IV we analyze the differences between our approach and the ones usually employed in literature, to stress why this approach is the most accurate to date.Sec.V discusses, in a very modelindependent fashion, the possible cosmological scenarios in which the dark sector temperature is determined.Finally, in Sec.VI we summarize our findings and conclude.
At finite temperature T , the Z N center of SU (N ) is a relevant global symmetry [47], making it possible to construct various gauge invariant operators charged under Z N .The Polyakov loop, which is charged with respect to the center Z N of the SU (N ) gauge group (it transforms as ℓ → zℓ with z ∈ Z N ), is an example where P denotes path ordering, A 0 is the time component of the vector potential associated with this gauge group, g is the SU (N ) coupling constant and (τ, x) are Euclidean spacetime coordinates.
The Polyakov loop serves as an order parameter for the confinement phase transition in Yang-Mills theory, which occurs at the energy scale Λ and is commonly used for this purpose [47].Below the critical temperature, T c , the expectation value of the Polyakov loop operator is zero, while it is non-zero at higher temperatures.The Polyakov Loop Model (PLM) is a mean field approach that models the phase transition in terms of Polyakov loops [39].This simplified model captures the essential characteristics of the confinement phase transition in SU (N ) theories with N ≥ 2 and has been applied to the study of heavy-ion collisions at the Relativistic Heavy Ion Collider [28,45].
The dark gluon-glueball dynamics can be effectively described by considering the dimension-4 glueball field H ∝ tr(G µν G µν ), with G µν QCD field strength tensor, and the Polyakov loop ℓ in an effective potential given by [25]: Here, the first term represents the zero-temperature glueball potential, determined by the trace anomaly constraint [48,49].The real polynomials V [ℓ] and P [ℓ] are invariant under Z N , with coefficients that are fitted to lattice data.The term V T [H] accounts for thermal corrections, which may involve non-analytic terms in H [46]. Remarkably, the potential in Eq. ( 2) reduces to the glueball dynamics and PLM model in the low and high temperature limits, respectively.Furthermore, the coupling between H and ℓ is the most general interaction term that can be constructed without violating the zerotemperature trace anomaly (see Eq.( 21) in Ref. [49]).This approach neglects heavier glueballs and pseudoscalar glueballs that are described by gauge-invariant operators with different charges under Z N .Despite its simplicity, this model captures the essential features of the Yang-Mills phase transition.
In the deconfined phase, T ≫ T c , the PLM term T 4 V[ℓ] dominates, i. e. dark gluons are the dominant component.The precise relation between the confinement scale Λ and the critical temperature of the phase transition T c depends mildly on the gauge group and is determined by lattice simulations.In this paper, we consider T c = (1.59+ 1.22/N 2 )Λ for N = {3, 4, 5} [15,50].
The Lagrangian that describes the glueball and Polyakov loop degrees of freedom is given by [25,51,52]: Here, c = (Λ/m gb ) 2 /2 √ e is a constant that depends on the glueball mass m gb , which we assume to be m gb = 6Λ [53].The Polyakov loop is a non-dynamical, homogeneous in space, order parameter that describes the average dynamics of the phase transition.It neglects bubble nucleation, which might have a significant impact on the formation of glueballs, as observed in presence of matter fields [54,55], but we leave this discussion for a future work.The kinetic term for the glueball field H is non-standard, due to its dimensionality.To canonically normalize this field, we redefine the glueball field as ϕ, where H = 2 −8 c −2 ϕ 4 , which evolves based on the following effective Lagrangian where c 1 is a free parameter relevant to the determination of the glueball relic abundance.Note that we keep only the lowest order in P[ℓ] satisfying the symmetries.The Polyakov loop potential V[ℓ] for a generic SU (N ) gauge group, with N = {3, 4, 5}, is determined from symmetry arguments to fit lattice thermodynamic quantities [22] V where and the parameters of this potential are shown in Tab.I, with the corresponding potentials shown in the left panel of Fig. 1.Here, we notice that the minima of ℓ do not differ strongly as a function of the chosen gauge group, albeit the overall potential being quite sensitive to this change.Since the Polyakov loop is a non-dynamical degree of freedom, its temperature evolution is determined by the location of the minimum in the effective potential.Being the order parameter of the phase transition, ℓ approaches 1 at high temperatures and vanishes for temperatures below the critical one.It is possible to numerically find the temperature evolution of ℓ by minimizing the potential in Eq. ( 4) with respect to this variable.The solution ℓ = 0 denotes the confined phase and it is a global minimum only for temperatures below the critical temperature.In the deconfined phase, the solution ℓ = 0 becomes metastable and a new solution ℓ = ℓ + becomes the new global minimum.The temperature evolution of the minimum of the Polyakov loop in these potentials is shown in the right panel of Fig. 1.The three gauge groups shown lead to a similar behavior for ℓ during the phase transition, with a slightly different critical temperature.Once the minimum of ℓ is determined, the Polyakov loop is "integrated out" using its equation of motion ℓ = ℓ(ϕ, T ), giving rise to a potential for the glueball field in the form V [ϕ, T ] = V [ϕ, ℓ(ϕ, T )].Moreover, we set the zero-point energy of the glueball field to zero in order to properly describe glueballs as matter.Fig. 2, in the left panel, shows the behavior of the glueball potential as function of the different gauge group.The deconfined phase (red lines) is the only one sensitive to the choice of gauge group, and the dependence is very mild.This is directly correlated with the observation that the minimum of ℓ does not change sensibly for different gauge groups with slightly larger ℓ when the color increases signaling the stronger jumping and thus a stronger first order phase transition.Connected with this potential, it is possible to calculate how the renormalized glueball mass evolves with temperature through the phase transition.This is defined as where ϕ min represents the minimum of the glueball field as function of the temperature.This quantity is shown in the right panel of Fig. 2, where we observe a minor difference between the different gauge groups in the deconfined phase.We note that, after confinement, the mass is fixed to be m gb = 6Λ by construction.The effect of the thermal potential V T will show up in a temperature dependence of the glueball mass in the confined phase.We estimated this effect by using the following potential [44] V We find consistent results with Ref. [1], since we verified that the impact of this term on the glueball potential is negligible.Namely, the glueball mass is affected less than 10%, with comparably small consequences on the considered phenomenology 1 .By comparing the temperature evolution of the glueball field to lattice simulations [56], it is possible to impose limitations on the value of the glueball-Polyakov loop coupling c 1 in Eq. ( 4), in the case of SU (3).We found this value to be c 1 = 1.225±0.19at 95% Confidence Level (CL) [1].Since less or no information from the lattice is available for other gauge groups, this will increase the uncertainties on the prediction of the relic abundance.To appropriately account for this, we increase the error associated to SU (3) by a factor of ≈ 4, such that c 1 = 1.225 ± 0.8 at 95% CL for SU (N ), N = {4, 5}, and we fix c 1 = 1.225 to generate the figures throughout the paper.This parameter significantly affects the location of the minimum of the Polyakov loop.Thus, it will play an important role in determining the initial conditions for the cosmological evolution of the glueball field and on the resulting DM abundance.

III. CALCULATION OF THE GLUEBALL RELIC DENSITY
As extensively discussed in Ref. [1], the glueball field evolves in a Friedmann-Lemaître-Robertson-Walker met-1 With the finite temperature contributions, the minimum of the thermal potential is shifted towards 0 with few percents.The larger effect is on changing the slope around the minimum.This small change cannot be calculated on the relic density because Eq.( 17) is already an oscillating function and all the uncertainties related to this average process are much bigger than the shift given by the inclusion of V T .ric as where the dot represents the derivative with respect to the cosmic time t and H = 1/2t is the Hubble parameter during a radiation-dominated era.The energy density of the glueball field is given by and this quantity is used to compute the glueball DM relic density.During a radiation-dominated era, where we expect the confinement to happen, there is a relation between the time and photon temperature where m P = 1.22 × 10 19 GeV is the Planck mass and g * ,ρ is the number of energy-related degrees of freedom.Thus, Eq. ( 9) can be rewritten as a function of the temperature as where the temperature of the dark sector T , such that T γ /T = ζ T , governs the moment of the phase transition and ζ T is model dependent, being determined by the interactions with the visible sector.Moreover, the second term can be neglected for a large range of temperatures as g * ,ρ is constant except at a few isolated events of entropy production (the QCD phase transition, for example).We consider it as a free parameter, depending on the moment in which the phase transition happens and, if g * ,ρ is constant, it can be reabsorbed in the definition of the temperature ratio by defining ξ ′ T = ξ T g The glueball evolution is analogous to a damped oscillator in a non-linear potential, and the energy stored in these oscillations around ϕ min ≈ 0.28Λ will determine the relic DM abundance since this energy density scales as ∼ T 3 , as Cold DM (CDM), when the harmonic approximation is valid.Using the following definitions 1/4 * ,ρ , T = Λ T and µ 2 = 4π 3 ξ ′4 T Λ 2 /45m 2 P , Eq. ( 12) can be written as which is solved from an arbitrary temperature T i > T c down to some final temperature T f in the confined phase, and from this temperature on the evolution is simply determined by the cosmological expansion.In Fig. 3 we show the results of Eq. ( 14).In the left panel we show the evolution of the glueball field for different gauge groups.
The amplitude of the oscillations in the confined phase is related to the relic energy density.Therefore, we expect that, for the choice of parameters shown in the figure, SU (5) gives an abundance slightly larger than SU (4) and significantly larger than SU (3).In the right panel of Fig. 3 we focus only on the case of SU (3), to highlight the dependence on the initial conditions.the evolution of the glueball field is shown for three different choices of initial conditions in the hot phase (on the right of the vertical dashed line denoting the critical temperature).After the phase transition the evolution of the three lines is qualitatively similar, suggesting that there is a weak de- pendence on the initial conditions.This was already observed in Ref. [1].Intuitively, in the hot phase, the glueball field is set to some arbitrary initial condition then starts to roll towards the minimum of its potential.This happens when the effective glueball mass m gb ≃ 3.5Λ in the deconfined phase becomes larger than the Hubble parameter, in the opposite case the glueball field evolution is frozen.Starting from the moment in which the Hubble parameter is comparable with the glueball mass, the glueball field efficiently converges to the minimum of its potential.Given the discontinuous nature of the phase transition, the minimum of the glueball potential jumps from the value immediately before the confinement to evolve in a temperature-independent potential.It means that regardless of the evolution of the glueball before the phase transition, only the initial condition set by the properties of the glueball potential V [ϕ, T ] in the deconfined phase is relevant to predict the glueball relic density.Moreover, the velocity of the field can be taken to be equal to zero since any velocity acquired immediately after the confinement is considerably larger than the velocity accumulated in the evolution in the deconfined phase.This feature can be interpreted as a washing out of the initial conditions due to the strong first order phase transition.
We discovered that, although the numerical solution of Eq. ( 14) is exact, the temperature-dependent potential makes its evaluation computationally expensive.Thanks to the weak sensitivity on the initial conditions, a good approximation is given by solving Eq. ( 14) only in the confined phase.The main advantage is using a temperature-independent potential to evolve the glueball field from the critical temperature T c down to a final temperature, T f , taking as initial conditions for the glueball field its minimum value φmin just before the phase transition (at a temperature Tc + ϵ with ϵ > 0) and a vanishing first derivative where the potential and the energy density are The energy density obtained from this calculation is shown in Fig. 4 for three different gauge groups.It is clear that, as the temperature drops, glueballs behave like CDM and their relic density is redshifted as ∼ a −3 , where a ∼ 1/T is the scale factor.Note that for T /T c ≳ 0.2 the glueball energy density does not redshift as CDM because in this intermediate regime non-linearities in the potential lead to deviation from the perfect fluid behavior for the glueball field, i.e. the glueball field does not behave like a pressureless fluid, p ̸ = 0. Indeed, only for T /T c ≲ 0.2, when the glueball field oscillate around its minimum under the influence of an effectively quadratic potential, the perfect fluid approximation is valid.Thus, the quantity shon in Fig. 4 approaches asymptotically a constant value.One can picture the intermediate phase 0.2 ≲ T /T c ≲ 1 as a continuation of the phase transition, when the non-relativistic glueball population is still being established.This result confirms the expectation from Fig. 3 that the relic density with SU (5) is larger than the one obtained with the two other gauge groups.
Since the energy density is an oscillating quantity, we evaluate an average during the last oscillations before Tf and this quantity saturates to a value independent on the phase transition scale.Then the relic density today is calculated by diluting the energy density in Eq. ( 17) with a factor (T γ,0 /ζ T T f ) 3 , to take into account the Universe expansion as where the critical density is ρ c /h2 = 1.05 × 10 4 eV cm −3 , h = 0.674 and the temperature of the photon bath today is T γ,0 = 0.235 meV [57].We defined Λ 0 as the phase transition scale that makes glueballs become the totality of DM.Naively, after combining T 3 γ,0 / ρ c /h 2 , we would roughly have Λ 0 ∼ ⟨ ρ T 3 ⟩ −1 f eV.Note that here we introduce the ζ T parameter, related to the glueball and photon temperatures, and not ζ ′ T which also includes the number of degrees of freedom in the Universe.This implies that the dependence on this parameter, equivalently µ, is weak in the limit µ ≪ 1, realized in the relevant case when the phase transition scale is much lower than the Planck scale 2 .This approximation is proven to be excellent, predicting the relic density with less than ∼ 10% uncertainty compared to the exact result.The good agreement of the two results is a proof that the detailed behavior of the glueball field in the deconfined phase has a minimal impact on the relic density prediction.

IV. RESULTS AND COMPARISON WITH LITERATURE
In Tab.II we summarize our findings for the glueball relic density.For each gauge group SU (N ) considered, labeled by N = {3, 4, 5} (first column), we recall the range of variability for the term c 1 at 95% CL (second column).Then we show the results of Eq. ( 17) (third column) and the corresponding Λ 0 (fourth column), as defined in Eq. ( 18), for a calculation running down to Tf = 0.1.
Note that the value of Λ 0 for N = 3 is 20% larger compared to the one presented in Ref. [1] because of the different critical temperatures considered.The values of the relic density found in Tab.II differ strongly from the estimates in literature, since they report a relic density one or two orders of magnitude higher [2,15,58].
The reason for this difference is due to a combination of several effects: inclusion of the higher-order interactions leading to n → m transitions; energy budget of the dark gluon field partially used for bubble formation; and different equation of state for the glueball field immediately after the phase transition.
In the following we expand each point individually.First, when solving the evolution equation, we are considering a non-trivial potential for the glueball field that, in the confined phase, can be expanded around the minimum φmin = 4e −1/4 √ c as where δ φ = φ − φmin .Only for δ φ ≪ 1 the perturbative concept of particle is valid.
However, our calculation is always valid, including all the interactions predicted by the glueball potential.It can be perturbatively understood as including all the possible interactions corresponding to different powers of the expansion in Eq. ( 19).In the literature, it is usually considered that glueballs interact only with a ϕ 5 interaction, which makes the 3 → 2 annihilation the only relevant process for DM formation.By contrast, in our case also the lower order terms are included.
We compared the glueball relic density in Tab.II with a calculation including only the ϕ 5 interaction, finding a factor ∼ 1.3 − 1.5 (depending on N ) of increase in this latter case.This shows that ϕ 3 and ϕ 4 interactions are The first column represents the gauge group, the second one is the value of c1 at 95% CL used in the calculation of the third column, which also shows the result of Eq. ( 17) at 95% CL, evaluated for µ = 10 −3 and Tf = 0.1.With these values it is possible to calculate the glueball relic density from Eq. ( 18) and the range of Λ0 is reported in the fourth column.important in the glueball thermalization process.Indeed, the 3 → 2 number changing process can happen both due to a ϕ 5 order vertex and because of a combination of ϕ 3 and ϕ 4 vertices as shown in Fig. 5. Without lower order interactions, just the ϕ 5 term induces a weaker interaction among the glueballs.Consequently, glueballs freeze-out earlier, when their number is higher and resulting in a larger relic density.In other words, the number-changing interactions have less time to reduce the number of glueballs.This picture is confirmed by the observation that including interactions up to the fourth order (ϕ 2 , ϕ 3 and ϕ 4 ) the relic density increases only of less than 1% compared to the exact calculation involving the log-potential.This reasoning also brings us to the conclusion that higher order number-changing processes, like 4 → 2 interactions, do not have a strong impact on the glueball thermalization.In conclusion, when comparing with results in literature, one should be careful in checking which potential is used.From the perspective of a complete model, not including ϕ 3 and ϕ 4 interactions is inconsistent, leading to a larger relic density.Moreover, note that the form of the glueball potential fixes uniquely, once expanded around the minimum, all the self-interaction couplings at any order.
The latter are usually taken to be O(1) in the literature, while our approach reveals that these numerical coefficients are rather different from 1 and any comparison with the literature must account for this important difference.As a final remark, Eq. 10 in Ref. [58] is obtained by setting the numerical factors in Eq. 8-9, involving the Lambert W-function, to 1.This is also causing a slight overestimation of the relic density.
The second important point can be understood on the basis of energy considerations.Starting from the effective Lagrangian, it is straightforward to compute the energy density of the gluon field at T → ∞, which corresponds to [20,59] with g = N 2 − 1 for SU (N ).As the temperature approaches the critical one at the phase transition, this energy reduces to match the fitted lattice data.The physical reason is that the dark gluons dissipate energy in the process of bubble nucleation and a smaller energy budget is actually available for the glueball formation.This effect is purely non-perturbative, and it has been taken into account in our analysis once that lattice fits are considered.Calculating the amount of energy stored in dark gluons at the phase transition we realize that it is 3 − 4 times smaller than Eq. ( 20) depending on N .Thus, the energy budget available to glueballs was overestimated if taken equal to the free dark gluon gas approximation.We remark the physical picture behind the energy exchanges within the dark sector.Most of the energy stored in the dark gluon plasma goes into formation of bound states effectively generating the "mass gap", i.e. converted into the glueball mass m gb .Some of that energy still remains as "heat" in the dark sector, essentially in the form of kinetic energy of glueballs, as well as through a contribution to the potential energy from glueball selfinteractions.Hence, due to the energy conservation and the absence of interactions with the SM sectors, the totality of energy stored in the dark gluon gas remains in the dark sector.Also, during the phase transition, some part of the energy of the transition (latent heat) is released in the form of gravitational waves [23,60], a dissipation effect that is neglected in this paper and postponed to a future work.The aforementioned reduction of the initial energy of the dark gluon gas is partially counterbalanced by a slower dilution of the glueball field compared to a pure cold DM case.Indeed, immediately after the phase transition, the glueball field is rolling fast in a potential which is much larger than its kinetic energy.This results into an equation of state p = wρ with −1 < w < 0, leading to a slower dilution of the glueball field.Only after this transient phase glueballs act like CDM.To summarize, here we list the various effects explaining the discrepancy with the literature: • we include ϕ 3 and ϕ 4 interactions, instead of only a ϕ 5 vertex.This makes the number changing processes more efficient, reducing the glueball relic density (note that in Ref. [15] the potential includes ϕ 3 and ϕ 4 terms, but not ϕ 5 ); • the glueball potential that we consider has 'large' couplings for self-interactions, leading to more than one order of magnitude of suppression in the relic abundance; • part of the energy stored in the dark gluon gas is dissipated in bubble nucleation, reducing the relic abundance of a factor 3 − 4; • the slower dilution of the glueball gas compared to CDM goes in the direction of increasing the DM relic abundance.
The interplay of all these effects is highly non-trivial, motivating the numerical analysis we developed in this work as a reliable method to compute the glueball DM relic density.

V. COSMOLOGICAL HISTORIES
In order to accurately determine the glueball relic density we must specify the temperature of the dark sector with respect to the photon one; i. e. the ζ T parameter.This quantity is vastly unconstrained because of the large number of models predicting different interactions between dark gluons and standard model particles.In a very model-independent fashion, we consider two possibilities for the dark gluon production in the Early Universe: freeze-out and the parent particle decay.
In the first case we assume that, at some point, dark gluons were in thermal equilibrium with the primeval plasma.This is possible due to feeble interactions between the dark and the visible sector.We prefer to keep a model-independent point of view in this work, without specifying the origin of this interaction, but just assuming that it is feeble enough that decoupling happens soon after the end of inflation and the DM is stable on cosmological timescales.A motivated example of feeble interaction is given by fermions that are charged under both the dark and SM gauge groups with a mass much larger than the confining temperature.In this case, all the interactions will be strongly suppressed by the mass scale of this mediator.In the following we do not assume any particular model of interaction and our considerations are completely model-independent.When their interaction rate becomes smaller than the Hubble parameter, dark gluons decouple from the thermal bath.In this case the temperature of dark gluons will trace the photon one up to the decoupling, then entropy production events will cause a cooling of the dark sector following where T d is the decoupling temperature, which determines the number of entropic degrees of freedom at the freeze-out.Without specifying the interaction between dark and visible sector, we know that the lowest temperature of the dark sector is obtained when g * ,s (T d ) = 106.75.This corresponds to a weak interaction between the two sectors that leads to a decoupling at T ≫ 100 GeV, where we assume only standard model particles to exist.This is a strong assumption given our ignorance of physics at very high-energy scales.For instance, in the minimal supersymmetric standard model extension g * ,s (T d ) = 228.75[61].This latter value will be used to fix the coldest dark gluon scenario, which gives ζ −1 T ≃ 0.26, where we have used g * ,s (T γ ) = 3.909.This consideration gives a sense of how the dark sector can be colder than the visible one, proving that in this scenario of freeze-out the two temperatures are never too different.
On the other extreme, a hot dark gluon sector is constrained by the measurements of the number of relativistic species.Indeed, dark gluons cannot contribute to the effective number of relativistic species N eff more than the constraint ∆N eff < 0.35 at the 95% CL [62].For a dark sector that is not in equilibrium with the thermal bath, this constraint translates into [63] ∆N eff = 4 7 11 4 where the number of degrees of freedom is g = N T .This constraint requires that dark gluons confine after the Big Bang Nucleosynthesis, which happens at T γ ≃ 1 MeV.
There are several realizations of the parent particle decay scenario, and we consider the one that we consider the simplest.In this case, a heavy field (that can be associated with the inflaton) decays into standard model particles and also in dark gluons.We take the branching ratio to decay into dark gluons to be f (and 1 − f is the branching ratio into standard model radiation).Compared to the photon energy density, the dark gluon one decreases because of entropy production events in the visible sector and the temperature evolves accordingly Compared to the freeze-out case, depending on the value of f , the dark sector can be extremely cold compared to the visible one (see also the recent Ref. [64] for a discussion on the temperature of confining dark sectors).We assume that the coldest dark sector case corresponds to a confinement happening soon after inflation, when the photon bath has a temperature T γ ≃ 10 16 GeV.This limit corresponds to the breakdown of our calculation, that is valid when the Universe is radiation dominated.This discussion reveals that dark gluons can be extremely cold compared to the thermal photon bath, leading to a very early phase transition, perhaps during inflation.In this case we expect a strong damping of the oscillations of the glueball field, suppressing the relic density.Therefore, we consider that, in order to produce glueball DM, the phase transition cannot happen before or during inflation.This discussion motivates us to use ζ −1 T as a free parameter: if ζ −1 T ∼ O(1), we are modeling a dark sector with interactions strong enough to establish thermal equilibrium soon after inflation; if ζ −1 T ≲ 1, the interactions are so feeble (or even absent) that there is no relation between the temperatures of the visible and dark sectors.Thus, casting this discussion in terms of ζ −1 T is an extremely powerful and general method describing in an unified way several possible models.
We remark that, despite glueballs undergo 3 → 2 number changing interactions (a cannibalistic phase) for a period of their evolution, this phase has to end before matter-radiation equality in the case that glueballs make up the majority of DM.This condition is verified in the allowed region of the parameter space.When cannibalism is finished, glueballs are mildly relativistic, with an average energy roughly 1.5 m gb .After that, they cool down quite rapidly because of the Universe expansion, effectively becoming CDM.As shown in Eq. ( 19), the glueball self-interactions are repulsive and one may wonder if this feature affects structure formation.A simple estimate of this effect is obtained by comparing the intensity of self-interactions, proportional to the glueball field amplitude ϕ 0 , with the particle mass m gb [65].The field amplitude is directly connected to the glueball number density n gb = Ωρ c /m gb ≃ 94 cm −3 (Λ/ eV) −1 through ϕ 0 = n gb /2m gb ≃ 2.5 × 10 −7 eV(Λ/ eV) −1 .Since m gb ≫ ϕ 0 , the structure formation is indistinguishable from collisionless DM.This discussion is motivated by the flexibility in modeling interactions between dark gluons and ordinary matter, determining the cosmological evolution of the hidden sector.This results in a broad parameter space available to glueballs to explain the nature of DM.

VI. CONCLUSIONS
In this work we explored in detail the formation of scalar glueballs in generic SU (N ) dark gauge sectors.We showed that this composite state is a good DM candidate, viable in several cosmological models.The delicate interplay between microphysics governing the phase transition (the confinement-deconfinement phase transition scale Λ) and the macroscopic cosmological evolution (the dark-tovisible temperature ratio ζ −1 T ) determines the relic density of the glueball DM.In Fig. 6 we summarized our findings in the Λ vs ζ −1 T parameter space, showing that a large portion of it is viable and unconstrained.The red band shows the parameter space where glueballs constitute the whole DM, depending on the gauge group considered.Above this line, glueballs would overclose the Universe and this gray region is excluded.Moreover, we require the phase transition to happen in a radiationdominated era, when the photon bath temperature is assumed to be approximately below 10 16 GeV, otherwise the glueball relic density would be strongly suppressed by inflation.Precisely, we require that the confining scale matches the glueball temperature after inflation, i.e.T = Λ = ζ −1 T T γ < ζ −1 T 10 16 GeV.As another consistency condition, we also mark the region where the confinement scale is super-Planckian.The blue region shows the range of ζ −1 T easily accommodated in a freezeout scenario for dark gluons.For example, we find that a window with 2 keV ≲ Λ ≲ 21 keV would explain the DM in a simple freeze-out scenario for the dark gluons forming glueballs.However, this region is in tension with constraints on DM self-interactions, requiring the crosssection to mass ratio be σ/m ≲ 0.19 cm 2 /g [66].In our model glueballs undergo 2 → 2 scatterings with a cross section σ ∼ m −2 Λ , in the non-relativistic limit, exceeding the current constraints if Λ ≲ 17 MeV and glueballs constitute the majority of DM.
Cosmological observations of ∆N eff constrain models in which dark gluons are extremely hot, excluding the possibility for a very low confinement scale Λ ≪ 100 eV.In conclusion, this study allowed us to delineate a useful parameter space for glueball DM.We expect future studies using observational constraints from cosmology and astrophysics, as well as those produced from upcoming laboratory experiments, will populate Fig. 6 to provide a clearer picture and better understanding of glueball DM properties and behaviors.
This work sets the necessary basis for future investigations of the glueball phenomenology as a promising candidate for unveiling the dark sector.Due to its very interdisciplinary nature, a line of research focused on glueballs also opens up the possibility of various synergies between cosmological surveys and collider searches [67][68][69].Besides phenomenological applications, future investigations should focus on comparing the results obtained in the presented formalism with alternative methods to describe the phase transition, in order to assess the robustness of these calculations.

FIG. 1 .
FIG. 1. Left panel: Polyakov loop potential V[ℓ] for different gauge groups: SU (3) solid lines, SU (4) dashed lines and SU (5) dotted lines.The colors correspond to the confined (black) or deconfined (red) phase.Note that the minimum of the potential is arbitrarily set to zero.Right panel: Polyakov loop evolution as function of the temperature for different gauge groups: SU (3) solid lines, SU (4) dashed lines and SU (5) dotted lines.

FIG. 2 .
FIG. 2. Left panel: glueball potential V [ϕ, T ] for different gauge groups: SU (3) solid lines, SU (4) and SU (5) dashed lines.The colors correspond to the confined (black) or deconfined (red) phase.In the confined phase the potential is independent on the gauge group, while it is weakly dependent in the deconfined phase.Note that in this case the potentials for SU (4) and SU (5) are indistinguishable.Right panel: effective glueball mass as function of the temperature for various gauge groups: SU (3) solid line, SU (4) dashed line and SU (5) dotted line.

FIG. 3 .
FIG. 3. Left panel: glueball evolution obtained by solving Eq. (14) for SU (3) (black line), SU (4) (red line) and SU (5) (blue line) with µ = 0.05 and initial condition φ( Ti) = 0.1.Right panel: similar to the left panel but only for SU (3), with µ = 0.05 and different initial conditions.The vertical dashed line marks the phase transition, and the red dashed line shows the evolution of the minimum of the glueball potential.

FIG. 4 .
FIG. 4. Evolution of the glueball energy density obtained by solving Eq. (14) for SU (3) (black line), SU (4) (red line) and SU (5) (blue line) with µ = 1.The lines start in the point where the phase transition happens, which is different for the three cases shown here.

FIG. 6 .
FIG. 6. Parameter space of glueball DM expressed in terms of the confinement scale Λ and the dark-to-visible sector temperature ratio ζ −1 T .Each point in this parameter space represents a possible DM scenario with Λ fixed by the properties of the gauge sector and ζ −1 T by its cosmological evolution.The gray regions are excluded by various arguments (see the text for more details): the confinement scale has to be sub-Planckian; the phase transition (PT) has to happen after inflation; dark gluons cannot exceed the observed ∆N eff ; glueballs cannot overclose the Universe and DM self-interactions are constrained by observations.The red region corresponds to DM fully composed of glueballs; the variability range takes into account the differences between the gauge groups.The blue region for ζ −1 T ≃ 0.26 is favored in freeze-out models, corresponding to a phase transition scale 2 keV ≲ Λ ≲ 21 keV (values corresponding to the largest variability shown in Tab.II, 30 eV ≲ Λ0 ≲ 372 eV) to have all the DM in the form of glueballs.However, this region is in tension with constraints.

TABLE II .
Results of the calculation of the relic density.