The Spectrum of Gravitational Waves, Their Overproduction in Quintessential Inﬂation and Its Inﬂuence in the Reheating Temperature

: One of the most important issues in an inﬂationary theory as standard or quintessential inﬂation is the mechanism to reheat the universe after the end of the inﬂationary period in order to match with the Hot Big Bang universe. In quintessential inﬂation two mechanisms are frequently used, namely the reheating via gravitational particle production which is, as we will see, very efﬁcient when the phase transition from the end of inﬂation to a kinetic regime (all the energy of the inﬂaton ﬁeld is kinetic) is very abrupt, and the so-called instant preheating which is used for a very smooth phase transition because in that case the gravitational particle production is very inefﬁcient. In the present work, a detailed study of these mechanisms is done, obtaining bounds for the reheating temperature and the range of the parameters involved in each reheating mechanism in order that the Gravitational Waves (GWs) produced at the beginning of kination do not disturb the Big Bang Nucleosynthesis (BBN) success.

However, an important difference occurs with respect to the standard inflationary paradigm, where the potential of the inflaton field has a local minimum and, thus, the inflaton field releases its energy while it oscillates, which allows particle production [24][25][26][27][28]. In contrast, for the "non oscillating" models, i.e., in quintessential inflation, where the inflaton field survives to be able to reproduce the current cosmic acceleration, a fast phase transition from the end of inflation to the beginning of kination (a regime where all the energy density of the inflaton field is kinetic) where the adiabatic regime is broken is needed in order to reheat the universe. This creates enough particles, which after decays and/or interactions with other fields, form a thermal relativistic plasma whose energy density will eventually become dominant. The mechanism of particle creation can be obtained in different ways, but the most used and the ones we will study in this work are the gravitational particle production [29][30][31][32][33][34][35][36][37][38][39][40] and the instant preheating [41][42][43][44] (see also [45] for a detailed description of both mechanisms).
Dealing with the mechanisms to reheat the universe, the question related to the bounds of the reheating temperature arises. Some works already considered the constraints for reheating in quintessential inflation models, both on instant preheating [46] and on gravitational particle production [47]. A lower bound is obtained recalling that the radiation dominated era is prior to the Big Bang Nucleosynthesis (BBN) epoch which occurs in the 1 MeV regime [48]. As a consequence, the reheating temperature has to be greater than 1 MeV (see also [49] where the authors obtain lower limits on the reheating temperature in the MeV regime assuming both radiative and hadronic decays of relic particles only gravitationally interacting and taking into account effects of neutrino self-interactions and oscillations in the neutrino thermalization calculations.) The upper bounds may depend on the theory we are dealing with; for instance, many supergravity and superstring theories contain particles such as the gravitino or a modulus field with only gravitational interactions and, thus, the late time decay of these relic products may disturb the success of the standard BBN [50], but this problem can be successfully removed if the reheating temperature is of the order of 10 9 GeV (see for instance [51]). This is the reason we will restrict the reheating temperature to remain, more or less, between 1 MeV and 10 9 GeV.
On the other hand, one has to take into account that a viable reheating mechanism has to deal with the affectation of the Gravitational Waves (GWs) in the BBN success by satisfying the observational bounds coming from the overproduction of the GWs [4] or related to the logarithmic spectrum of its energy density [52]. As we will see throughout this work, the overproduction of GWs constrains very much the value of the parameters involved in the different reheating mechanisms and also impose hard bounds in the reheating temperature.
In addition, another issue related to quintessential inflation is the possibility to explain the present abundance of dark matter. Effectively, assuming that dark matter is made of non-decaying superheavy particles only coupled to gravity which are gravitationally created during the abrupt phase transition, one can show that a certain range of mass values of the dark matter leads to a viable model overpassing all the bounds coming from the overproduction of GWs [53,54].
The manuscript is organized as follows: In Section 2 we deduce the initial condition to apply the WKB approximation and ensure that the vacuum fluctuation of a massive field coupled to gravity does not affect the classical evolution of the inflaton field. Section 3 is devoted to the presentation of our quintessential inflation model, inspired in the well-known Peebles-Vilenkin one [4], i.e., depending on two parameters and containing an abrupt phase transition from the end of inflation to the beginning of kination, and the subsequent study of its dynamical evolution. Next, in Section 4 we study both reheating mechanisms in quintessential inflation, namely via gravitational particle production and via instant preheating, obtaining bounds for the reheating temperature. In Section 5 we deal with the constraints to preserve the BBN success coming from the logarithmic spectrum of GWs and also from its overproduction during the phase transition from the end of inflation to the beginning of kination, obtaining the range of values of the parameters involved in each reheating mechanism and also more restrictive bounds for the reheating temperature. In Section 6 we consider the present abundance of dark matter, assuming that it is composed by superheavy particles conformally coupled to gravity, which are also produced during the abrupt phase transition from the end of inflation to the beginning of kination, obtaining bounds for its mass. In Section 7 we consider another quintessential inflation model with a more abrupt phase transition and we show the importance of this fact and the differences with the previous model. Finally, in the conclusions we discuss the obtained results.

