Hidden Dark Matter from Starobinsky Inflation

The Starobinsky inflation model is one of the simplest inflation models that is consistent with the cosmic microwave background observations. In order to explain dark matter of the universe, we consider a minimal extension of the Starobinsky inflation model with introducing the dark sector which communicates with the visible sector only via the gravitational interaction. In Starobinsky inflation model, a sizable amount of dark-sector particle may be produced by the inflaton decay. Thus, a scalar, a fermion or a vector boson in the dark sector may become dark matter. We pay particular attention to the case with dark non-Abelian gauge interaction to make a dark glueball a dark matter candidate. In the minimal setup, we show that it is difficult to explain the observed dark matter abundance without conflicting observational constraints on the coldness and the self-interaction of dark matter. We propose scenarios in which the dark glueball, as well as other dark-sector particles, from the inflaton decay become viable dark matter candidates. We also discuss possibilities to test such scenarios.


Introduction
A simple model of inflation consistent with the cosmic microwave background (CMB) observation by the Planck satellite [1] is the so-called Starobinsky inflation model [2]. In the Starobinsky model, the R 2 term is added to the Einstein-Hilbert action where R denotes the Ricci scalar. This model contains a scalar degree of freedom, called scalaron, and it can be regarded as the inflaton if the coefficient of the R 2 term takes an appropriate value. One interesting feature of the Starobinsky model is that the inflaton universally couples to the Standard Model sector as well as other sectors through the trace of the energy-momentum tensor, and hence the reheating temperature is predictable. In order to explain the dark matter (DM) of the universe, we may consider models with the dark sector containing DM, which is sequestered from the visible sector and communicates with the visible sector only via the gravitational interaction. Then, it is also possible to predict the inflaton decay rate into the dark sector.
In this paper we consider a simple dark sector extension to the Starobinsky inflation in order to explain DM. Such a scenario has been discussed in [3,4,5]. It has been shown that an addition of massive scalar or fermion, as well as SU (N ) gauge boson may be enough to explain the correct DM abundance through the inflaton decay. We revisit this issue. In particular, we pay particular attention to the case of dark (pure) non-Abelian gauge sector. A characteristic property of the dark glueball DM is that it has a sizable self-interaction cross section. In addition, we also consider the cases that the DM candidate is a scalar boson, a fermion, and a massive Abelian gauge boson which was not discussed in the previous studies.
Let us mention several relevant works. In the past, the self-interacting DM was proposed to alter the density perturbation [6] where the entropy of the dark sector is given by hand (see also [7,8]). Recently such kinds of scenarios have become attracting attentions again [9,10,11,12] since they may explain the small-scale crisis and alleviate the S 8 tension [13,14,15,16,17]. The DM production in this context has also been popularly studied. For instance the DM can be produced via the coupling to SM particles [18,19,20], or via the inflaton decay by assuming the vanishing coupling to the SM particles [21]. (See also the DM production from inflaton decay in a general context for heavy DM candidates [22,23,24,25,26,27] and light sub-keV DM candidates [28,29].) Our scenario may be regarded as a concrete model of self-interacting DM produced by the inflaton decay. In particular, we consider the Starobinsky inflation, which is motivated from the current CMB observation, and assume a dark sector containing a DM candidate, which automatically suppresses the coupling to the SM particles. We pay particular attention to the case with a dark SU (N ) gauge group whose dark glueball becomes the DM candidate [30,31,32,33,34]. Interestingly, the inflaton coupling to the dark sector is determined, and hence we have rigid predictions. We find that it is difficult to explain the present DM in the minimal setup since the constraint on the DM self-interaction cross section is in tension with the prediction DM abundance. We also discuss several loopholes to avoid such constraints. This paper is organized as follows. In Section 2, we give an overview of the model of our interest. In Section 3, we study the possibility of the dark glueball DM in the minimal setup. We will see that the dark glueball DM is hardly realized in the minimal setup. In Section 4, we propose scenarios to make the dark glueball a viable candidate of the DM. In Section 5, we discuss the cases that the dark sector contains a scalar boson, a fermion, or a massive Abelian gauge boson as a candidate of the DM. In Section 6, we summarize our results.

