Production of Purely Gravitational Dark Matter

In the purely gravitational dark matter scenario, the dark matter particle does not have any interaction except for gravitational one. We study the gravitational particle production of dark matter particle in such a minimal setup and show that correct amount of dark matter can be produced depending on the inflation model and the dark matter mass. In particular, we carefully evaluate the particle production rate from the transition epoch to the inflaton oscillation epoch in a realistic inflation model and point out that the gravitational particle production is efficient even if dark matter mass is much larger than the Hubble scale during inflation as long as it is smaller than the inflaton mass.


Introduction
There are many particle physics models to explain cosmological dark matter (DM). Most models introduce (small) DM interactions with some other fields to obtain the desired DM abundance. In this respect, the simplest possibility is that DM interacts with other sectors only through gravity. We call it as the purely gravitational DM (PGDM). Interestingly, even in the case where DM interacts only through gravity, there are several processes which contribute to the DM particle production.
The so-called gravitational particle production often refers to the particle creation due to the expansion of the Universe, or the time dependence of the cosmic scale factor [1,2]. The classic paper [3] studied the scalar particle production during the transition from inflation to the matter dominated (MD) or radiation dominated (RD) Universe and found that number density of n χ ∼ H 3 inf is produced, where H inf denotes the Hubble scale during inflation, if the scalar particle χ is effectively massless or m χ H inf with m χ being the mass of χ. Refs. [4,5] considered the gravitational particle creation at the end of inflation as a mechanism for the production of supermassive DM. More recently, Refs. [6,7] pointed out that gravitational particle production always happens during the inflaton oscillation era after inflation, since the cosmic scale factor is also oscillating. There it is found that the production is most efficient at early epoch and the produced number density is n χ ∼ H 3 inf , similar to that found in Ref. [3]. It becomes clear that this estimation applies to m χ as large as the inflaton mass m φ , even if it is much larger than H inf . It may be easily understood once one regards the gravitational particle production due to the oscillating background as the "gravitational annihilation" of inflaton. #1 On the other hand, Refs. [11][12][13][14] studied the production of PGDM through the annihilation of standard model particles in thermal bath with s-channel graviton exchange. It may also be regarded as a kind of gravitational particle creation. The efficiency of DM production is dominated at earlier epoch or high cosmic temperature. However, in the very early epoch the Universe is dominated by the (coherently oscillating) inflaton field, hence the DM production is expected to be dominated by the annihilation of inflaton rather than particles in thermal bath in most cases [13].
In this paper we revisit the gravitational particle production in a realistic situation. It is often misunderstood that the gravitational particle production is not efficient if m χ H inf . It is not always true, however. In many inflation models, there is a large hierarchy between the inflaton mass scale and Hubble scale: m φ H inf . Gravitational particle production is efficient even for m φ > m χ H inf . This point is not stressed in literature except for a few works [6,7]. Therefore we want to study the gravitational production in a comprehensive manner.
This paper is organized as follows. In Sec. 2 some basics of a scalar field in the expanding background are described and equations to estimate the gravitational particle production are derived. In Sec. 3 we (semi)-analytically and numerically evaluate the gravitational particle #1 Refs. [8,9] pointed out that there is a tachyonically resonant production if χ has a (large) non-minimal coupling to gravity. See also Refs. [10]. production rate in a realistic inflationary cosmology. Sec. 4 is devoted to summary and discussion.
2 Scalar field in cosmological background

