Gravitational production of scalar dark matter

We investigate the gravitational production of scalar dark matter particles during the inflationary and reheating epochs. The oscillatory behavior of the curvature scalar R during the reheating phase generates two different enhancement mechanisms in the particle production. On the one hand, as it has been already discussed in previous works, it induces tachyonic instabilities in the field which are the dominant enhancement mechanism for light masses. On the other hand, we have found that it also provokes a resonant effect in the ultraviolet region of the spectrum which becomes dominant for masses in the range 109 GeV to 1013 GeV. We have developed an analytical approximation to describe this resonance effect and its consequences on the ultraviolet regime. Once we have calculated the theoretical gravitational production, we constrain the possible values of the phenomenological field parameters to be considered as a dark matter candidate. We do so by comparing the theoretically predicted abundance with the observed one and ensuring that the theoretical prediction does not lead to overproduction. In particular, we find that there is a region of intermediate masses that is forbidden as they would lead to overproduction.


Introduction
The nature of dark matter has been an open question since F. Zwicky first proposed its existence to explain the dynamics of the Coma cluster galaxies [1]. In the last decades, dark matter has become a key ingredient to explain cosmological observations. However, we are still lacking a fundamental description of its nature. The community has followed several approaches to understand it with special effort in looking for extensions of the Standard Model (SM) of particles as many beyond SM proposals include new different fields that could be potential candidates to explain dark matter [2]. However, without any conclusive experiment so far, there is great uncertainty in the intrinsic properties of the dark matter candidates as each theoretical proposal predicts different values for its parameters. The only experimental certainty we have so far is that the cross-section of dark matter with the SM fields must be very small [3,4]. Hence, the question of how it was produced arises naturally. Within this work we study the production of particles due to gravitational effects during both the inflationary and post-inflationary epochs. We focus on discerning if this mechanism can provide enough dark matter density to explain the observed abundance. In particular, we study the case of a scalar field with a non-minimal coupling to the Ricci curvature scalar.
Many different authors have studied the particle production due to gravitational effects in detail. The first works in this arena focused on the production of elementary particles due to the expansion of the Universe [5,6]. The classic article [7] analysed the gravitational production of particles after the phase transition to a radiation dominated cosmology after inflation, obtaining a general result for the number density of created particles which is of the order of the energy scale of inflation cube (n ∼ H 3 0 ) if the particles can be considered effectively massless during inflation, i.e., if the mass is much smaller than the JHEP06(2020)084 energy scale during this phase (m H 0 ). More recent papers [8][9][10] have considered the gravitational production as a possible mechanism to produce enough supermassive dark matter particles (WIMPZillas), obtaining then a lower bound for the possible dark matter mass in this scenario. All these works have neglected the importance of the oscillations of the background quantities due to the inflaton dynamics and their impact on the evolution of the quantum field. However, the effects of these oscillatory behaviors on the gravitational production can be important. In [11,12], they analyse the impact of the scale factor oscillations on the gravitational production, and the works [13,14] study the tachyonic instability induced on the field by the oscillations of the scalar curvature. There are some recent works which deals with the gravitational production of self-interacting dark matter during the early Universe, imposing some tight constraints on the possible values of the self-coupling parameter [15,16].
In this work we have focused on a scalar field with negligible interactions with the SM but with a direct coupling to the scalar curvature through a term in the action ξRϕ 2 . We start with an initial de Sitter phase which mimics the inflationary dynamics and where the initial vacuum state for the dark matter field is well defined for values of the coupling to the scalar curvature ξ ≥ 1/6. For couplings ξ < 1/6 it is known that in de Sitter geometry the vacuum state is unstable [6]. On the other hand, ref. [17] suggests that ξ should not exceed the value 10. Furthermore, analyses based on quantum cosmology considerations suggest that it should be of order 1 (between 1/6 and 1/3, to be more precise) [18]. For these reasons together with the numerical difficulties encountered when computing the production for higher values, we will restrict the range of the parameter ξ to lie between the values 1/6 and 1. Then we have taken into account the oscillatory behavior of the inflaton during the reheating phase and its implications for the curvature scalar. We have seen that besides the already studied tachyonic instability, there is a resonance between the frequency of the field and the frequency of the curvature scalar oscillations that significantly enhances the production of ultraviolet modes for the scalar field. We study this phenomenon using an analytical approximation which allows us to describe only the ultraviolet modes. Comparing the obtained produced gravitational abundance with the observed dark matter abundance, we have set constraints on the possible parameters of the dark matter field in order not to be overproduced. It is important to note that dark matter produced by this mechanism induces adiabatic perturbations instead of isocurvature perturbations as is discussed in [14,19]. Therefore the stringent constraints coming from the CMB observations on the amplitude of isocurvature perturbations [20] do not apply to our study.
This manuscript is organized as follows: in section 2 we introduce the model we have considered throughout this work. It consists on a massive scalar field non-minimally coupled to gravity through the Ricci scalar. We are neglecting the back-reaction of the field on the background as its energy density will be small as compared to the one of SM degrees of freedom. Hence, dark matter is considered as a test field on a given background spacetime described by the Friedman-Lemaître-Robertson-Walker metric. We are interested in the cosmological production during the early Universe, specifically during the inflationary and reheating epochs. In order to have a well-defined vacuum, we have mimicked the infla-JHEP06(2020)084 tionary epoch with the de Sitter solution. Then we consider the reheating dynamics from the oscillations of a massive inflaton which leads to a matter dominated cosmology. We also define what we understand by gravitational production of particles. This production mechanism is based on the time evolution of the vacuum state and the fact that the vacuum defined by an observer at each instant of time does not necessarily match with the evolved vacuum. However, these two states are related by a Bogolyubov transformation and the coefficients of this transformation are related to the produced particle density. In general the mode equation can not be solved analytically so, in section 3 we have performed a numerical study of a representative section of the field space of parameters, computing the spectral production and the particle density, which will be used later to set constraints on the dark matter parameters based on the predicted abundance. In section 4 an analytical approach is followed to understand the resonance mechanism for the ultraviolet modes that we have observed in the numerical results. Finally, in section 5 we show and discuss the numerical results for the production spectra of particles at the end of reheating, and the dark matter abundance that is observed nowadays. We discuss the constraints on the parameter space of the scalar field in order to have a viable dark matter candidate.