Starobinsky inflation model
In this section, we summarize basic properties of the Starobinsky inflation which are necessary for our analysis. We start with the action in the Jordan frame for the Starobinsky inflation [2]: where M Pl 2.4 × 10 18 GeV is the reduced Planck scale, R denotes the Ricci scalar, and S vis and S dark are actions for the visible sector (which contains the SM particles) and dark sector (which contains DM), respectively. #1 Here and hereafter, the "hat" is for quantities in the Jordan frame.
The action S grav is equivalent to the following one (for the review of so-called f (R) gravity, see, for example, [35]): Regarding ϕ as an auxiliary field, the Eular-Lagrange equation gives ϕ = 1 − 1 3µ 2R , with which the above action becomes equal to the original one given in Eq. (2.2). In the following, we use the action given in Eq. (2.3), which is more convenient for our discussion.
With the action S grav given above, there exists a physical scalar degree of freedom which is often called scalaron; the scalaron plays the role of the inflaton in the Starobinsky inflation model. To see this, it is convenient to work in the Einstein frame. The Einstein-Hilbert action can be obtained with the following Weyl transformation: with whichR and the Ricci scalar in the Einstein frame R are related aŝ We can see that there exists a scalar degree of freedom; the canonically normalized scalar field φ is related to ϕ as In terms of φ, the Einstein frame action for the gravity sector is given by where the potential of the scalaron is given by (2.8) The minimum of the potential of φ is chosen to be φ = 0; expanding the potential around the minimum, we can find We adopt the sign convention of the metric of g µν ∼ (+, −, −, −). The Ricci scalar during inflation is given by R −12H 2 , with H being the expansion rate. and hence µ can be regarded as the mass of φ in the vacuum: (2.10) As one can see from Eq. (2.8), the potential of φ becomes flat when φ → ∞. Thus, if the initial amplitude of φ is large enough, the slow-roll inflation can occur.
Regarding V (φ) as the inflaton potential, the curvature perturbation amplitude A s , the scalar spectral index n s , and the tensor-to-scalar ratio r are evaluated as 13) where N e is the e-folding number for the epoch at which the mode k (with k being the present physical wave number) exits the horizon, and is related to the inflaton amplitude during inflation as (2.14) Taking N e (k = 0.05 Mpc −1 ) 50 − 60, the observed value of A s 2.1 × 10 −9 [36] is realized with µ 3 × 10 13 GeV. The scalar spectral index becomes consistent with the observed value while the tensor-to-scalar ratio is well below the current upper bound [1].

Reheating
A brief overview the cosmic evolution is as follows. Assuming that the initial amplitude of φ is large enough, the Starobinsky inflation occurs, during which φ slowly moves towards the minimum of the potential. When the inflaton amplitude becomes ∼ M Pl , the slow-roll condition breaks down and the inflation ends. Then, the inflaton starts to oscillate around the potential minimum; soon after the start of the oscillation and well before the decay, the energy density of the inflaton scales as a −3 (with a being the cosmic scale factor). When the expansion rate of the universe becomes comparable to the total decay rate of the inflaton, denoted as Γ φ , the inflaton decays effectively and the universe is reheated.
In order to discuss the detail of the reheating process, we consider the interaction of the inflaton with the matter particles. The interaction terms of the inflaton in the Einstein frame show up after the Weyl transformation, and are expressed in the following form: where T µν is the energy momentum tensor: With Eq. (2.15), we can calculate the decay rate of the inflaton. The partial decay rates depend strongly on the properties of the final-state particles. Below we list them for the cases of a final state scalar, a fermion and a vector boson.
• For a scalar field χ, we consider the following Lagrangian: #2 where we introduce the non-minimal coupling constant ξ χ . Here we assume that the scalar DM is charged under a Z 2 symmetry to forbid the interactionRχ. We will come back to remove this assumption later. Then, the decay rate is given by where we assume that χ is much lighter than φ.
• For a Dirac fermion ψ with the decay rate is given by (2.20) • We may also consider the decay of φ into a pair of gauge bosons in association with a gauge group G. At the tree level, the partial decay rate Γ φ→GµGµ (with G µ denoting the gauge boson) vanishes because the energy momentum tensor of the gauge boson is traceless. With taking into account quantum effects, however, this is not the case because where α G is the gauge coupling constant, β G is the β-function of α G , and G (a) µν is the field strength tensor. We use the 1-loop β-function:

22)
#2 Here, we neglect the interaction of the scalar field because we are interested in the two body decay process φ → χχ.
where b G is a constant with which the coupling constant at two different renormalization scales, Q and Q 0 , are related as Then, the decay rate is given by [4] where D G is the number of gauge bosons (i.e., dimension of the group G).
• One may consider a massive Abelian gauge boson V µ (with its field strength denoted as V µν ) with the Lagrangian When the mass of V µ is much lighter than the inflaton mass, the production of the longitudinal mode dominates and the decay rate is given by (2.26) The mass of V µ may arise from the spontaneous breaking of the gauge symmetry due to the vacuum expectation value of a scalar field charged under the U (1) symmetry. For simplicity, here and hereafter we consider the case that the Higgs mode in association with the breaking of the U (1) symmetry is much heavier than the inflaton. Such a case is realized when v → ∞ and g U (1) → 0 (where v is the vacuum expectation value of the scalar field responsible for the U (1) symmetry breaking and g U (1) the gauge coupling constant) while keeping the product g U (1) v finite. If g U (1) is exactly equal to 0, the decay process should be understood as the production of the Nambu-Goldstone (NG) mode. When v → ∞, the Higgs mode should decouple from the low energy spectrum and we can consider an effective field theory containing only the inflaton and the NG mode. In such a case, using the fact that the non-minimal coupling of the NG scalar is forbidden by the shift symmetry, we can see that the decay rate given in Eq. (2.26) parametrically agrees with the decay rate to the NG mode (see Eq. (2.18)).
In order for a viable thermal history, the inflaton should dominantly decay into visible sector fields. The inflaton dominantly decays into scalar fields unless the scalars are all conformally coupled (i.e., ξ χ 1 6 ). If the decay processes into the scalars are suppressed, the inflaton may also decay into gauge fields. #3 Here, we evaluate the decay rate of the inflaton #3 There are also decay processes to SM fermion pairs and a Higgs boson via the trace anomaly. The three-body decay rates are suppressed by the phase spaces and/or loop factors.
into the visible-sector field as by taking account of the SM Higgs field and gauge bosons as final states, where ξ h is the SM Higgs non-minimal coupling to the Ricci scalar in the form of (2.17). In addition, A = 1, 2, and 3 are contributions of U (1) Y , SU (2) L , and SU (3) C gauge groups, respectively, with D 1 = 1, D 2 = 3, and D 3 = 8. Hereafter, for simplicity, we assume that the β-functions of the SM gauge coupling constants are determined only by the SM fields (at least at the 1-loop level), and take b 1 = 41 6 , b 2 = − 19 6 , and b 3 = −7. We can relate Γ φ→vis with the reheating temperature after the inflaton decay. In our discussion, we consider the case that the inflaton dominantly decays into visible-sector fields. Then, the reheating temperature, which is the visible-sector temperature at the end of the inflaton decay, is evaluated as where t R is the cosmic time at the end of the reheating, and g is the relativistic degrees of freedom of the SM. Assuming that φ → h † h is the dominant decay channel of φ, with h being the SM Higgs field, g (t R ) = 106.75 and µ = 3 × 10 13 GeV, we obtain Here let us comment on the "gravitational particle production" in the Starobinsky model. The gravitational particle production refers to the creation of fields that do not have any direct interaction to the inflaton at the Lagrangian level. In this sense, the inflaton decay in the Starobinsky model may be regarded as a result of gravitational production, since the inflaton sector and the matter sector are decoupled in the Jordan frame (2.1). The gravitational particle production during the reheating phase #4 has been discussed in Refs. [45,46,47,42]. It is shown that the gravitational particle production picture in the Jordan frame gives the same result as the inflaton decay picture in the Einstein frame [48,46].