Model and equations of motion
Let us consider an action where M P is the reduced Planck scale, R is the Ricci scalar, φ denotes the inflaton field with V (φ) being its potential and χ denotes a real scalar field. It has a Z 2 symmetry under which χ changes its sign, and hence χ is stable and is a candidate of DM. We assume that χ does not have a direct coupling to the inflaton and other standard model fields. It interacts only through the metric or the gravity. The coupling strength to the gravity is controlled by the non-minimal coupling ξ. Pure Einstein gravity corresponds to ξ = 0 and the conformal coupling corresponds to ξ = 1/6. We use the Friedmann-Robertson-Walker metric: where a(τ ) denotes the cosmic scale factor with τ being the conformal time, which is related to the physical time as dt = adτ . Defining χ ≡ aχ, the action of χ is given by where the prime denotes the derivative with respect to τ . Thus χ satisfies the equation of motion Treated as classical background, φ has the following equation of motion, where the dot denotes the derivative with respect to the physical time t and the Hubble parameter H is given by the Friedmann equation, with the χ contribution to the energy density neglected. Thus for any given inflation model, we can calculate the production rate of χ through the time dependence of the scale factor a in (3). These equations are written in terms of the conformal time as where the conformal Hubble parameter H is The Friedmann equation of the second kind is given by

Quantization and adiabatic vacuum
Since m (eff)2 χ is time-dependent in the expanding Universe, we should be careful about the choice of mode function and vacuum state. Let us define the Fourier mode as where χ k = χ −k must be satisfied from the reality of χ, and the ladder operator satisfies The Fourier mode satisfies the equation of motion: From the canonical commutation relation we obtain the normalization condition The vacuum state |0 is defined as a k |0 = 0 for some mode function χ k at some initial time τ = τ i . In the Heisenberg picture, the state does not evolve once we fix it at the initial time. Instead, the mode function develops with time, which may be interpreted as the particle production as will be shown later.
Given initial conditions, one can directly solve (12), but here we use a different technique [18]. Let us rewrite χ k (τ ) as follows: Since both α k (τ ) and β k (τ ) are time-dependent, one can always write χ k in this form. Still there is a degree of freedom for the choice of α k (τ ) and β k (τ ). One can impose the following condition, which is consistent with the equation of motion (12): Instead of solving (12), one may solve (16), which is often much easier for the purpose of evaluating the particle production numerically. Note that (14) ensures In order to extract the number density from α k and β k , we must first specify the initial condition for the mode function (which corresponds to the choice of the vacuum state) and the observer state which counts the number density. In general, there is no preferred choice for the states of the vacuum and the observer in curved spacetime. In our case, however, we may formally assume that the spacetime is asymptotically static in the far past τ → −∞ (deep in the inflationary era) as well as in the far future τ → +∞ (deep in the MD or RD era). In such a case, it is natural to take the vacuum/observer as the negative frequency modes in the limit τ → −∞/∞, respectively. The negative frequency mode approaches to in the limit τ → −∞. Thus we take α(τ i ) = 1 and β(τ i ) = 0 in the limit τ i → −∞ as the initial condition. Note that it represents the adiabatic vacuum of the infinite order [2], since the spacetime is assumed to be static in the far past/future regions. #2 In our numerical calculation, however, it is of course impossible to run from τ i = −∞, and hence we start our numerical calculation with α(τ i ) = 1 and β(τ i ) = 0 at some large but finite τ i . We will discuss how to infer the result with τ i → −∞ from our numerical results with finite τ i in the next section.
Here is a comment on the size of m 2 χ and ξ. The exact solution to (12) during the (pure) de Sitter era with the initial condition (18) is given by #2 One can show that our definition of f χ is equivalent to the number density defined by the Bogoliubov coefficient where the vacuum and the observer states are taken as the zero-th order adiabatic vacuum as long as ω k /ω k = 0 at τ = τ i and τ f , where τ f is the conformal time at which the number density is evaluated. Since the adiabatic expansion is exact in the limit in the far past/future regions, our f χ coincides with the number density defined by the adiabatic vacuum of the infinite order for −τ i , τ f → ∞ as well.
where H (1) ν denotes the Hankel function of the first kind. In this paper we mainly concentrate on the case of ν 2 < 1/4. Otherwise long wave superhorizon fluctuation develops and χ obtains a homogeneous classical field value, #3 and we must take account of the coherent oscillation of the homogeneous χ field for estimation of the final DM abundance. For the small nonminimal coupling |ξ| ≤ 1/6, ν 2 < 1/4 roughly means that χ is heavier than the Hubble scale during and after inflation. For the nonminimal coupling ξ ≥ 1/6, ν 2 < 1/4 always holds even if m 2 χ H 2 as far as m 2 χ is positive.

