Slightly Ultra-violet Freeze-in a Hidden Gluonic Sector

TL;DR

This paper proposes a mechanism involving higher dimensional operators to produce a viable dark glueball dark matter candidate with a wide mass range, addressing over-density issues and potential cosmic ray signals.

Contribution

It introduces higher dimensional operators coupling the dark gauge sector to the Standard Model, naturally generating the necessary temperature hierarchy for dark glueballs to be viable dark matter.

Findings

Dark glueballs can be viable dark matter over a wide mass range.

Operators enable decay channels, potentially observable in cosmic rays.

The model addresses over-density issues in dark glueball scenarios.

Abstract

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.

Equations

$$m_0=\left(\alpha +\beta /N_d^2\right){\Lambda }_d,$$

$${\mathcal{L}}_{\mathrm{DBG}}=\frac{1}{2}{\left(\partial s\right)}^2-\frac{1}{2!}{m}_s^2{s}^2-\frac{c_3}{3!}\frac{4\pi }{N_d}{m}_s{s}^3-\frac{c_4}{4!}{\left(\frac{4\pi }{N_d}\right)}^2{s}^4-\frac{c_5}{5!}{\left(\frac{4\pi }{N_d}\right)}^3\frac{s^5}{m_s}-...$$

$$m_s\gtrsim 100\times (3/N_d{)}^{4/3}\mathrm{MeV}.$$

$${\Omega }_s{h}^2\simeq \frac{g_S^d{\xi }_T^3}{g_S}\frac{{\mathcal{T}}_d}{3.6\mathrm{eV}},$$

$${\mathcal{L}}_{int}=\frac{g_d^2}{M_6^2}\mathrm{tr}{\mathcal{F}}^{\mu \nu }{\mathcal{F}}_{\mu \nu }|H{|}^2+\frac{g_d^2}{M_8^4}\mathrm{tr}{\mathcal{F}}^{\mu \nu }{\mathcal{F}}_{\mu \nu }\mathrm{tr}{G}^{\mu \nu }{G}_{\mu \nu }$$

$$⟨\mathrm{SM}|{\mathcal{O}}_{\mathrm{SM}}|0⟩⟨0|{\mathcal{O}}_d|J^{CP}⟩.$$

$$g_d^2{F}_s\simeq 3.06m_0^3.$$

$${\Gamma }_{s\to \mathrm{SM}+\mathrm{SM}}={\left(\frac{v{g}_d^2{F}^s}{M^2(m_s^2-m_h^2)}\right)}^2\Gamma (h\to \mathrm{SM}+\mathrm{SM}){|}_{m_h\to m_s},$$

$${\tau }_s\simeq 1.1\times 10^{30}{\left(\frac{M_6}{10^{17}\mathrm{GeV}}\right)}^4{\left(\frac{100\mathrm{GeV}}{{\Lambda }_d}\right)}^{5}s.$$

$${\tau }_s\simeq 5.3\times 10^{27}{\left(\frac{M_6}{10^{10}\mathrm{GeV}}\right)}^4{\left(\frac{0.1\mathrm{GeV}}{{\Lambda }_d}\right)}^7{\left(\frac{0.1\mathrm{GeV}}{m_f}\right)}^{2}s.$$

$${\tau }_s\simeq 2.3\times 10^{26}{\left(\frac{M_6}{10^4\mathrm{GeV}}\right)}^4{\left(\frac{1\mathrm{MeV}}{{\Lambda }_d}\right)}^{9}s.$$

$${\Gamma }_{s\to gg}=\frac{1}{2\pi }\frac{1}{M_8^8}{m}_s^3(g_d^2{F}_s{)}^2.$$

$${\tau }_s\simeq 1.1\times 10^{27}{\left(\frac{M_8}{10^{8.5}\mathrm{GeV}}\right)}^8{\left(\frac{10\mathrm{GeV}}{{\Lambda }_d}\right)}^{9}s.$$

$$k{\left(m_{\varphi }/M_{\mathrm{Pl}}\right)}^4\ll {\alpha }^8,$$

