Slightly Ultra-violet Freeze-in a Hidden Gluonic Sector

The dark glueball (DGB) from a hidden Yang-Mills sector is a simple non-WIMP dark matter candidate characterized by very few parameters. However, it suffers the over dense issue. To overcome it, in general the dark sector is required to be hierarchically cooler than the visible sector. To naturally generate the desired hierarchy, in this paper we introduce higher dimensional operators coupling the dark gauge field strength tensor to the standard model (SM) Higgs doublets or gauge field strength tensors. By tracking the different phases of the universe from the end of inflation, prethermalization, reheating to the radiation dominance era, we show that these operators can make the DGB be a viable dark matter candidate over a wide mass region, from the sub-GeV to multi-PeV or even beyond. At the same time, these operators open decay channels of DGB to the SM species, and partial of the parameter space could leave hints in the cosmic ray.


I. INTRODUCTION
It becomes more and more channeling to naturally reconcile the weakly interacting massive particle (WIMP) paradigm of dark matter (DM) and a bunch of strong exclusions limits on the strength of DM-SM (standard model) interactions. The dark glueballs (DGBs), predicted by a hidden SU(N d ) pure Yang-Mills theory in the confining phase, may furnish one of the simplest non-WIMP DM candidate; its nature is specified by very few parameters. Not a WIMP, it can easily resolve the puzzle of null DM searches. In addition to that, the SU(N d ) non-Abelian gauge sector is well expected in various contexts of new physics beyond the SM such as dark matter/radiation [1][2][3], origin of baryon asymmetry [4], naturalness problem [5], and especially in the string theory [6].
One attractive merit of the DGB DM is that its mass is generated dynamically. There is a tower of DGBs classified by their quantum numbers labeled as J P C , and they are absolutely stable without a connection to the SM. The one lying at the bottom, namely J P C = 0 ++ , is supposed to be the dominant DM components. We will concentrate on this state throughout this paper. The mass of DGBs are determined by the confinement scale Λ d . According to the lattice simulation without introducing the topological θ term [7], the mass of 0 ++ is where α and β are order one parameters. In the large N d limit, the correction from color number is suppressed by 1/N 2 d and therefore the mass approaches a constant universal to all SU(N d ). From the well-studied N d = 3 case, α ≈ 7 [7]. One may expect that the mass spectra of N d > 3 assembles that of SU (3).
How to get the correct relic density of the dark glueball DM is an issue. Usually, it is assumed that the dark gluons form a thermal bath at the early universe, experiencing the highest temperature T max far above the confining scale Λ d . When T cools down below Λ d , the confining phase transition happens and the dark gluons turn into the non-relativistic dark glueballs. In terms of the effective theory for 0 ++ , its number density freezes-out after the 3 → 2 process incapable of catching the Hubble expansion rate. It is found that the DGB overcloses the universe, except for T (Λ d ) ≪ T (Λ d ), namely a much cooler dark sector than the SM sector. The above conclusion holds no matter whether the two sectors used to reach thermal equilibrium or not.
In order to produce DGBs in the early universe, we should introduce proper interactions for the SU(N d ) dark sector. A popular way to link it with the SM sector is adding higher dimensional operators in the form of O SM O d , suppressed by very high cut-off scale in order to guarantee a sufficiently long lifetime of 0 ++ . In this setup, the dark gluons can be slowly produced through the SM species scatterings with rates which are very suppressed by the cut-off scale appearing in the higher dimensional operators. Consequently the dark gluons fail to keep in thermal equilibrium with the SM sector, and thereby, as desired, they naturally have a much lower temperature than the SM bath. This way of producing the dark sector is nothing but an application of the UV freeze-in mechanism [8,9]. Different than the usual freeze-in via renormalizable operators which is not senstive to UV [8,10,13,14], in UV freeze-in special attention should be paid to the very early universe before reheat, when the universe went through eras with an even higher temperature than the reheat temperature [11,12,[15][16][17][18].
The goal of this paper is to detailedly investigate if the d = 6 and d = 8 operators could appropriately reheat the dark gluonic sector via the UV freeze-in mechanism. We take into account the contributions from different phases of the universe: 1) Phase I just after inflation during which the radiation from inflaton decay has not reached thermalizaiton yet; DM production in this phase receives attention just very recently [17]. 2) Phase II after thermalization of radiation but before the completion of reheat. 3) Phase III, the usual radiation era after reheat. We find that the contribution from phase I can be the dominant one in the d = 8 case, provided that a very heavy inflaton decaying slowly, for instance via Planck suppressed operators. Whereas in the d = 6 case dark gluons tend to be thermally produced.
The paper is organized as the following: In Section II we present the effective model of DGB and introduce higher dimensional operators bridging the dark and visible sectors, and we also present the resulting decay widths of DGB. In Section III we focus on how the higher dimensional operators slightly reheat the dark sector. Section IV contains the conclusions and discussions.