Energy and number density
In order to evaluate the amount of created particles in the cosmological background, one must define the energy or number density. Of course, they are actually UV divergent and hence we need to carefully renormalize them. The energy momentum tensor of χ is given by : where G µν denotes the Einstein tensor and partial derivative is taken with respect to physical time t. The energy density may be defined through ρ χ = T χ 00 . Using (10), we obtain If we naively introduce a physical UV cutoff scale Λ, we find terms like ρ χ ∝ Λ 4 , H 2 Λ 2 and H 4 log Λ. These divergences are renormalized by the constant term, Planck constant and the coefficient of the R 2 term in the action, respectively [2]. #4 Since we want to calculate the energy density of χ in the far future in which H → 0 and H → 0, one can safely define the renormalized energy density as #3 The long wavelength limit is H > 0 always holds if ν 2 < 1/4 during inflation and hence ω k is real for all the wavenumber k.
#4 As for the quadratic divergence, one can explicitly check that a term likeḢΛ 2 is cancelled out.
where f χ (k) denotes the phase space density of the χ particle with physical momentum k/a. Similarly, the number density is given by Therefore, given an cosmological model, we can calculate the evolution β k by using (16) and then estimate the number of produced particles. In particular, in the most cases of our interest β k is so small that we can safely take α k 1. Then we have As we discussed in the previous subsection, we take the limit τ i → −∞, and evaluate f χ at (23) is rigorous, assuming that the spacetime is static in these regions.
3 Production of purely gravitational dark matter 3.1 Analytic estimation

Conformal coupling
Now let us evaluate β k through (25). In the following we always assume m χ m φ . First, let us consider the case of conformal coupling ξ = 1/6 for simplicity. Of course, in the massless limit m χ = 0 the χ field does not feel the Universe expanding due to the conformality and no particle production occurs. However, for finite m χ , particle production happens. We can approximate Ω k as where we have ignored additional irrelevant constant. The integrand is given by The integrand approaches to zero both in the limit τ → −∞ (deep in the inflationary era) and τ → +∞ (deep in the MD or RD era). Thus the integration (25) is well-defined, where τ k is defined as the conformal time at which a(τ k )m χ = k. Now we need the time dependence of the scale factor a(t) for evaluating this integral. As explicitly noticed in Ref. [6], the scale factor a(t) contains rapidly oscillating part during the inflaton coherent oscillation, as where v φ is the final VEV of inflaton and a(t) is the time-averaged value over the inflaton oscillation period. The particle production may be caused by both the time-dependence of the averaged value a(t) and also by the rapidly oscillating part. We stress that such a clear separation is possible only in the harmonic inflaton oscillation regime and distinction between "slowly changing" part and "rapidly changing" part is not well defined in general. Nevertheless, this picture helps the analytic understandings of the particle production. First let us focus on the effect by the time-averaged part a(t) . We discuss this contribution in detail in App. A and the resultant phase space distribution is given by (56) and (57): where c ∼ 5 and a end denotes the scale factor at the end of inflation and k c is defined by the momentum so that m χ t k = 1.
Next we consider the effect from the rapidly oscillating part in (29). The particle production due to such a rapidly oscillating effective mass term is analogous to the process known as the narrow resonance due to the oscillating scalar field [15][16][17][18]. In our situation the parametric resonance effect is neglected since the Hubble expansion is fast enough to expel the created particles from the resonance band in the momentum space [6] and in such a case the production rate is the same as the perturbative decay/annihilation rate. The effective inflaton annihilation rate in this case is if m χ m φ , where Φ denotes the oscillation amplitude of the inflaton field and C is a numerical factor. Thus the created number density during one Hubble time is given by One immediately notices that the dominant contribution comes from the earliest epoch, i.e., H ∼ H inf . The typical physical momentum of χ is k/a ∼ m φ for each epoch. The final momentum distribution of χ is not exponentially suppressed for large k since the production continues until the inflaton decays. We obtain the following approximated form: and there will be an exponential cutoff at k ∼ a(t = Γ −1 inf )m φ where Γ inf denotes the total decay width of the inflaton. #5 This contribution is not suppressed even for m χ H inf . Note also that the low momentum behavior of (33) may not be so simple because of the nontrivial time scale of the inflaton dynamics during the transition from inflation to the reheating era. The final momentum distribution is roughly the sum of (33) and (30). In realistic situation, as we will see later, these two contribution should be smoothly connected because the "end" of inflation is rather vague and the inflaton oscillation is far from harmonic at the early stage of the reheating and its typical time scale changes from H inf to m φ gradually. The number density is then given by where we will numerically find C ∼ 10 −3 − 10 −2 and with A ∼ 6 × 10 −4 for m χ H inf . Thus the contribution from the pure Hubble expansion likely to dominate for m χ < H inf while that from the oscillation effect (33) becomes important for m χ > H inf . In the latter case, we can evaluate the present DM energy density from the gravitational production divided by the entropy density as for m χ < m φ , where T R denotes the reheating temperature after inflation and we assumed that the inflaton coherent oscillation behaves as non-relativistic matter.
#5 For the completion of the reheating, we need to introduce inflaton coupling to some other sector.