Mass Density of Dark Glueball: Minimal Case
First, we discuss the possibility that the dark glueball in the dark sector becomes DM [30,31,32,33,34]. In this Section, we consider the minimal case that the visible sector #4 Gravitational particle production often means the particle creation during inflation or at the transition epoch from de Sitter to the reheating phase [37,38,39,40]. In particular, the vector DM production in this context has been discussed in Refs. [41,42,43,44,34] contains the SM particles while the dark sector contains an SU (N ) pure Yang-Mills gauge sector. In this case, after the inflation, the inflaton decays into the SM particles as well as the dark gluon. Then, after the confinement of SU (N ), which is expected to occur when the temperature in the dark sector becomes lower than the dynamical scale, the dark glueballs are formed. We denote the lightest dark glueball as X. We are interested in the production of the dark glueball from the inflaton decay, and we assume that dark glueball mass m X is smaller than the inflaton mass.

Abundance of dark glueball
With the cosmic expansion, the universe evolves as follows: • After inflation, the dark gluon as well as the SM particles are produced by the inflaton decay. The partial decay rate of the inflaton into the dark gluon pair is given by Eq.
(2.24) with the β-function coefficient: The gauge coupling constant at the renormalization scale Q is evaluated as where Λ is the dynamical scale of the dark SU (N ) gauge theory. In addition, • When ρ dark Λ 4 , the dark sector is described as a relativistic matter consisting of dark gluons. In such an era, assuming that the SM particles are also relativistic, the evolutions of the energy densities of φ, the visible sector, and the dark sector, denoted as ρ φ , ρ vis , and ρ dark , respectively, are governed bẏ where the "dot" denotes the derivative with respect to time, and the expansion rate H is given by with ρ tot ≡ ρ φ + ρ vis + ρ dark . The temperature of the dark sector in such an era is related to ρ dark as while the entropy density in the dark sector is given by After the completion of the inflaton decay, s vis (which is the entropy density of the visible sector) and s dark are both proportional to a −3 . The dark-sector temperature at the time of the reheating is estimated as where we have assumed that the inflaton dominantly decays into the Higgs pair. Notice that, even after the confinement of dark SU (N ), s dark ∝ a −3 , assuming that there is no extra entropy production like a first order phase transition to the confinement phase.
• When the dark-sector temperature T dark becomes lower than ∼ Λ, we expect that the confinement of the dark SU (N ) occurs. Then, the dark sector contains only the dark glueball degrees of freedom. Because there is no symmetry forbidding 2 ↔ 3 (and other number changing) processes of dark glueballs, the dark glueballs are in chemical equilibrium as far as the scattering rate of the number changing processes is larger than the expansion rate of the universe. Assuming that the thermal bath of the dark sector is dominated by the lightest dark glueball, and that the dark glueball is in chemical equilibrium, the energy density, pressure, and entropy density of the dark sector are given by with m X being the mass of dark glueball. The number density of the dark glueball is given by (3.14) Because of the conservation of the entropy in the comoving volume, the dark-sector temperature does not decrease rapidly after the confinement of the dark SU (N ) as far as the chemical equilibrium is maintained in the dark sector [6]. When the chemical equilibrium is maintained in the dark sector, n X decreases slightly slower than a −3 , while the temperature scales as T dark ∝ 1/ ln a. In Fig. 1, we plot T dark (normalized by Λ), s dark (normalized by Λ 3 ), and n dark (normalized by Λ 3 ) as functions of the scale factor, assuming the chemical equilibrium; in the figure, the scale factor is normalized as a| T dark =Λ = 1. • When the rate of the number changing processes becomes smaller then the expansion rate, the number of the glueball in the comoving volume is fixed. Even in such an era, the 2 to 2 elastic scattering may be effective; we will consider the effect of the elastic scattering later.
In our analysis, we assume the adiabatic expansion of the universe after the reheating. This is the case for N = 2; the confinement phase transition is second order for N = 2, while it is first order for N ≥ 3 [49]. #5 In the following discussion, we take N = 2. A Comment on the case of N ≥ 3 will be given at the end of this Section. In order to estimate the present mass density of the dark glueball, we first solve Eqs.
In addition, we use the fact that the entropy ratio s dark /s vis is a constant of time; in our analysis, we calculate the entropy ratio by solving Eqs. For the calculation of the mass density of the dark glueball, we need to estimate when the chemical equilibrium of the dark sector breaks down and the number of the dark gluon in the comoving volume freezes out; the cosmic time at the freeze out is denoted as t F . We evaluate t F by solving the following equation: #6 where the right-hand side is the scattering rate for the 3X ↔ 2X scattering processes in the chemical equilibrium. The accurate calculation of the scattering rate is hardly performed because of the strongness of the interaction. Here, from a dimensional consideration, we parameterize where c 32 is a constant; we assume c 32 ∼ O(1). In fact, the resultant mass density of the dark glueball is insensitive to the parameter c 32 because the number density of the dark glueball approximately scales as a −3 (although, as mentioned before, it decreases slightly faster than a −3 ). Because Γ 3X↔2X has an exponential dependence on T dark , the dark-sector temperature at the time of the freeze out is insensitive to ξ h . Let us define Then, with numerically solving Eq. (3.16), we found x F ∼ 20 for m X = 100 MeV. We estimate the present mass density of the dark glueball as where, in the second equality, we have used the relation n X /s dark T dark /m X in the nonrelativistic limit. The ratio of the entropy densities s dark /s vis , which is independent of time, #6 We use Eq. (3.16) to estimate the freeze-out temperature of the dark glueball. After the freeze out, the inverse process becomes irrelevant so that the Boltzmann equation becomeṡ Assuming the radiation dominance at the time of the freeze out, the evolution of the number density after the freeze out is given by Thus, the number of the glueball in the comoving volume is almost conserved if H Γ 3X↔2X . In [6], the freeze-out condition is slightly different; we have checked that the mass density of the glueball does not change significantly even if we use the freeze-out condition given in [6]. is calculated by numerically solving Eqs. (3.4) -(3.6), as we have mentioned. The ratio is insensitive to Λ because it is determined by the relative size of Γ φ→vis and Γ φ→dark ; for Λ ∼ 100MeV, we found Comparing Eq. (3.19) with ρ crit /s now 3.6h 2 × 10 −9 GeV (where ρ crit and s now are critical density and the present entropy density, respectively, and h is the present expansion rate in units of 100 km/sec/Mpc), we can obtain the density parameter of the dark gluon: Ω X = m X n X (t now )/ρ crit (with t now being the present cosmic time). In Fig. 2, we show Ω X h 2 as a function of Λ for several values of ξ h ; the density parameter is approximately obtained as where we have assumed that the inflaton dominantly decays into the SM Higgs pair.