Background, field dynamics and gravitational production
We consider a massive scalar field ϕ as a dark matter candidate. Its dynamics in a given curved spacetime is encoded in the action where g is the determinant of the spacetime metric (with signature (−+++)), R is the Ricci scalar curvature, ξ is the non-minimal dimensionless coupling constant to gravity, and m is the mass of the dark matter field ϕ. In this action we have not included any direct coupling of the dark matter field to the SM sector due to the experimental constraints on the corresponding cross-sections. However, as long as these interactions exist, no matter how weak they are, the renormalization group flow will generate a non-minimal coupling to spacetime curvature [21,22] and therefore we have included it in the action. Note that this action shares some common terms with the action for Higgs inflationary models due to the presence of the non-minimal coupling ξ. However, there is an important difference between them as the Higgs field has quartic self-interactions in its potential while our action is just quadratic, i.e., it only includes the mass and non-minimal coupling terms. This quartic interaction term is definitively a game changer in the sense that makes the whole dynamics entirely different and hence the role that these two different fields play in the cosmological models. In cosmology, spacetime geometry is well described in average by the FLRW metric. For later convenience, we write it in conformal time η. Furthermore, both for simplicity and in accordance with observations [23], we are considering flat spatial sections. Therefore the spacetime geometry will be described by

JHEP06(2020)084
where the functional form of the scale factor a(η) is determined by the energetic content of the Universe during each cosmological epoch. From the point of view of cosmological particle production, the relevant epochs are inflation and reheating. We model the inflationary epoch using the de Sitter solution which is a good first approximation to single-field inflation except for the transition to reheating. On the other hand, during reheating, the behavior of the scale factor is determined by the potential of the inflaton field φ near its minimum value. In the case at hand the potential under consideration is with m φ being the mass of the inflaton field. With this potential, it is well known that the background geometry will behave in average as a matter dominated Universe [24]. For our study, it suffices to deal only with the average behavior of both the scale factor and the Hubble parameter, for which we keep the symbols a(η) and H(η). On the other hand, we will take into account the oscillating terms in the curvature scalar sourced by the dynamics of the inflaton. As we will discuss, the gravitational particle production is enhanced by the instabilities provoked by these oscillations in the effective mass term of the field. If we set the end of inflation at η = 0, the averaged scale factor for these two epochs can be modelled as [25] a(η) = where H 0 = √ 12πm φ is the value of the Hubble parameter at the onset of the single-field inflationary epoch for the model we are considering (2.3).
Although this parameterization provides both a continuous scale factor and a continuous Hubble parameter H(η) for all values of η, it fails to provide a continuous Ricci scalar at η = 0. For the approximations we are going to use to solve the mode equations for the field we will need R to be at least a C 4 function. To construct such a function we used an eight order polynomial to interpolate between the value of R at the end of the de Sitter phase and the one coming from the dynamics of the inflaton field during reheating. We start the interpolation then at η = 0 and finish it at η osc , which is defined as the time in which the Hubble parameter from the background parameterization (2.4) coincides with the averaged one from the considered dynamics of the inflaton.
The dynamics of the inflaton is governed, in cosmological time (related to conformal time via dt = a(η)dη), by the equation of motion Inflation ends when H(t) ∼ m φ and hence during reheating the friction term will be subdominant with respect to the mass term. Henceforth, the solution for the inflationary field can be approximated at first consistent order in a Wentzel-Kramers-Brillouin (WKB) expansion as