Minimal coupling
Now let us consider the case of minimal coupling, ξ = 0. In this case, we have and we obtain ω k 2ω k 1 2 where we have assumed m χ H inf to avoid the growth superhorizon modes as already mentioned. The first term is the same as the conformal coupling case and hence we focus on the second term hereafter. In the following we also assume m χ m φ .
The Ricci scalar R is expressed in terms of φ as (9) and hence it is a rapidly oscillating function when the inflaton oscillates coherently. Therefore we can divide R into its oscillation averaged part and the oscillating part, schematically as First let us consider the effect due to the averaged part. After the partial integration, the integral (25) is expressed as It may be dominated around the transition τ ∼ τ end where a 2 R/ω k takes its maximum value. Thus we roughly have #6 Next we consider a particle production due to the oscillating part in R. As already mentioned in the case of conformal coupling, such a rapidly oscillating term induces particle production similar to the standard narrow resonance. On these grounds, the effective inflaton annihilation rate during the inflaton coherent oscillation is given by #6 Even for m χ H inf , the result for the high momentum mode k > a end H inf looks similar to that for the high momentum limit of the case of m χ H inf . These high momentum contribution goes as n χ ∼ H 3 inf (a end /a(t)) 3 , as given in Ref. [3], but the long wavelength contribution can be dominant in such a case. See Sec. 4. if m χ m φ where C is a numerical factor and Φ denotes the inflaton oscillation amplitude. Thus the created number density during one Hubble time is given by Again, the dominant contribution comes from the earliest epoch, i.e., H ∼ H inf . Similarly to the case studied in the conformal coupling case, the typical physical momentum of χ is k/a ∼ m φ for each epoch and the final momentum distribution of χ is not exponentially suppressed for large k since the production continues until the inflaton decays. We expect that the momentum distribution looks like and there is an exponential cutoff at k ∼ a(t = Γ −1 inf )m φ . It is not suppressed even for m χ H inf . Note again that the low momentum behavior of (43) may not be so simple because of the nontrivial time scale of the inflaton dynamics during the transition from inflation to the reheating era. As mentioned above, there are another contribution as (40), (33) and (30). Again we stress that, in realistic situation, these contributions should be smoothly connected because we cannot strictly define the "end" of inflation and the typical time scale of the inflaton motion changes from H inf to m φ gradually around the transition epoch. The number density is dominated by that from the oscillation effect (43) for m χ H inf and two contributions are comparable for m χ H inf . In any case, the number density is given by for H inf < m χ < m φ where we numerically find C ∼ 10 −2 . The present DM energy density from gravitational production divided by the entropy density is then given by for m χ < m φ and T R denotes the reheating temperature after inflation and we assumed that the inflaton coherent oscillation behaves as non-relativistic matter.