Constraints
For the dark glueball DM, we should consider several cosmological and astrophysical constraints. Hereafter, the constraints are discussed in order.
The dark glueball may have sizable velocity at the time of the freeze out, and it may affect the structure formation. In particular, if the glueball mass is too small, it may affect the formation of the damped Lyman-α system; according to [50], the dark glueball mass as the dominant component of DM is required to be heavier than ∼ O(10) keV, and hence Λ O(10) keV. (3.22) From the observations of bullet clusters, the self interaction of the DM is known to be weak. The self-interaction cross section σ 2X↔2X is bounded from above [51,52,53,54,55]; here we adopt σ 2X↔2X m X 1 cm 2 /g, Then, it requires ξ h 70 to make the dark gluon density consistent with the observed DM density of Ω DM h 2 0.12 [36]. Such a value of ξ h is, however, disfavored from the point of view of the stability of the electroweak vacuum (see below).
We should also consider the stability of the electroweak vacuum during and after inflation. It has been well known that the electroweak vacuum is unstable if the SM is valid up to a scale much higher than the electroweak scale (like the Planck scale). With the central values of the SM parameters inferred from experiments, the lifetime of the electroweak vacuum for the flat spacetime has been calculated to be much longer than the present cosmic time [56,57,58,59]. However, the stability of the electroweak vacuum is not guaranteed during [60,61,62,63,64,65,66,67] and after inflation [68,69,70,71,72,73,74,75,76] in particular with the non-minimal coupling of the SM Higgs doublet. In the present scenario, the inflaton is assumed to decay into the SM Higgs field via the following interaction: In the inflaton oscillating era after inflation, the non-minimal coupling of the Higgs doublet is dangerous because the oscillating behaviors of R (as well as that of the inflaton) may cause an instability [69]. The Einstein frame action of the Higgs is given by By defining the canonical Higgs field as ϕ −1/2 h, the effective mass of the canonical Higgs field is calculated as #7 In conventional inflation models like the chaotic inflation model, the expression of the Higgs effective mass is different and is given by When the universe is dominated by the inflaton,R is given bŷ During inflation, we obtain ∆m 2 h 12ξ h H 2 . Thus, if ξ h is positive and is larger than ∼ 0.1, the Higgs field is forced to stay at the origin during the inflation even with taking into account the quantum fluctuation; it can help stabilizing the electroweak vacuum. In [67], a lower bound ξ h 0.05 was obtained in the Starobinsky inflation. In the inflatonoscillating era after inflation, on the contrary, the non-minimal coupling of the Higgs doublet may destabilize the electroweak vacuum because ∆m 2 h oscillates and it takes negative value during the oscillation. The constraint on the Higgs non-minimal coupling has been studied in [69] in the chaotic inflation model. We translate the result of [69] to constrain the nonminimal coupling of the Higgs in the Starobinsky inflation model. (More detailed study about the stability of the electroweak vacuum in the Starobinsky inflation is left as a future project [77].) For the chaotic inflation, the smallest effective mass allowed by the stability of the electroweak vacuum is found to be about −5 × 10 27 GeV 2 . Using the fact that the maximal value of ϕ −2R in the Starobinsky inflation model is about 9 × 10 26 GeV 2 , we find the upper bound on the non-minimal coupling to be #8 ξ h 6. (3.30) The constraint (3.30) from the stability of the electroweak vacuum is in tension with the lower bound of ξ h 70 from the bullet cluster constraint (3.24). We have come to this conclusion with the study of the case of N = 2. For the case of N ≥ 3, as we have mentioned, the confinement phase transition is expected to be first order [49]. Thus, in such a case, an extra entropy production in the dark sector is expected as a result of the release of the latent heat, which enhances the energy density of the dark sector compared to the case of second-order phase transition. Thus, for N ≥ 3, the value of ξ h to realize Ω X = Ω DM becomes larger for a fixed value of Λ. As a consequence, the tension in two constraints is expected to become more severe for N ≥ 3. Thus, in the minimal setup, the dark glueball can hardly play the role of DM with satisfying astrophysical and cosmological constraints. In the next Section, we discuss possibilities to make the dark glueball a viable candidate of DM.
Before closing this section, we comment on the stability of the dark glueball. In the present scenario, the stability is not guaranteed by any symmetry [34]. The dark glueball #8 Note that the second term of (3.28), which is linear in φ, is the dominant term to derive this bound. The same term also induces the inflaton perturbative decay into the Higgs pair. may decay into a pair of the SM particles if there exists a direct interaction between the dark sector and the SM sector. Allowing Planck-suppressed higher dimensional operators connecting two sectors, the following dimension-6 operator determines the lifetime of the dark glueball: µν is the field strength of the dark gluon and c hhGG is a coefficient of O(1). Assuming that the decay mode X →f f (with f being an SM fermion) is kinematically allowed, the decay rate of the dark glueball is estimated as Γ X ∼ m 2 f Λ 7 /M 4 Pl m 4 h , where m f denotes the mass of f and m h 125 GeV is the Higgs boson mass. Then, the lifetime is given by (3.32) Thus, for the mass range of the dark glueball considered in our analysis, the lifetime is much longer than the present cosmic time and we can safely neglect the constraints on decaying DM. Notice that, if the dark glueball is heavier, its lifetime becomes shorter. In particular, if the dark glueball is heavier than the electroweak scale, it is severely constrained by astrophysical and cosmological observations. Unstable DM decaying into the SM particles with 1 sec Γ −1 X 10 26 sec affects the big-bang nucleosynthesis, the spectrum of the CMB, and/or the flux of high energy cosmic rays and hence is severely constrained. Such constraints excludes the dynamical scale of 10 5 GeV Λ 10 10 GeV [34].

