Smooth reheating and dark matter via non-Abelian gauge theory

We demonstrate how a dark sector can provide a crucial link between a vacuum-dominated early universe and the later radiation-dominated 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.


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][2][3].Certainly, some mechanism had to be active and efficient 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 efficiently constrain inflationary models [5][6][7], as well as dark sectors [8][9][10].
It 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][12][13][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][16][17][18][19][20][21][22][23][24][25][26][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) Email addresses: simone.biondini@unibas.ch(S.Biondini), helena.kolesova@uis.no(H.Kolešová), simona.procacci@unige.ch(S.Procacci) epoch.In this work, we restrict to an inflaton solely coupled to a dark SU(N) gauge theory.
As 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][30][31][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 sufficiently hot SM bath could be realized.
The 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.
Glueball states are often classified in terms of their quantum numbers under angular momentum (J), parity (P), and charge conjugation (C), which ultimately determine the effective 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-date [22].This is the option that will be explored in our work.
We 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.
The 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 offer some discussion and concluding remarks in section 5.

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.

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 The 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 and fix decay constant, f a , 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.
The 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 suffices 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. 1 The values we use are f a ≈ 1.25 m pl , m φ ≈ 1.09×10 −6 m pl .This potential is known as natural since it might be generated in a UV completion of the theory.

Dark vectors
In the deconfined phase, the non-Abelian dark sector is expected to thermalize2 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.
The 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]. 4A 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 effective description of the interactions can be employed and leading operators have mass dimension six and eight [34,35].
In particular, the dimension-six operator that involves the SM Higgs doublet, H, reads where c 6 is the matching coefficient, while color and Lorentz indices are suppressed.Whenever M ≪ m pl the corresponding effective operators are expected to be more efficient 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, The 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 → SM SM, are suppressed as (v h /T DS ) 2 ≪ 1. Dimension-eight operators are instead suppressed by (T DS /M) 4 .
We verify a posteriori that the portal interactions are never efficient 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 × 106 GeV ≪ T DS ≈ 8.4 × 10 10 GeV.

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][44][45].Since the minimal YM action respects angular momentum, charge, and parity, glueballs are systematized in terms of the relevant quantum numbers J PC .
We present 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. 6This 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.
Dark glueballs interact with the SM particles via the same set of effective 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][44][45].Moreover, the qualitative difference 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 ++ ≫ m h ,7 Although 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 x 2 − z 2 extracted as in footnote 5. Conversely, the lightest C-odd state can only decay via dimension-8 operators of the schematic form [34] where B µν is the SM hypercharge field strength and c 8 is again a matching coefficient.This operator induces the decays 1 +− → 0 ++ γ and 1 +− → 0 ++ Z, with widths scaling as Λ9 DS /M8 .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 ++ ≫ m Z approximation, the full decay width is where 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 large-N limit suggests [22].In 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 large-N expansion (see Appendix A for details).

Evolution equations
In order to track the evolution of the different species, we rely on a network of evolution equations, that comprise the interaction rates introduced in sec.2, and the Hubble rate, We use the approximate initial value H 2 ref ≡ 8πV(φ 0 )/(3m 2 pl ) as an inverse-time unit. 8The 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.

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. 9If 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, ėφ + 3He φ = −Υe φ (9) ėDS + 3H(e DS + p DS ) = Υe φ − Γ e DS (10) ėSM The friction coefficient Υ 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].As 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 efficiently populated.Relying on their fast thermalization [38,50], the SM is in turn produced from frequent 2→2 scatterings of vectors à la freeze-in.The corresponding rate is Γ = Γ G as given in eq. ( 4).
During 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 different decay rates to SM for the two phases.
Once 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 where i runs over different glueball states.The modeling in eq. ( 12) is supported by lattice simulation of a pure SU(3) gauge theory [51]. 11Eqs.( 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).
In the confined phase, the inflaton can still lose energy in the form of dark vectors, which hadronize quickly and effectively 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 affect the freeze-out dynamics.
At 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.

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. 10For SU (3) we have m 0 ++ ∼ 5.4T c DS , hence the non-relativistic approximation is not appropriate shortly after the phase transition.However, we verify the negligible effect of this error on our conclusions, mainly because the freeze out of stable glueballs happens at lower temperatures T/T c ∼ 0.5, see Fig. 2. 11 More specifically, for T DS /T c DS ≲ 0.7, the lattice results agree very well with a system modeled by the two lightest glueballs 0 ++ and 2 ++ .
As the dark sector is made of glueballs with different masses and interactions, the composite states fall out of chemical equilibrium at different 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.
The relative abundance of different 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.

Numerical solutions and constraints
The main result of our study is a scan of numerical solutions for different values of the parameters Λ DS and M.An example of full solution for the energy densities of different sectors throughout the different stages described in sec.3 is shown in 12 To estimate this threshold value Λ φ DS , we observe that, when the inflaton decays, the Hubble rate satisfies H ∼ Υ.Moreover, the energy density of the dark sector around the phase transition can be approximated as e DS ∼ 0.1Λ 4  DS .If we assume that e DS dominates the Hubble rate we find Λ   fig. 1. 13 We rescale dimensionful quantities by powers of H ref , defined below eq.( 8). 13 For the parameters choice Λ DS ≈ 1.5 × 10 6 GeV, and M ≈ 8.5 × 10 13 GeV we obtain T rh SM ≈ 4.3 MeV, and Ω DM h 2 ≈ 0.11.
The 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 → SM SM 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.
When 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.

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 sufficiently 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 ) Eq. ( 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 different 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 affected also by the uncertainties in the non-perturbative input and matching coefficient, c 8 , as presented in eq. ( 7).Keeping these opposite effects in mind, in this work we resort ourselves to In 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 /(3m 2 pl ), we find an analytical estimate for the minimal decay rate required by eq. ( 16), Combining eqs.( 18) and ( 5) sets a constraint on a further combination of Λ DS and M. For the results illustrated in fig.3, we extract 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.
Finally, 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).

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 efficient, 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 inefficient.
Since 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.
The 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 different for the three-glueballs system (solid lines) and for the two-glueballs system (dashed lines).
The 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 inefficient 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) differ considerably for the 2-glueball and 3-glueball scenarios as expected based on the discussion of the interaction strengths above.
Our 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 efficient, and radiation domination appears later, leading to more dilution of the DM density. 14Finally, exclusion limits from indirect detection (blue), c.f. eq. ( 15), are identical for both cases , but affected by the systematic uncertainty discussed above (see the difference between the dotted and solid lines showing the possibility of enlarging the white region).
Let us stress the large impact of adding the next-to-leading glueball contribution to the dark energy density e DS .A realistic Figure 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).determination 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 effect 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.

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 non-Abelian 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 different 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 the dark matter relic density, as well as the SM reheating temperature.
Our 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 different gauge groups, e.g.SO(N) with N ≥ 8 [25,27,35].The decay rate of stable glueballs scales then as Λ 2N−3 DS /M 2N−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 offer further observational constraints, as discussed in [33,[61][62][63][64].In particular, the confinement phase transition in different 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].
Other 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 coefficient, Υ, 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 affect the freeze-out dynamics.
Finally, 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.
Note Added: During the revision of our manuscript the study of ref. [69] appeared that investigates a closely related scenario of dark glueballs.

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 ṅ1 +− + 3Hn 1 +− = B The 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 different 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

Figure 1 :
Figure 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.

Figure 2 :
Figure 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.