II. A DECAYING DARK GLUEBALL
A. The self-interaction & relic density issues of DGB dark matter At the energy scale below the scale Λ d where the SU(N d ) gauge coupling becomes strong, the degrees of freedom of the theory are DGBs instead of gluons. But how to construct the low energy effective theory for these new ingredients is out of the comfort zone of QFT. The following effective model for the lowest lying state 0 ++ ≡ s [19] is proposed in the spirit of non-perturbative methods, where c i are expected to be order one numbers in the light of the large-N limit and naive dimension assumption (NDA). However, it is definitely not beyond dispute. Different viewpoints lead to significant differences in the estimation of couplings. In Ref. [20] the estimation is just based on NDA without considering the large-N limit, then the 1/N d suppression factor is absent. On the contrary, in Ref. [21] which just follows the large-N limit, the 4π factor is not presented. Thus, we may keep an open mind on the strength of these couplings. In this paper we follow Eq. (2.4) unless explicitly stated otherwise. The values of effective couplings have direct implications to the self-interaction and relic density of the dark matter candidate s. First, the self-interacting term s 4 gives rise to the elastic scattering ss → ss, which has a large cross section σ 2→2 ∼ (4π/N d ) 4 /m 2 s . The Bullet Cluster imposes an upper bound on σ 2→2 /m s 1cm 2 /g [22], which in turn gives the lower bound on the DGB mass [19,23], In the large-N limit, DGB is supposed to be weakly interacting and thus could relieve this bound, but not much given a normally large N d . In addition to leave hints in the large scale structure today, this fast elastic scatterings ss → ss could keep the DGBs in thermal equilibrium within themselves just after the dark confining phase transition. Second, the s 5 term leads to a sizable particle number depletion process 3s → 2s, which could maintain the chemical equilibrium among the DGBs. After the 3s → 2s depletion rate falls below the Hubble expansion rate at T d , the chemical equilibrium is lost and the DGB number density freezes out. As a consequence of entropy conservation in both sectors, the relic density is determined to be [19,23,24], where g S and g d S are the entropy degrees of freedom in the visible and dark sectors densities, respectively; ξ T ≡ T /T is the ratio between the temperatures of the two sectors. All of them are calculated at a temperature T i Λ d , as the initial conditions for the DGB freeze-out era. And they are constant up to any high temperature until entropy conservation is violated.
It is expected that T d Λ d . In particular, T d ≃ Λ d in the limit c 5 → 0 [20,25]. Otherwise, one may parameterize the deviation by T d = w f Λ d with w f 1. Anyway, from Eq. (2.3) we immediately see that, as stated in the introduction, ξ T ≪ 1 is indispensable to reduce the large ratio T d /(3.6eV) ≃ Λ d /(3.6eV) 10 8 indicated by Eq. (2.2). The main goal of this paper is to specify the UV freeze-in as the natural mechanism for generating ξ T ≪ 1 .