Dark Glueball as Dark Matter
In the previous section, we have seen that the mass density of the dark glueball hardly becomes consistent with the present DM density in the minimal case, adopting the astrophysical and cosmological constraints. In the parameter region consistent with the constraints, the relic abundance of the dark glueball becomes so large that the dark glueball cannot play the role of DM. In non-minimal cases, on the contrary, this conclusion may be altered. One possibility to realize the dark glueball DM is to introduce an extra entropy production in the SM sector after the reheating due to the inflaton decay. With an extra entropy production, the dark glueball density is diluted so that the difficulty due to the overproduction may be avoided. Another possibility is to add extra fields to make the electroweak vacuum (absolutely) stable or to enhance the inflaton decay rate into the visible sector. In the following, we consider such possibilities.

Entropy production
The extra entropy production often occurs when there exists a long-lived particle which dominates the universe before it decays. In the minimal setup discussed in the previous section, there is no candidate of such a long-lived particle. Here, we consider a simple extension of the setup; we introduce a new sector containing another non-Abelian gauge interaction. For simplicity, we assume that the new sector, called dark sector, consists only of SU (N ) gauge boson with its dynamical scale Λ . (Here and hereafter, we put "prime" on quantities in the dark sector.) Then, because the inflaton of the Starobinsky inflation couples universally to the trace of the energy momentum tensor, it decays also into the dark gauge bosons. With the confinement of the SU (N ) gauge interaction, glueball in the dark sector (which we call X ) should be formed and it behaves as non-relativistic matter. Assuming an interaction connecting dark sector and the SM sector, like the one given in Eq. (3.31), X decays into SM particles. If the dynamical scale Λ is high enough, the lifetime of X becomes shorter than the present cosmic time. In particular, with a proper choice of Λ , X dominates the universe before its decay, and an extra entropy production may occur which reduces the mass density of X.
In the following, we quantitatively estimate the dilution factor due to the decay of X and see if the mass density of X can become consistent with the present DM density. For simplicity, we concentrate on the case with N = 2 and Λ ∼ 10 12 GeV. Notice that the dynamical scale of SU (N ) is assumed to be lower than the inflaton mass so that the inflaton can decay into the dark gauge bosons; the decay rate is given in Eq. (2.24). In addition, in order to make X decay into SM particles, we assume that there exists the following operator: where G (a) µν is the field strength of the SU (N ) gauge field. The decay rate of X is estimated as The "cut-off" scale M is model-dependent. The interaction may arise from the effect of quantum gravity with M ∼ O(M Pl ). A smaller value of the cut-off scale is also possible in particular if there exists a heavy particle, which couples to the SM Higgs doublet, with an SU (N ) quantum number. An example is an SU (N ) non-singlet heavy scalar particle (which we call η) with the quartic interaction of ∼ h † hη † η; in such a case, M is given by the mass scale of η (with a loop factor). We do not specify the origin of the interaction given in Eq. (4.1) and leave M as a free parameter. Let us estimate the dilution factor due to the decay of X . As in the case discussed in the previous section, the decay process into the dark sector occurs via the trace anomaly. Because we take Λ Λ, α dark (m φ ) (which is the SU (N ) gauge coupling constant at Q = m φ ) is larger than α dark (m φ ). Consequently, the partial decay rate Γ φ→dark becomes significantly larger than Γ φ→dark , and hence the production of the SU (N ) gauge boson is more effective than that of SU (N ) gauge boson because With our choice of parameters, Γ φ→dark /Γ φ→dark ∼ O(10 −2 ). For Λ ∼ 10 12 GeV, (i) the energy density of the dark sector is much smaller than Λ 4 at the time of the reheating and hence the dark sector is always in the confinement phase after the reheating, and (ii) the scattering rate for the number changing processes of X is smaller than the expansion rate. As a consequence, X behaves as a non-relativistic matter after the reheating (but before the epoch of H ∼ Γ X →h † h ), and its energy density scales as ∼ a −3 . Assuming that X decays after dominating the universe, the energy density of X at the time of the decay (denoted as t dec ) is given by ρ X (t dec ) ∼ M 2 Pl Γ 2 X →h † h , and the visible sector temperature after the decay is estimated as Numerically, (4.5) Let us define the dilution factor: where s (would-be) vis is the visible-sector entropy density without the effect of the entropy production and is given by s and ρ X both scale as a −3 for t R t t dec , the dilution factor is given by Here, it is assumed that Γ φ→dark Γ φ→vis , and also that g (t dec ) = g (t R ). With our choice of model parameters, (4.8) Thus, if M and Λ are properly chosen, the mass density of X can be reduced down to the present DM density in the parameter region consistent with the astrophysical and cosmological constraints. So far, we have considered the entropy production due to the decay of heavy dark glueball. We comment that the entropy production is also possible by the decay of visible sector particles. If there exists an extra boson or fermion in the visible sector and also if it dominates the universe at some epoch, its decay into the SM particles dilutes the dark glueball. Then, the dynamical scale relevant for the dark glueball DM can be larger. One candidate may be the right-handed neutrino that couples very weakly to the SM leptons; it does not contribute too much to the active neutrino masses as far as its Yukawa interaction is weak enough. If its mass is around 10 12 GeV, the production due to the inflaton decay can be efficient. Then, if the Yukawa coupling is smaller than O(10 −10 ), the mass of the glueball DM can be heavier than O(100) MeV.