JHEP06(2020)084
where Φ 0 is the amplitude of the field at the end of inflation and depends on the specific inflationary model. In the case of massive chaotic inflation that we are analysing, the observations of the CMB fluctuations [23] set the following values for both the mass and the amplitude of the inflaton field at the end of inflation: where κ 2 = 8πM −2 P , and M P is the Planck mass M P 1.22 × 10 19 GeV. The dynamics of spacetime during the considered epochs is dominated by the energy density of the inflaton field. Hence, the behavior of the curvature scalar during reheating is given by the trace of the Einstein equations sourced by the energy-momentum tensor of the inflaton: Using the solution for the inflaton field (2.6), and keeping all the terms up to the first order WKB approximation, the curvature scalar can be written as (2.9) The resulting curvature scalar from the above approximations is shown in figure 1.
In cosmological scenarios, the quantization procedure suffers from an ambiguity which can be traced back to the choice of the canonical pair of variables that we are going to quantize [26]. In our case, there is a preferred choice for the canonical pair of variables if one demand that the quantum theory is invariant under the spatial symmetries of the background and that its quantum dynamics admits an unitary implementation [27]. This criterion entails that the field to be quantized must be χ = a(η)ϕ. (2.10) From the action, we can see that the rescaled field satisfies the equation of motion Exploiting the fact that the spatial sections of the background are flat, we can expand the rescaled field in Fourier modes that satisfy the equation where we have defined the time-dependent frequency as (2.14) Figure 1. Curvature scalar R(η) as a function of conformal time. During de Sitter inflation (η < 0) the Ricci scalar remains constant, which is in good agreement with all the inflationary models. The most appreciable difference between our approach and a realistic inflationary model should occur near η = 0. We have modeled the transition between the end of inflation and the onset of oscillations using a polynomial function so the scalar curvature is a C 4 function. For η ≥ η osc we have considered the curvature sourced by the oscillations of the inflaton field (2.9).
Due to the isotropy of the background, we can always choose bases of solutions of this equation with elements v k and v * k that only depend on k = |k|. The Fourier mode amplitude χ k is a linear combination of the mode functions v k and v * k with creation and annihilation variables a k and a * −k : Moreover, to preserve the standard Poisson bracket structure for the field and hence for the creation and annihilation variables, the mode functions must be normalized (we follow the convention of [28] The field quantization is carried out by promoting the field variable an its canonically conjugate momentum to operators defined through the creation and annihilation operatorŝ a k andâ † k acting on a Fock space. These operators satisfy the following commutation relations:

JHEP06(2020)084
The Fock space is generated by the action of the creation operators on the vacuum state, defined as the state that is annihilated by all the annihilation operators: In curved spacetimes the selection of the vacuum state is not unique but depends on the choice of the set of mode functions. The definition of the creation an annihilation operators depends on the choice of basis and hence it introduces an ambiguity in the quantization procedure. In general, different choices will define inequivalent quantizations for the same field. We have followed the criterion presented in [26,27] to select the field variable to be quantized, requiring the unitary implementation of its quantum dynamics.
In cosmology, the gravitational production occurs because of the time evolution of the background geometry. In these scenarios, the vacuum state for a quantum field is not stationary. Particles are produced because the evolved vacuum state at a given instant of time is in general an excited state with respect to the instantaneous vacuum defined at that instant of time.
We can define the vacuum state of the scalar field through the initial conditions given in the de Sitter phase. These initial conditions define a set of modes v k and v * k for the solution of (2.13). The general isotropic solution for the mode equation in the de Sitter phase (η < 0) can be written as: ν are the Hankel functions [29] of order We define the vacuum state of the field as the one determined by the set of mode functions that in the asymptotic past (η → −∞) behave as the positive-frequency plane waves in Minkowski spacetime. This behavior is obtained if we take c 2 = 0 in (2.19) and normalize the asymptotic expansion of v k according to (2.16). Then, modes defining the initial vacuum have the form (in the de Sitter phase, i.e. for and evolve to the reheating epoch (i.e. for η > 0) according to (2.13).
On the other hand, a comoving observer will have a different instantaneous notion of vacuum. We can use the so-called adiabatic prescription to define this instantaneous vacuum, commonly used in cosmology. This prescription is appropriate and well defined if the geometry evolves slowly enough as compared with the characteristic time scales of the field, i.e., whenever the so-called adiabatic condition is satisfied:

JHEP06(2020)084
This condition will not be fulfilled in general during the reheating phase for all the modes of the field. However, as we are interested in the production during the whole reheating phase, we only need to have a good prescription for the adiabatic vacuum after this phase has ended. Once the Ricci scalar is no longer oscillating, condition (2.22) holds and we can define the adiabatic vacuum. The modes u k (η) that define the adiabatic vacuum at the time η 0 instantaneously behave as the positive-frequency plane wave modes for the vacuum in Minkowski spacetime at η 0 : where ω k (η) is given in (2.14).
Once the adiabatic vacuum has been defined for each instant of time, we can compare the mode expansion associated with each Fock space by computing the Bogolyubov coefficients. These coefficients define the linear transformation between two different mode expansions evaluated at the same time. Comparing the evolved modes v k until the time η, defined through the initial vacuum state, with the set of modes u k associated with the adiabatic vacuum at each time η, we obtain a time-dependent Bogolyubov transformation. For computing the particle production, we are only interested in the coefficient β k relating the negative-frequency modes of one expansion with the positive-frequency modes of the other one. This coefficient can be easily calculated and turns out to be [6]: Through the computation of the expectation value for the number operator defined with the instantaneous vacuum modes in the original vacuum, the gravitational produced density of particles can be obtained directly from the β k Bogolyubov coefficient. Indeed the total number of particles with momentum k created at the time η is given by |β k | 2 /[2(2π) 3 ], which integrated over all possible momenta and divided by the spatial volume a 3 , gives where we have already used cosmological time instead of the conformal one.
As stated before, we will focus on the computation of the particle density produced at the end of reheating. Afterwards, we will evolve the computed number density taking into account that the evolution of the Universe is isoentropic after this moment. Therefore, the gravitational production of particles is negligible after reheating has ended. With the evolved number density we can compute the predicted abundance for dark matter from our model as it is explained in section 5.

JHEP06(2020)084 3 Numerical analysis
In general, the computation of the mode functions during the reheating epoch cannot be performed analytically, although as we discuss in section 4, there are certain regions of the parameter space in which we can carry out an analytical approximation.
In this section we solve (2.13) numerically for a wide range of the scalar field parameters. We then introduce the numerical solutions for the mode functions v k in (2.25), and via (2.26) we obtain the total particle production. We present and discuss the numerical results for the particle production spectra obtained in this way for different parameters of the field, namely for different values of the mass m and the coupling to curvature ξ.
Throughout the numerical exploration of the space of parameters for the field, we have found that there are four different spectral behaviors as the dynamics of the curvature scalar produces three different phenomena on the particle production as we will discuss.
On the one hand, its constant value during the de Sitter phase contributes to the effective mass of the field m 2 +12H 2 0 (ξ −1/6). If this effective mass is below the energy scale of inflation H 2 0 , the production grows but if it is above it, then the production decreases. On the other hand, the oscillations of the curvature scalar during reheating produces two new phenomena. Since the frequency of the field is oscillating, there will be resonant effects whenever the two frequencies (the frequency ω k itself and the frequency at which it oscillates) are equal. These resonances will be important mostly in the ultraviolet regions of the spectra. Furthermore, these oscillations make the field suffer from tachyonic instabilities, because the effective mass for the field becomes tachyonic for certain times. This instability affects different k-bands as the instability amplitude decreases in each oscillation. Let us now discuss how these effects affect the different regions in parameter space.
For masses much smaller than the interaction term in the sense that m 2 (ξ − 1/6)|min(R osc )|, the most important production enhancement mechanism is the tachyonic instability already studied in [14]. This instability affects mostly the infrared region of the spectrum as it is shown in figure 2. As the coupling constant ξ grows, the first infrared peak is damped because the interaction term during the de Sitter epoch makes the effective mass larger and the enhancement occurs on top of the Planckian spectral production at the end of inflation. Furthermore, the larger the interaction term, the larger the ultraviolet peaks, as the instability regions during the oscillatory regime contributes with a larger imaginary effective mass.
If the mass term is larger than the energy scale of inflation, i.e, m 2 > H 2 0 and hence larger than the interaction term during the oscillatory phase, then, the resulting spectrum does not differ much from the one obtained in the de Sitter phase. In this scenario, the coupling term increases the effective mass during the inflationary epoch, giving rise to a different Planckian spectrum with no instability-enhancements as it is shown in figure 3. Therefore, in this regime, as the coupling increases, the particle production decreases.
There is a third regime, shown in figure 4, in which the mass is below the energy scale of inflation and it is of the order of the coupling term. Within this regime, increasing the coupling constant increases the spectral production as long as the effective mass during inflation is still below the energy scale of inflation. Furthermore, there are resonances due JHEP06(2020)084   . Representative spectral particle production for masses in the regime 10 9 GeV < m < H 0 . We are showing the particular case for m = 7 × 10 11 GeV. In these cases the spectral production increases with the coupling term as the resonance effect becomes more important.
to the fact that the frequency at which the field oscillates coincides at some times with the frequency at which the frequency itself oscillates. This occurs in the ultraviolet region of the spectrum.
Finally, there is a fourth regime for masses m ∼ 10 9 GeV in which there is a transition between the region dominated by the tachyonic instability and the resonance dominated one. The typical spectra for this transition regime is shown in figure 5.
Once we have obtained the spectral production density, we can integrate it to obtain the produced particle density at the end of reheating following equation (2.26). In figure 6 this production is shown as a function of the field parameters m and ξ. In this plot, we can distinguish the four regions discussed above. The first one is situated in the lightest masses, where the particle production depends only on the value of the coupling to the scalar curvature except for coupling constants in the neighbourhood or the conformal one. The dominant production mechanism in this region is the tachyonic instability induced by the oscillations of the scalar curvature. Around m ∼ 10 9 GeV there is a transition to the resonant enhancement dominated region. In the range from m ∼ 10 10 GeV to 10 13 GeV there is a dependence of the production on both field parameters, with a monotonic growth in ξ and a monotonic decrease with m. Finally, once the mass is above the energy scale of inflation and the coupling to gravity is not too large, we recover the usual behaviour for the production, as it decays rapidly with both the mass and the coupling term for the considered range. It is worth noting that for larger couplings in this regime a resonance enhanced region would be reached eventually and the production would raise again monotonically with the coupling constant value. JHEP06(2020)084

Analytical approximations
The differential equation for the modes (2.13) can be approximated in two different scenarios: either for masses well above the energy scale of inflation m H 0 or in the ultraviolet regime of the spectrum for all the parameter space.

Large mass
In the case in which the mass of the field is well above the energy scale of inflation (m H 0 ), we can approximate the solution using a second order WKB expansion. Within this approximation, we assume an ansatz for the mode equation of the form Here, W k is determined by introducing this ansatz in the equation of motion for the modes (2.13). If we can neglect the fourth and higher order derivatives of the field frequency, i.e., we can keep the approximation up to second order, this function turns out to be (4.2) Figure 6. In this figure we are showing the particle density in logarithmic scale at the end of reheating, i.e., the contours of log 10 n(t rh ) as a function of the mass m and the coupling to curvature ξ. We can see clearly the different behaviors of the production. For the lightest masses, the enhancement mechanism is the tachyonic instability. There is a transition regime for masses of about 10 9 GeV in which there are contributions from the tachyonic instability and the resonant effect. For masses lower than the energy scale of inflation the only enhancement mechanism is the resonance between the frequency of the field and the one of the Ricci scalar. For masses above the energy scale of inflation, there are no enhancement mechanisms and the typical production from a de Sitter phase is recovered.
On the other hand, A k and B k are coefficients obtained demanding that the modes are C 1 functions in the transition from the de Sitter phase to the reheating one: with v dS k being the mode solution during the de Sitter era, given by eq. (2.21). The Bogolyubov coefficient encoding the spectral production can be straightforwardly computed by introducing the WKB approximation for the modes in (2.25), and it turns to be We should note that had we considered the lowest order of the WKB expansion, i.e., W k = ω k , the expression (4.5) would has been trivial. Indeed, the mode expansion describing the field in this case would be the same one which defines the adiabatic vacuum at each instant of time and hence the only spectral production would be the one due to the change of vacuum at the end of inflation, encoded in the B k coefficients of the mode expansion (4.1). Hence, the effects of the evolution of the background during reheating would be neglected.
We can obtain an analytic approximation to the number density of particles in this regime by integrating the spectral particle production for all the modes, and dividing by the spatial volume (see (2.26)): (4.6) In this expression it is simple to see that in order to include the non-trivial effects of the background evolution in the number density, we need to approximate the solution to (4.1) with a higher order WKB approximation that the one used to define the adiabatic vacuum (otherwise W k = ω k ). The integration in modes can be performed with any numerical method and the agreement with number density obtained from the numerical solutions to (4.1) is exceptionally good.

Ultraviolet regime
Let us now study the ultraviolet regime (k m φ ). In this case the numerical calculation shows a significant increase in the particle production at the time when the time dependent frequency ω k becomes comparable with the characteristic frequency at which it oscillates. The bare WKB approximation does not capture this resonance as shown in figure 7.
In order to analytically deal with this resonance, we have followed a different approximation scheme for this case. Let us start by defining an auxiliary set of mode functions x k defined as the difference between the exact solution and the WKB approximation: (4.7) These modes are also normalized according to (2.16) and satisfy a forced time-dependent harmonic oscillator equation where the source term turns out to be: The reason to define this new set of modes is to extract the resonant behavior explicitly in the equations of motion for the modes and convert it into a classical parametric resonance JHEP06(2020)084 induced by the source term. This source term has two different frequencies: the one from the mode function oscillatory behavior and the one of the oscillation frequency of the scalar curvature. The general solution to the equation (4.8) can be written as where x h k is a solution to the homogeneous equation and x p k is a particular solution of the full equation, b 1 and b 2 being arbitrary constants. This particular solution can be written as Before the resonance is triggered, the approximate WKB solution v WKB k approximates very well the exact solution v exact k , i.e. x k = 0. Therefore, we set b 1 and b 2 to zero. For the solution x h k to the homogeneous equation, necessary to calculate the particular solution (4.11), we can simply use the WKB approximation we were using before because we have now isolated the parametric resonant frequency of interest.
To summarise, the approximate mode functions can be written, in this refined WKB approximation, in the form

JHEP06(2020)084
where x p k is given by (4.11) and, in this expression, The Bogolyubov coefficients are linear in the mode functions, so the corrections to the β k coefficients obtained through this refined WKB approximations can be added straightforwardly: (4.14) As it is shown in figure 7 this approximation method gives analytical estimations that agree with the numerical results and capture the resonant enhancement perfectly. Using this analytical approximation we can obtain also the number density of particles by replacing (4.14) into (2.26).

Constraints for dark matter
Once we have computed the gravitational particle production at the end of the reheating phase, we can use this result to set constraints on the field parameters if it was to be a dark matter candidate. The constraints are obtained imposing that the predicted abundance equals the one obtained from the observations. In order to make this comparison, we need to evolve the computed production from the end of reheating to the present day. This evolution is computed taking into account that the scalar field is decoupled from the rest of the constituents of the Universe as we have neglected its interactions. Therefore, the evolution of its particle density is solely due to the isoentropic expansion of the Universe. We can rewrite the time dependence of the number density as a dependence on the temperature of the radiation that fills the Universe. This relation arises from the statistical interpretation of the temperature and the evolution of the particle geodesics in a cosmological background [30]. Hence, the evolved number density can be rewritten in terms of the temperature of the background radiation as where g rh * S denotes the effective entropic relativistic degrees of freedom of the Standard Model of particles at the end of reheating, g 0 * S denotes the same quantity today, T 0 is the temperature today and T rh is the temperature at the end of reheating. The actual temperature at which reheating took place is still unknown and depends heavily on the particular model considered. Hence, our constraints will have an uncertainty coming from the indetermination on the specific reheating temperature.
Note that we have evolved the number density instead of the energy density of the field. There are two advantages in doing so. On the one hand, we do not need to take into account the effective mass which depends on the curvature term and, on the other hand, we do not need to take into account when each region of the parameter space becomes non-relativistic. Hence, it is simpler to compute its abundance as the energy density today JHEP06(2020)084 Figure 8. We represent the predicted abundance for the scalar field in logarithmic scale, i.e., the contours of log 10 Ω considering the maximum allowed reheating temperature, T rh = 10 15 GeV. The white contours represents the observed abundance of dark matter Ω DM = 0.268. The allowed parameter space of the scalar field is to the left of the first white contour and to right of the second one.
is obtained using just the mass of the particles as they are non-relativistic. The predicted abundance is then obtained through the following expression: where the nought quantities are evaluated at present time and we have emphasized that the density number of particles at the end of reheating depends on the model parameters. The predicted abundance from gravitational production during the early stages of the Universe is shown in figure 8 for the maximum allowed value of the reheating temperature: T rh = 10 15 GeV. This upper bound for the reheating temperature comes from the lack of observation of primordial gravitational waves [31]. Although our constraints depend on the reheating temperature, we can set solid constraints if we consider this maximum possible value. Also in figure 8, we are showing the bounds on the parameter space for the maximum allowed reheating temperature. Note that in this figure we are not showing the region near the conformal coupling to gravity, as the computed abundance depends highly on the coupling value within this region and we have not explored it numerically in sufficient detail. In figure 9 we show how the constraints on the parameter space depend heavily on the specific reheating temperature considered. We see that there are two allowed regions to describe dark matter: one for light candidates and another one for heavy candidates. It is a general result that there is a forbidden mass gap between these two regions for physically meaningful reheating temperatures. This forbidden region is larger as the considered temperature of reheating is lower.
JHEP06(2020)084 Figure 9. The constraints on the field parameter space are represented through the contours corresponding to the observed abundance: Ω = 0.268 for different reheating temperatures. The ones obtained for a specific reheating temperature T rh = 10 15 GeV are displayed in solid line. In dot-dashed line we plot the same constraints but for T rh = 10 12 GeV and finally they are displayed in dotted line for T rh = 10 9 GeV.
We have treated the conformal case in its own and the bounds for the mass in this case are shown in figure 10 for different reheating temperatures. It is important to note that, as the production is minimal for the conformal case, the less stringent bounds will come from it. For the maximum allowed reheating temperature, we can see that for the conformal coupling case, the maximum allowed mass for the light sector would be m ∼ 10 11 GeV while for couplings far from the conformal the typical maximum allowed mass would be m ∼ 10 5 GeV. So in conclusion, the maximum allowed mass for the light candidates sector depends heavily on the coupling constant to curvature ξ due to all the enhancement mechanisms we have discussed throughout this work. Note that the maximum mass coming from the conformal case in this scenario is the upper bound on all the possible masses for the light sector while the minimum mass for this same case is the minimum mass for a WIMPzilla like particle. Furthermore, we can give an empirical law to determine which mass will be the upper bound for the created dark matter particle as a function of the reheating temperature. Fitting the particle production n to an empirical law, we found that it is proportional to the mass in the conformal case. Hence, the upper bound on the mass as a function of the reheating temperature can be written as: For large enough masses we see that the predicted abundance becomes independent of the coupling and has a maximum value for a mass around the energy scale of inflation. Depending on the temperature at which reheating ended, for high enough masses the abundance can again match the observed one in accordance with [8], giving a possible allowed superheavy region. From our computations, this minimum mass would be m 5 × 10 13 GeV for the maximum allowed reheating temperature.

Conclusions
In this work we have studied the gravitational particle production of a scalar field in the early Universe. We have considered a de Sitter phase prior to the reheating scenario, mimicking the behavior of a geometry sourced by a massive inflationary field. Using this approach we can mimic any inflationary model by adjusting the de Sitter phase parameters and inflaton mass accordingly. We have focused in mimicking a chaotic massive inflation with potential (2.3). Through this work, we have realized that the oscillations of the curvature scalar R affect the gravitational particle production of a scalar field in different ways. In addition to the already studied tachyonic instability [13,14] produced by making the effective mass term imaginary, there is a resonance mechanism that has not been studied before. This resonance affects the ultraviolet region of the spectra and dominates the particle production enhancement for m > 10 9 GeV. This resonant behavior is related to the micro-oscillations in the frequency of the field due to the Ricci oscillations, in analogy to what happens in the Mathieu equation [29]. This result is in disagreement with previous JHEP06(2020)084 works as [14] where they ignored the effect of having many oscillations of the scalar curvature and they have not taken into account the resonant effect induced by these oscillations. One of the conclusions of this work is that gravitational production can account for the observed dark matter abundance. Furthermore, this gravitational production has allowed us to set constraints on the possible values of the mass of the dark matter particle if any possible direct coupling with the SM is negligible. In particular we have found that there is a forbidden band of masses which split the dark matter into either light or supermassive candidates. The width of this band is very sensitive to the reheating temperature. We have seen that the largest possible mass for the light candidate is obtained in the conformal coupling case and that for masses above the inflationary energy scale, the production become independent of the coupling constant ξ. Hence, analysing the conformal case is sufficient to set the constraints on the possible masses for the scalar dark matter candidate. However, the most stringent constrains come from the cases in which ξ > 1/6 as the possible maximum mass is much smaller than in the conformal case. For values of the coupling constant ξ not covered in this work, i.e., ξ > 1, the gravitational production would be further enhanced, widening then the range of forbidden masses for a dark matter candidate. The most affected region would be the large mass one because the field would suffer from tachyonic instabilities for large enough values of ξ, enhancing then heavily the production and dragging the allowed lightest heavy mass to larger values. For the small and intermediate mass regions, considering larger values of ξ would increase monotonically the production, dragging the allowed values of the mass for the dark matter candidate to even lighter ones.