B. Higher dimensional operators and the decaying dark glueball
Although it is possible that the dark SU(N d ) Yang-Mills sector thus dark glueballs interact with the SM sector purely gravitationally, many people introduce messengers communicating interactions between the two sectors. These messengers are supposed to be charged under SU(N d ) and at the same time communicating with the SM particles via gauge or Yukawa gauge interactions. Moreover, they are very heavy, having masses m Q ≫ Λ d . Alternatively, the messengers may be a heavy moduli-like field which is neutral under any gauge groups [26]. After integrating out the messengers one obtains the higher dimensional operators O SM O d that describe the interactions between the dark glueball and SM species. The concrete effective Lagrangian is model dependent, on the choice of messengers. In this article, to demonstrate the idea, we just consider two representative examples at the d = 6 and d = 8 level, respectively, where M 6,8 are the cut-off scales, which in a concrete UV model are the combinations of the heavy messenger mass scales and various couplings, but for later use here we leave g 2 d with g d the dark QCD gauge coupling. The dark QCD gluon field strength tensor is written as F µν = F a µν T a , with T a the generators of SU(N d ) satisfying the normalization conditions tr(T a T b ) = δ ab /2. The same convention applies to other non-Abelian gauge groups. We only consider the scalar operator S ≡ trF µν F µν , and likely other operators accompany, e.g., trF µν F µν trG µν G µν and trF µα F ν α trG α µ G αν , and so on [27,28]. But their presence does not change the main line of our discussion, though a possible sizable numerical correction [30].
Several comments on UV models are in orders. First, the Higgs-portal operator selects a class of UV complete models, for instance those with scalar messengers which are neutral under the SM gauge groups but have Higgs portal interactions [29]. Otherwise, one may have to arrange the proper Yukawa interactions between the fermionic messengers and the Higgs doublet; see Ref. [28]. Second, at the d = 8 level if the messengers are color neutral, the portal does not contain gluons; on the contrary, the portal might be pure gluon-portal provided that the messengers do not carry electroweak quantum numbers. In the context of grand unification, it is expected that all kinds of SM vector bosons are present. (2.5) In this matrix element, the SM part can be calculated perturbatively, but the other part must fall back on non-perturbative methods and is parameterized as F O d J CP , known as the decay constant of the state J CP . Hereafter we will specify J CP = 0 ++ ≡ s and take F S 0 ++ as F s . From the lattice result [31], one may take For a given O d , the accessible annihilation decay modes of s depend on its mass. In the following we will discuss the decay pattern of DGB separately in the d = 6 and d = 8 cases.

C. Higgs-portal DGB (d = 6)
For the Higgs portal case, the decay pattern of s is well-understood as an additional Higgs state mixing with the SM Higgs boson, which for a wide mass region of m h has been widely studied [32]. The partial decay width of s into a pair of SM states is written as where the factor in the bracket could be identified as the mixing angle between h and s; v = 246 GeV and m h = 125 GeV is the SM-like Higgs mass. If the DGB is heavy having a mass m s = m 0 ≫ m W , then s dominantly decays into the vector bosons and the resulting lifetime of s is estimated to be In order to evade the constraints on decaying DM into vector bosons [33], most stringently by the diffuse gamma ray searches, the lifetime must be much longer than the age of the universe. For the decaying DGB DM around the TeV scale, it suggests that the interactions between two sectors actually must be suppressed by the near-Planck scale.
In the opposite limit, the DGB is very light and lies much below the weak scale but above 100 MeV due to the self-interaction bound Eq. (2.2). Then, s dominantly decays into the heaviest SM fermion pair ff kinematically accessible. Probably the most suppressed case 1 Transitions between dark blueballs are at the radiative level and thus may play significant roles in the models where the DGBs decay very fast without protection by quantum numbers. But in our scope, where J CP = 0 ++ is the dominant DM component, the leading order consideration is sufficient.
of DGB decay is that it dominantly decays into a pair of muons for 1GeV m s ≫ 2m µ . In this case one has Even if the region 2m e ≪ m s 2m µ is marginally realized after taking into account the uncertainty of the lower bound on m s , the above estimation will not be dramatically changed. However, if we do not insist on this bound and allow a much lighter DGB then the scale of M 6 can be fairly low. Now since m s ≪ 2m e , only a pair of photon is accessible in s decay and the resulting DGB lifetime is It is seen that in this case the cutoff scale can be as low as the interesting TeV scale. Therefore, the lower bound on Λ d by virtue of self-interaction places a strong lower bound on M 6 .

D. Vector-boson-portal DGB (d = 8)
Now we move to the next case, DGBs communicating with the SM sector via the vectorboson-portal at the dimension-8 level. s can annihilate decay into a pair of vector bosons such as a pair of free QCD gluons, and the decay width is [28] Γ s→gg = 1 2π Compared to the d = 6 case, it has a much stronger dependence on the confining scale Λ d , which then leads to a much longer lifetime of DGB. Note that if the DGB mass is below the GeV scale, the hadronic mode S → gg is closed today. The lifetime of the decaying DGB in terms of Eq. (2.11) is estimated to be Lowering Λ d down to the sub-GeV scale, M 8 can be down to the PeV scale. In our analysis, we focus on the DGB mass regions which admit a simple analytical expressions of decay width,that is convenient to our final numerical display.