Other loopholes
Before closing this section, we discuss other possibilities to realize the dark glueball DM without conflicting the astrophysical and cosmological constraints.
The main reason for the minimal dark glueball DM scenario does not work was that the large non-minimal coupling ξ h , required for correct cold dark glueball DM abundance, would lead to the collapse of the electroweak vacuum during reheating, as discussed in the previous section. However, the stability of the electroweak vacuum crucially depends on the existence of new particles.
For example, one may introduce an additional scalar field Ψ, which develops a vacuum expectation value Ψ and has a portal coupling to the SM Higgs as ∼ |h| 2 |Ψ| 2 , can make the electroweak vacuum absolutely stable due to the scalar threshold effect depending on the mass and coupling constant of the scalar [78,79]. In this case, large non-minimal coupling ξ h does not lead to any cosmological problem. Such a scalar field Ψ may be identified with the B − L Higgs field or the Peccei-Quinn scalar, and hence is theoretically motivated well. See also Ref. [80] for other models to make the electroweak vacuum stable.
If there exist extra scalar fields in association with, for example, the B − L symmetry or the Peccei-Quinn symmetry, the inflaton may also decay into those scalars in particular when their non-minimal coupling ξ extra is sizable. This fact motivates us to consider another possibility to realize the dark glueball DM. If ξ extra is sizable, the inflaton may dominantly decay into such an extra scalar field. Assuming that the potential of such a scalar field does not have instability, Γ φ→vis can be safely enhanced by taking large ξ extra . With large enough ξ extra , the branching ratio into the dark sector is suppressed and Λ consistent with Ω X = Ω DM becomes so large that the bullet cluster constraint can be avoided.

Other Dark Matter Candidates
So far, we have studied the possibility of dark glueball DM. In this section, we discuss other possibilities, i.e., the cases that the dark sector contains a scalar field, a fermion, or a massive Abelian gauge boson. (They are collectively denoted by X in this Section.) We assume that they are so weakly interacting that they freely propagate after being produced.