$$t_{th}-t_{end}\simeq {\Gamma }_{\varphi }^{-1}{\alpha }^{-16/5}{\left(\frac{{\Gamma }_{\varphi }{m}_{\varphi }^2}{M_{\mathrm{Pl}}^3}\right)}^{2/5},$$

$$\frac{d}{dt}{\rho }_{\varphi }+3H{\rho }_{\varphi }+{\Gamma }_{\varphi }{\rho }_{\varphi }$$

$$\frac{d}{dt}{\rho }_R+4H{\rho }_R-{\Gamma }_{\varphi }{\rho }_{\varphi }+{\tilde{\gamma }}_{g_d}$$

$$\frac{d}{dt}{\rho }_{g_d}+4H{\rho }_{g_d}-{\tilde{\gamma }}_{g_d}$$

$$H\equiv \dot{a}/a=\sqrt{({\rho }_{\varphi }+{\rho }_R+{\rho }_{g_d})/3M_{\mathrm{Pl}}^2}.$$

$${\rho }_{end}\sim m_{\varphi }^2{M}_{\mathrm{Pl}}^2,$$

$$t=\frac{3}{2}\left({\left(\frac{a}{a_{end}}\right)}^{3/2}-1\right)H_{end}^{-1}\approx \frac{3}{2}{H}^{-1}.$$

$$t=\frac{1}{2}\left({\left(\frac{a}{a_{re}}\right)}^2-1\right)H_{re}^{-1}\approx \frac{1}{2}{H}^{-1}.$$

$${\Phi }^{\prime }$$

$$R^{\prime }$$

$$G_d^{\prime }$$

$$R^{\prime }$$

$$G_d^{\prime }$$

$$\Phi (\hat{a})\approx \frac{4}{3}\frac{M_{\mathrm{Pl}}^2}{m_{\varphi }^4}{\hat{a}}_{end}^3{H}_{end}^2,$$

$$R(\hat{a})\approx \frac{8}{15}\frac{{\Gamma }_{\varphi }{M}_{\mathrm{Pl}}^2{H}_{end}}{m_{\varphi }^4}{\hat{a}}_{end}^{3/2}{\hat{a}}^{5/2}.$$

$$T=T_{th}{\hat{a}}_{th}^{3/8}{\hat{a}}^{-3/8},$$

$$T_{max}=T_{th}\simeq {\alpha }^{4/5}{m}_{\varphi }{\left(\frac{24}{{\pi }^2{g}_{re}}\right)}^{1/4}{\left(\frac{{\Gamma }_{\varphi }{M}_{\mathrm{Pl}}^2}{m_{\varphi }^3}\right)}^{2/5}.$$

$$T=T_{re}{a}_{re}{a}^{-1}.$$

$$\sigma (s)=\frac{1}{4\sqrt{(p_1\cdot p_2{)}^2-m_1^2{m}_2^2}}\int \frac{d^3{p}_3}{(2\pi {)}^{3}2E_3}\frac{d^3{p}_4}{(2\pi {)}^{3}2E_4}(2\pi {)}^4{\delta }^4(p_1+p_2-p_3-p_4)|\mathcal{M}{|}^2,$$

Full text

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 DM1 ; DM2 ; Ko:2017uyb , origin of baryon asymmetry baryon , naturalness problem twin , and especially in the string theory string .

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^{PC}$, and they are absolutely stable without a connection to the SM. The one lying at the bottom, namely $J^{PC}=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 $\Lambda_{d}$. According to the lattice simulation without introducing the topological $\theta$ term lattice , the mass of $0^{++}$ is

[TABLE]