III. HIGHER DIMENSIONAL OPERATORS: REHEAT THE DARK SECTOR
The UV freeze-in mechanism is sensitive to the scattering process characterized by hard momentum, so we have to backdate the universe to the extremely hot state just after inflation. It is usually known as the reheating stage and gives the broadly quoted "highest" temperature of the universe T re ∼ Γ φ M Pl , with Γ φ the width of the inflaton perturbative decay to radiation which thermalizes instantaneously at t re ≃ 1/Γ φ . But there is an even hotter phase than T re during reheating [34], whose possible impacts on UV freeze-in are investigated by several groups [11,12,15], to find that the impacts are small given d < 9.
These studies assume that the radiation is always instantaneously thermalized during reheating. However, instantaneous thermalization may be a bad approximation. Actually, during reheating there exists a phase prior to the thermal equilibrium of radiation, namely prethermalization when the radiation has not reached full thermal equilibrium yet. Radiation in this phase is dominated by the extremely hard primary modes with momentum around the inflaton mass m φ , and thus the UV freeze-in process, whose cross section is greatly enhanced, may be very effective in this phase. Ref. [17] carefully studied the production of FIMP DM during the prethermlization era and found that it indeed could be the dominant contribution to the final FIMP relic density.
Before proceeding, we briefly review the thermalization process studied in Ref. [35]. The inflaton decay width is parameterized as Γ φ = km 2 φ /M Pl , from which one can see that k 1 leads to T re m φ and thus thermal effects on the inflaton decay become important [36,37]. To avoid this complication, in this paper we just focus on the slow decay limit with k ≪ 1, for instance, inflaton decay via a Plank suppressed dimension five operator. Furthermore, it is reasonable to consider that the radiation is charged under a non-Abelian gauge group with normal gauge coupling α ∼ O(0.01). The energy of primary hard radiation could be effectively dissipated away via the fast collinear soft gauge boson emissions, eventually reaching thermalization. It is shown that as long as thermalization is completed at a time scale t th well before the reheating time scale t re .
Concretely, the thermalization time scale is determined by with t end denoting the end of inflation.
In what follows, we will follow the thermal history of the universe, to describe how UV freeze-in via the higher dimensional operators just provides a good way to slightly reheat the dark gluonic sector.
A. UV freeze-in dark gluons