Dark matter from inflaton decay
The number density of X (denoted as n X ) evolves aṡ Here X is assumed to be pair produced by the decay of φ. We solve Eqs. (3.4) and (3.5) as well as Eq. (5.1) to calculate the relic density of X. Because X survives until today once produced, the number density of X at the end of the reheating can be approximated as For the cases that X is a real scalar, a Dirac Fermion, and a massive Abelian gauge boson, we calculate the density parameter Ω X as follows.
• When X is a real scalar field, we can use Eq. (2.18) to calculate the decay rate.
Assuming that the inflaton dominantly decays into the Higgs doublet pair, we find [3,4] where ξ X is the non-minimal coupling of X. #9 • When X is a Dirac fermion, the partial decay rate of the inflaton into the X pair is given in Eq. (2.20), and we obtain [3,4] (5.4) • When X is a massive Abelian gauge boson, we use Eq. (2.26). Then the relic density is given by For the case that X is a fermion, the relic density consistent with the present DM density is realized with the mass of m X ∼ O(10 7 ) GeV. On the contrary, for the bosonic case, the mass of X is required to be as light as m X ∼ O(10 − 100) keV for Ω X = Ω DM . For light DM candidates, we should consider the constraint from the Lyman-α forest data [81,82]. Such a constraint is important for the cases of scalar and massive Abelian gauge #9 In Ref. [4] the possibility of scalar DM heavier than the inflaton has also been considered. In such a case, the gravitational production is exponentially suppressed but still there may be a parameter region that is consistent with the observed DM abundance. We do not pursue this case further in this paper.
boson. DM models with sizable DM velocity at the time of the structure formation can be characterized by the free streaming length. Here, we calculate the averaged free streaming length: where t eq is the time of radiation-matter equality and a now is the scale factor at present. In addition, v X (t) is the averaged velocity of X produced by the inflaton decay; in the present case, the energy of X is equal to m φ /2 when produced, and the velocity of X is determined solely by the redshift after the production. In Fig. 3, we plot L FS as a function of ξ h for the case that X is a scalar boson, taking the mass of X realizing Ω X = Ω DM (see Eq. (5.3)). The non-minimal coupling is taken to be ξ X = 0 and 0.1. (Notice that the result for ξ X = 0 also applies to the case of massive Abelian gauge boson.) The free streaming length becomes smaller in the limit of large |1 − 6ξ h | because, in such a limit, (i) the reheating temperature is enhanced and hence the velocity of X at the time of the structure formation is more suppressed due to the redshift, and (ii) the branching ratio of φ into X is suppressed so that the mass of X realizing the correct DM density becomes larger. The free streaming length is of O(10) Mpc for ξ h ∼ 0, which is strongly disfavored, while it can be O(0.1) Mpc (or smaller) when ξ h is larger than a few. In order to discuss the consistency with the structure formation, we simply compare L FS in the present model with the free streaming length of warm DM (WDM). It has been pointed out that the Lyman-α forest constraint excludes the WDM with its mass lighter than 5.3 keV [82], for which the free streaming length defined in Eq. (5.6) is ∼ 0.06 Mpc. Requiring L FS 0.06 Mpc, #10 the dark scalar field with ξ X = 0 hardly becomes a viable candidate of DM because such a small free streaming length requires ξ h 10 which conflicts with the constraint from the stability of the electroweak vacuum (see (3.30)). (The same conclusion holds for the case of massive Abelian gauge boson.) With non-minimal coupling of ξ X ∼ 0.1, such a difficulty can be avoided, as indicated by Fig. 3. In addition, the minimally coupled dark scalar field and the massive Abelian gauge boson can be DM without conflicting the astrophysical and cosmological constraints if there exists a sizable entropy production or if the electroweak vacuum is somehow stabilized (as we discussed in the case of dark glueball DM). We also note here that future 21cm line observations will improve the sensitivity to the free streaming length of DM [83,84]. The Square Kilometer Array, for example, will have sensitivity to WDM lighter than ∼ 24 keV (for which L FS 0.01 Mpc. Thus, the future observations of 21cm line signals will provide crucial test of the present scenario with the dark scalar DM. #10 Rigorously speaking, the velocity distribution of the present DM model is different from that in the case of the WDM. We neglect the possible error in the estimation of the upper bound on L FS due to such a difference. A more accurate study of the Lyman-α forest constraint requires a detailed study with a hydrodynamical simulations, which is beyond the scope of our analysis.  Figure 3: L FS as a function of ξ h for the case that X is a scalar boson, taking ξ X = 0 (blue) and 0.1 (red). The result for ξ X = 0 also applies to the case of massive Abelian gauge boson.