Numerical simulation in realistic inflation model
Let us now calculate the phase space density numerically. As a concrete inflation model, let us consider the hilltop inflation or new inflation models [19][20][21][22][23][24][25]. The inflaton potential is given by where n ≥ 6 is favored from the recent cosmological data [26]. In a numerical calculation we take n = 6. This model is of our interest because it in general predicts a large hierarchy between m φ and H inf . Actually we have The field value at the end of inflation is defined by where the slow-roll parameter ≡ (M P ∂V /∂φ) 2 /( √ 2V ) 2 becomes equal to unity. The dimensionless power spectrum of the large scale curvature perturbation is given by where N e 50-60 denotes the e-folding number at which the perturbation with the present horizon scale exits the horizon. The Planck observation suggests P ζ 2.2 × 10 −9 [28]. We numerically solved (16) along with the background (7) and (8). Fig. 1 shows the resulting phase space distribution f χ (k) for the conformal coupling case ξ = 1/6 (left) and the minimal coupling case ξ = 0 (right). We have taken v φ = 0.5M P , which means m φ /H inf 29, and M is taken to satisfy the observed density perturbation. Three lines correspond to m χ = (0.2, 0.5, 1, 2) × m φ , respectively. The initial condition is set to be φ(τ i ) 0.4φ end with where a(τ i ) can be arbitrary value. The final evaluation time is taken to be τ f = −τ i . The wave number in the horizontal axis is normalized by a end m φ where a end denotes the scale factor at the end of slow-roll regime φ = φ end . The results are consistent with (33) for the conformal coupling case ξ = 1/6 and (43) for the minimal coupling ξ = 0 for the high momentum mode k a end m φ . In particular, the overall size of f χ is almost the same as long as m φ m χ for ξ = 0, and hence we conclude that the particles with m φ m χ are indeed produced even if m χ H inf . On the other hand, the overall size of f χ decreases as m χ increases beyond m φ even for ξ = 0 (see the m χ = 2m φ line in the right panel). The lower momentum behavior comes mainly from the transition epoch from inflation to the oscillation epoch, which is difficult to discuss analytically, but the estimations (36) and (45) remain valid because the number or energy density is dominated by the mode around k ∼ a end m φ . #7 Now we comment on subtlety in evaluating f χ (k). As we discussed in the previous section, we formally define the particle number density by taking −τ i , τ f → ∞. In our numerical calculation, however, we take some finite τ i and τ f due to the limitation of computational cost, and hence we should carefully subtract effects of finite τ i and τ f . More precisely, we want to take the limit −kτ i , kτ f → ∞, since τ i itself can be taken arbitrarily. For   fixed physical momentum k/a end , such a limit can be achieved by taking the duration of inflation long in a numerical calculation. In order to identify such effects, we plot f χ (k) for different initial condition in Fig. 2. Three lines correspond to different initial condition, φ i = (0.38, 0.41, 0.44) × φ end for "long," "mid" and "short," respectively. The DM mass is taken to be m χ = m φ (left) and m χ = 2m φ (right). The wave number in the horizontal axis is normalized by a end m φ . As seen from the left panel, the flat part at large k is initial condition dependent and becomes smaller as the initial time is taken to be earlier. The following simple example may be helpful to understand this behavior. Let us consider the integral: where τ 0 is a real number. In the limit −τ i = τ f = ∞, we can exactly solve it by using the residue theorem and obtain exponential form: I(k) = πe −kτ 0 /τ 0 . If we take large but finite −τ i = τ f ( τ 0 ), we instead have power law tail as I(k) ∼ πe −kτ 0 /τ 0 +1/(τ 2 f k). The integrand is more complicated in a realistic situation, but we expect that a similar phenomena occur in our numerical calculation. Therefore the flat part at large k is interpreted as an effect of finite τ i and τ f , and we expect that it disappears for −τ i , τ f → −∞. On the other hand, the modes with smaller k are not affected by the change of the initial condition, and hence are expected to be intact in the limit −τ i , τ f → ∞. We have checked that the result in Fig. 1 is not affected by the change of the initial condition, and hence we expect that it provides good estimation of f χ in the limit −τ i , τ f → ∞. This issue is related to the choice of the initial condition at τ = τ i . If we could carefully choose the initial conditions of α k and β k at τ = τ i so that they match the solution with α k = 1 and β k = 0 at τ = −∞, the τ i dependence would be gone. The phase space density of χ after the gravitational particle production for the minimal coupling ξ = 0. The three lines "short", "mid" and "long" correspond to the duration of inflation in a numerical calculation. We have taken m χ = m φ (left) and m χ = 2m φ (right).
It is worth noting here that we do not actually know the initial condition of the Universe. Inflation should last at least for ∼ 60 e-foldings, but the dynamics before the "observable" inflation is unclear. Therefore, strictly speaking, taking the limit τ i → −∞ may not be justified in the real Universe. Similarly, the present Universe has of course finite age, and we cannot formally take τ f → ∞. We expect that taking these limits will give good estimation practically, but it would be interesting that information of the initial condition of the Universe may be contained in the momentum distribution.

