We demonstrate how a dark sector can provide a crucial link between a vacuumdominated early universe and the later radiationdominated epoch. As an example we consider a feeble axionlike coupling of dark nonAbelian gauge fields with a single inflaton. In this scenario a dark heat bath emerges at the end of inflation. As the dark sector cools down its gauge fields confine into composite glueball states. The reheating of the visible sector is realized through portal interactions. After confinement discrete symmetries protect part of the glueball states that form relic dark matter while the others keep decaying into Standard Model degrees of freedom. The dark relic abundance depends mostly on the departure from equilibrium dynamics after the phase transition constraining the confinement scale of the dark sector. Moreover indirect detection and BigBang nucleosynthesis set bounds on the scale of portal interactions. We show that in a narrow but nonvanishing region of the parameter space the correct dark matter abundance and sufficient reheating temperature are simultaneously reached.
['Dark matter', 'Glueballs', 'Inflation', 'Phase transition', 'Reheating', 'High Energy Physics - Phenomenology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology']
Smooth reheating and dark matter via nonAbelian gauge theory
{'Smooth reheating and dark matter via non-Abelian gauge theory': "S. B iondini a , H. K ole ˇ sov ' a b , S. P rocacci c  \na Department of Physics, University of Basel, Klingelbergstr. 82, CH-4056 Basel, Switzerland  \nb Department of Mathematics and Physics, University of Stavanger, 4036 Stavanger, Norway  \nc Departement of Theoretical Physics, University of Geneva, Quai Ernest Ansermet 24, CH-1211 Geneva, Switzerland", 'Abstract': 'We demonstrate how a dark sector can provide a crucial link between a vacuum-dominated early universe and the later radiationdominated epoch. As an example, we consider a feeble axion-like coupling of dark non-Abelian gauge fields with a single inflaton. In this scenario, a dark heat bath emerges at the end of inflation. As the dark sector cools down, its gauge fields confine into composite glueball states. The reheating of the visible sector is realized through portal interactions. After confinement, discrete symmetries protect part of the glueball states that form relic dark matter, while the others keep decaying into Standard Model degrees of freedom. The dark relic abundance depends mostly on the departure from equilibrium dynamics after the phase transition, constraining the confinement scale of the dark sector. Moreover, indirect detection and Big-Bang nucleosynthesis set bounds on the scale of portal interactions. We show that in a narrow but non-vanishing region of the parameter space, the correct dark matter abundance and su ffi cient reheating temperature are simultaneously reached.  \nKeywords: dark matter, glueballs, inflation, phase transition, reheating', '1. Introduction': "Understanding the universe in terms of a few fundamental fields is perhaps the holy grail of particle cosmology. Assuming inflation to explain the origin of primordial perturbations implies the existence of a transition period, when the cold and inflated universe starts being populated by particles originating from the conversion of the inflaton energy. A successful nucleosynthesis demands for the onset of a Standard Model (SM) thermal bath at temperatures T SM ≳ 4 MeV [1-3]. Certainly, some mechanism had to be active and e ffi cient in producing a cosmological abundance of dark matter, which is about five times larger than that of the ordinary baryonic matter [4]. Current precision measurements e ffi ciently constrain inflationary models [5-7], as well as dark sectors [8-10].  \nIt is then quite appealing and nowadays accepted as part of the speculations in particle cosmology, that a dark sector may be coupled to the inflaton. This entails the possibility of a dark sector being populated first from inflaton decays, and the SM particles being in turn generated via portal interactions with the dark sector, see e.g. [11-14]. The main scope of our work is to link the various stages of the cosmological evolution, from inflation to the onset of the SM nucleosynthesis, for a quite minimal hidden sector: a non-Abelian pure Yang-Mills (YM) theory [15-27]. Dreaming of a complete framework, here we want to keep track of the continuous transition from the vacuumdominated (inflation) to the radiation-dominated (hot big bang)  \nepoch. In this work, we restrict to an inflaton solely coupled to a dark SU(N) gauge theory.  \nAs far as the inflationary stage is concerned, we consider axion-like (or natural) inflation [28] as the minimal inflation model, where a weakly-coupled thermal bath does not imply large back-reactions to the slow-roll regime [29-32]. The setup for the evolution equations that connect the inflaton field with the non-Abelian sector, together with a detailed study of the temperature of the YM plasma, have been put forward in ref. [33]. We follow up on the same scenario of a warm end of inflation by investigating the fate of the dark sector, asking whether the observed dark matter energy density and the generation of a su ffi ciently hot SM bath could be realized.  \nThe salient features of the dark sector are dictated by the behaviour of the non-Abelian gauge group. As the gauge coupling becomes large at low energies, dimensional transmutation implies a scale similar to Λ QCD for the SM strong interactions. We refer to this scale as Λ DS in the following. During the universe's thermal history, whenever the temperature of the dark sector crosses the scale Λ DS (from above), the physical degrees of freedom become a tower of hidden glueballs, as a result of a strongly-coupled dynamics.  \nGlueball states are often classified in terms of their quantum numbers under angular momentum ( J ), parity ( P ), and charge conjugation ( C ), which ultimately determine the e ff ective interactions with the visible sector [34, 35]. Glueballs can then generically decay into SM particles and provide a way to reheat the SM. However, some of the glueball states can be rendered cosmologically stable. More specifically, C-odd glueballs feature a suppressed decay rate with respect to C-even states, hence making the lightest C-odd glueball a viable dark matter candi-  \ndate [22]. This is the option that will be explored in our work.  \nWe show that it is challenging but possible to account for the observed dark matter relic density and to reheat the SM at temperatures larger than 4 MeV, while keeping the C-odd glueball dark matter cosmologically stable. Our study is complementary to earlier investigations [21, 22, 25, 26, 36], where the dark non-Abelian plasma is absent or very sub-dominant during the early stages of the universe thermal history, namely the inflaton is assumed to preferably reheat the SM.  \nThe structure of the paper is as follows. In section 2 we discuss the model for inflation and the dark sector, together with the portals with the SM. We present the evolution equations for the various components in section 3 and their numerical solution, together with the viable parameter space, in section 4. We o ff er some discussion and concluding remarks in section 5.", '2. Setup': 'This section introduces the benchmark models we use to describe inflation, the emergence of a dark sector, and the interactions of the latter with the SM.', '2.1. Minimal non-Abelian axion-like inflation': "We consider a variant of axion-like inflation [28] where the inflaton couples feebly to a non-Abelian dark sector, resulting in the production of YM fields already during the slow-roll inflationary phase [29, 30, 33]. Denoting φ the inflaton field, the Lagrangian reads  \nL ⊃ 1 2 ∂µφ ∂ µ φ -V ( φ ) -φ fa αϵ µνρσ F c µν F c ρσ 16 π + L dark YM . (1)  \nThe dark YM sector is described by the field strength, F c µν , the coupling, α , a colour index, c ∈{ 1 , . . . , N } , and L dark YM , which comprises the SU( N ) kinetic term. Notice that we set N = 3 in this work. For quantitative studies, we set the inflaton potential V ( φ ) = m 2 φ f 2 a [1 -cos( φ/ fa )], and fix decay constant, fa , and mass, m φ , using Planck's constraints [5], as in ref. [32]. 1 However, our final results do not depend on the details of the slow-roll phase, as discussed in sec. 5.  \nThe heating-up dynamics of the system in eq. (1) has been studied in ref. [33], showing a crucial dependence on the assumed confinement scale of the dark non-Abelian bath, Λ DS. In particular, if Λ DS is small with respect to the scale of inflation, the energy released by the weakly coupled inflaton su ffi ces to maintain the bath at a temperature above the critical scale, T c DS ≃ 1 . 24 Λ DS [37]. A phase transition must then occur in a later stage, while the bath is cooling down after inflation. In the present work, we show that adding a SM bath to this scenario can lead to successful reheating, and provide a viable dark matter candidate. Portal interactions between the SM and the dark sector might be active in both stages: first, when the dark vectors are weakly coupled, and later, when they are confined in composite glueball states.", '2.2. Dark vectors': 'In the deconfined phase, the non-Abelian dark sector is expected to thermalize 2 and for a given confinement scale Λ DS, its state can be described solely by its temperature, T DS. In order to capture the onset of non-perturbative dynamics approaching the critical point, we use a lattice-fit approach [33, 39]. Also for the SM we assume that it is in a thermal state within itself, 3 meaning that the dynamics of the system in this regime is governed by equilibrium attractors.  \nThe universality of gravitational interactions provides an unavoidable source for the visible sector from a dark thermal bath. SM particles are produced through gravitational interactions at tree level, through the exchange of a graviton in the s-channel [40, 41]. 4 A second option, which is often considered in the literature, envisages massive mediator states that couple to both the dark and visible sectors. If the characteristic mass scale of the mediators satisfies M ≫ Λ DS , T DS, an e ff ective description of the interactions can be employed and leading operators have mass dimension six and eight [34, 35].  \nIn particular, the dimension-six operator that involves the SM Higgs doublet, H , reads  \nL portal 6 ⊃ c 6 α tr( FF ) H † H M 2 , (2)  \nwhere c 6 is the matching coe ffi cient, while color and Lorentz indices are suppressed. Whenever M ≪ m pl the corresponding e ff ective operators are expected to be more e ffi cient for the production of the visible sector through scatterings. We adopt M as a free parameter of the model. The corresponding rates for the SM production are estimated as Γ G ≃ n eq G ⟨ σ GG→ SMSM v rel ⟩ with n eq G being the equilibrium number density for the dark gluons. In particular,  \nΓ grav G = T 5 DS m 4 pl 283 40 π 3 ζ (3) , (3)  \nΓ G = T 5 DS M 4 12 α 2 c 2 6 π 3 ζ (3) . (4)  \nThe graviton-mediated thermally averaged cross-section in eq. (3) is extracted from refs. [25, 42], bearing in mind that one needs dark vectors annihilation into SM particles; we carry out a similar calculation for the operator in (2). 5 The rate in eq. (4) is the leading contribution for the broken and unbroken phases of the SM and it corresponds to the processes GG → hh and GG → HH † respectively. In the broken phase, the 2 → 2 processes that are mediated by a Higgs-boson exchange, GG → h → SMSM, are suppressed as ( vh / T DS) 2 ≪ 1. Dimension-eight operators are instead suppressed by ( T DS / M ) 4 .  \nWe verify a posteriori that the portal interactions are never e ffi cient enough to set equal the temperatures in the visible and dark sectors. The maximal SM temperature reached in the benchmark viable scenario illustrated in fig. 1 is T SM ≈ 3 . 5 × 10 6 GeV ≪ T DS ≈ 8 . 4 × 10 10 GeV.', '2.3. Dark glueballs': 'Relevant cosmological information is contained in the departure from equilibrium of the dark glueballs after the phase transition. The glueball spectrum for SU(3) YM theories has been studied on the lattice [43-45]. Since the minimal YM action respects angular momentum, charge, and parity, glueballs are systematized in terms of the relevant quantum numbers J PC .  \nWepresent our analysis illustrating the dynamics of the lightest C -even and C -odd glueball, with masses m 0 ++ ≃ 6 . 8 Λ DS and m 1 + -≃ 11 . 5 Λ DS, respectively. 6 This minimal system has been shown to capture qualitatively the dynamics of the glueball ensemble [21, 22]. However, we need to model the thermodynamics of the dark glueball system also at temperatures close T c DS , where other glueballs are considerably abundant. To verify the two-glueballs approximation, we thus add the next lightest glueball, 2 ++ with m 2 ++ ≃ 9 . 4 Λ DS, to our numerical analysis.  \nDark glueballs interact with the SM particles via the same set of e ff ective operators as dark vectors [25, 34, 35], originating from graviton-mediation or a new-physics scales. However, such operators are now projected onto the composite states, hence requiring non-perturbative inputs [43-45]. Moreover, the qualitative di ff erence between C -even and C -odd states shapes the phenomenology. The lightest C-even glueball can decay through the dimension-6 operator (2), with the decay rate scaling as Λ 5 DS / M 4 . Going beyond the parametric estimation, we adapt the expression for the dominant decay process 0 ++ → hh from refs. [22, 34] for m 0 ++ ≫ mh , 7  \nΓ GB ≡ ⟨ Γ 0 ++ ⟩ ≃ c 2 6 32 π M 4 m 0 ++ ( α F 0 ++ S ) 2 K 1( z ) K 2( z ) GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> z = m 0 ++ / T DS . (5)  \nAlthough non-perturbative inputs for the matrix element exists also from lattice studies, 4 πα F 0 ++ S = 2 . 3(5) m 3 0 ++ [44, 47], we rather relay to the large N-limit to align eq. (5) with the approximation in eq. (7), F 0 ++ S ∼ m 3 0 ++ /α [22]. The thermal average is represented by the factor K 1 / K 2 ∈ [0 . 8 , 1 . 0], with Kn = z -n R ∞ z d x x n -1 e -x √ x 2 -z 2 extracted as in footnote 5.  \nConversely, the lightest C-odd state can only decay via dimension-8 operators of the schematic form [34]  \nL portal 8 ⊃ c 8 α 3 / 2 α 1 / 2 QED tr( FFF ) µν B µν M 4 , (6)  \nwhere B µν is the SM hypercharge field strength and c 8 is again a matching coe ffi cient. This operator induces the decays  \n1 + -→ 0 ++ γ and 1 + -→ 0 ++ Z , with widths scaling as Λ 9 DS / M 8 . The 1 + -state is cosmologically stable, if its lifetime exceeds the age of the universe, yielding a dark matter candidate, even if C-even glueballs decay before BBN. In the m 1 + -, m 0 ++ ≫ mZ approximation, the full decay width is  \nΓ DM ≡ ⟨ Γ 1 + -⟩ ≃ c 2 8 α QED 24 π c 2 W m 3 1 + -GLYPH<16> α 3 / 2 M 1 + -0 ++ GLYPH<17> 2 M 8       1 -m 2 0 ++ m 2 1 + -      3 , (7)  \nwhere c W is the cosine of Weinberg angle. We omit the thermal average suppression, since K 1 / K 2 → 1 at the temperatures relevant for this decay. For the matrix element M 1 + -0 ++ , lattice results are not available. The largeN limit suggests M 1 + -0 ++ α 3 / 2 ∼ √ 4 π/ N m 3 0 ++ [22].  \nIn addition to glueball-SM interactions, the dark composite states will interact with one another as a remnant of the strong dynamics that bind them. The corresponding processes, most notably 3 → 2 and 2 → 2 reactions, are key to determining the cosmological evolution of their densities (see sec. 3). Despite the interactions between glueballs being hard to determine from first principles, one expects the multi-glueball processes to be consistent with J , P , and C conservation, and we again borrow parametric estimates for interaction strengths from results in the largeN expansion (see Appendix A for details).', '3. Evolution equations': 'In order to track the evolution of the di ff erent species, we rely on a network of evolution equations, that comprise the interaction rates introduced in sec. 2, and the Hubble rate,  \nH 2 = 8 π 3 m 2 pl GLYPH<16> e φ + e DS + e SM GLYPH<17> . (8)  \nWe use the approximate initial value H 2 ref ≡ 8 π V ( φ 0) / (3 m 2 pl ) as an inverse-time unit. 8 The energy density of inflaton, e φ , dark sector, e DS, and SM, e SM, follow coupled evolution equations, which we now turn to study. A benchmark solution is shown in fig. 1.', '3.1. Energy densities': "The evolution of the inflaton embedded within a medium is taken over from refs. [30, 33]. In the post-inflationary era, the inflaton energy density is dominated by the kinetic term, while its pressure vanishes. 9 If we consider also the energy transfer between dark and visible sectors, the evolution of the system can be described by the following set of coupled equations,  \n˙ e φ + 3 He φ = -Υ e φ (9)  \n˙ e DS + 3 H ( e DS + p DS) = Υ e φ -Γ e DS (10)  \n˙ e SM + 4 He SM = Γ e DS . (11)  \nThe friction coe ffi cient Υ in eqs. (9) and (10), transfers energy from the inflaton to the dark sector degrees of freedom. For temperatures T ≪ m φ relevant in our setup, it is dominated by the vacuum decay rate Υ ≃ α 2 m 3 φ / (32 π 3 f 2 a ) [30, 33, 49].  \nAs long as the dark sector stays in the deconfined phase, we parameterize its pressure p DS and energy density e DS by a temperature, T DS, using SU(3) lattice fits [39]. Even if the friction term is suppressed, Υ ≪ H , in this first stage, dark vectors are e ffi ciently populated. Relying on their fast thermalization [38, 50], the SM is in turn produced from frequent 2 → 2 scatterings of vectors ' a la freeze-in. The corresponding rate is Γ = Γ G as given in eq. (4).  \nDuring the phase transition at T DS = T c DS , the dark sector is in a mixed state, where dark vectors and glueballs are present at the same time [33]. In this regime, within the adiabatic approximation, we assume equilibrium number densities for the glueball states and take into account the di ff erent decay rates to SM for the two phases.  \nOnce all vectors have been confined into composite states, they are expected to remain in thermal equilibrium until a certain freeze-out temperature T fo DS < T c DS . The thermodynamic quantities describing the full dark sector, e DS and p DS, are then functions of glueball number densities and T DS, according to the non-relativistic relations 10  \ne DS = X i mini 1 + 3 2 T DS mi ! , p DS = T DS X i ni (12)  \nwhere i runs over di ff erent glueball states. The modeling in eq. (12) is supported by lattice simulation of a pure SU(3) gauge theory [51]. 11 Eqs. (9)-(11) for the energy densities need hence to be complemented by Boltzmann equations for the glueball number densities, capturing their departure from equilibrium (c.f. sec. 3.2).  \nIn the confined phase, the inflaton can still lose energy in the form of dark vectors, which hadronize quickly and e ff ectively contribute to the glueball energy density. However, in our parameter space of interest, the inflaton field usually vanishes around the time of the phase transition and does not a ff ect the freeze-out dynamics.  \nAt the same time, C-even glueballs decay into SM particles, hence Γ e DS in eqs. (10) and (11) has to be replaced by Γ GB e UG, where Γ GB was introduced in eq. (5) and e UG is the energy density in unstable glueballs. In turn, C-odd glueballs remain cosmologically stable.", '3.2. Boltzmann equations for dark glueball states': 'This is how we model the freeze-out dynamics of the glueballs, that eventually sets the dark matter relic abundance.  \nAs the dark sector is made of glueballs with di ff erent masses and interactions, the composite states fall out of chemical equilibrium at di ff erent temperatures and develop individual chemical potentials. In order to track the glueball number densities, we employ a network of Boltzmann equations that is based on earlier studies [21], augmented by the feed-down of glueballs from the decaying inflaton, as well as the glueballs decay into SM. It was shown in [21] that a simplified scenario for the densities of the 0 ++ and 1 + -glueballs captures to a very good approximation the freeze-out dynamics of the whole ensemble. For simplicity, we adopt this choice when presenting the coupled Boltzmann equations in the main body of the paper. Hence, e UG ≃ e 0 ++ and the stable glueball energy density is stored in e 1 + -. We stress that, in order to improve the thermodynamic description of the YM sector [51], we shall extend the treatment by including the next-to-lightest glueball state 2 ++ in our numerical analysis (see section 4). We collect the corresponding equations in Appendix B.  \nThe coupled Boltzmann equations for the two-glueballs system read  \n˙ n 1 + -+ 3 Hn 1 + -= B 1 + -φ Υ e φ m 1 + -- ⟨ σ 2 → 2 v ⟩        n 2 1 + - -      n 0 ++ n eq 0 ++       2 ( n eq 1 + -) 2        , (13)  \n˙ n 0 ++ + 3 Hn 0 ++ = B 0 ++ φ Υ e φ m 0 ++ + ⟨ σ 2 → 2 v ⟩        n 2 1 + - -      n 0 ++ n eq 0 ++       2 ( n eq 1 + -) 2        -n 0 ++ Γ GB - ⟨ σ 3 → 2 v 2 ⟩ n 2 0 ++ ( n 0 ++ -n eq 0 ++ ) . (14)  \nThermally averaged cross-sections ⟨ σ 2 → 2 v ⟩ (here for 1 + -1 + -→ 0 ++ 0 ++ ) and ⟨ σ 3 → 2 v 2 ⟩ (here changing the number of 0 ++ states) are extracted from ref. [21], as given in Appendix A.  \nThe relative abundance of di ff erent glueball species sourced (non-thermally) from the inflaton decays may be inferred phenomenologically from statistical hadronization models. Most notably, the branching ratio of a given glueball is proportional to its equilibrium distribution evaluated at a hadronization temperature , T had ≃ T c DS [52, 53]. This is how we extract the branching ratios B J PC φ that enter eqs. (13) and (14). These terms are of little importance below a certain threshold value of Λ DS, 12 since then the phase transition in the dark sector happens after e φ → 0.', '4. Numerical solutions and constraints': "The main result of our study is a scan of numerical solutions for di ff erent values of the parameters Λ DS and M . An example of full solution for the energy densities of di ff erent sectors throughout the di ff erent stages described in sec. 3 is shown in  \nFigure 1: Numerical solution for the energy density for the various components. The dotted line stands for the inflaton, the dashed line for the dark YM, and the solid curve for the SM. The vertical gray band indicates the time at which occurrence of the phase transition in the YM sector.  \n<!-- image -->  \nFigure 2: Number densities for the C-even 0 ++ and C-odd 1 + -glueballs respectively in blue (upper) and red (lower) solid lines. By n eq i we denote the equilibrium density with zero chemical potential. The normalization n 0 ≡ n 0 ++ ( T c DS ) corresponds to the number density of 0 ++ at the end of the phase transition.  \n<!-- image -->  \nfig. 1. 13 We rescale dimensionful quantities by powers of H ref, defined below eq. (8).  \nThe period when e φ ≈ const corresponds to the slow-roll phase of inflation. When φ approaches the minimum of the potential, e φ starts decreasing and the expansion of the universe becomes decelerated, translating in a milder dilution of e DS and e SM. At the same time, both dark and visible sectors are still sourced by the highly oscillating inflaton, so that, at this stage, they reach their maximum temperatures. This turning point is followed by an early matter-dominated era, when e φ ≫ e DS ≫ e SM. For viable parameter choices, e φ vanishes close to the phase transition in the YM sector. Although SM fields are sourced both before and after the phase transition, the decay of a glueball is more likely than the GG → SMSM scattering of two dark gluons (c.f. eqs. (5) and (4) respectively). The consequence is a 'kink' on the e SM line at the beginning of the phase transition in the dark sector.  \nWhen the phase transition is completed, e DS contains both stable and unstable glueballs. Their respective fractions vary according to the Boltzmann equations, until finally the unstable glueballs decay. This happens when H ≲ Γ GB. The DS consists then of the stable glueballs only, which evolve decoupled from the SM sector, and the standard radiation-dominated era starts.", '4.1. Cosmological and astrophysical constraints': 'Let us now turn to the question of what parameter choices provide a viable scenario. By viable we mean, in particular, capable of providing the correct relic abundance of su ffi ciently long-lived dark matter, while achieving a successful SM reheating. In practice, we keep the inflationary dynamics intact and inspect for which dark sector confinement scale, Λ DS, and which portal scale, M , the following conditions hold  \nτ DM = 1 / Γ DM ≳ 10 26 s , (15)  \nT rh SM ≥ 4 MeV , (16)  \nΩ DM h 2 ≤ 0 . 1200 ± 0 . 0012 . (17)  \nEq. (15) corresponds to an approximate bound coming from indirect detection experiments, widely used in literature [54, 55]. This constraint is expected to become more stringent, up to τ DM ∼ 10 28 s [56, 57]. On the other hand, the precise theoretical estimate of Γ DM crucially depends on the branching ratios of the decaying dark matter to di ff erent SM states, and we choose to leave this analysis to future work. Moreover, the interpretation of eq. (15) in terms of constraints on Λ DS and M is a ff ected also by the uncertainties in the non-perturbative input and matching coe ffi cient, c 8, as presented in eq. (7). Keeping these opposite e ff ects in mind, in this work we resort ourselves to τ DM c 2 8 ∼ 10 26 -δ s, δ ∈ [0 , 1].  \nIn eq. (16), T rh SM is the SM temperature when its energy density starts dominating, i.e., after the decay of unstable glueballs. The Hubble rate at this time satisfies H ≈ Γ GB, and, approximating H 2 ≈ 8 π e SM / (3 m 2 pl ), we find an analytical estimate for the minimal decay rate required by eq. (16),  \nΓ GB ≥ Γ min GB ≈ 5 . 8 × 10 -43 m pl . (18)  \nCombining eqs. (18) and (5) sets a constraint on a further combination of Λ DS and M . For the results illustrated in fig. 3, we  \nextract T rh SM from the numerical solution for e SM of the evolution equations, using the SM equation of state provided in [58], finding a slightly larger value, Γ min GB | num ≈ 1 . 4 × 10 -42 m pl. Qualitatively, however, the SM reheating temperature indeed depends solely on the value of Γ GB and not on other combination of the parameters Λ DS and M .  \nFinally, for eq. (17), the dark matter relic abundance is dictated by the freeze-out dynamics, which requires a numerical solution of the full system, as presented in eqs. (9)-(14).', '4.2. Freeze-out dynamics and parameter space': 'At the core of the extraction of the dark matter energy density one finds the solution of coupled Boltzmann equations. We solved the equations for a two-glueballs system, described by eqs. (13)-(14), and for a three-glueballs system, as in eqs. (B.1)-(B.3). The result shown in fig. 2 corresponds to the latter case. The observed dynamics is qualitatively the same: the freeze-out of even states happens when the 3 → 2 interactions cease to be e ffi cient, and a non-zero chemical potential develops. After this point, the 2 → 2 reactions maintain the value of the chemical potentials such that n 1 + -/ n eq 1 + - ≈ n 0 ++ / n eq 0 ++ . The relic abundance of the stable glueballs is set at the point when also the 2 → 2 interactions become ine ffi cient.  \nSince the 2 → 2 cross section scales as Λ -2 DS , (c.f. eq. (A.4)), for smaller Λ DS, the interaction strength is larger, the freeze-out happens later, and the relic abundance is suppressed. Quantitatively, the addition of the 2 ++ state acts as a doubling of the interaction channels and thus of the cross-sections, which translates into n 1 + -being half as large at freeze-out completion.  \nThe resulting parameter space is shown in fig. 3, where the viable region is unshaded (white). The interplay of the conditions in eqs. (15)-(17) is di ff erent for the three-glueballs system (solid lines) and for the two-glueballs system (dashed lines).  \nThe numerical results for T rh SM = 4 MeV (purple) confirm the estimate in eq. (18): a too small Λ DS or a too large M suppress the decay of the even glueballs, resulting in an ine ffi cient energy transfer to the visible sector. Adding the 2 ++ state has little impact on this line. On the other hand, the curves depicting Ω DM h 2 = 0 . 12 (green) di ff er considerably for the 2-glueball and 3-glueball scenarios as expected based on the discussion of the interaction strengths above.  \nOur numerical results also confirm that larger Λ DS leads to larger DM relic abundance. Larger M , in turn, implies that the decay of the glueballs to SM is less e ffi cient, and radiation domination appears later, leading to more dilution of the DM density. 14 Finally, exclusion limits from indirect detection (blue), c.f. eq. (15), are identical for both cases , but a ff ected by the systematic uncertainty discussed above (see the di ff erence between the dotted and solid lines showing the possibility of enlarging the white region).  \nLet us stress the large impact of adding the next-to-leading glueball contribution to the dark energy density e DS. A realistic  \nFigure 3: Parameter space of the model in the plane ( M , Λ DS). The white unshaded triangle is the viable region that complies with the cosmological and astropysical conditions (15)-(17). The solid lines correspond to the three-glueball case described by eqs. (B.1)-(B.3) whereas the dashed lines correspond to the two-glueball one, eqs. (13)-(14). For the latter case, the open window is in the bottom-left corner, shown magnified in the inset (top-left corner).  \n<!-- image -->  \ndetermination of the parameter space must therefore take into account higher excitations when the system falls out of equilibrium. On the other hand, we expect that adding even heavier glueballs will have less dramatic e ff ect due to more significant Boltzmann suppression of their number densities. We remark that O (1) uncertainties are ubiquitous in the modeling of glueball dynamics, which range from the complexity of the dark phase transition to the proper inclusion of glueball interactions.', '5. Conclusions and outlook': 'We have presented a model providing successful inflation and reheating of SM fields and, at the same time, a viable dark matter candidate with correct relic abundance. In summary, the thermal history that is expected to occur in our framework can be summarized as follows. The inflaton decays into nonAbelian dark vectors that form a thermal bath of dark vectors at high temperatures. Already at this stage, portal interactions between the dark and the visible sectors are active, sourcing a non-negligible visible bath. After inflation, the expansion of the universe induces a cooling of the dark plasma that undergoes a phase transition from a deconfined to a confined phase, at a temperature T c DS ∼ Λ DS. The interactions among glueballs trigger various processes that control the number density of the di ff erent dark composite states [18, 21]. Together with the glueball decays into visible particles, 3 → 2 and 2 → 2 scatterings enter a network of evolution equations that ultimately determine  \nthe dark matter relic density, as well as the SM reheating temperature.  \nOur framework is relatively general and can be applied to different choices in the dark and inflationary sectors. In particular, although the viable region in the parameter space appears small for our particular implementation (see fig. 3), generalizations of our scenario might enlarge it. For example, one may consider di ff erent gauge groups, e.g. SO( N ) with N ≥ 8 [25, 27, 35]. The decay rate of stable glueballs scales then as Λ 2 N -3 DS / M 2 N -4 , such that the blue line in fig. 3 would be shifted upwards correspondingly. However, apart from few pioneering studies like [59, 60], the non-perturbative input on the glueball spectrum and interactions is mostly missing in those theories. For these groups, the stable glueball would likely be too long-lived to be probed by indirect detection experiments. On the other hand, gravitational wave signatures might o ff er further observational constraints, as discussed in [33, 61-64]. In particular, the confinement phase transition in di ff erent pure Yang-Mills theories is known to be of first-order [65, 66], for the case at hand leading to a peak frequency in the LISA range, f 0 ∼ O (1) Hz [62]. Per contra, the spectral amplitude might be suppressed by the long-lasting matter-dominated period [67].  \nOther generalizations may also consider direct decays of the inflaton into SM degrees of freedom, e.g. via a similar axionlike coupling. This has been considered in the literature of warm inflation, e.g. [68]. Since Λ QCD ≪ Λ DS, the couplings satisfy α QCD <α , inducing a corresponding hierarchy in the friction terms. Hence, the dark sector could still dominate the total energy density for a certain amount of time, while the DM relic abundance would be reduced and the parameter space enlarged. Further freedom is also in the details of inflation itself, as these turn out to barely enter the cosmological observables in fig. 3. In scenarios with a larger friction coe ffi cient, Υ , the inflaton decays earlier, and our results are unchanged. On the other hand, for feebler interaction between inflaton and the dark sector, the inflaton might a ff ect the freeze-out dynamics.  \nFinally, UV completions of our scenario, leading to the operators in eqs. (2) and (6), could be related to other puzzles in particle cosmology as neutrino masses or baryon asymmetry.  \nNote Added : During the revision of our manuscript the study of ref. [69] appeared that investigates a closely related scenario of dark glueballs.', 'Acknowledgements': 'We thank Mikko Laine for comments on the manuscript and Torsten Bringmann, Jacopo Ghiglieri, Admir Greljo, Giacomo Landini, Harvey Meyer, Marco Panero and Roman Pasechnik for useful discussions. The work of S.B. and S.P. is supported by the Swiss National Science Foundation, under the Ambizione grant 2185783 and under grant 212125 respectively. H.K. is supported by the Research Council of Norway under the FRIPRO Young Research Talent grant no. 335388.', 'Appendix A. Glueball Lagrangian and cross sections': "Based on naive dimensional analysis [70-72] and largeN scaling [73, 74], we consider the following Lagrangian density for a given glueball ϕ GB [17, 21, 22]  \nL ⊃ 1 2 ( ∂µϕ GB) 2 -a 2 2! m GB ϕ 2 GB -a 3 3! 4 π N ! m GB ϕ 3 GB -a 4 4! 4 π N ! 2 ϕ 4 GB . (A.1)  \nIn the first line, the glueball kinetic and mass terms are accompanied by three- and four-glueball self-interactions. Terms with higher powers of ϕ GB within the e ff ective potential approach are considered in other works, e.g., [23, 24]. Such additional interactions postpone the freeze-out and reduce the dark matter relic abundance. However, they also introduce additional unknown coe ffi cients within our approach and we choose not to include them. The coe ffi cients ai are expected to be of order unity.  \nThe glueball interaction terms lack a detailed nonperturbative treatment at the moment. Even if the schematic form of the e ff ective Lagrangian in eq.(A.1) may hold for interactions among di ff erent glueballs, it is strictly appropriate only for scalar glueball states, such as 0 ++ . For the 3 → 2 process of the lightest glueball, 0 ++ , we improve upon the results of ref. [21] by computing the diagrams at order (4 π/ N ) 3 and find  \n⟨ σ 000;00 v 2 ⟩ ≃ 1805 π 5 a 2 3 388962 N 6 m 5 1 GLYPH<16> 11 a 2 3 + 7 a 4 GLYPH<17> 2 , (A.2)  \nat leading order in the velocity expansion. The cross-section in eq. (A.2) is about five times smaller than the estimation in refs. [21, 22] for ai = 1, which reads ⟨ σ 000;00 v 2 ⟩≃ (4 π ) 3 / ( N 6 m 5 0 ++ ).  \nProceeding with analogous computations for the decays and interactions involving glueballs of higher angular momenta would require a precise formulation of the corresponding effective Lagrangians, which goes beyond the goal of this paper. In order to align the approximation on the inputs for our equations, we eventually adopt the estimation in terms of the power counting and use [21, 22]. The corresponding cross-section for the next-to-lightest C-even glueball, 2 ++ is then  \n⟨ σ 222;22 v 2 ⟩ ≃ (4 π ) 3 N 6 m 5 2 ++ . (A.3)  \nAlso for the 2 → 2 processes, we simply estimate of the squared amplitude in terms of the power counting of the e ff ective Lagrangian, and we obtain at leading order in the velocity expansion and with ai = 1 the expression  \n⟨ σ ab ; cdv ⟩ ≃ (4 π/ N ) 4 gagbScd 64 π 2 s λ 1 / 2 GLYPH<16> 1 , m 2 c / s , m 2 d / s GLYPH<17> λ 1 / 2 GLYPH<16> 1 , m 2 a / s , m 2 b / s GLYPH<17> . (A.4)  \nHere, λ ( x , y , z ) = ( x -y -z ) 2 -4 yz is the Kall'en function, gi = (2 Ji + 1), s = ( ma + mb ) 2 , and Scd is the symmetry factor accounting for identical particles in the final state.", 'Appendix B. Three-glueballs evolution equations': 'The evolution equations describing the freeze-out dynamics of two unstable states, 0 ++ and 2 ++ , and one stable, 1 + -, read  \n˙ n 1 + -+ 3 Hn 1 + -= B 1 + -φ Υ e φ m 1 + -- ⟨ σ 11;00 v ⟩        n 2 1 + - -      n 0 ++ n eq 0 ++       2 ( n eq 1 + -) 2         \n- ⟨ σ 11;22 v ⟩        n 2 1 + - -      n 2 ++ n eq 2 ++       2 ( n eq 1 + -) 2         \n- ⟨ σ 11;02 v ⟩       n 2 1 + - -n 0 ++ n 2 ++ n eq 0 ++ n eq 2 ++ ( n eq 1 + -) 2       , (B.1)  \n˙ n 0 ++ + 3 Hn 0 ++ = B 0 ++ φ Υ e φ m 0 ++ -n 0 ++ Γ 0 ++ + n 2 ++ Γ 2 ++  \n+ ⟨ σ 11;00 v ⟩        n 2 1 + - -      n 0 ++ n eq 0 ++       2 ( n eq 1 + -) 2         \n- ⟨ σ 000;00 v 2 ⟩ n 2 0 ++ ( n 0 ++ -n eq 0 ++ )  \n+ ⟨ σ 22;00 v ⟩        n 2 2 ++ -      n 0 ++ n eq 0 ++       2 ( n eq 2 ++ ) 2         \n+ ⟨ σ 11;02 v ⟩ 2       n 2 1 + - -n 0 ++ n 2 ++ n eq 0 ++ n eq 2 ++ ( n eq 1 + -) 2        \n+ " n 0 ++ ⟨ σ 20;00 v ⟩ + n 2 ++ ⟨ σ 22;20 v ⟩ 2  \n+ n 1 + -⟨ σ 12;10 v ⟩ #       n 2 ++ -n eq 2 ++ n eq 0 ++ n 0 ++       , (B.2)  \n˙ n 2 ++ + 3 Hn 2 ++ = B 2 ++ φ Υ e φ m 2 ++ -n 2 ++ Γ 2 ++  \n+ ⟨ σ 11;22 v ⟩        n 2 1 + - -      n 2 ++ n eq 2 ++       2 ( n eq 1 + -) 2         \n- ⟨ σ 222;22 v 2 ⟩ n 2 2 ++ ( n 2 ++ -n eq 2 ++ )  \n- ⟨ σ 22;00 v ⟩        n 2 2 ++ -      n 0 ++ n eq 0 ++       2 ( n eq 2 ++ ) 2         \n+ ⟨ σ 11;02 v ⟩ 2       n 2 1 + - -n 0 ++ n 2 ++ n eq 0 ++ n eq 2 ++ ( n eq 1 + -) 2        \n-" n 0 ++ ⟨ σ 20;00 v ⟩ + n 2 ++ ⟨ σ 22;20 v ⟩ 2  \n+ n 1 + -⟨ σ 12;10 v ⟩ #       n 2 ++ -n eq 2 ++ n eq 0 ++ n 0 ++       . (B.3)  \nThe lightest state, 0 ++ , decays to the SM via Γ 0 ++ = Γ GB in eq. (5). Since the next-to-lightest state, 2 ++ , has spin two, it decays through di ff erent operators. In particular, the direct decay to SM fields in analogy to Γ GB implies a dimension-eight operator. However, a dimension-six contribution is allowed for our setup via the process 2 ++ → 0 ++ h [34]. Its decay rate is  \n⟨ Γ 2 ++ ⟩ ≃ c 2 6 480 π m 3 2 ++        α M 2 ++ , 0 ++ S M 2        2 1 -m 0 ++ m 2 ++ ! 5 K 1( z ) K 2( z ) GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> z = m 2 ++ / T DS , (B.4)  \nwith the matrix element M 2 ++ , 0 ++ S ≈ 0 . 03 m 0 ++ [35].', 'References': '- [1] M. Kawasaki, K. Kohri, N. Sugiyama, MeV scale reheating temperature and thermalization of neutrino background, Phys. Rev. D 62 (2000) 023506. arXiv:astro-ph/0002127 , doi:10.1103/PhysRevD.62. 023506 .\n- [2] S. Hannestad, What is the lowest possible reheating temperature?, Phys. Rev. D 70 (2004) 043506. arXiv:astro-ph/0403291 , doi:10.1103/ PhysRevD.70.043506 .\n- [3] F. De Bernardis, L. Pagano, A. Melchiorri, New constraints on the reheating temperature of the universe after WMAP-5, Astropart. Phys. 30 (2008) 192-195. doi:10.1016/j.astropartphys.2008.09.005 .\n- [4] N. Aghanim, et al., Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641 (2020) A1. arXiv: 1807.06205 , doi:10.1051/0004-6361/201833880 .\n- [5] Y. Akrami, et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10. arXiv:1807.06211 , doi:10.1051/ 0004-6361/201833887 .\n- [6] P. A. R. Ade, et al., BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2 / Keck Observations through the 2015 Season, Phys. Rev. Lett. 121 (2018) 221301. arXiv:1810.05216 , doi:10.1103/PhysRevLett.121.221301 .\n- [7] P. A. R. Ade, et al., Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP / Keck Observations through the 2018 Observing Season, Phys. Rev. Lett. 127 (15) (2021) 151301. arXiv:2110.00483 , doi:10.1103/PhysRevLett.127.151301 .\n- [8] P. Harris, V. V. Khoze, M. Spannowsky, C. Williams, Constraining Dark Sectors at Colliders: Beyond the E ff ective Theory Approach, Phys. Rev. D 91 (2015) 055009. arXiv:1411.0535 , doi:10.1103/PhysRevD. 91.055009 .\n- [9] G. Elor, N. L. Rodd, T. R. Slatyer, W. Xue, Model-Independent Indirect Detection Constraints on Hidden Sector Dark Matter, JCAP 06 (2016) 024. arXiv:1511.08787 , doi:10.1088/1475-7516/2016/06/024 .\n- [10] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, F. S. Queiroz, The waning of the WIMP? A review of models, searches, and constraints, Eur. Phys. J. C 78 (3) (2018) 203. arXiv:1703.07364 , doi:10.1140/epjc/s10052-018-5662-y .\n- [11] S. Biondini, K. Sravan Kumar, Dark matter and Standard Model reheating from conformal GUT inflation, JHEP 07 (2020) 039. arXiv: 2004.02921 , doi:10.1007/JHEP07(2020)039 .\n- [12] A. Berlin, D. Hooper, G. Krnjaic, PeV-Scale Dark Matter as a Thermal Relic of a Decoupled Sector, Phys. Lett. B 760 (2016) 106-111. arXiv: 1602.08490 , doi:10.1016/j.physletb.2016.06.037 .\n- [13] T. Tenkanen, V. Vaskonen, Reheating the Standard Model from a hidden sector, Phys. Rev. D 94 (8) (2016) 083516. arXiv:1606.00192 , doi: 10.1103/PhysRevD.94.083516 .\n- [14] A. Paul, A. Ghoshal, A. Chatterjee, S. Pal, Inflation, (P)reheating and Neutrino Anomalies: Production of Sterile Neutrinos with Secret Interactions, Eur. Phys. J. C 79 (10) (2019) 818. arXiv:1808.09706 , doi:10.1140/epjc/s10052-019-7348-5 .\n- [15] E. D. Carlson, M. E. Machacek, L. J. Hall, Self-interacting dark matter, Astrophys. J. 398 (1992) 43-52. doi:10.1086/171833 .\n- [16] A. E. Faraggi, M. Pospelov, Selfinteracting dark matter from the hidden heterotic string sector, Astropart. Phys. 16 (2002) 451-461. arXiv: hep-ph/0008223 , doi:10.1016/S0927-6505(01)00121-9 .\n- [17] K. K. Boddy, J. L. Feng, M. Kaplinghat, T. M. P. Tait, Self-Interacting Dark Matter from a Non-Abelian Hidden Sector, Phys. Rev. D 89 (11) (2014) 115017. arXiv:1402.3629 , doi:10.1103/PhysRevD.89. 115017 .\n- [18] A. Soni, Y. Zhang, Hidden SU(N) Glueball Dark Matter, Phys. Rev. D 93 (11) (2016) 115025. arXiv:1602.00714 , doi:10.1103/ PhysRevD.93.115025 .\n- [19] G. D. Kribs, E. T. Neil, Review of strongly-coupled composite dark matter models and lattice simulations, Int. J. Mod. Phys. A 31 (22) (2016) 1643004. arXiv:1604.04627 , doi:10.1142/S0217751X16430041 .\n- [20] J. Halverson, B. D. Nelson, F. Ruehle, String Theory and the Dark Glueball Problem, Phys. Rev. D 95 (4) (2017) 043527. arXiv:1609.02151 , doi:10.1103/PhysRevD.95.043527 .\n- [21] L. Forestell, D. E. Morrissey, K. Sigurdson, Non-Abelian Dark Forces and the Relic Densities of Dark Glueballs, Phys. Rev. D 95 (1) (2017) 015032. arXiv:1605.08048 , doi:10.1103/PhysRevD.95.015032 .\n- [22] L. Forestell, D. E. Morrissey, K. Sigurdson, Cosmological Bounds on Non-Abelian Dark Forces, Phys. Rev. D 97 (7) (2018) 075029. arXiv: 1710.06447 , doi:10.1103/PhysRevD.97.075029 .\n- [23] P. Carenza, R. Pasechnik, G. Salinas, Z.-W. Wang, Glueball Dark Matter Revisited, Phys. Rev. Lett. 129 (26) (2022) 261302. arXiv:2207. 13716 , doi:10.1103/PhysRevLett.129.261302 .\n- [24] P. Carenza, T. Ferreira, R. Pasechnik, Z.-W. Wang, Glueball dark matter, precisely, Phys. Rev. D 108 (12) (2023) 123027. arXiv:2306.09510 , doi:10.1103/PhysRevD.108.123027 .\n- [25] C. Gross, S. Karamitsos, G. Landini, A. Strumia, Gravitational Vector Dark Matter, JHEP 03 (2021) 174. arXiv:2012.12087 , doi:10.1007/ JHEP03(2021)174 .\n- [26] M. Redi, A. Tesi, H. Tillim, Gravitational Production of a Conformal Dark Sector, JHEP 05 (2021) 010. arXiv:2011.10565 , doi:10.1007/ JHEP05(2021)010 .\n- [27] M. Yamada, K. Yonekura, Dark baryon from pure Yang-Mills theory and its GW signature from cosmic strings, JHEP 09 (2023) 197. arXiv: 2307.06586 , doi:10.1007/JHEP09(2023)197 .\n- [28] K. Freese, J. A. Frieman, A. V. Olinto, Natural inflation with pseudo Nambu-Goldstone bosons, Phys. Rev. Lett. 65 (1990) 3233-3236. doi: 10.1103/PhysRevLett.65.3233 .\n- [29] K. V. Berghaus, P. W. Graham, D. E. Kaplan, Minimal Warm Inflation, JCAP 03 (2020) 034, [Erratum: JCAP 10, E02 (2023)]. arXiv:1910. 07525 , doi:10.1088/1475-7516/2020/03/034 .\n- [30] M. Laine, S. Procacci, Minimal warm inflation with complete medium response, JCAP 06 (2021) 031. arXiv:2102.09913 , doi:10.1088/ 1475-7516/2021/06/031 .\n- [31] P. Klose, M. Laine, S. Procacci, Gravitational wave background from nonAbelian reheating after axion-like inflation, JCAP 05 (2022) 021. arXiv: 2201.02317 , doi:10.1088/1475-7516/2022/05/021 .\n- [32] P. Klose, M. Laine, S. Procacci, Gravitational wave background from vacuum and thermal fluctuations during axion-like inflation, JCAP 12 (2022) 020. arXiv:2210.11710 , doi:10.1088/1475-7516/2022/12/020 .\n- [33] H. Kolesova, M. Laine, S. Procacci, Maximal temperature of stronglycoupled dark sectors, JHEP 05 (2023) 239. arXiv:2303.17973 , doi: 10.1007/JHEP05(2023)239 .\n- [34] J. E. Juknevich, Pure-glue hidden valleys through the Higgs portal, JHEP 08 (2010) 121. arXiv:0911.5616 , doi:10.1007/JHEP08(2010)121 .\n- [35] J. E. Juknevich, D. Melnikov, M. J. Strassler, A Pure-Glue Hidden Valley I. States and Decays, JHEP 07 (2009) 055. arXiv:0903.0883 , doi: 10.1088/1126-6708/2009/07/055 .\n- [36] Q. Li, T. Moroi, K. Nakayama, W. Yin, Hidden dark matter from Starobinsky inflation, JHEP 09 (2021) 179. arXiv:2105.13358 , doi: 10.1007/JHEP09(2021)179 .\n- [37] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, H. Ohno, Critical point and scale setting in SU(3) plasma: An update, Phys. Rev. D 91 (9) (2015) 096002. arXiv:1503.05652 , doi:10.1103/PhysRevD.91.096002 .\n- [38] Y. Fu, J. Ghiglieri, S. Iqbal, A. Kurkela, Thermalization of non-Abelian gauge theories at next-to-leading order, Phys. Rev. D 105 (5) (2022) 054031. arXiv:2110.01540 , doi:10.1103/PhysRevD.105.054031 .\n- [39] L. Giusti, M. Pepe, Equation of state of the SU(3) Yang-Mills theory: A precise determination from a moving frame, Phys. Lett. B 769 (2017) 385-390. arXiv:1612.00265 , doi:10.1016/j.physletb. 2017.04.001 .\n- [40] M. Garny, M. Sandora, M. S. Sloth, Planckian Interacting Massive Particles as Dark Matter, Phys. Rev. Lett. 116 (10) (2016) 101302. arXiv: 1511.03278 , doi:10.1103/PhysRevLett.116.101302 .\n- [41] Y. Tang, Y.-L. Wu, Pure Gravitational Dark Matter, Its Mass and Signatures, Phys. Lett. B 758 (2016) 402-406. arXiv:1604.04701 , doi: 10.1016/j.physletb.2016.05.045 .\n- [42] Y. Tang, Y.-L. Wu, On Thermal Gravitational Contribution to Particle Production and Dark Matter, Phys. Lett. B 774 (2017) 676-681. arXiv: 1708.05138 , doi:10.1016/j.physletb.2017.10.034 .\n- [43] C. J. Morningstar, M. J. Peardon, The Glueball spectrum from an anisotropic lattice study, Phys. Rev. D 60 (1999) 034509. arXiv: hep-lat/9901004 , doi:10.1103/PhysRevD.60.034509 .\n- [44] Y. Chen, et al., Glueball spectrum and matrix elements on anisotropic lattices, Phys. Rev. D 73 (2006) 014516. arXiv:hep-lat/0510074 , doi:10.1103/PhysRevD.73.014516 .\n- [45] V. Mathieu, N. Kochelev, V. Vento, The Physics of Glueballs, Int. J. Mod. Phys. E 18 (2009) 1-49. arXiv:0810.4453 , doi:10.1142/', 'S0218301309012124 .': "- [46] S. Aoki, et al., Review of Lattice Results Concerning Low-Energy Particle Physics, Eur. Phys. J. C 74 (2014) 2890. arXiv:1310.8555 , doi:10. 1140/epjc/s10052-014-2890-7 .\n- [47] H. B. Meyer, Glueball matrix elements: A Lattice calculation and applications, JHEP 01 (2009) 071. arXiv:0808.3151 , doi:10.1088/ 1126-6708/2009/01/071 .\n- [48] A. R. Liddle, P. Parsons, J. D. Barrow, Formalizing the slow roll approximation in inflation, Phys. Rev. D 50 (1994) 7222-7232. arXiv: astro-ph/9408015 , doi:10.1103/PhysRevD.50.7222 .\n- [49] M. Laine, A. Vuorinen, Y. Zhu, Next-to-leading order thermal spectral functions in the perturbative domain, JHEP 09 (2011) 084. arXiv:1108. 1259 , doi:10.1007/JHEP09(2011)084 .\n- [50] W. DeRocco, P. W. Graham, S. Kalia, Warming up cold inflation, JCAP 11 (2021) 011. arXiv:2107.07517 , doi:10.1088/1475-7516/ 2021/11/011 .\n- [51] H. B. Meyer, High-Precision Thermodynamics and Hagedorn Density of States, Phys. Rev. D 80 (2009) 051502. arXiv:0905.4229 , doi:10. 1103/PhysRevD.80.051502 .\n- [52] A. Andronic, P. Braun-Munzinger, K. Redlich, J. Stachel, Decoding the phase structure of QCD via particle production at high energy, Nature 561 (7723) (2018) 321-330. arXiv:1710.09425 , doi:10.1038/ s41586-018-0491-6 .\n- [53] A. Batz, T. Cohen, D. Curtin, C. Gemmell, G. D. Kribs, Dark sector glueballs at the LHC, JHEP 04 (2024) 070. arXiv:2310.13731 , doi:10.1007/JHEP04(2024)070 .\n- [54] D. Eichler, TeV particles as weakly unstable dark matter, Phys. Rev. Lett. 63 (1989) 2440. doi:10.1103/PhysRevLett.63.2440 .\n- [55] A. Ibarra, D. Tran, C. Weniger, Indirect Searches for Decaying Dark Matter, Int. J. Mod. Phys. A 28 (2013) 1330040. arXiv:1307.6434 , doi:10.1142/S0217751X13300408 .\n- [56] O. K. Kalashev, M. Y. Kuznetsov, Constraining heavy decaying dark matter with the high energy gamma-ray limits, Phys. Rev. D 94 (6) (2016) 063535. arXiv:1606.07354 , doi:10.1103/PhysRevD.94.063535 .\n- [57] T. Cohen, K. Murase, N. L. Rodd, B. R. Safdi, Y. Soreq, γ -ray Constraints on Decaying Dark Matter and Implications for IceCube, Phys. Rev. Lett. 119 (2) (2017) 021102. arXiv:1612.05638 , doi:10.1103/ PhysRevLett.119.021102 .\n- [58] http: // www.laine.itp.unibe.ch / eos15 / .\n- [59] R. Lau, M. Teper, SO(N) gauge theories in 2 + 1 dimensions: glueball spectra and confinement, JHEP 10 (2017) 022. arXiv:1701.06941 , doi:10.1007/JHEP10(2017)022 .\n- [60] M. Teper, Pfa ffi an particles and strings in SO(2N) gauge theories, JHEP 01 (2019) 136. arXiv:1810.04546 , doi:10.1007/JHEP01(2019) 136 .\n- [61] J. Halverson, C. Long, A. Maiti, B. Nelson, G. Salinas, Gravitational waves from dark Yang-Mills sectors, JHEP 05 (2021) 154. arXiv:2012. 04071 , doi:10.1007/JHEP05(2021)154 .\n- [62] W.-C. Huang, M. Reichert, F. Sannino, Z.-W. Wang, Testing the dark SU(N) Yang-Mills theory confined landscape: From the lattice to gravitational waves, Phys. Rev. D 104 (3) (2021) 035005. arXiv:2012.11614 , doi:10.1103/PhysRevD.104.035005 .\n- [63] E. Morgante, N. Ramberg, P. Schwaller, Gravitational waves from dark SU(3) Yang-Mills theory, Phys. Rev. D 107 (3) (2023) 036010. arXiv: 2210.11821 , doi:10.1103/PhysRevD.107.036010 .\n- [64] M. Yamada, K. Yonekura, Cosmic strings from pure Yang-Mills theory, Phys. Rev. D 106 (12) (2022) 123515. arXiv:2204.13123 , doi:10. 1103/PhysRevD.106.123515 .\n- [65] K. Holland, M. Pepe, U. J. Wiese, The Deconfinement phase transition of Sp(2) and Sp(3) Yang-Mills theories in (2 + 1)-dimensions and (3 + 1)-dimensions, Nucl. Phys. B 694 (2004) 35-58. arXiv:hep-lat/ 0312022 , doi:10.1016/j.nuclphysb.2004.06.026 .\n- [66] B. Lucini, M. Teper, U. Wenger, Properties of the deconfining phase transition in SU(N) gauge theories, JHEP 02 (2005) 033. arXiv:hep-lat/ 0502003 , doi:10.1088/1126-6708/2005/02/033 .\n- [67] H. Kolesova, M. Laine, Update on gravitational wave signals from postinflationary phase transitions, Phys. Lett. B 851 (2024) 138553. arXiv: 2311.03718 , doi:10.1016/j.physletb.2024.138553 .\n- [68] K. V. Berghaus, M. Forslund, M. V. Guevarra, Minimal warm inflation with a heavy QCD axion (2 2024). arXiv:2402.13535 .\n- [69] D. McKeen, R. Mizuta, D. E. Morrissey, M. Shamma, Dark Matter from  \nDark Glueball Dominance (6 2024). arXiv:2406.18635 .  \n- [70] A. Manohar, H. Georgi, Chiral Quarks and the Nonrelativistic Quark Model, Nucl. Phys. B 234 (1984) 189-212. doi:10.1016/ 0550-3213(84)90231-1 .\n- [71] A. G. Cohen, D. B. Kaplan, A. E. Nelson, Counting 4 π 's in strongly coupled supersymmetry, Phys. Lett. B 412 (1997) 301-308. arXiv: hep-ph/9706275 , doi:10.1016/S0370-2693(97)00995-7 .\n- [72] M. A. Luty, Naive dimensional analysis and supersymmetry, Phys. Rev. D 57 (1998) 1531-1538. arXiv:hep-ph/9706235 , doi:10.1103/ PhysRevD.57.1531 .\n- [73] G. 't Hooft, A Planar Diagram Theory for Strong Interactions, Nucl. Phys. B 72 (1974) 461. doi:10.1016/0550-3213(74)90154-0 .\n- [74] E. Witten, Baryons in the 1 / n Expansion, Nucl. Phys. B 160 (1979) 57115. doi:10.1016/0550-3213(79)90232-3 ."}