Dark matter from scalar coherent oscillation
For the dark scalar DM X, we can consider another possibility than the production through the inflaton decay, i.e., the misalignment mechanism [85,86,87]. The scalar amplitude may be non-vanishing during inflation because of the long wavelength quantum fluctuation or Hubble-induced mass. The former possibility may be dangerous due to the isocurvature problem [47]. In the latter case, the scalar may acquire a large Hubble-induced mass if non-minimal coupling ξ X is |ξ X | 1.
Let us consider the case of ξ X < 0. Then X feels the negative Hubble mass during inflation and its potential minimum is away from that in the present universe. This may lead to the misalignment production of DM. To be concrete, let us assume with ξ X < 0. Here we have introduced a quartic term of X. Then during inflation, X has the amplitude of 8) due to the valance between the (negative) Hubble induced mass and the quartic term (see Eq. (3.28)). Soon (but not very soon) after the inflation, the Hubble induced mass becomes negligible since it is redshifted away faster than the quartic term. Then the X oscillation, driven mainly by the quartic term, starts. At the end of the inflation, the oscillation energy density is given by where t end is the cosmic time at the end of the inflation. The X energy density decreases as that of the radiation as far as the quartic term is dominant for the coherent oscillation. Eventually the mass term becomes dominant and then the X energy density behaves like that of non-relativistic matter thereafter. This transition happens when m 2 X X 2 /2 ∼ λ X X 4 /4; the cosmic time of the transition is denoted as t tr . Then, After this epoch, the ratio ρ X /s becomes a constant and is estimated as (c.f. [88,89,90,91]) As a consequence, we obtain the DM abundance as Ω (osc) We note that X is also produced by the direct decay of the inflaton, but it is sub-dominant and tends to be hot DM if m X O(1) MeV. The transition temperature T vis (t tr ) should be larger than ∼ O(1) keV. It gives a lower bound of m X as m X 1 eV [89,90,91]. There is also a lower bound of m X O(1) keV if |X inf | M Pl . Note also that soon after inflation there may be a tachyonic/parametric resonance if ξ X −1. On the other hand, we may need ξ X −1 to suppress the isocurvature bound. #11 To clarify the resonance effect on this scenario as well as the detailed parameter region requires a further study [77].
Another simple possibility of realizing the misalignment production is to (softly) break the Z 2 symmetry under which X changes its sign. Even without the Z 2 symmetry, the scalar field is long-lived if it is light enough while we may introduce the interaction ofRX. Neglecting the X self-coupling or the direct coupling to other SM particles, the Lagrangian is given by #12 where c X is the Z 2 breaking coupling while the non-minimal coupling ξ X is assumed to be positive for simplicity. We choose the minimum of the potential at the present universe as X = 0 without loss of generality. Assuming that the bare scalar mass m X is much smaller than the expansion rate during inflation, the X amplitude during inflation is given by X −c X M Pl due to the Hubble induced mass term (with ξ X > 0). After the reheating the non-minimal coupling becomes irrelevant. The scalar field X starts to oscillate when the expansion rate becomes comparable to m X . The present energy density of the coherent oscillation is estimated as (5.14) Thus the scalar field can be a fuzzy DM. Assuming ξ X 1, the Hubble-induced mass becomes comparable to or larger than the expansion rate during inflation. Then, the quantum fluctuation of X during inflation is suppressed and hence the isocurvature problem can be avoided. The scalar particle X produced by the inflaton decay is still relativistic in the present universe if m X 1 eV (for which c X 10 7 to make Ω (osc) X = Ω DM ). In such a case, X produced by the decay contributes to the dark radiation. The increase of the effective neutrino number is estimated as Because ξ h 6 due to the vacuum stability bound and also because ξ X 1 for the isocurvature bound, ∆N eff is bounded from below as ∆N eff 0.02. (5.16) Thus, interestingly, this scenario can be tested by the future measurement of dark radiation [92,93,94].

Summary
In this paper, we have studied a production mechanism of glueball DM from a confining dark (pure) SU (N ) Yang-Mills theory. Due to the gauge symmetry the dark gauge particles are rarely interacting with the SM particles, and the DM stability is guaranteed if the confinement scale is much lower than PeV. This also implies that it is difficult to be produced thermally if the reheating temperature is much below the Planck scale. In the Starobinsky inflation, on the other hand, it has been argued that the DM is naturally and easily produced from inflaton decay. Carefully studying the Higgs potential stability during and after inflation, the cosmic evolution, and constraints from bullet clusters as well as structure formations, we have found that this possibility is disfavored in the minimal setup. On the contrary, if there is a further entropy dilution to the SM sector or alternatively the Higgs instability scales happen to be much higher than central values suggested by the current experimental data, we can evade the constraints and a dark glueball becomes a viable DM candidate. As a concrete example, we build a simple model with an additional dark pure Yang-Mills gauge sector, whose confinement scale is close to the inflaton mass. The decay of this additional glueball dilutes the DM density to evade the constraints found by us. To have the reheating temperature higher than the electroweak scale, the DM mass is around hundred MeVs and thus this scenario may explain the small scale crises and alleviate the S 8 tension.
We have also discussed the possibility that a scalar, a fermion, or a massive Abelian gauge boson in the dark sector becomes DM. As well as the production of the DM particle from the inflaton decay, we have discussed the possibility that a scalar field in the dark sector becomes DM. If there exists a scalar field in the dark sector, its coherent oscillation can be induced by the coupling with the inflaton, which can become a viable candidate of DM. We show explicit examples of such scenarios.
The DM candidates in the dark sector interact with SM particles extremely weakly. Thus, the direct detection of such a DM is hopeless with the current technology. However, with cosmological and astrophysical observations, we may acquire some hints of the hidden dark matter in the dark sector. For example, a signal of the dark glueball DM may be seen if the observational constraints on the self interaction of DM can be improved. In addition, for the light scalar (or massive Abelian) DM from the inflaton decay, the future 21cm line observations, which is expected to improve the sensitivity to the free streaming length of DM, can provide an interesting test.
We have also discussed that the non-minimal coupling of the SM Higgs plays important roles in the present scenario. The non-minimal coupling may enhance the partial decay rate of the inflaton into the visible sector, while it may be also important to stabilize the electroweak vacuum during and after inflation. More detail about this issue will be discussed in elsewhere [77].