Discussion
We studied the gravitational production of scalar PGDM by tracing the evolution of its wave function from inflation to the reheating era. We confirmed that PGDM mass with m χ H inf can be efficiently produced as long as m χ m φ , where H inf and m φ denote the Hubble scale during inflation and the inflaton mass, respectively. Roughly, PGDM number density of ∼ H 3 inf is produced even for H inf m χ ( m φ ) at the end of inflation and the beginning of the inflaton oscillation. However, we note that our results depend on the assumption that the inflaton field remains homogeneous over the horizon scale until the inflaton coherent oscillation becomes purely harmonic. In actual situation, there can be some instability that tend to make the inflaton field inhomogeneous. Although our results may not be affected much unless the inflaton field becomes highly relativistic, further investigation is needed in order to correctly estimate the PGDM abundance in various inflation models.
So far we have focused on the case of m χ H inf in order to avoid the growth of large scale fluctuation during inflation. If m χ H inf , the superhorizon modes of the χ field are generated during inflation, #8 and it will be effectively regarded as the homogeneous mode, which results in the coherent oscillation of χ field after inflation. This is another source of the DM production. Typical field variance averaged over the superhorizon mode is derived from the solution (19) as [29,30] In a horizon patch, the χ field is almost homogeneous with its typical value given by χ 2 . The χ field then begins a coherent oscillation when H ∼ m χ . The abundance of the coherent oscillation is estimated as Here we assumed that χ begins to oscillate before the reheating is completed: Γ inf < m χ . Otherwise, the abundance is suppressed by the factor ∼ m χ /Γ inf . This also contributes to the DM density. One crucial difference from the gravitationally produced contribution ρ (GP) χ is that the large scale fluctuation of χ has an (uncorrelated) isocurvature perturbation, which is severely constrained from observation. #9 The magnitude of DM isocurvature perturbation is estimated as where R χ ≡ ρ (CO) χ /ρ DM denotes the fraction of χ coherent oscillation energy density in the total DM energy density. The observational constraint is S DM 9 × 10 −6 [28]. Note however that (50) is an asymptotic averaged value when the inflation lasts long enough and it may be possible to have larger/smaller field value in the actual Universe if m χ H inf . It is also affected by other terms like L ∼ −λχ 4 . Scattering of standard model particles in thermal bath also produces PGDM through the graviton exchange [11][12][13][14]. The cross section for the process like ψψ → χχ, where ψ collectively denotes the standard model fields that are in thermal bath with temperature T , is #8 The condition of the growth of long wavelength mode in the dS era is ν 2 > 1/4 in (19), or m 2 χ < (2 − 12ξ)H 2 inf . #9 Note on terminology: coherent oscillation may also be regarded as gravitational production, since (50) is induced gravitationally. We just call such long wavelength (superhorizon modes) contribution as coherent oscillation and short wavelength (subhorizon modes) one as gravitational production. The long dashed red line shows the isocurvature perturbation limit and region above this line is excluded. The dominant production mechanism is coherent oscillation for m χ < √ 2H inf and gravitational production for for T m χ , A 1/24π, 1/48π and 1/12π for complex scalar, Dirac fermion and massless vector, respectively. The final DM abundance from thermal production is dominated by those produced in the earliest epoch and estimated as where T max ∼ (T 2 R H inf M P ) 1/4 is the maximum cosmic temperature of the dilute plasma after inflation. Compared with the inflaton-induced gravitational contribution (45), the thermal production contribution is suppressed by T 4 R /(H 2 inf M 2 P ). Thus thermal contribution is subdominant in most cases, and some DM scenarios implementing thermal production solely were discussed in Ref. [31][32][33][34][35], for examples. However, it is possible that χ is much heavier than the inflaton, m χ m φ , so that all the gravitational production is not effective, while T max > m χ and hence thermal production is active. In such a case, thermal production can give a dominant contribution to the PGDM abundance.
In Fig. 3, we show the total contribution for PGDM's relic abundance Ω χ as functions of H inf and DM mass m χ . We have chosen the reheating temperature T R = 10 9 GeV (10 7 GeV) and inflaton mass m φ = 10 12 GeV (10 13 GeV) for minimal coupling ξ = 0 in the left (right) panel. From top to down, different lines (dot-dashed, solid, dashed, and dotted) correspond to Ω χ = (4, 1, 0.1, 0.01) × Ω DM . The long dashed line marks the isocurvature perturbation limit (52) and region above this line is excluded. Note that the dominant production mechanism is coherent oscillation for m χ < √ 2H inf and gravitational production for m χ > √ 2H inf . In the shown parameter space, the thermal production can always be neglected. DM in the shown region might be probed by cosmic ray and CMB experiments if there is (small) Z 2 breaking term and it decays, see Ref. [12,36].
Finally we briefly comment on the case of large nonminimal coupling ξ. Compared with the case of minimal coupling ξ = 0, the PGDM abundance is expected to scale as (ξ − 1/6) 2 . However, for very large ξ m 2 φ /H 2 inf , the χ field feels large tachyonic mass because the Ricci scalar is oscillating between positive and negative value during inflaton oscillation and low momentum modes are exponentially enhanced due to the tachyonic instability [8][9][10]. In such a case the final PGDM abundance may be greatly enhanced.
around t t k and t t end are comparable. As a result, f χ (k) = |β k | 2 is given as where we have introduced a numerical factor c. We will see that c ∼ 5 fits the numerical results well. Next let us consider the high momentum mode k > a end m χ . In this case we have t k > t end and we obtain where the first (second) term corresponds to the integral around t ∼ t end (t ∼ t k ), and k c ≡ a end H inf (m χ /H inf ) 1/3 is defined by the momentum so that m χ t k = 1. If m χ < H inf , the second term mostly dominates. We have omitted the mixing term since it is subdominant in most cases. Intuitively, modes with k > k c are always adiabatic against the expansion of the Universe. If m χ > H inf , the first term dominates and all the modes are suppressed by at least the factor ∼ e −mχ/H inf , since they are all adiabatic throughout the whole history of the Universe and no significant excitations are expected. As a result, the number density due to the gravitational production in this case is approximated as where we find that the numerical factor A is weakly dependent on m χ and A ∼ 6 × 10 −4 for m χ H inf . We solved (16) with the Friedmann equation H 2 = ρ φ /(3M 2 P ) by using the toy model (55) for the conformal coupling ξ = 1/6 to estimate the gravitational particle production rate. The result of f χ (k) is shown in Fig. 4 for m χ = (1, 0.1, 0.01) × H inf . The results are well fitted by the formula (56) and (57) with c = 5.