where $\alpha$ and $\beta$ 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, $\alpha\approx 7$ lattice . 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 ${\cal T}_{max}$ far above the confining scale $\Lambda_{d}$. When ${\cal T}$ cools down below $\Lambda_{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\rightarrow 2$ process incapable of catching the Hubble expansion rate. It is found that the DGB overcloses the universe, except for ${\cal T}(\Lambda_{d})\ll T(\Lambda_{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 ${\cal O}_{\rm SM}{\cal 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 Hall:2009bx ; Elahi:2014fsa . Different than the usual freeze-in via renormalizable operators which is not senstive to UV McDonald:2001vt ; Hall:2009bx ; Bernal:2017kxu ; Kang:2014cia , 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 Co:2015pka ; McDonald:2015ljz ; Garcia:2017tuj ; Chen:2017kvz ; Garcia:2018wtq ; Forestell:2018dnu .

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 Garcia:2018wtq . 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

II.1 The self-interaction & relic density issues of DGB dark matter

At the energy scale below the scale $\Lambda_{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^{++}\equiv s$ Forestell:2017wov is proposed in the spirit of non-perturbative methods,

[TABLE]

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. Boddy:2014yra 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. Soni:2016gzf which just follows the large-$N$ limit, the $4\pi$ factor is not presented. Thus, we may keep an open mind on the strength of these couplings. In this paper we follow Eq. (5) 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\to ss$, which has a large cross section $\sigma_{2\rightarrow 2}\sim(4\pi/N_{d})^{4}/m_{s}^{2}$. The Bullet Cluster imposes an upper bound on $\sigma_{2\rightarrow 2}/m_{s}\lesssim 1{\rm cm}^{2}/g$ bullet , which in turn gives the lower bound on the DGB mass DGB ; Forestell:2017wov ,

[TABLE]

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\to 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\rightarrow 2s$, which could maintain the chemical equilibrium among the DGBs. After the $3s\rightarrow 2s$ depletion rate falls below the Hubble expansion rate at ${\cal 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 DGB ; DGB:relic ; Forestell:2017wov ,

[TABLE]

where $g_{S}$ and $g_{S}^{d}$ are the entropy degrees of freedom in the visible and dark sectors densities, respectively; $\xi_{T}\equiv{\cal T}/T$ is the ratio between the temperatures of the two sectors. All of them are calculated at a temperature ${\cal T}_{i}\gtrsim\Lambda_{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 ${\cal T}_{d}\lesssim\Lambda_{d}$. In particular, ${\cal T}_{d}\simeq\Lambda_{d}$ in the limit $c_{5}\rightarrow 0$ Boddy:2014yra ; Soni:2017nlm . Otherwise, one may parameterize the deviation by ${\cal T}_{d}=w_{f}\Lambda_{d}$ with $w_{f}\lesssim 1$. Anyway, from Eq. (4) we immediately see that, as stated in the introduction, $\xi_{T}\ll 1$ is indispensable to reduce the large ratio ${\cal T}_{d}/(3.6{\rm eV})\simeq\Lambda_{d}/(3.6{\rm eV})\gtrsim 10^{8}$ indicated by Eq. (3). The main goal of this paper is to specify the UV freeze-in as the natural mechanism for generating $\xi_{T}\ll 1$ .

II.2 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}\gg\Lambda_{d}$. Alternatively, the messengers may be a heavy moduli-like field which is neutral under any gauge groups Chowdhury:2018tzw . After integrating out the messengers one obtains the higher dimensional operators ${\cal O}_{\rm SM}{\cal 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,

[TABLE]

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_{d}^{2}$ with $g_{d}$ the dark QCD gauge coupling. The dark QCD gluon field strength tensor is written as ${\cal F}_{\mu\nu}={\cal F}^{a}_{\mu\nu}T^{a}$, with $T^{a}$ the generators of $SU(N_{d})$ satisfying the normalization conditions ${\rm tr}(T^{a}T^{b})=\delta^{ab}/2$. The same convention applies to other non-Abelian gauge groups. We only consider the scalar operator $S\equiv{\rm tr}{\cal F}^{\mu\nu}{\cal F}_{\mu\nu}$, and likely other operators accompany, e.g., ${\rm tr}{\cal F}^{\mu\nu}\widetilde{\cal F}_{\mu\nu}{\rm tr}{G}^{\mu\nu}\widetilde{G}_{\mu\nu}$ and ${\rm tr}{\cal F}^{\mu\alpha}{\cal F}_{\alpha}^{\nu}{\rm tr}{G}_{\mu}^{\alpha}{G}_{\alpha\nu}$, and so on DG:eff1 ; DG:eff2 . But their presence does not change the main line of our discussion, though a possible sizable numerical correction Ellis:2018cos .

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 Cline:2013zca . Otherwise, one may have to arrange the proper Yukawa interactions between the fermionic messengers and the Higgs doublet; see Ref. DG:eff2 . 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. Third, $M_{6,8}$ do not directly correspond to the mass scale of loop particles in the UV model, but usually the difference is not big if there are no very small/large extra Yukawa couplings involved in the loop.

The operators ${\cal O}_{\rm SM}{\cal O}_{d}$ may open decay channels for the DGBs into the SM sector. The calculation of their decay width is straightforward given the factorized matrix elements, which for a certain DGB state $J^{CP}$ decaying into the SM states without involving another DGB in the final state is 111Transitions 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.

[TABLE]

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_{J^{CP}}^{{\cal O}_{d}}$, known as the decay constant of the state $J^{CP}$. Hereafter we will specify $J^{CP}=0^{++}\equiv{s}$ and take $F_{0^{++}}^{S}$ as $F_{{s}}$. From the lattice result FS , one may take

[TABLE]

For a given ${\cal 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.

II.3 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 Hdecay . The partial decay width of ${s}$ into a pair of SM states is written as

[TABLE]

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}\gg m_{W}$, then $s$ dominantly decays into the vector bosons and the resulting lifetime of $s$ is estimated to be

[TABLE]

In order to evade the constraints on decaying DM into vector bosons Cohen:2016uyg , 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. (3). Then, ${s}$ dominantly decays into the heaviest SM fermion pair $f\bar{f}$ kinematically accessible. Probably the most suppressed case of DGB decay is that it dominantly decays into a pair of muons for $1{\rm GeV}\gtrsim m_{s}\gg 2m_{\mu}$. In this case one has

[TABLE]

Even if the region $2m_{e}\ll m_{s}\lesssim 2m_{\mu}$ 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}\ll 2m_{e}$, only a pair of photon is accessible in $s$ decay and the resulting DGB lifetime is

[TABLE]

It is seen that in this case the cutoff scale can be as low as the interesting TeV scale. Therefore, the lower bound on $\Lambda_{d}$ by virtue of self-interaction places a strong lower bound on $M_{6}$.

II.4 Vector-boson-portal DGB ($d=8$)

Now we move to the next case, DGBs communicating with the SM sector via the vector-boson-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 DG:eff2

[TABLE]

Compared to the $d=6$ case, it has a much stronger dependence on the confining scale $\Lambda_{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\rightarrow gg$ is closed today. The lifetime of the decaying DGB in terms of Eq. (12) is estimated to be

[TABLE]

Lowering $\Lambda_{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}\sim\sqrt{\Gamma_{\phi}M_{\rm Pl}}$, with $\Gamma_{\phi}$ the width of the inflaton perturbative decay to radiation which thermalizes instantaneously at $t_{re}\simeq 1/\Gamma_{\phi}$. But there is an even hotter phase than $T_{re}$ during reheating reheating , whose possible impacts on UV freeze-in are investigated by several groups Co:2015pka ; McDonald:2015ljz ; Chen:2017kvz , 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_{\phi}$, and thus the UV freeze-in process, whose cross section is greatly enhanced, may be very effective in this phase. Ref. Garcia:2018wtq 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. Harigaya:2013vwa . The inflaton decay width is parameterized as $\Gamma_{\phi}=km_{\phi}^{2}/M_{\rm Pl}$, from which one can see that $k\gtrsim 1$ leads to $T_{re}\gtrsim m_{\phi}$ and thus thermal effects on the inflaton decay become important thermal ; k>1 . To avoid this complication, in this paper we just focus on the slow decay limit with $k\ll 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 $\alpha\sim{\cal 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 $\Gamma_{\phi}\ll\alpha^{8}M_{\rm Pl}^{3}/m_{\phi}^{2}$, or equivalently,

[TABLE]

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

[TABLE]

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.

III.1 UV freeze-in dark gluons

III.1.1 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-radiation-dark gluon system are

[TABLE]

where the collision term $\widetilde{\gamma}_{d}$ will be specified later. The Hubble parameter is determined by the Friedmann equation:

[TABLE]

In this set of BE describing the energy density transfer, inflaton dominates the system throughout the reheating stage thus $H\approx\sqrt{\rho_{\phi}/3M_{\rm Pl}^{2}}$. Its energy density is simply red shifting, and therefore $H\propto a^{-3/2}$ or $H(a)=H_{end}\left(a_{end}/a\right)^{3/2}$, where $H_{end}=\sqrt{\rho_{end}/3M_{\rm Pl}^{2}}$ with $\rho_{end}$ the energy density of inflaton at the end of inflation. In a class of popular inflation models $\rho_{end}$ is estimated to be

[TABLE]

From the Friedmann equation one can estimate the age of the universe $t=\int_{a_{end}}^{a(t)}da^{\prime}/(a^{\prime}H(a^{\prime}))$, arriving at the $t-a-H$ relation

[TABLE]

The last relation holds for $a\gg a_{end}$. This $t-a$ scaling rule is different than the one in the radiation dominant era Eq. (22) because the reheating era is matter dominant.

Reheating is completed when the inflatons decay away at $t_{re}\simeq 1/\Gamma_{\phi}$, after which the inflaton domination gives way to a period of radiation domination. Radiation energy density contains two components, the ordinary $\rho_{R}$ produced by inflaton decay and the dark radiation $\rho_{g_{d}}$ produced by the ordinary radiation. Such a system is described by Eq. (17) and Eq. (18) with $\Gamma_{\phi}\rightarrow 0$. 222It 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. During this era $\rho_{R}$ dominates over $\rho_{g_{d}}$ and $H\approx\sqrt{\rho_{R}/3M_{\rm Pl}^{2}}\propto a^{-2}$. Then similar to the previous derivation one gets the $t-a-H$ relation

[TABLE]

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 $\Phi=\rho_{\phi}a^{3}/m_{\phi}$, $R=\rho_{R}a^{4}$ and $G_{d}=\rho_{g_{d}}a^{4}$. We also introduce the dimensionless variable $\hat{a}\equiv am_{\phi}$. The radiation 333It is could be the SM gluons or any other members which are charged under the SM non-Abelian gauge groups, such as the Higgs doublet. is gradually built up from the inflaton perturbative decay $\phi\rightarrow RR$, while the dark sector is only slightly reheated by the radiation scattering $R+R^{\prime}\rightarrow g_{d}+g_{d}^{\prime}$. 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. (23) the $\Gamma_{\phi}\rho_{\phi}$ term can be dropped; similarly, in Eq. (17) the $\widetilde{\gamma}_{g_{d}}$ term is removed. Then, the BEs take the simple form during reheating,

[TABLE]

where the prime denotes the derivative with respect to ${\hat{a}}$. While during the radiation dominated era, the BEs are reduced to

[TABLE]

In the no back reaction limit, it is ready to solve the BEs one by one. First, the solution to Eq. (23) is

[TABLE]

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. (23), one may multiply the above solution by the decaying factor $e^{-\hat{a}}$. Next, with Eq. (28), the energy density of radiation background in the phase II can be obtained by directly integrating over $\hat{a}$,

[TABLE]

Whereas in the phase III, the solution to Eq. (26) is a trivial redshiting, and then in Eq. (27) $R$ takes a constant value $R_{re}\equiv R(\hat{a}_{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 $\rho_{R}=\frac{\pi^{2}}{30}g^{re}_{*}T^{4}$ with $g_{re}$ the degree of freedoms of the relativistic particles during reheating. Along with Eq. (29), one gets the scaling rule $T\propto a^{-3/8}$ during the reheating era after thermalization:

[TABLE]

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. (15), Eq. (21), one can determine

[TABLE]

So $T_{max}\sim\alpha^{4/5}\left(kM_{\rm Pl}/m_{\phi}\right)^{2/5}m_{\phi}$. Its magnitude relative to $m_{\phi}$ depends on the value of $k$. The referred value $k\sim k_{0}\equiv m_{\phi}/M_{\rm Pl}$ is of special interest since it makes $T_{re}\ll T_{max}\ll m_{\phi}$, which implies that the thermal effects from radiation is fully under control. However, in general $T_{max}$ can exceed $m_{\phi}$ as $k\gg k_{0}$ and $m_{\phi}$ 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 $\phi\bar{\psi}\psi/M_{\rm Pl}$. After the inflaton decays away, entropy is conserved and thus temperature drops much faster, following the well-known rule

[TABLE]

With Eq. (30) and Eq. (32) one can express the collision terms in the light of $a$.

III.1.2 The cross sections

We now move to the calculation of the collision terms in Eq. (18), which involve the cross sections for dark gluons prodcution. In our setup, the dark gluons are produced via the $2\rightarrow 2$ processes. In general, the Lorentze invariant cross section of the process $1+2\rightarrow 3+4$ is defined as average

[TABLE]

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/(2{s})$. Because $\sigma(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. (5), the squared amplitudes for the two representative cases are given by

[TABLE]

with $N_{c}=3$ for QCD. The process $H+H^{*}\rightarrow 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

[TABLE]

Averages of the initial state over the spin and internal degrees of freedom give rise to the factors $1/D_{H}^{2}=1/2^{2}$ and $1/D_{g}^{2}=1/(2^{2}(N_{c}^{2}-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 $\sim M_{6,8}\gg m_{\phi}$. Otherwise one should take into account the production of mediator Chen:2017kvz .

III.1.3 The collision term

With the Lorenz invariant cross sections Eq. (36) and Eq. (37), the collision term in Eq. (18) can be calculated using the standard approach average .

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_{\phi}/2$; due to the expansion rule Eq. (21), this spectrum is universally redshifted to the later time $t$, located at the momentum $(m_{\phi}/2E)^{3/2}$. Using this argument, the spectrum of the radiation is derived to be Garcia:2018wtq ; Harigaya:2013vwa

[TABLE]

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),

[TABLE]

where $d_{\phi}\equiv\Gamma_{\phi}(t-t_{end})\ll 1$ and $n_{\phi,end}={\rho_{end}}/{m_{\phi}}$, the inflaton number density at the end of inflation. With the above distribution function, one can calculate the collision term

[TABLE]

with $x\equiv m_{\phi}/\sqrt{s}>1$ and $D_{R}$ the internal degrees of freedom of radiation $R$. The Lorentz invariant cross section $\sigma(x)\propto x^{-2}$ and $x^{-6}$ for the $d=6$ and $d=8$ operators, respectively. And the corresponding integral of $x$ are $\frac{1}{30}+\frac{\pi}{192}\approx 0.05$ and $\frac{1}{90}+\frac{7\pi}{5120}\approx 0.015$, giving rise to the following collision terms

[TABLE]

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, $\propto 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

[TABLE]

with $z=\sqrt{s}/T$, while $\sigma(z)\propto 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 $\sqrt{s}$ can be way larger than the given $T$. Actually, the dominant integrating range is near $z\sim 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\times 10^{7}$ for $d=6$ and $d=8$, respectively; consequently, $\widetilde{\gamma}_{g_{d}}\propto T^{2d-3}$. The larger $d$ case benefits from more significant UV enhancement in the $z\sim 10$ region, so it enjoys an impressive numerical enhancement. Note that these big numbers do not appear in $\widetilde{\gamma}_{g_{d}}^{I}$ in phase I where temperature does not exist and $\sqrt{s}<m_{\phi}$.

Eq. (44) 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. (30) and Eq. (32), thus leading to $\widetilde{\gamma}^{II}_{g_{d}}(\hat{a})$ and $\widetilde{\gamma}^{III}_{g_{d}}(\hat{a})$, respectively. Concretely, in the two phases they are given by

[TABLE]

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. (30): For $d=8$, $\widetilde{\gamma}^{III}_{g_{d}}\propto T^{13}$, whereas it is proportional to $T^{12}$ in Ref. Chen:2017kvz .

III.1.4 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 $\hat{a}_{f}>\hat{a}_{re}$ is obtained by integrating over the scaled scale factor $\hat{a}$ 444Different than Ref. Chen:2017kvz , here we adopt $a$ instead of $T$ as the integration variable. The reason is understood by nothing but that there is no temperature in the prethermal era. in the prethermalization, thermalization and radiation dominant era successively until $\hat{a}_{f}$,

[TABLE]

As a feature of UV-freeze in production, it is soon frozen as the scale factor $\hat{a}_{f}$ becomes significantly larger than $\hat{a}_{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}\gg a_{end}$, one has

[TABLE]

Thereby, increasing the dimension of operators leads to a suppressed production yield in the phase I. This contribution has a strong dependence on $\alpha$, stemming from the dependence on $\hat{a}_{th}$. But $\alpha$ will not appear in the contributions from other phases.

Phase II

In the phase II, freeze-in production for the $d=6$ case is again IR dominated, and the contribution is

[TABLE]

where $a_{re}\gg a_{th}$ is assumed. The pure numerical factor including the $\pi$ is $\approx 0.046$. For $d=8$, the integral is,

[TABLE]

The pure numerical factor is about $1.9\times 10^{5}$ taking $N_{c}=3$, which is an impressive enhancement. Note that the one more power of $T$ of $\widetilde{\gamma}_{g_{d}}^{II}$ commented below Eq. (45) is not sufficient to make $G^{II,8}_{d}$ turn into UV dominance, and thus we have the relation shown in Eq. (51).

Phase III

The contribution from Phase III is always UV dominated. This leads to the same expression of $G^{III}_{d}$ as in $G^{II}_{d}$, up to the pure numerical factors: $G^{III,6}_{d}\approx 1.37G^{II,6}_{d}$ and $G^{III,8}_{d}\approx 0.17G^{II,8}_{d}$.

The total comoving energy density of dark gluons is $G^{6/8}_{d}(\hat{a}_{f})=G_{d}^{I,6/8}(\hat{a}_{th})+G_{d}^{II,6/8}(\hat{a}_{re})+G_{d}^{III,6/8}(\hat{a}_{f})$, which as mentioned before is frozen soon after reheating. Therefore, one may take $G^{6/8}_{d}(\hat{a}_{f})$ as a constant $G^{6/8}_{d}$ as long as $\hat{a}_{f}$ is sufficiently large.

It is of interest to estimate the ratio of the contributions from thermal and nonthermal productions,

[TABLE]

It is consistent with the argument at the beginning of this section: The production in the prethermal phase may become dominant for a very heavy $m_{\phi}$, in particular for the case with a larger $d$, which likely belongs to this situation except for a very light $m_{\phi}$ even below the weak scale. So, Our conclusion is consistent with the one drawn in Ref. Garcia:2018wtq .

III.2 DGB Relic density

III.2.1 A much cooler dark sector after reheating

To estimate the relic density of DGB, we adopt the working assumption that the dark gluons are thermalized and establish their own plasma temperature ${\cal T}$ not later than $t_{re}$. Moreover, we approximate that the energy transfer between the dark and visible sectors completely ceased just at $t_{re}$ hence $G_{d}({\hat{a}=\hat{a}_{re}})\approx G_{d}({\hat{a}\gg\hat{a}_{re}})$, amounting to dropping the part from phase III. The resulting error is negligible if UV freeze-in is dominated by phase 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 $\rho_{g_{d}}=\pi^{2}(N_{d}^{2}-1){\cal T}_{re}^{4}/15$, which gives

[TABLE]

Since after reheating both sectors keep cooling down, respectively following the rules ${\cal T}(a)={\cal T}_{re}a_{re}/a$ and Eq. (32), from which we immediately get ${\cal T}=\frac{{\cal T}_{re}}{T_{re}}T=\xi_{T}T$. As a natural consequence of our mechanism reheating the dark sector, $\xi_{T}\ll 1$, which just provides the desired hierarchy to accommodate the correct DGB relic density. Let us show the expression of $\xi_{T}$ for $d=6$:

[TABLE]

Increasing $m_{\phi}$ helps to enhance $\xi_{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.

• •

For a relatively light $m_{\phi}\ll 10^{13}{\rm GeV}(1.0/k)^{1/16}\left(100/g_{re}\right)^{\frac{27}{32}}\left(\alpha/0.01\right)^{2}$, the phase I contribution is absolutely negligible. Note that this condition is insensitive to $k$ because of the power suppression $k^{1/16}$. Now one has the estimation

[TABLE]

with $r_{\phi,6}\equiv m_{\phi}/M_{6}\ll 1$. 555As we have commented before, $M_{6}$ is not exactly corresponding to the mass scale of loop particles, but in the concrete UV model the actual mass scale usually is even lower than $M_{6}$, except for the presence of large couplings in the loop which could overcome the loop suppression. The smallness of $\xi_{T,6}$ is readily obtained by a sufficiently small $r_{\phi,6}$ or $k$. For the reference value $k=k_{0}$, $\xi_{T,6}\sim r_{\phi,6}\left({m_{\phi}}/{M_{\rm Pl}}\right)^{1/8}$ and then how to get a large enough $\xi_{T,6}\gtrsim 10^{-3}$ is of concern; hence instead $r_{\phi,6}$ should have a moderate size.

• •

To the contrary, if $m_{\phi}\gg 10^{13}{\rm GeV}$, then the phase I is the dominant contribution to reheat the dark sector and $\xi_{T,6}\sim k^{5/12}\alpha^{-\frac{4}{3}}r_{\phi,6}\left(\frac{m_{\phi}}{M_{\rm Pl}}\right)^{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_{\phi}\lesssim 10^{15}$ GeV chaotic . This means that there is barely room for this scenario. It is tempting to fix $m_{\phi}=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

[TABLE]

where $N_{c}=3$ has been used. The first term denotes for the contribution from phase I, and it dominates the other contributions when

[TABLE]

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_{\phi}\gtrsim 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\ll 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 $\xi_{T,6}$ is estimated to be

[TABLE]

Taking $k=k_{0}$, to make $\xi_{T,8}\gtrsim 10^{-4}$ the inflaton must be as heavy as around $10^{13}$GeV. But if we work in the intermediate region $1\gg k\gg k_{0}$, the required $m_{\phi}$ can be accordingly lighter.

III.2.2 Dark confining phase transition and relic density of the dark gluon ball

As ${\cal T}$ drops below $\Lambda$, 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$,

[TABLE]

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 $\propto a^{-3}$ until today, when the visible sector temperature $T_{0}=2.37\times 10^{-13}$ GeV. Hence, its fraction in the total energy budget $\Omega_{s}h^{2}\equiv\frac{\rho_{s}(a_{0})}{\rho_{c}}=\frac{\rho_{s}(a_{\Lambda})}{\rho_{c}}\left(\frac{a_{\Lambda}}{a_{0}}\right)^{3}$ is

[TABLE]

where we have used the critical energy density $\rho_{c}=8.1\times 10^{-47}h^{2}{\rm 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. (5). The two approaches give the same scaling behavior $\Omega_{s}h^{2}\propto\xi_{T}^{3}\Lambda_{d}$.

As we have stressed, a small $\xi_{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 of DGB with correct relic density on the $(m_{s},\tau_{s})$ plane in Fig. 1, choosing several values of $r_{\phi}$. This is done by trading $M_{6,8}$ with $\tau_{s}$ in Eq. (9) or Eq. (10), $\Lambda_{d}$ with $\Omega_{s}h^{2}$ via Eq. (58). Several observations are made:

• •

In general, for a given DGB mass, a larger $r_{\phi}$ leads to a significantly shorter DGB lifetime. This is because, to maintain a constant $\xi_{T}$, accordingly $M_{6,8}=m_{\phi}/r_{\phi}$ becomes smaller thus much faster DGB decay. In other words, the parameter space having a larger $r_{\phi}$ 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. Cohen:2016uyg , displaying the FERMI-LAT gamma ray data FERMILAT and IceCube neutrino data IceCube 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 ${(G^{II,8}_{d}+G^{III,8}_{d})}/{G^{I,8}_{d}}>5$. In the heavy DGB region $\sim{\cal O}(\rm TeV)$, referring to the first and second panels, it is clearly seen that the larger $r_{\phi}$ 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_{\phi}$, 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 $\gtrsim 10^{28}s$ for a PeV scale DGB Cohen:2016uyg .

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 $\sim 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_{\phi}$, 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.

Acknowledgements

I would like to thank Institute of Theoretical Physics, Chinese Academy of Sciences for the hospitality during my visiting there to finalize this manuscript. This work is supported in part by the National Science Foundation of China (11775086).