Initial Conditions for Inflation and the Application of the WKB Approximation
We want to know when one can apply the WKB solution in the early universe (see for instance [55,56] in order to approximately find the modes of a field coupled to gravity. This is very important because it allows us to compute analytically important quantities such as the vacuum polarization and the energy density of the produced particles after an abrupt phase transition. To do all the analytic calculations we will consider a potential such as the one used by Peebles and Vilenkin in [4] with a discontinuity in some derivative and, thus, we can obtain an analytic expression of the reheating temperature depending on the parameters involved in the reheating mechanism (the mass of the produced particles, the decay rate, the coupling constant between the quantum field which produces the particles, the inflaton field, ...).
Therefore, first at all it is well-known that at temperatures of the order of the Planck's mass quantum effects become very important and the classical picture of the universe is not possible. However, at temperatures below M pl , for example at GUT scales (i.e., when the temperature is of the order of T ∼ 4 × 10 −3 M pl ∼ 10 16 GeV), the beginning of the Hot Big Bang (HBB) scenario is possible. For the flat FLRW universe the energy density of the universe, namely ρ, and the Hubble parameter H are related through ρ = 3H 2 M 2 pl , and, for a universe filled with radiation, the temperature of the universe is related to the energy density via ρ = (π 2 /30)g * T 4 , where the degrees of freedom for the Standard Model are g * = 106.75 (see for instance [57]). Thus, one can conclude that a classical picture of the universe would be possible when H ∼ = 5 × 10 −5 M pl ∼ = 10 14 GeV. Now we consider that inflation starts at this scale, i.e., we take the value of the Hubble parameter at the beginning of inflation (denoted by H beg ) as H beg = 5 × 10 −5 M pl , and we assume that a quantum χ-field coupled to gravity and/or to the inflaton field, which will be the responsible to reheat the universe, is in the vacuum at the beginning of inflation. If we choose the mass of the χ-field at least one order greater than this value of the Hubble parameter (m χ ≥ H beg ∼ = 5 × 10 −4 M pl ∼ = 10 15 GeV, which is a mass of the same order as those of the vector mesons responsible for transforming quarks into leptons in simple theories with SU(5) symmetry [58]), one can apply the WKB approximation to calculate the re-normalized energy density of the vacuum. After subtracting the adiabatic modes up to order four, we obtain an energy density of the order H 6 /m 2 χ [59], which is subdominant compared to the energy density of the background 3H 2 M 2 pl and, thus, does not affect the classical evolution of the inflation up to an abrupt phase transition where the adiabatic regime is broken, the χ-field stops being in the vacuum and particles are copiously produced with an energy density which decays slower than the one of the inflation, thus becoming eventually dominant.
The dynamical evolution of the vacuum modes could be understood as follows: the k-vacuum mode during the adiabatic regime can be approximated by χ (n) k,WKB , where n is the order of the WKB approximation, but, when the adiabatic regime breaks down during a period of time, the WKB approximation cannot be used and only at the end of this period one can again use it. However, now the vacuum mode is a combination of positive and negative frequency modes which can be approximated by a linear combination of χ (n) k,WKB and its conjugate of the form α k,n χ (n) k,WKB + β k,n (χ (n) k,WKB ) * , where α and β are the so-called Bogoliubov coefficients, and it is the manifestation of the gravitational particle production. Basically this is the viewpoint of particle creation in curved space-times [29][30][31], where the β-Bogoliubov coefficient, which is calculated matching the modes before and after the discontinuity for models with a discontinuity in some derivative of the potential as the one introduced by Peebles-Vilenkin in [4]. This is the key ingredient to calculate the energy density of the produced particles.
In fact, the energy density of the produced particles after the end of the phase transition evolves as [60] where ω k (τ) is the time dependent frequency of the k-mode, and when β k is known we have an analytic expression of this energy density that allows us to calculate the reheating temperature and deduce its bounds.

The Peebles-Vilenkin Model
To deal with an analytically solvable problem, i.e., having an analytic expression of the β-Bogoliubov coefficient, we consider a sudden phase transition where the third derivative of the Hubble parameter is discontinuous, which happens for the following improved version of the well-known Peebles-Vilenkin quintessential inflationary potential [4], where λ is a dimensionless parameter and M is a very small mass compared with the Planck one.
Here, it is important to point out that the inflationary part of the original Peebles-Vilenkin potential is a quartic potential and, thus, the theoretical values of the spectral index and the ratio of tensor to scalar perturbations do not enter in the marginalized joint confidence contour in the plane (n s , r) at 2σ CL [61] without the presence of the running [12]. This is the reason one has to change the quartic part by a Starobinsky-type potential, whose spectral values do actually enter in this contour.
The value of the parameter λ is calculated as follows: we use the theoretical and observational values of the power spectrum of the curvature fluctuation in a co-moving coordinate system when the pivot scale leaves the Hubble radius [62], P ζ ∼ = is the main slow-roll parameter and the star " * " means that the quantity is evaluated when the pivot scale leaves the Hubble radius, obtaining where we have used that for our model one has * ∼ = 3 16 (1 − n s ) 2 , where n s denotes the spectral index and during inflation H 2 * ∼ = λ 3 M 2 pl . From the recent observations by Planck [61] the value of the spectral index is constrained to be n s = 0.968 ± 0.006. Thus, taking its central value one gets λ ∼ = 9 × 10 −11 , which means that H * ∼ = 5.48 × 10 −6 M pl . The tensor-to-scalar perturbation ration r for this model yields r = 16 * ≈ 3(1 − n s ) 2 , which leads for this range of values of n s to a small enough quantity (r ≤ 0.00581 at 2σ C.L.) which agrees with the observational constraints.
On the other hand, note that for our toy model the second derivative of the potential is discontinuous at ϕ = 0, nearly at the beginning of the kination phase (To simplify, we will assume that kination starts when ϕ = 0 because, as is shown in Figure 1, the maximum value of the kinetic energy is very close to ϕ = 0). In addition, using Raychaudhuri equation, one can see that the third derivative of the Hubble rate is discontinuous at the beginning of kination, hence allowing particle production because the adiabatic evolution is broken. For example, if one considers a massive χ-field coupled to gravity, the fourth derivative of the frequency ω k (τ) = k 2 + a 2 (τ)m 2 χ is discontinuous for any k-mode.
In fact, this kind of potentials with discontinuities was studied by Starobinsky and others in [63,64], who showed that the discontinuity of the effective potential could be obtained introducing a second scalar field coupled to the inflaton that experiences a cosmological second order phase transition (see for instance the introduction of Linde's book [58] for some simple examples of first and second order phase transitions), as is explained in Section 4 of [64] considering the standard toy model used many times in the hybrid inflationary scenario [65].
What is important is that we have to understand the breakdown of the adiabatic behavior, at least for a more smooth potential, as follows: 1 in a region close to the beginning of kination with a characteristic time less than (H(t)) −1 and, thus, in this region the adiabatic regime is broken, allowing the production of particles. Unfortunately, in this situation the analytic calculation of the energy density of the produced particles is not possible. This is the reason we consider our toy model (2), where one can get an analytic expression of this energy density (Please note that the second derivative of (2) is discontinuous, meaning that the third derivative of the Hubble rate is discontinuous at the beginning of kination, i.e., the fourth derivative of ω k (τ) is discontinuous at that moment and, thus, the non-adiabatic condition (4) is met). Finally, numerical calculations (namely event-driven integration with an ode RK78 integrator) [53] show that at the beginning of kination one has H kin ∼ = 1.44 × 10 −6 M pl and, thus, the energy density of the background at the beginning of kination is given by ρ ϕ,kin ∼ = 6.26 × 10 −12 M 4 pl .

The Dynamics of the Model
To deal with the evolution of the system we need to consider the back-reaction of particle production. Effectively, after particle production one has the so-called semi-classical Friedmann equation H 2 = ρ ϕ +ρ χ 3M 2 pl , where ρ ϕ stands for the energy density of produced particles and ρ χ is the energy density of produced particles. Then, when the inflaton dominates, i.e., during inflation, kination and at late times, one has H 2 = ρ ϕ 3M 2 pl , but from the end of kination up to beginning of quintessence one has H 2 = ρ χ 3M 2 pl , i.e., the background is driven by the energy density of created particles, meaning that during radiation one has H = 1 3t and in the matter-domination era H = 2 3t , which influences the evolution of the infation field, which is given by the conservation equation Taking into account this fact, we start with the analytic analysis and after it we do the numerics.

Analytic Results
We start with the initial conditions at the beginning of kination for our improved version of the Peebles-Vilenkin model: During kination, the scale factor and the Hubble rate evolve as a ∝ t 1/3 =⇒ H = 1 3t and, from the Friedmann equation, the evolution in this phase will be Here two different situations can occur: The superheavy χ-particles created during the phase transition from the end of inflation to the beginning of kination could decay 1.
After the end of kination.

2.
Before the end of kination.
In the first case, at the end of kination one has where we used the relation H end = √ 2H kin Θ (see Section 4.2 for the deduction), being Θ ≡ ρ χ,kin ρ ϕ,kin (the ratio of the energy density of the χ-field to the one of the inflaton at the beginning of the kination phase) the so-called heating efficiency [66]. During the period between t end and t R (t R denotes the reheating time, i.e., when the universe starts to be radiation-dominated), in the case that the χ-particles were superheavy, the universe is matter-dominated and, thus, the Hubble parameter becomes H = 2 3t . During this epoch, the gradient of the potential could also be disregarded, hence the equation of the scalar field becomesφ + 2 tφ = 0 and, thus, where we used thatφ(t) = −φ end t t end 2 withφ end = 2 3 M pl t end . Then, one gets having employed that H 2 R = 2ρ ϕ,R 3M pl with ρ ϕ,R = π 2 30 g * T 4 R and we also have thaṫ During the radiation period one can continue disregarding the potential and the dynamical equation becomesφ + 3 2tφ = 0, whose solution is given by and, thus, sinceφ R t R = π 6 g * 30 H kin Θ (being T R the reheating temperature) at the matter-radiation equality, one has where g eq ∼ = 3.36 are the degrees of freedom at the matter-radiation equality and T eq is the temperature of the radiation at the matter-radiation equality, which is related to the energy density via the relation ρ eq = π 2 15 g eq T 4 eq ∼ = 8.8 × 10 −1 eV 4 and, thus, given by T eq ∼ = 7.9 × 10 −10 GeV T R . In the same way, Remark 1. To obtain the value of ρ eq , we chose as the value of the cosmic red-shift at the matter-radiation equality z eq ≡ −1 + a 0 a eq = 3365, the value of the ratio of the energy density of the matter to the critical energy density at the present time equal to Ω matt,0 = 0.308 and the value of the Hubble rate at the present time equal to H 0 = 1.42 × 10 −33 eV. Then, since ρ matt,0 = 3H 2 0 M 3 pl Ω matt,0 , one finally gets In the second case, i.e., when the decay of the χ-particles is before the end of kination, which always happens when reheating is via instant preheating, the beginning of the radiation era coincides with the end of kination. Thus, and, taking into account that During the radiation era, disregarding once again the potential, we will have but nowφ R t R = 2 3 M pl , meaning that given that T eq T R , andφ After the matter-radiation equality the dynamical equations cannot be solved analytically and, thus, one needs to use numerical methods to compute them. To do that we need to use a "time" variable that we choose to be minus the number of e-folds up to the present epoch, namely N ≡ − ln(1 + z) = ln a a 0 . Now, using the variable N, one can recast the energy density of radiation (the energy density of the decay products of the χ-field which we continue denoting by ρ χ ) and matter respectively as where N eq = − ln(1 + z eq ) ∼ = −8.121 is the value of N at the matter-radiation equality.
To obtain the dynamical system for our model, we introduce the following dimensionless variables, where once again H 0 ∼ = 1.42 × 10 −33 eV denotes the current value of the Hubble parameter. Now, using the variable N = − ln(1 + z) defined above and also the conservation equationφ + 3Hφ + V ϕ = 0, we will have the following non-autonomous dynamical system [67]: where the prime represents the derivative with respect to N, the Friedmann equation now looks as where we introduced the following dimensionless energy densitiesρ χ = Then, we have to integrate the dynamical system, starting at N eq = −8.121, with initial conditions x eq and y eq , and the value of the parameterM is obtained equaling at N = 0 Equation (23) to 1, i.e., imposingH(0) = 1. For the first case (the decay after the end of kination), the initial conditions are obtained analytically in Equations (12) and (13). Effectively, from Formula (35) in Section 4.2, where c.c. means that the χ-field is conformally coupled to gravity and n.c. non-conformally coupled. Then, for viable reheating temperatures T rh ≤ 10 9 GeV and as we will see in Section 4.2 for m χ ∼ = 10 15 GeV, one has y eq 1. In addition, for x eq , after a simple calculation, Last, the initial conditions for the second case (the decay before the end of kination) are , y eq ∼ = 4π 9 g eq 5

Numerical Results
In all cases compatible with the constraints found in this manuscript, which are summarized in Table 1 (see Conclusions), we obtained that M/M pl ∼ = 10 −13 , which coincides with the result obtained in [4]. Please note that for this model, the energy scale of inflation V 1/4 (ϕ −M pl ) ∼ λ 1/4 M pl ∼ 10 15 GeV is close to the GUT scale, while the energy scale for dark energy V 1/4 (ϕ ∼ = 0) ∼ λ 1/4 M ∼ 10 2 GeV is near the electroweak scale.
Next we show in Figures 2 and 3 the reduced densities {ρ i } i=χ,m,ϕ , the density parameters {Ω i } i=χ,m,ϕ and the effective Equation of State (EoS) parameter ω e f f for all of them, showing that at the present time w e f f ∼ = −0.6 < −1/3, which proves the current cosmic acceleration, and at late time w e f f goes to −1, meaning that this model leads to an eternal acceleration. As clearly seen in the figures, the results remain almost unchanged for the different considered cases, corresponding to different values of the reheating temperature, given that for all of them the value of parameterM yields almost the same value, namelyM ∼ 10 13 GeV ∼ 10 −5 M pl .

Compatibility of the Model with the Cosmological Perturbations
After having studied the dynamics of the model, in order to verify its compatibility with the cosmological perturbations, we are going to compare the number of e-folds for our considered potentials, namely N = 2 1−n s , with the one obtained from [68] where M symbolizes the beginning of the matter domination era. Analogously as in [17], it leads to where g R = 107, 90 and 11 respectively for T R ≥ 175 GeV, 175 GeV ≥ T R ≥ 200 MeV and 200 MeV ≥ T R ≥ 1 MeV; ln a end a E = H end H kin H(t)dt, which was numerically calculated for both considered potentials when we take the spectrum index to be the central value n s = 0.968, and Therefore, we obtain the value of the reheating temperature in function of n s for both potentials, which is represented in Figure 4. We observe that all the important bounds for our model, namely the BBN ones and the ones summarized in Table 1 (see conclusions) lay within the allowed values for the spectral index, namely n s = 0.968 ± 0.006 [61].

Reheating in Quintessential Inflation
In this section, we will discuss the most common ways to reheat the universe: Reheating via gravitational production of light or superheavy particles and instant preheating.

Gravitational Production of Light Particles
When the produced particles during the phase transition are very light, the energy density of the relativistic plasma formed by these light particles is given by [4,36,56,60,69,70] where R ∼ 10 −2 N s , being N s the number of scalar fields, which for the minimal GUT is 4 (the electro-weak Higgs doublet) [4]. Therefore, we will use that R ∼ = 10 −1 .

Remark 2.
Here it is important to recall that this formula is only obtained for toy models (see for instance [36,70]) and we understand that it will also work for more realistic models in which the adiabatic evolution is broken near the beginning of the kination phase.
Since as we immediately show the thermalization process of the plasma is an instantaneous process, the universe will become reheated at the end of the kination epoch, i.e., when the energy densities of the scalar field and that of the relativistic plasma were of the same order. This occurs when a kin a R 2 = Θ, where we used, once again, the so-called heating efficiency defined in Section 3.1.1 as Θ ≡ ρ χ,kin ρ ϕ,kin and the fact that the energy density of the produced particles decays as a −4 while the one of the inflaton field decays as a −6 during kination. Thus, the reheating temperature is given by where g * = 106.75 are the degrees of freedom for the Standard Model. Then, and, since a simple calculation leads to the value Θ ∼ = H 2 kin 30M 2 pl ∼ = 6.9 × 10 −14 , we can conclude that the reheating temperature when the reheating is via the gravitational production of light particles is T R ∼ = 213 TeV, which is basically the same result as the one obtained by Peebles and Vilenkin in their paper [4]. Finally, we will show that the thermalization is nearly an instantaneous process compared with the duration of the kinetic era. Following the reasoning of [3] and [4], the decay products have a typical energy of the form¯ ∼ H kin a kin a(τ) and their number density is n ∼ R¯ 3 ∼ = 10 −1¯ 3 . Now, we take into account that if the particles interact by the exchange of gauge bosons and establish thermal equilibrium among the fermions and gauge bosons, the interaction rate will be nσ, where the cross section is given by σ ∼ α 2 2 , with the coupling constant satisfying the inequality 10 −2 ≤ α ≤ 10 −1 . Therefore, the thermal equilibrium will be accomplished when the interaction rate becomes comparable to the , which happens when a kin a th 2 = 10 −1 α 2 , where the subscript "th" attached to any quantity refers to its value at the time when the thermal equilibrium was established. On the other hand, one can calculate the scale factor at the reheating time t = t R , which occurs at the end of kination, as follows: Since a kin a R 2 = Θ, then we will have and, thus, where we used that during kination the scale factor evolves as t 1/3 . This result means that the thermal equilibrium occurs well before the equality between the energy density of the scalar field and the one of the decay products, i.e., well before to the end of kination which in this case coincides with the beginning of the radiation era. Hence, one can safely assume an instantaneous thermalization.

Gravitational Production of Superheavy Particles
In this subsection, we will assume that the χ-field has a mass m χ greater than 10 15 GeV. Therefore, since this mass is greater than the Hubble rate, we can apply the WKB solution when we calculate the evolution of the modes. In the appendix of [15] it was shown that the leading term of the re-normalized energy density of the produced particles after the phase transition is given by and, for our model, in order to obtain the β-Bogoliubov coefficient we use the WKB approximation.
Please note that the second iteration W k including temporal derivatives up to order four was obtained in [71], which is enough for our calculations because for the model (2) the third derivative of the Hubble rate is discontinuous and the term responsible for the leading contribution to the k . As we show in Appendix A, for the conformally coupled case, i.e., when ξ = 1/6, the term leading to the main contribution is given by for the conformally coupled case for the nonconformally coupled case, (33) where in the nonconformally coupled case we took ξ − 1 6 ∼ = 1, which is its maximum value because the WKB approximation is a perturbative one that only holds when m 2 χ ξ − 1 6 R and, thus, since at GUT scales R ∼ = 12H 2 ∼ = 10 29 GeV 2 and M pl m χ ≥ 10 15 GeV, we can conclude that ξ − 1 6 ≤ 1.

Remark 3.
Please note that in the nonconformally coupled case in (A11) we corrected by a factor of 2 the result obtained in [53].
As we will see in next section dealing with the overproduction of GW, when considering reheating via gravitational production of superheavy particles, in order to prevent the BBN success we have to impose the decay of the χ-field to be after the end of kination, i.e., after the equality between the energy density of the field and the one of the produced particles.
Since the decay is after t end (t end denotes the instant when kination ends), one has to impose Γ ≤ H end , where Γ is the decay rate of χ-particles. Taking this into account, one has where we used that for a superheavy field the heating efficiency satisfies Θ = a kin a end

3
. On the other hand, a simple calculation leads to the result Consequently, from Equation (34) one can easily find H end = √ 2H kin Θ and, thus, one obtains that the decay rate has to satisfy Γ ≤ √ 2H kin Θ, which means that GeV for the conformally coupled case GeV for the nonconformally coupled case.
Since as we already showed that the thermalization is nearly instantaneous, the reheating temperature (i.e., the temperature of the universe when the thermalized plasma starts to dominate) will be T R = 30 where we used that after t end the energy density of the produced particles dominates the energy density of the inflaton field. Then, we will have that T R ∼ = 0.54 In addition, finally, since we are assuming that the mass of the χ-field is greater than 10 15 GeV, we obtain the following upper bound for the reheating temperature, for the conformally coupled case 178 TeV for the nonconformally coupled case.
We end this subsection noting that these bounds are obtained without taking into account the production of GWs, which leads to more restrictive bounds as we will see in next section.

Instant Preheating
In this subsection, we will assume that the bare mass of the χ-field, which we impose to be conformally coupled to gravity, is zero and we also consider an interaction between the inflaton field ϕ and the quantum χ-field, whose interacting Lagrangian is given by L int = − 1 2 g 2 ϕ 2 χ 2 , where g is a dimensionless coupling constant. The enhanced symmetry point was chosen ϕ = 0 because at this point the velocity of the scalar field is nearly maximum as one can see in Figure 1. In this situation the χ-particles, which have an effective mass m χ,e f f (t) = g|ϕ(τ)|, are created via a mechanism named instant preheating, which was introduced in [42] in the framework of standard inflation and was applied for the first time to quintessential inflation in [43].

Remark 4.
The reheating via instant preheating is usually used in models with very smooth potentials because in these models the gravitational production of particles is completely inefficient due to the adiabatic regime during all the evolution (see for instance [20,44]). On the contrary, the introduction of the interacting Lagrangian term depicted above breaks down the adiabatic evolution at the beginning of the kination phase, which as we will see allows the production of enough particles to reheat the universe in a viable way. However, as we will see in Section 6.2, if we assume that dark matter is also created during the phase transition from inflation to kination, then this dark matter cannot be created via instant preheating and one needs another mechanism to create it, which could be the gravitational particle production. Therefore, in this hypothetical situation the potential cannot be so smooth.
Then, if the reheating is via instant preheating, soon after the beginning of kination the χ-field acquires an effective mass equal to m χ,e f f = gM pl and the energy density of the χ-field is given by [43] ρ χ (τ) = gM pl n χ (τ) = gM pl n χ,kin a kin a(τ) where the number density of particles at the beginning of kination is calculated as follows: Near the beginning of kination, i.e., when ϕ = 0, one has ϕ(τ) ∼ = ϕ kin (τ − τ kin ) and the frequency of the k-mode of the field χ is ω k (τ) = k 2 + g 2 a 2 kin ϕ 2 kin (τ − τ kin ) 2 , where the expansion of the universe is not considered and for this reason we approximated the scale factor by its value at the beginning of kination. Then, the k-mode of the χ-field satisfies the equation of a time dependent harmonic oscillator obtaining an over-barrier problem in scattering theory, whose β-Bogoliubov coefficient is related to the reflexion coefficient via the formula [72][73][74][75] where γ denotes a closed path that wraps around the turning points τ ± = τ kin ± i k ga kin ϕ kin and the average number of produced particles in the k-mode is given by Thus, the average number density of χ-particles at the beginning of kination is given by Remark 5. Since the effective mass of the χ-field is g|ϕ(τ)|, in order to prevent the vacuum polarization effects from affecting the evolution of the inflaton field during inflation, one has to impose the effective mass of the χ-field to be greater than the Hubble rate, which leads to the condition which always holds if we assume that g ≥ λ 3 M pl |ϕ END | (because |ϕ(τ)| is a decreasing function during inflation), where ϕ END denotes the value of the ϕ-field at the end of inflation. Since inflation ends when the slow-roll parameter is equal to one, one easily gets that ϕ END = 3 2 ln( 94M pl , which leads to the constraint for the parameter g of g ≥ 5.83 × 10 −6 .
These particles are very massive and, in order to avoid a second inflationary epoch due to the χ-field, one has to assume that the decay is well before the end of the kination regime [43]. Since the thermalization is nearly instantaneous as we have already seen, in this case the reheating is completed at the end of kination and, thus, the reheating temperature is calculated as follows: Using that H dec H kin = Γ H kin = a kin a dec 3 we have ρ ϕ,dec = 3Γ 2 M 2 pl , and ρ χ,dec = On the other hand, from the condition ρ χ,dec ≤ ρ ϕ,dec (the decay is before the end of the kination phase), one gets Γ ≥ 6.18 × 10 −6 g 5/2 M pl (48) and, since we showed that g ≥ 5.83 × 10 −6 , we see that the decay rate is greater than 5.07 × 10 −19 M pl ∼ = 1.2 GeV. We note that the evolution of the energy density of the created particles and the background are respectively ρ χ (t) = ρ χ,dec a dec a(t) 4 , ρ ϕ (t) = ρ ϕ,dec a dec a(t) GeV.
This reheating temperature [i.e., Equation (50)] could be bounded using (48), obtaining On the other hand, since the decay is after the beginning of kination, we have that Γ ≤ H kin ∼ = 1.44 × 10 −6 M pl ∼ = 3.51 × 10 12 GeV, getting the bound which means that in order to preserve the BBN success, a reheating temperature approximately between 1 MeV and 10 6 TeV is required and, thus, the constraint g ≤ 2.76 × 10 −4 has to be satisfied, which restricts the value of g in the following narrow band, Finally, if for example we choose g ∼ = 10 −5 and Γ ∼ = 10 −10 M pl , which satisfy the constraint (48), one obtains a reheating temperature equal to T R ∼ = 2.18 × 10 4 TeV.

BBN Constraints Coming from the Production of Gravitational Waves
This section is devoted to present the bounds of the proposed improved version of the quintessential inflationary model using the Big Bang Nucleosynthesis (BBN), where we explicitly use the BBN constraints from the logarithmic spectrum of GWs and consequently the BBN bounds from the overproduction of GWs.

BBN Constraints from the Logarithmic Spectrum of Gws
It is well-known that during inflation GWs are produced (known as primordial GWs, in short PGWs) and in the post-inflationary period, i.e., during kination, the logarithmic spectrum of GWs, namely Ω GW defined as Ω GW ≡ 1 ρ c dρ GW (k) d ln k (where ρ GW (k) is the energy density spectrum of the produced GWs; ρ c = 3H 2 0 M 2 pl , where H 0 is the present value of the Hubble parameter, is the so-called critical density) scales as k 2 [66], producing a spike in the spectrum of GWs at high frequencies. Then, so that GWs do not destabilize the BBN, the following bound must be imposed (see Section 7.1 of [52]), where h 0 ∼ = 0.678 parametrizes the experimental uncertainty to determine the current value of the Hubble constant and k BBN , k end are the momenta associated with the horizon scale at the BBN and at the end of inflation respectively. As was shown in [69], the main contribution of the integral (55) comes from the modes that leave the Hubble radius before the end of the inflationary epoch and finally re-enter during the kination, that means, for k end ≤ k ≤ k kin , where k end = a end H end and k kin = a kin H kin .
For these modes one can calculate the logarithmic spectrum of GWs as in [76] (see also [66,[77][78][79] where the graviton spectra in quintessential models were reassessed, in a model-independent way, using numerical techniques), where is the amplitude of the GWs; Ω γ ∼ = 2.6 × 10 −5 h −2 0 is the present density fraction of radiation, and the quantity˜ , which is approximately equal to 0.05 for the Standard Model of particle physics, takes into account the variation of massless degrees of freedom between decoupling and thermalization (see [66,69] for more details). As was derived in [69], the specific form of the expression above comes from the behavior of the Hankel functions for small arguments. Now, plugging expression (56) into (55) and disregarding the sub-leading logarithmic terms, one finds Remark 6. A further bound on primordial gravitational waves is imposed by the CMB constraint on additional massless degrees of freedom. As GWs with frequencies larger than the corresponding horizon at CMB decoupling contribute to the radiation density of the Universe, constraints on their total energy density can be phrased in terms of the effective number of massless neutrino species N e f f , which is bounded by N e f f = 3.04 ± 0.17 (see Section 5 of [80]), namely and, thus, turning to the following constraint, which is practically the same constraint obtained above. Please note that this constraint is more restrictive than the one obtained with the effective number of massless neutrino species from BBN, namely N e f f = 3.28 ± 0.28 [80], which leads to the same constraint as in Equation (57).
To calculate the ratio k kin /k end , we will have to study the following three different situations: 1.
When the produced particles are very light and its energy density decays as a −4 . In this case, as was shown in [66], one will have where Θ is once again the heating efficiency introduced previously. Thus, the constraint (57) eventually leads to Taking into account that when the reheating is due to the creation of very light particles during the phase transition the reheating temperature is (see Section 4.1) one has the following lower bound, Finally, note that we showed that when reheating is due to the gravitational production of light particles the reheating temperature is T R ∼ = 213 TeV, which means that the reheating via the gravitational production of light particles satisfies the bound (55).

2.
When the reheating is due to the production of superheavy particles which decay after the end of kination, as we showed in Section 4.2, we have that Θ = a kin a end 3 and H end = √ 2H kin Θ and, thus, Therefore, the constraint (57) leads to Now, since in Section 4.2 we obtained the following value of the heating efficiency, we deduce that in the conformally coupled case the mass of the χ-field has to satisfy m χ ≤ 3.52 × 10 14 GeV, which is incompatible with our assumption m χ ≥ 10 15 GeV. This shows that the gravitational production of superheavy particles conformally coupled to gravity is not viable. On the contrary, when the χ-field is not conformally coupled to gravity, one gets the bound m χ ≤ 1.5 × 10 15 GeV, which means that the viability of our model requires its mass to be m χ ∼ = 10 15 GeV when reheating is due to the gravitational production of superheavy particles noncoformally coupled to gravity.

3.
When the decay happens before the end of kination, as the case of instant reheating, a simple calculation leads to and, consequently, the constraint (57) becomes which applied to our model finally leads to another lower bound of the reheating temperature, which is obtained via instant preheating (see Section 4.3), T R ≥ 1.97 × 10 5 g 15/8 Θ 3/4 GeV ∼ = 1.84 × 10 6 g 3/8 GeV, (69) where we used that in the case of instant preheating one has Θ = g 2 2π 2 . Now, since the reheating temperature has to be less than 10 6 TeV, one gets the bound g ≤ 1.97 × 10 7 , which is less restrictive than the one obtained in Section 4.3, meaning that when reheating is due to the production of particles via instant preheating, the bound (55) is clearly overpassed.

BBN Bounds from the Overproduction of Gws
The success of the BBN demands that [44] where ρ GW (t) is the energy density of the GWs produced at the phase transition and both quantities are evaluated at the reheating time. The value of the energy density of the GWs is ρ GW (t) ∼ = (see for example [4,70]). Then, when the reheating is via the gravitational production of light particles, we have which as pointed out by Peebles and Vilenkin, results 1 N s = 0.25 for a minimal GUT. Therefore, this bound is never reached, meaning that in this case the overproduction of GWs could affect the BBN process. However, if one goes beyond a minimal GUT and accepts minimal supersymmetric (SUSY) theories, then N s = 104 and, thus, the bound (70) is overpassed.
On the other hand, in the case in which superheavy particles (which could decay in lighter ones to match with the HBB) are gravitationally created during the phase transition, we firstly see that the decay can never be before the end of kination because, if so, at the decay time, which occurs when H dec = Γ, we would have  (74) and, thus, Therefore, imposing the constraint (70), one gets that m χ ≤ 3.5 × 10 14 GeV, which contradicts our assumption that m χ ≥ 10 15 GeV.
Hence, the decay must occur after t end and, assuming once again the instantaneous thermalization, the reheating time will coincide with the decay one. Then, we will have ρ χ,dec = 3Γ 2 M 2 pl and, since we obtain and, thus, Now, we use that for the nonconformally coupled case we already showed that m χ ∼ = 10 which reduces the maximum reheating temperature T R ∼ = 0.54 Finally, in the case of instant preheating, when the decay is before the end of kination as we already explained, using Formula (72) and the fact that Θ = g 2 2π 2 we arrive at and, taking into account the bound Γ M pl ≤ 1.44 × 10 −6 , we get Therefore, since at the end of Section 4.3 we already showed that g ≥ 5.83 × 10 −6 , we reach which assures that the constraint (70) is fulfilled when the reheating is via instant preheating. Summing up, we showed that a viable reheating in the case of gravitational reheating requires the creation of superheavy nonconformally coupled particles with mass nearly 10 15 GeV, which must decay after the end of kination, obtaining a maximum reheating temperature around 37 TeV. In addition, when the reheating is via instant preheating, the coupling constant g has to be close to 10 −5 , obtaining a reheating temperature greater than 20 TeV (where we used the bound (69)).

Abundance of Dark Matter
In this section, we will explore the possibility that the breakdown of the adiabatic regime leads to the possibility to explain the abundance of dark matter through the gravitational production of superheavy particles [40,81], although gravitational production of dark matter could also occur in standard inflation during the oscillations of the inflaton field [82][83][84][85] (see also the earlier papers [9,86,87]).

Reheating Via Production of Superheavy Nonconformally Coupled to Gravity
As we already showed, the reheating via gravitational production of superheavy particles conformally coupled is not viable and, when these superheavy particles are nonconformally coupled to gravity, its mass must be very close to 10 15 GeV. Now we also assume that there is another kind of superheavy particles conformally coupled to gravity, named Y-particles, which do not decay and only interact gravitationally, and could be the responsible for the current abundance of the dark matter. In Section 4.2, we saw that the energy density of the Y-particles will be ρ Y (τ) ∼ = . Therefore, since in order to preserve the BBN success the decay of the χ-particles has to be after the end of kination, as we saw in Section 3.1, the thermalization is instantaneous. Therefore, at the reheating time we will have Now, taking into account that m χ ∼ = 10 15 GeV and the energy density of the χ-particles must be greater than the one of the Y-particles so that the universe reheats, we deduce the following bound for the mass of Y-particles, After reheating the evolution of the corresponding energy densities will be meaning that at the matter-radiation equality On the other hand, taking the following observational data at present time (H 0 ∼ = 1.42 × 10 −33 eV, Ω Y,0 = 0.262, Ω b,0 = 0.048 and Ω matt,0 = 0.31, where b denotes the baryonic matter and matt the total matter (dark+baryonic)), it is satisfied that where z eq denotes once again the cosmic red-shift at the matter-radiation equality. Then, since ρ χ,eq = ρ matt,eq at the matter-radiation equality, we will have ρ χ,eq = ρ matt,0 (1 + z eq ) 3 and, thus, which combined with (87), leads to the relation and, consequently, where we used that m χ ∼ = 10 15 GeV. Now, taking into account that ρ Y,eq = 3H 2 0 M 2 pl Ω Y,0 (1 + z eq ) 3 and choosing 3365 as the value of the cosmic red-shift at matter-radiation equality, we obtain ρ Y,eq ∼ = 3.598 × 10 −1 eV 4 and, inserting this expression in (90), one gets the following relation between the reheating temperature and the mass of dark matter, GeV.
Then, since the decay is after the end of kination and the thermalization is nearly instantaneous, the reheating temperature is which gives us the following relation between the mass of dark matter and the decay rate of the χ-particles, Finally, we have to use the bounds that must satisfy the decay rate to bound the mass of dark matter. For example the overproduction of gravitational waves leads to Γ M pl ≤ 7.90 × 10 −28 and, thus, m Y ≤ 7.94 × 10 17 GeV. In addition, taking into account that the reheating temperature must be greater than 1 MeV because the BBN occurs approximately at 1 MeV when the universe is already reheated, from Equation (92) one gets that is, a viable quintessential inflation model where dark matter is created gravitationally requires the mass of dark matter to be bounded as follows, and a maximum reheating temperature around 37 TeV.

Reheating via Instant Preheating
First of all it is important to note that when reheating is via instant preheating, the production of dark matter cannot be via the same mechanism because the coupling constant g is restricted to be close to 10 −5 and, thus, the energy of the dark matter and the one of the particles that reheat the universe after its decay would be of the same order, which would forbid a radiation phase, which is essential for correctly depicting the evolution of our universe.
Therefore, in that case we also have to consider the possibility that dark matter was created gravitationally. Therefore, we consider once again superheavy Y-particles only conformally coupled to gravity, which would be the responsible for the present abundance of dark matter in the universe in our model.
As we already discussed, in order to avoid a second inflationary period it is mandatory that unlike the superheavy particles created gravitationally studied in the previous section, these χ-particles decay well before the end of kination. Then, at the matter-radiation equality we will have ρ χ,eq = ρ χ,dec a dec a eq 4 , ρ Y,eq = ρ Y,dec a dec a eq 3 , and, thus, a dec a eq = gM pl n χ,kin where we used that at the decay time ρ χ,dec = gM pl n χ,kin On the other hand, as we already saw in the previous subsection, meaning that a dec a eq = Ω matt,0 Ω Y,0 ρ Y,kin gM pl n χ,kin .
Therefore, we will have ρ Y,eq = ρ χ,eq Ω Y,0 Ω matt,0 = gM pl n χ,kin Γ H kin we get that but the energy density of the dark matter at the matter-radiation equality is ρ Y,eq ∼ = 3.598 × 10 −1 eV 4 . Then, we have the following relation between the mass of the dark matter and the decay rate of the χ-particles where we have used that g ∼ = 10 −5 . Finally, using that when the reheating is via instant preheating the decay of the χ-particles will be during the kination phase, we have the bounds obtained in Section 4.3, which bound the mass of the dark matter to be in the domain that is, the mass of dark matter must be very close to 10 17 GeV when reheating is via instant preheating.

Other Kind of Potentials
In this section, we would like to check the importance of the breakdown of the adiabatic evolution. For this reason we will consider a more abrupt phase transition than the one given by the potential (2). For example, we choose the following potential, where now α denotes a positive dimensionless parameter. Please note that in this case, the inflationary piece is an Exponential SUSY inflation-type potential.
Here, when the field vanishes the potential has a discontinuity in its first derivative, which was pointed out in [63,64], and could be obtained introducing a second scalar field experiencing a first order phase transition. This means, using Raychaudhuri equation, that the second derivative of the field and, thus, the second derivative of the Hubble rate are discontinuous at the beginning of kination. Therefore, once again, we have to understand this model as a toy model, which allows us to perform analytically all the calculations, belonging to the class of potentials satisfying 1 near the beginning of kination, where ω k (τ) = k 2 + m 2 χ a 2 (τ) is the frequency of the k-mode of the χ-field. Another important thing is that since the phase transition is abrupter than in the previous case, this means that now the gravitational production will be greater, thus obtaining a greater reheating temperature than for the potential (2).
As proved in [67], both potentials lead to equivalent expressions for the spectrum index, number of e-folds and ratio of tensor to scalar perturbations, that is, the spectral index is n s ∼ = 1 − 2 N and the tensor/scalar ratio is given by r = 8 α 2 N 2 , where the number of e-folds depends on the reheating temperature as follows, Hence, the computations will be done for the same value of α as used in the previous potential, namely α = 2 3 , for which we obtain the corresponding numerical value ofφ kin = 5.75 × 10 −6 M 2 pl . If we first consider the case in which the reheating takes place via the gravitational production of light particles, the energy density of produced particles is, once again, given by ρ χ (τ) ∼ = a kin a(τ) 4 and, following step by step the reasoning done in Section 4.1, we get that the reheating temperature becomes greater than before, namely T R ∼ = 568 TeV.
With regards to the case of gravitational production of superheavy particles, here we first need to recalculate the β-Bogoliubov coefficients both when the χ field is conformally coupled to gravity or not.
In this case we only need the first order WKB approximation, i.e., W (1) k , which is obtained in Appendix A, and the value of the β-Bogoliubov coefficient (see Formula (A14) of Appendix B), to obtain the energy density of the produced particles at the beginning of kination which leads for the decay being after the end of kination to a reheating temperature of which is respectively 3.77 × 10 3 TeV and 4.56 × 10 4 TeV when restricting m χ ≥ 10 15 GeV, so we obtain minimum reheating temperatures considerably greater than the ones obtained for the potential with a discontinuity in the second derivative. Now, when considering instant preheating, the results are very similar to the other potential, obtaining as well a very narrow band for g corresponding to g ∼ = 10 −5 . Then, by taking into account the BBN constraints, we first consider the ones coming from the logarithmic spectrum of GW. When the reheating is produced via the production of light particles, we get the constraint T R ≥ 260 TeV, which is fulfilled by the computed value T R ∼ = 568 TeV. With regards to the production of superheavy particles decaying after the end of kination we obtain the bounds m χ ≤ 1.80 × 10 16 GeV for the conformally coupled case and m χ ≤ 2.17 × 10 17 GeV when non-conformally coupled. Please note that differently from the other potential, both constraints are compatible with m χ ≥ 10 15 GeV.
As we will immediately see, when dealing with the overproduction of GWs, for this potential the decay of superheavy particles is possible before the end of kination. In this case, the reheating temperature is given by (50), namely T R = 30 with which leads to the following reheating temperature.
As for the case of instant preheating. we get a minimum reheating temperature around 55 TeV. Thus, for this kind of potentials, if the particles decay before the end of kination the reheating via gravitational production of superheavy particles could lead to a reheating temperature greater than the one obtained when the reheating is via instant preheating.
Regarding the constraints from the overproduction of GW, reheating via production of light particles is also forbidden and, in the case of the production of superheavy particles decaying before the end of kination, we obtain that which restricts m χ to be m χ ≤ 9.94 × 10 14 GeV and m χ ≤ 1.2 × 10 16 GeV for the conformal and non-conformal cases respectively. Hence, both cases are compatible with the minimum mass for m χ . In addition, when superheavy particles decay after the end of kination, the upper bound for the reheating temperature is T R ≤ 6.26 × 10 3 TeV for ξ = 1/6 9.20 × 10 5 TeV for ξ − 1 which does not reduce the bounds in (111) for m χ ≥ 10 15 GeV. Therefore, the constraints coming from the production of GWs lead to completely different results depending on how abrupt the phase transition is. For instance, for the first potential reheating via gravitational production of superheavy conformally coupled particles is forbidden, which does not happen with the second potential if these particles decay after the end of kination. In addition, for the second potential the decay of superheavy particles, both conformally and non-conformally coupled to gravity, could be produced before or after the end of the reheating, obtaining a very efficient reheating mechanism which could lead to reheating temperatures greater than the one obtained using instant preheating as a reheating mechanism.
Dealing with the present abundance of dark matter, we are going to consider two types of particles: χ-particles, which can be now both conformally and non-conformally coupled to gravity, and Y-particles, conformally coupled and responsible for the abundance of dark matter. In contrast with the former potential, the decay of χ-particles can be both before and after the end of kination. If we first proceed analogously as in Section 6.1, i.e., by considering that the decay is produced after the end of kination when χ-particles are non-conformally coupled, we obtain that reaching finally the following bounds, where we used that a dec a eq = , given that ρ a,dec = ρ a,kin Γ H kin for a = χ, Y. Now, using Equation (68) and Equation (118), we find that bounds for the dark matter mass are considerably over the Planck's mass. The same happens when χ-particles are conformally coupled. Hence, the presence of dark matter gravitationally created during the phase transition from the end of inflation to the beginning of kination cannot be explained with a reheating via the gravitational creation of superheavy particles decaying before the end of kination for this potential.
Finally, when the reheating is produced via instant preheating, the energy density of the dark matter particles at the matter-radiation equality results ρ Y,eq ∼ = 3.14 × leading, as well, to bounds for m Y over the Planck's mass. Therefore, instant preheating cannot either be used for this potential in order to explain the presence of gravitationally produced dark matter.

Concluding Remarks
In the present work, we studied with all details the reheating of the universe via gravitational particle production and via instant preheating in quintessential inflation, taking into account the bound imposed by the production of GWs during the phase transition between the end of inflation and the beginning of kination.
To perform analytically all the calculations, we considered a toy model inspired in the well-known Peebles-Vilenkin model in which the discontinuity occurs in the second derivative of the potential.
Our study shows that the reheating via gravitational production of light particles is forbidden due to the overproduction of GWs, i.e., the bounds imposed to prevent the success of the BBN are not overpassed. A similar situation occurs when the reheating is via the gravitational production of superheavy particles conformally coupled to gravity, in this case the bound imposed by the spectrum of the GWs is not accomplished. Therefore, only two situations are acceptable to have a viable reheating that does not affect the success of the BBN: 1.
Reheating via graviational particle production of superheavy particles not conformally coupled to gravity.
However, several restrictions must be imposed to the parameters appearing in the theory: In the case of gravitational production of superheavy particles nonconformally coupled to gravity, the decay of theses particles in lighter ones in order to obtain a relativistic plasma has to be after the end of kination obtaining a maximum reheating temperature around 37 TeV. In addition, the mass of these superheavy particles has to be approximately equal to 10 15 GeV.
On the contrary, when reheating is via instant preheating, the produced particles have to decay before the end of the kination phase, obtaining a minimum temperature around 10 TeV. Moreover, the dimensionless coupling constant between the inflaton field and these particles has to be of the order of 10 −5 .
On the other hand, when one assumes that dark matter could be created via gravitational particle production of conformally coupled particles during the phase transition from the end of inflation to the beginning of kination, its mass has to range between 10 16 GeV and 10 18 GeV, when the reheating is via gravitational production of superheavy particles nonconformally coupled to gravity. In addition, when the reheating is via instant preheating the mass of the dark matter would only be of the order of 10 17 GeV.
In last section, we have considered another toy model inspired in the Peebles-Vilenkin model in which the discontinuity occurs in the first derivative of the potential and we have shown the differences with respect to the first potential considered, i.e., with the one with the discontinuity in the second derivative of the potential. Basically, in that case, if the reheating is via gravitational production of superheavy particles the reheating temperature is considerably increased being able to be even greater than the one obtained when the reheating mechanism is via the instant preheating.
Moreover, the constraints coming from the production of GWs allow in this case the decay of superheavy particles before and after the end of kination. However, if one assumes that the abundance of dark matter is due to its gravitational production during the phase transition, then neither the reheating via gravitational production of superheavy particles decaying before the end of kination nor via instant preheating could be able to explain this abundance.
Finally, we show the allowed cases with the corresponding values of the parameters in the following table, where c.c. and n.c. stand for χ-particles being respectively conformally and non-conformally coupled to gravity, V 1 (ϕ) and V 2 (ϕ) refer to the two potentials that were considered, and a line was drawn where we have achieved no constraints because the corresponding process was proved as forbidden.  Acknowledgments: This investigation has been supported by MINECO (Spain) grant MTM2017-84214-C2-1-P, and also in part by the Catalan Government 2017-SGR-247.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The WKB Approximation in Cosmology
For a non-conformally coupled with gravity χ-field, its k-mode satisfies the equation [60] where Ω 2 k = ω 2 k + ξ − 1 6 a 2 R, being ξ the coupling constant to gravity, R = 6(Ḣ + 2H 2 ) the Ricci scalar and ω k = k 2 + a 2 m 2 χ the time dependent frequency of the k-mode.
and the second iteration W (2) k including temporal derivatives up to order four was obtained in [71]. The result of this calculation in terms of the cosmic time is given by As a consequence, for the conformally coupled case, i.e., when ξ = 1/6, the term leading to the main contribution is given by
On the other hand, from Raychaudhuri equation ...