The BEs for the inflaton-radiation-dark gluon system
The dark gluons are radiation, so it is more useful to derive BE for the energy density instead of number density. During the reheating region, the BEs for the inflaton-radiationdark gluon system are d dt where the collision term γ d will be specified later. The Hubble parameter is determined by the Friedmann equation: In this set of BE describing the energy density transfer, inflaton dominates the system throughout the reheating stage thus H ≈ ρ φ /3M 2 Pl . Its energy density is simply red shifting, and therefore H ∝ a −3/2 or H(a) = H end (a end /a) 3/2 , where H end = ρ end /3M 2 Pl with ρ end the energy density of inflaton at the end of inflation. In a class of popular inflation models ρ end is estimated to be From the Friedmann equation one can estimate the age of the universe t = a(t) The last relation holds for a ≫ a end . This t − a scaling rule is different than the one in the radiation dominant era Eq. (3.9) because the reheating era is matter dominant.
Reheating is completed when the inflatons decay away at t re ≃ 1/Γ φ , after which the inflaton domination gives way to a period of radiation domination. Radiation energy density contains two components, the ordinary ρ R produced by inflaton decay and the dark radiation ρ g d produced by the ordinary radiation. Such a system is described by Eq. (3.4) and Eq. (3.5) with Γ φ → 0. 2 During this era ρ R dominates over ρ g d and H ≈ ρ R /3M 2 Pl ∝ a −2 . Then similar to the previous derivation one gets the t − a − H relation It is ready to show that in any phase dt/da = H −1 a −1 .
Simplifications to the BEs can be made. To scale out the effect of the universe expansion, we consider evolution of the dimensionless quantities Φ = ρ φ a 3 /m φ , R = ρ R a 4 and G d = ρ g d a 4 . We also introduce the dimensionless variableâ ≡ am φ . The radiation 3 is gradually built up from the inflaton perturbative decay φ → RR, while the dark sector is only slightly reheated by the radiation scattering R + R ′ → g d + g ′ d . These considerations motivate the no back reaction approximation: In the inflaton dominance era, the energy leaking to the radiation is negligible and thus in Eq. (3.10) the Γ φ ρ φ term can be dropped; similarly, in Eq. (3.4) the γ g d term is removed. Then, the BEs take the simple form during reheating, 2 It is assumed that the transformation from phase II to phase III is prompt. A more accurate treatment should take into account the decaying term of inflaton. where the prime denotes the derivative with respect toâ. While during the radiation dominated era, the BEs are reduced to In the no back reaction limit, it is ready to solve the BEs one by one. First, the solution to Eq. (3.10) is which is a constant as expected, because the inflaton energy density is (approximately) purely red-shifting as a −3 . To account for the perturbative decay of inflaton in Eq. (3.10), one may multiply the above solution by the decaying factor e −â . Next, with Eq. (3.15), the energy density of radiation background in the phase II can be obtained by directly integrating overâ, Whereas in the phase III, the solution to Eq. (3.13) is a trivial redshiting, and then in Eq. (3.14) R takes a constant value R re ≡ R(â re ). Finally, the energy density of dark gluons will be presented at the end of this subsection. We need the relation between a and T to solve the BEs. After the radiation reaching thermal equilibrium at t th , the visible sector temperature can be defined through ρ R = π 2 30 g re * T 4 with g re the degree of freedoms of the relativistic particles during reheating. Along with Eq. (3.16), one gets the scaling rule T ∝ a −3/8 during the reheating era after thermalization: with T th the thermalization temperature, corresponding to the time scale t th . It is also the maximum temperature of the thermalized universe, T max . From the above scaling rule and Eq. (3.2), Eq. (3.8), one can determine Its magnitude relative to m φ depends on the value of k. The referred value k ∼ k 0 ≡ m φ /M Pl is of special interest since it makes T re ≪ T max ≪ m φ , which implies that the thermal effects from radiation is fully under control. However, in general T max can exceed m φ as k ≫ k 0 and m φ in the relatively low mass region. We have to remind the readers that if T max > M 6,8 , again the on-shell production of mediators may be important. Another merit of using k 0 is the greatly narrowing the many possibilities in properly reheating the dark gluonic sector. Hereafter we will quote k 0 frequently in the following analysis. Maybe k 0 relates to the inflaton decay via the Planck scale suppressed d = 5 operator like φψψ/M Pl . After the inflaton decays away, entropy is conserved and thus temperature drops much faster, following the well-known rule T = T re a re a −1 .

The cross sections
We now move to the calculation of the collision terms in Eq. (3.5), which involve the cross sections for dark gluons prodcution. In our setup, the dark gluons are produced via the 2 → 2 processes. In general, the Lorentze invariant cross section of the process 1 + 2 → 3 + 4 is defined as [38] σ(s) = 1 where the spin summation for the final state and average of the initial states are implied. The prefactor is known as the Moller velocity. In the massless limit, it is 1/(2s). Because σ(s) is Lorenz invariant, we will calculate it in the CM frame of the incident partiles for simplicity. From the effective operators given in Eq. (2.4), the squared amplitudes for the two representative cases are given by with N c = 3 for QCD. The process H + H * → g d + g d happens before electroweak spontaneously breaking, so there is a factor 2 to account for the doublet. Then, the invariant cross sections are .

(3.24)
Averages of the initial state over the spin and internal degrees of freedom give rise to the factors 1/D 2 H = 1/2 2 and 1/D 2 g = 1/(2 2 (N 2 c −1) 2 ) for the d = 6 and d = 8 cases, respectively. Additionally, we have taken into account the symmetry factors 1/2 (for d = 6) and 1/2 2 (for d = 8). To guarantee the validation of the effective theory throughout the universe evolution after inflation, it is required that the heavy states integrated out should have masses ∼ M 6,8 ≫ m φ . Otherwise one should take into account the production of mediator [15].

The collision term
With the Lorenz invariant cross sections Eq. (3.23) and Eq. (3.24), the collision term in Eq. (3.5) can be calculated using the standard approach [38].
Let us start from the phase I, the prethermal phase. We need the distribution function of radiation, f R (E, t), which can be derived by such a fact: At a certain time t i , the radiation from inflaton decay at rest carries a fixed momentum m φ /2; due to the expansion rule Eq. (3.8), this spectrum is universally redshifted to the later time t, located at the momentum (m φ /2E) 3/2 . Using this argument, the spectrum of the radiation is derived to be [17,35] where n R is the number density of the radiation, obtained by means of directly counting the number of inflatons that have decayed away (assuming that one radiation is produced for per inflaton decay), where d φ ≡ Γ φ (t − t end ) ≪ 1 and n φ,end = ρ end /m φ , the inflaton number density at the end of inflation. With the above distribution function, one can calculate the collision term Stressed again, in the prethermal era there is no conceptual of temperature, and thus the collision term has no dependence on T . Moreover, those two kinds of operators demonstrate a common behavior as the scale factor increases, ∝ a −3 . The reason is traced back to the fact that the distribution function is the only source of a. After the thermalization is completed at t th , the distribution function takes the well known Bose-Einstein or Fermi-Dirac distributions, but in order to derive a simple analytical expression, both of them are approximated by the Maxiwell-Boltzman distribution. After thermalization, the collision term becomes with z = √ s/T , while σ(z) ∝ T 2 z 2 and T 6 z 6 for the d = 6 and d = 8 operators. We would like to pause to make a comment on the upper limit for z, which is simply set to infinity in the sense that √ s can be way larger than the given T . Actually, the dominant integrating range is near z ∼ 10. This treatment results in no T -dependence after integrating with respect to z, which just contributes a pure numerical factor, 3074 and 1.77 × 10 7 for d = 6 and d = 8, respectively; consequently, γ g d ∝ T 2d−3 . The larger d case benefits from more significant UV enhancement in the z ∼ 10 region, so it enjoys an impressive numerical enhancement. Note that these big numbers do not appear in γ I g d in phase I where temperature does not exist and √ s < m φ .
Eq. (3.30) holds as long as the plasma has been thermalized, both in the reheating and the radiation domination eras, which have different a−T relations Eq. (3.17) and Eq. (3.19), thus leading to γ II g d (â) and γ III g d (â), respectively. Concretely, in the two phases they are given by The higher dimensional operator leads to a higher negative power of a, rather than common. We are considering the evolution of energy density, so the collision terms gain one more power of T , compared with the collision terms in the BEs for number density. This is easily seen by using the T − a relation Eq. (3.17): For d = 8, γ III g d ∝ T 13 , whereas it is proportional to T 12 in Ref. [15].

The energy relic density of dark gluons
With all of those ingredients, we are now at the position to calculate the relic density of dark gluons. The comoving energy density of dark gluons atâ f >â re is obtained by integrating over the scaled scale factorâ 4 in the prethermalization, thermalization and radiation dominant era successively untilâ f , As a feature of UV-freeze in production, it is soon frozen as the scale factorâ f becomes significantly larger thanâ re . The frozen maybe happen even much earlier, and we will come back to this point later.
To reveal the mechanism of production, we present the contributions phase by phase. The integrals are trivial and has simple but lengthy analytical expressions: Phase I: In the phase I, dark gluons production is "IR" dominated both for the d = 6 and d = 8 operators (or any other d in the more general context), because their collision terms share the same power of a, −3. Taking a th ≫ a end , one has d = 6 : Thereby, increasing the dimension of operators leads to a suppressed production yield in the phase I. This contribution has a strong dependence on α, stemming from the dependence onâ th . But α will not appear in the contributions from other phases.
I; the error is still insignificant even if phase I is not important, because the contributions from phase II and III are comparable. Now we are capable of computing the ratio of the temperatures of the two sectors at the reheating time. The temperature of the dark gluon plasma at this time can be derived via Since after reheating both sectors keep cooling down, respectively following the rules T (a) = T re a re /a and Eq. (3.19), from which we immediately get T = Tre Tre T = ξ T T . As a natural consequence of our mechanism reheating the dark sector, ξ T ≪ 1, which just provides the desired hierarchy to accommodate the correct DGB relic density. Let us show the expression of ξ T for d = 6: Increasing m φ helps to enhance ξ T . The first and second terms in the square bracket originate from phase I and phase II plus phase III, respectively. To have a more direct impression on this ratio, let us analyze two limits.
• To the contrary, if m φ ≫ 10 13 GeV, then the phase I is the dominant contribution to reheat the dark sector and ξ T,6 ∼ k 5/12 α − 4 3 r φ,6 m φ M Pl

5/12
. On the other hand, in cosmology the inflation models having convex potentials give the largest inflaton mass, e.g., in the popular chaotic model m φ 10 15 GeV [39]. This means that there is barely room for this scenario. It is tempting to fix m φ = 10 15 GeV for this case, dubbed the convex limit. But will not expand discussions on this special case.
A similar analysis can be made for the case with d = 8, yielding where N c = 3 has been used. The first term denotes for the contribution from phase I, and it dominates the other contributions when m φ ≫ 0.4 × 10 13 GeV (α/0.01) 2 (100/g re ) 39/32 (k/1.0) 11/16 , almost coinciding with the counterpart in the d = 6 case. But here taking a very small k could substantially relax the condition. For instance, consider k = k 0 then the above condition is weakened to m φ 10 3 GeV, which is a loose condition. Compared to the d = 6 case, where it is hard to realize the phase I production of dark gluons, in the d = 8 case this phase tends to play a much more prominent role, in particular for a very slowly decaying inflaton (namely k ≪ 1). The cause is traced back to the enhanced UV-sensitive by virtue of the larger d. Choosing k = k 0 has one more reason, though not physics related: The discussions on prethermalization production of dark species is rare, and here we present one good example. Anyway, in the d = 8 case we will just focus on the scenario of prethermal production dark gluons and then ξ T,6 is estimated to be Taking k = k 0 , to make ξ T, 8 10 −4 the inflaton must be as heavy as around 10 13 GeV. But if we work in the intermediate region 1 ≫ k ≫ k 0 , the required m φ can be accordingly lighter.

Dark confining phase transition and relic density of the dark gluon ball
As T drops below Λ, dark confining phase transition happens, and the relativistic dark gluons are confined to dark glue balls. The energy density of dark gluons is assumed to be transferred to that of the lowest state of DGB, s, This is a reasonable assumption to estimate the DGB relic density. But the excited states may be also produced and moreover the DGBs still carry some momentum (the corresponding kinetic energy does not contribute to the final DM relic density) just after the phase transition, although they soon become nonrelativistic. Estimation of the resulting uncertainty is beyond the scope of this article. After the confining phase transition, the energy density of the nonrelativistic DGB scales as matter ∝ a −3 until today, when the visible sector temperature T 0 = 2.37 × 10 −13 GeV.
Hence, its fraction in the total energy budget Ω s h 2 ≡ ρs(a 0 ) where we have used the critical energy density ρ c = 8.1 × 10 −47 h 2 GeV 4 . Note that here we estimate the DGB relic density directly from energy conservation, but the precise relic density of DGB relies on the s 5 term in the low energy effective model of DGB, namely Eq. (2.4). The two approaches give the same scaling behavior Ω s h 2 ∝ ξ 3 T Λ d . As we have stressed, a small ξ T is a built-in feature of our way to reheat the dark sector and thus the relic density problem can be naturally solved. It is still of concern that if the DGB could leave observable via its decay. To that end, we demonstrate the parameter space We also require m φ < 10 15 GeV. Upper panels: The first (k = k 0 ) and second (k = 1000k 0 )panels are for the relatively heavier DGB much above the weak scale which dominantly decays into a pair of vector bosons, while the third (k = k 0 ) and fourth panel are for the relatively light DGB which dominantly decays into µμ. Lower panels: first (k = k 0 ) and second (k = 10k 0 ) for prethermal UV freeze-in; third (k = 10 3 k 0 ) and fourth (k = 10 4 k 0 ) for thermal UV freeze-in. The thick green line and dahsed red line denote the FERMI-LAT and IceCube constraints, respectively.
of DGB with correct relic density on the (m s , τ s ) plane in Fig. 1 • In general, for a given DGB mass, a larger r φ leads to a significantly shorter DGB lifetime. This is because, to maintain a constant ξ T , accordingly M 6,8 = m φ /r φ becomes smaller thus much faster DGB decay. In other words, the parameter space having a larger r φ demonstrate a more promising detect prospect. Partial parameter space which gives DGB lifetime significantly shorter than 10 28 s has already been ruled out. We quote the results from Ref. [33], displaying the FERMI-LAT gamma ray data [40] and IceCube neutrino data [41] constraints in thick green and dashed red lines.
• In the d = 6 namely Higgs portal case (upper panels), as explained before, it is hard to realize the scenario of dark gluon production mainly in the phase I, so all four panels In the heavy DGB region ∼ O(TeV), referring to the first and second panels, it is clearly seen that the larger r φ cases, in particular given a larger k, have been ruled out by FERMI-LAT already. Whereas the low mass region is difficult to probe, owing to both its soft final states and much longer lifetime.
• In the d = 8 case (lower panels), the decay rate is highly suppressed by (m s /M 8 ) 9 , while this suppression does not appear in the relic density where m s is replaced by the much higher scale m φ , so this case stands little chance of observation in the light DGB region. To get a larger decay rates, we merely focus on the region above PeV scale; increasing k could help to enhance the decay rates and thus some of the parameter space is accessible in the experiments which are sensitive to super energetic cosmic ray, for instance the gamma ray data (up tp 2TeV) from FERMI-LAT and high energy neutrino data from IceCube. We are considering the gluon-portal and its contributions to gamma ray and neutrinos require simulations. But if one replaces gluons with the electroweak gauge bosons, then one may interpret the data similar to the Higgs portal, to find that it is sensitive to lifetime 10 28 s for a PeV scale DGB [33].
Anyway, in this paper we just schematically show the features of the parameter space which may admit a future discovery, and the systematic discussions on the indirect detection bounds deserves a specific study elsewhere.

IV. CONCLUSION AND DISCUSSIONS
Many new physics give rise to a pure SU(N d ) gauge sector surviving at low energy, and dark glueball is a prediction of such a sector in the confining phase. DGB provides a simple non-WIMP DM candidate characterized by very few parameters. However, its correct relic density needs careful treatment because it in general leaves too many relics to overcome the universe. The way out this problem is that the gluonic sector should be much cooler than the visible sector after the two sectors are reheated. We pointed out that if the two sectors are linked via the higher dimensional operators, and the SM species freezes in the gluonic sector in the very early universe, a very cool gluonic sector is a natural consequence. For concreteness, two kinds operators are introduced: • At the d = 6 level the dark gauge field strength tensor couples to the SM Higgs doublet; • At the d = 8 level the dark gauge field strength tensor couples to the SM vector bosons.
We carefully studied the yields of dark gluons in the different phases of the universe, from the end of inflation, reheating in the prethermalization and thermalization stages to radiation dominance, to find that the temperature and the production of the dark gluonic sector is sensitive to the inflaton mass and decay width; they determine two key temperatures T max and T re . For instance, considering the well-motivated Planck suppressed decay of the inflaton, it is found it is difficult to sufficiently reheat the dark gluonic sector; moreover, it is difficult to realize prethermalization production of dark gluons (it is true even for the much faster inflaton decay). By contrast, in the d = 8 case prethermalization production tends to be the main production mechanism. Anyway, the DGB is a viable dark matter candidate over a wide mass region, from sub-GeV to multi-PeV. The higher dimensional operators at the same time open decay channels for the DGB, and some of the parameter space leave hints in the cosmic ray: The d = 6 case is hopeful to be detected in the TeV DGB mass region, while in the d = 8 case a PeV scale DGB with lifetime ∼ 10 28 s might be probed by FERMI-LAT and IceCube. Interestingly, the IceCube PeV events can be interpreted by a decaying DGB with electroweak gauge bosons portal at the d = 8 level.
There are some open questions. In this paper we are confined to the assumption that the messengers are very heavy, having mass much greater than m φ , and thus they can not be produced on-shell during reheat. But this is questionable if one has the interest in a low cut-off scale and thus the messengers can be produced even if the reheat temperature is low. In this case we have to carefully take into account the production of messengers, to see if their own could be thermalized. Another open question is the thermalization of the gluonic sector. If it never establishes thermal equilibrium, there will be no confining phase transition and the dark gluons will leave as dark radition. Actually, thermalization of a generic FIMP sector is also of interest because it may has cosmological implications. Last but not least, a detailed study on the indirect detection bounds in the wide DBG mass region should be done in the coming paper. The possible gravitationals wave associated with the confining phase transition may also furnish a window to probe this simple but fairly hidden DM candidate.