# Gravitational particle production of superheavy massive particles in
Quintessential Inflation: A numerical analysis

Llibert Aresté Saló<EMAIL_ADDRESS>School of Mathematical Sciences,
Queen Mary University of London, Mile End Road, London, E1 4NS, United Kingdom
Jaume de Haro<EMAIL_ADDRESS>Departament de Matemàtiques, Universitat
Politècnica de Catalunya, Diagonal 647, 08028 Barcelona, Spain

###### Abstract

We compute numerically the reheating temperature due to the gravitational
production of conformally coupled superheavy particles during the phase
transition from the end of inflation to the beginning of kination in two
different Quintessential Inflation (QI) scenarios, namely Lorentzian
Quintessential Inflation (LQI) and $\alpha$-attractors in the context of
Quintessential Inflation ($\alpha$-QI). Once these superheavy particles have
been created, they must decay into lighter ones to form a relativistic plasma,
whose energy density will eventually dominate the one of the inflaton field in
order to reheat after inflation our universe with a very high temperature, in
both cases greater than $10^{7}$ GeV, contrary to the usual belief that heavy
masses suppress the particle production and, thus, lead to an inefficient
reheating temperature. Finally, we will show that the over-production of
Gravitational Waves (GWs) during this phase transition, when one deals with
our models, does not disturb the Big Bang Nucleosynthesis (BBN) success.

Gravitational Particle production; Quintessential Inflation; Reheating;
Gravitational Waves; Numerical Calculations.

###### pacs:

04.20.-q, 98.80.Jk, 98.80.Bp

##  I Introduction

Today, the inflationary paradigm is the most accepted implementation to the
Big Bang (BB) theory in order to solve a number of shortcomings associated
with the standard BB cosmology, such as the horizon, the flatness or the
primordial monopole problems guth ; linde , and to explain correctly the early
universe at background level. More remarkable is the fact that inflation is
also able to explain the origin of inhomogeneities in the universe as quantum
fluctuations chibisov ; starobinsky ; pi ; bardeen ; Linde:1982uu , leading to
theoretical results that match greatly with the recent observational data
provided by the Planck’s team Planck .

Once having a viable theory explaining the early universe, one can extend it
in order to deal with the whole evolution of the universe, and thus, unifying
its early- and late-time accelerated expansions (see odintsov ; li ;
Copeland:2006wr  for a review of the current dark energy models). One of the
most attractive scenarios able to do it is the so-called Quintessential
Inflation (QI), introduced for the first time by Peebles and Vilenkin in his
seminal paper pv  (see hap19  for a review of the Peebles-Vilenkin model),
where the idea behind their proposal comes through the introduction of a
single scalar field, also named inflaton, that at early times is the
responsible for inflation while at late times it allows the current cosmic
acceleration via quintessence.

Due to the simplicity of this proposal and since the behavior of the slow-roll
regime is the one of an attractor -the dynamics of the model are simply
obtained with the initial value of the scalar field and its derivative at some
moment during this regime-, the models of QI deHaro:2016hpl ; deHaro:2016hsh ;
deHaro:2016ftq  -which generally only depend on two parameters- caught the
attention of some researchers who wanted to confront QI with the observational
data A ; B ; C ; D ; E ; F ; G ; H ; I ; J ; K ; K1 ; K2 ; K3 ; K4 ; K5 ; K6 ;
K7 ; K8 , thus becoming a popular topic in some reduced circles.

In this way, dealing with QI, it is well-known that all the scenarios
containing a period of inflation need a reheating mechanism to match with the
hot BB universe guth  because the particles existing before the beginning of
this period were completely diluted at the end of inflation resulting in a
very cold universe. Here, we will choose as a reheating mechanism the so-
called gravitational particle production Parker ; fmm ; glm ; gmm ; ford ;
Zeldovich ; Damour ; Giovannini  of superheavy particles conformally coupled
with gravity, which was applied to standard inflation (potentials with a deep
well) in kolb ; kolb1 ; Birrell1 ; hashiba ; hashiba1 . However, the
gravitational reheating in QI is normally applied to very light fields
Spokoiny ; pv ; A ; vardayan  and only in few papers, which deal with toy
discontinuous models as the Peebles-Vilenkin one, it is applied to massive
particles H ; ha ; hap18 ; F ; J  (see also Hashiba  where the authors also
deals with the gravitational production of dark matter). In fact, regarding
smooth QI potentials, particle creation and in particular gravitational
particle production is sometimes associated to an overbarrier problem (see for
instance Section VII of kofman ), which has to be analytically studied using
the complex WKB approximation, that is, the Stokes phenomenon hashiba , whose
real application is limited to the creation of particles by parabolic
potentials kofman .

For this reason, in the present work we deal numerically with smooth QI
potentials such as the ones provided by Lorentzian Quintessential Inflation
(LQI) benisty ; benisty1 ; benisty2  and $\alpha$-attractors in Quintessential
Inflation ($\alpha$-QI) benisty3 . Then, in order to get the reheating
temperature, we will calculate numerically the energy density of the
superheavy particles produced during the phase transition from the end of
inflation to the beginning of kination (the period where all the energy of the
field becomes kinetic Joyce ). To do it, we will use the well-known
Hamiltonian diagonalization method (see gmmbook  for a review), showing that
the time dependent $\beta$-Bogoliubov coefficient encodes the polarization
effects associated to the creation and annihilation of the so-called
quasiparticles gmmbook  and also the real particles created during the phase
transition. However, through a toy model inspired in the Peebles-Vilenkin one,
and presenting a discontinuity of the first derivative of the potential at the
beginning of kination, we will show that these polarization effects disappear
when the universe evolves adiabatically, which happens soon after the
beginning of kination. Thus, in order to calculate the energy density of the
produced particles (the real particles), one can safely use the value of the
$\beta$-Bogoliubov coefficient after the beginning of kination. This is the
key point of our investigation, and we have numerically checked that this also
happens for our QI smooth potentials.

Another important point is the overproduction of Gravitational Waves (GWs) in
QI, which are also produced during the phase transition from the end of
inflation to the beginning of kination. This overproduction in many QI models
may disturb the success of the Big Bang Nucleosynthesis (BBN), but, as we will
show, the reheating via gravitational particle production of superheavy
particles in the QI scenarios studied in this work prevents the
incompatibilities of the BBN with the overproduction of GWs, and the reason
why this happens is the fact that the gravitational production of superheavy
particles, contrary to the standard belief that in analytic calculations only
ultraviolet modes are taken into account, is very efficient for long
wavelength modes, leading to a high reheating temperature able to overcome all
the constraints ensuring the BBN success.

Finally, a few words about the viable values of the reheating temperature are
in order. A lower bound for the reheating temperature comes from the fact that
the radiation-dominated era occurs before the BBN epoch, which takes place in
the $1$ MeV regime gkr , and thus, the reheating temperature should naturally
be greater than $1$ MeV. On the contrary, the upper bound of this temperature
is dependent on the theory we are concerned with. In fact, in some
supergravity theories such as $\alpha$-attractors containing particles with
only gravitational interactions, the late time decay of these relics may
jeopardize the success of the standard BBN lindley . To solve this problem one
has to consider sufficiently low reheating temperature (of the order of
$10^{9}$ GeV or less) eln .

The units used in the manuscript are $\hbar=c=1$ and the reduced Planck’s mass
is denoted by $M_{pl}\equiv\frac{1}{\sqrt{8\pi G}}\cong 2.44\times 10^{18}$
GeV.

##  II The diagonalization method

This short section is a review of our previous work hpa  (see also the
pioneering works fmm ; glm ; gmm ; zs  for a more detailed vision of the
topic). The idea of the method goes as follows: given a quantum scalar field
of superheavy particles conformally coupled to gravity, namely $\chi$, the
Klein-Gordon (K-G) equation in the Fourier space, which is satisfied by the
modes in the flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime, is
given by Birrell

$\displaystyle\chi_{k}^{\prime\prime}(\tau)+\omega^{2}_{k}(\tau)\chi_{k}(\tau)=0,$
(1)

where the prime denotes the derivative with respect to the conformal time
$\tau$ and $\omega_{k}(\tau)=\sqrt{k^{2}+m_{\chi}^{2}a^{2}(\tau)}$ is the
time-dependent frequency being $m_{\chi}$ the mass of the quantum field
$\chi$.

As usual, the modes that define the vacuum state at a given initial time
$\tau_{i}$ are the ones that minimize the energy density, so they must satisfy
the conditions

$\displaystyle\chi_{k}(\tau_{i})=\frac{1}{\sqrt{2\omega_{k}(\tau_{i})}}e^{-i\int^{\tau_{i}}\omega_{k}(\bar{\eta})d\bar{\eta}},\quad\chi_{k}^{\prime}(\tau_{i})=-i\omega_{k}(\tau_{i})\chi_{k}(\tau_{i}),$
(2)

and thus, the vacuum expectation value of the energy density will be given by
Bunch

$\displaystyle\langle\rho(\tau)\rangle\equiv\langle
0|\hat{\rho}(\tau)|0\rangle=\frac{1}{4\pi^{2}a^{4}(\tau)}\int_{0}^{\infty}k^{2}dk\left(|\chi_{k}^{\prime}(\tau)|^{2}+\omega^{2}_{k}(\tau)|\chi_{k}(\tau)|^{2}-\omega_{k}(\tau)\right),$
(3)

where, in order to obtain a finite energy density gmmbook , we have subtracted
the energy density of the zero-point oscillations of the vacuum
$\frac{1}{(2\pi)^{3}a^{4}(\tau)}\int d^{3}k\frac{1}{2}\omega_{k}(\tau)$.

Following the method developed in zs  (see also Section $9.2$ of gmmbook ), we
will write the modes as follows,

$\displaystyle\chi_{k}(\tau)=\alpha_{k}(\tau)\frac{e^{-i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}}{\sqrt{2\omega_{k}(\tau)}}+\beta_{k}(\tau)\frac{e^{i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}}{\sqrt{2\omega_{k}(\tau)}},$
(4)

where $\alpha_{k}(\tau)$ and $\beta_{k}(\tau)$ are the time-dependent
Bogoliubov coefficients. Now, imposing that the modes satisfy the condition

$\displaystyle\chi_{k}^{\prime}(\tau)=-i\omega_{k}(\tau)\left(\alpha_{k}(\tau)\frac{e^{-i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}}{\sqrt{2\omega_{k}(\tau)}}-\beta_{k}(\tau)\frac{e^{i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}}{\sqrt{2\omega_{k}(\tau)}}\right),$
(5)

one can show that the Bogoliubov coefficients must satisfy the system

$\displaystyle\left\\{\begin{array}[]{ccc}\alpha_{k}^{\prime}(\tau)&=&\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}e^{2i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}\beta_{k}(\tau)\\\
\beta_{k}^{\prime}(\tau)&=&\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}e^{-2i\int^{\tau}\omega_{k}(\bar{\tau})d\bar{\tau}}\alpha_{k}(\tau),\end{array}\right.$
(8)

in order for the expression (4) to be a solution of the equation (1).

Finally, inserting (4) into the expression for the vacuum energy density (3),
and taking into account that the Bogoliubov coefficients satisfy the equation
$|\alpha_{k}(\tau)|^{2}-|\beta_{k}(\tau)|^{2}=1$, one finds that

$\displaystyle\langle\rho(\tau)\rangle=\frac{1}{2\pi^{2}a^{4}(\tau)}\int_{0}^{\infty}k^{2}\omega_{k}(\tau)|\beta_{k}(\tau)|^{2}dk,$
(9)

where it is important to notice that $|\beta_{k}(\tau)|^{2}$ encodes the
vacuum polarization effects and also the production of particles, which only
happens when the adiabatic evolution breaks. In fact, the quantity

$\displaystyle\langle
N(\tau)\rangle=\frac{1}{2\pi^{2}a^{3}(\tau)}\int_{0}^{\infty}k^{2}|\beta_{k}(\tau)|^{2}dk$
(10)

was named, in the Russian literature, as the number density of quasiparticles
gmmbook , which, as we will see in next section dealing with a toy model, is
very different from the number density of the produced particles because it
also contains the vacuum polarization effects, that is, the creation and
annihilation of pairs.

##  III Particle creation of superheavy particles conformally coupled to
gravity

This section is devoted to the numerical calculation of the energy density of
the gravitationally produced particles in two different QI scenarios with
smooth potentials, namely:

  1. 1.

Lorentzian Quintessential Inflation.-

Based on the well-known Lorentzian distribution, the authors of benisty ;
benisty1  considered following the ansatz

$\displaystyle\epsilon(N)=\frac{\xi}{\pi}\frac{\Gamma/2}{N^{2}+\Gamma^{2}/4},$
(11)

where $\epsilon$ is the main slow-roll parameter, $N$ denotes the number of
e-folds, $\xi$ is the amplitude of the Lorentzian distribution and $\Gamma$ is
its width. From this ansatz, one can find the exact corresponding potential of
the scalar field, namely

$\displaystyle V(\varphi)=\lambda
M_{pl}^{4}\exp\left[-\frac{2\xi}{\pi}\arctan\left(\sinh\left(\gamma\varphi/M_{pl}\right)\right)\right]\bm{\cdot}\left(1-\frac{2\gamma^{2}\xi^{2}}{3\pi^{3}}\frac{1}{\cosh\left(\gamma\varphi/M_{pl}\right)}\right),$
(12)

where $\lambda$ is a dimensionless parameter and the parameter $\gamma$ is
defined by

$\gamma\equiv\sqrt{\frac{\pi}{\Gamma\xi}}.$

Here, to simplify the structure of the potential although without modifying
its properties, we set $\xi=\gamma$ and disregard the term
$\left(1-\frac{2\gamma^{2}\xi^{2}}{3\pi^{3}}\frac{1}{\cosh\left(\gamma\varphi/M_{pl}\right)}\right)$,
which has no influence in the dynamics. Hence, we obtain the simplified
version

$\displaystyle V(\varphi)=\lambda
M_{pl}^{4}\exp\left[-\frac{2\gamma}{\pi}\arctan\left(\sinh\left(\gamma\varphi/M_{pl}\right)\right)\right],$
(13)

plotted in Figure 1, and where in order to match with the current
observational data one has to choose $\lambda\sim 10^{-69}$ and $\gamma\cong
122$ (see for details benisty2 ).

Summing up, in this scenario kination starts when $H\sim H_{kin}\cong 4\times
10^{-8}M_{pl}$ with $\varphi_{kin}\cong-0.03M_{pl}$, inflation ends when
$\varphi_{END}\cong-0.078M_{pl}$ and the pivot scale leaves the Hubble radius
when $\varphi_{*}\cong-0.154M_{pl}$. Finally, at very late times the effective
Equation of State (EoS) parameter goes to $-1$, which leads to an eternal
acceleration.

Figure 1: Plot of the Lorentzian Quintessential Inflation potential.

  2. 2.

Exponential $\alpha$-attractor in Quintessential Inflation.-

The corresponding potential is obtained, combined with a standard exponential
potential, from the following Lagrangian motivated by supergravity and
corresponding to a non-trivial Kähler manifold (see for instance K  and the
references therein),

$\displaystyle\mathcal{L}=\frac{1}{2}\frac{\dot{\phi}^{2}}{(1-\frac{\phi^{2}}{6\alpha})^{2}}M_{pl}^{2}-\lambda
M_{pl}^{4}e^{-\kappa\phi},$ (14)

where $\phi$ is a dimensionless scalar field, and $\kappa$ and $\lambda$ are
positive dimensionless constants.

In order that the kinetic term has the canonical form, one can redefine the
scalar field as follows,

$\displaystyle\phi=\sqrt{6\alpha}\tanh\left(\frac{\varphi}{\sqrt{6\alpha}M_{pl}}\right),$
(15)

obtaining the following potential plotted in Figure 2,

$\displaystyle V(\varphi)=\lambda
M_{pl}^{4}e^{-n\tanh\left(\frac{\varphi}{\sqrt{6\alpha}M_{pl}}\right)},$ (16)

where we have introduced the notation $n\equiv\kappa\sqrt{6\alpha}$, and by
taking $\alpha\sim 10^{-2}$ one has to choose $n\sim 10^{2}$ and $\lambda\sim
10^{-66}$ in order to match with the observational data (see benisty3  for
details).

Finally, it is interesting to note that for this case kination starts later
than in LQI, more precisely, when $\varphi_{kin}\cong-0.5M_{pl}$ with
$H_{kin}\sim 4\times 10^{-7}M_{pl}$. However, the end of inflation and the
horizon crossing occurs earlier than in LQI, $\varphi_{END}\cong-0.89M_{pl}$
and $\varphi_{*}\cong-1.7M_{pl}$ respectively.

Figure 2: Plot of the Exponential $\alpha$-attractor potential.

###  III.1 A toy model

Before dealing with our models and in order to understand better the
gravitational particle production, we warm up revewing a toy potential
inspired in the Peebles-Vilenkin model hpa ,

$\displaystyle
V(\varphi)=\left\\{\begin{array}[]{ccc}\frac{1}{2}m^{2}\left(\varphi^{2}-M_{pl}^{2}+{M}^{2}\right)&\mbox{for}&\varphi\leq-
M_{pl}\\\
\frac{1}{2}m^{2}\frac{M^{6}}{(\varphi+M_{pl})^{4}+M^{4}}&\mbox{for}&\varphi\geq-
M_{pl},\end{array}\right.$ (19)

where $m\sim 5\times 10^{-6}M_{pl}$ is the mass of the inflaton field and
$M\cong 20$ GeV is a very small mass needed to match the theoretical results
provided by the model with the current cosmic acceleration hap19 .

The toy model contains a discontinuity of the first derivative of the
potential at the beginning of kination, i.e., when $\varphi=0$. Then,
according to the conservation equation, the second temporal derivative of the
scalar field is discontinuous at the beginning of kination, as well as the
second temporal derivative of the Hubble parameter. Consequently, the third
derivative of the frequency $\omega_{k}(\tau)$ (which depends on the scale
factor) is discontinuous at the beginning of kination, namely $\tau_{kin}$,
which is the moment when particles are gravitationally created because it is
when the adiabatic evolution is broken.

On the other hand, before performing analytic calculations with this toy
model, the following remark, which helps us to understand the more realistic
ones, is needed: we have to make some assumptions to be sure that the
polarization effects do not affect the evolution of the inflaton field during
the slow-roll period, that is, during the slow-roll regime we will demand that
the polarization effects will be sub-dominant and do not affect the dynamics
of the inflaton field. Analytically, a condition that ensures that
polarization effects are sub-dominant is $H/m_{\chi}\ll 1$, which implies that
$\omega_{k}^{\prime}/\omega_{k}^{2}\ll 1$, i.e., an adiabatic evolution, and
where once again $m_{\chi}$ is the mass of superheavy field $\chi$. Taking
into account that a classical picture of the universe (quantum gravitational
effects could be neglected) appears at GUT scales with $H_{GUT}\cong 5\times
10^{-5}M_{pl}\cong 10^{14}$ GeV, we have to choose superheavy massive fields
with a mass greater or equal than $m_{\chi}\sim 5\times 10^{-4}M_{pl}\cong
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 lindebook .

Therefore, once we have chosen the mass of the $\chi$-field, in order to
obtain the value of the $\beta$-Bogoliubov coefficient we come back to the
equation (8) and, in the first approximation, we take $\alpha_{k}(\tau)=1$,
getting

$\displaystyle\beta_{k}(\tau)=\int^{\tau}\frac{\omega_{k}^{\prime}(\eta)}{2\omega_{k}(\eta)}e^{-2i\int^{\eta}\omega_{k}(\bar{\eta})d\bar{\eta}}d\eta,$
(20)

which, after integration by parts, yields before the beginning of kination

$\displaystyle\beta_{k}(\tau)=\left(-\frac{\omega^{\prime}_{k}(\tau)}{4i\omega_{k}^{2}(\tau)}+\frac{1}{8\omega_{k}(\tau)}\left(\frac{\omega^{\prime}_{k}(\tau)}{\omega_{k}^{2}(\tau)}\right)^{\prime}+\frac{1}{16i\omega_{k}(\tau)}\left(\frac{1}{\omega_{k}(\tau)}\left(\frac{\omega^{\prime}_{k}(\tau)}{\omega_{k}^{2}(\tau)}\right)^{\prime}\right)^{\prime}+....\right)e^{-2i\int^{\tau}\omega_{k}(\bar{\eta})d\bar{\eta}}.$
(21)

However, after the beginning of kination the $\beta$-Bogoliubov coefficient
must be given by

$\displaystyle\beta_{k}(\tau)=\left(-\frac{\omega^{\prime}_{k}(\tau)}{4i\omega_{k}^{2}(\tau)}+\frac{1}{8\omega_{k}(\tau)}\left(\frac{\omega^{\prime}_{k}(\tau)}{\omega_{k}^{2}(\tau)}\right)^{\prime}+\frac{1}{16i\omega_{k}(\tau)}\left(\frac{1}{\omega_{k}(\tau)}\left(\frac{\omega^{\prime}_{k}(\tau)}{\omega_{k}^{2}(\tau)}\right)^{\prime}\right)^{\prime}+....\right)e^{-2i\int^{\tau}\omega_{k}(\bar{\eta})d\bar{\eta}}+C,$
(22)

where the constant $C$ has to be chosen in order that the $\beta$-Bogoliubov
coefficient becomes continuous at $\tau_{kin}$ because the equation (8) is a
first order differential equation, and mathematically it is necessary to
demand the solution to be continuous. Thus, after some cumbersome calculations
(see hpa  for details) one has

$\displaystyle
C=\left(\frac{m_{\chi}^{2}m^{3}a^{5}_{kin}}{16i\omega^{5}_{k}(\tau_{kin})}+....\right)e^{-2i\int^{\tau_{kin}}\omega_{k}(\bar{\eta})d\bar{\eta}},$
(23)

where we have introduced the definition $a_{kin}\equiv a(\tau_{kin})$.

The terms of the $\beta$-Bogoliubov coefficient different from $C$ lead to
sub-leading geometric quantities in the energy density. Effectively, the term
$-\frac{\omega^{\prime}_{k}(\tau)}{4i\omega_{k}^{2}(\tau)}$ leads to the
following contribution to the energy density,
$\frac{m_{\chi}^{2}H^{2}}{96\pi}$, which is negligible compared with
$H^{2}M_{pl}^{2}$. The same happens with
$\frac{1}{8\omega_{k}(\tau)}\left(\omega^{\prime}_{k}(\tau)/\omega_{k}^{2}(\tau)\right)^{\prime}$
leading to a term of order $H^{4}$, which satisfies ${H^{4}}\ll
H^{2}M_{pl}^{2}$. The product of the first and second term generates in the
right-hand side of the modified semi-classical Friedmann equation a term of
the order ${H^{3}m_{\chi}}$, which is also sub-leading compared with
$H^{2}M_{pl}^{2}$. Finally, the third term of (22) leads in the right-hand
side of the semi-classical Friedmann equation to the sub-leading term
$\frac{H^{6}}{m_{\chi}^{2}}$.

Fortunately, this does not happen with $C$, whose leading term gives the main
contribution of the vacuum energy density due to the gravitational particle
production. In fact, the time dependent terms, which as we have already shown
are always sub-leading, are vacuum polarization effects, and they rapidly
disappear in the adiabatic regime, that is, soon after the beginning of
kination $|\beta_{k}(\tau)|^{2}$ approaches to $|C|^{2}$, obtaining

$\displaystyle\langle\rho(\tau)\rangle\cong\left\\{\begin{array}[]{ccc}0&\mbox{
when}&\tau<\tau_{kin}\\\
10^{-5}\left(\frac{m}{m_{\chi}}\right)^{2}m^{4}\left(\frac{a_{kin}}{a(\tau)}\right)^{3}&\mbox{
when}&\tau\geq\tau_{kin},\end{array}\right.$ (26)

which at the beginning of kination is sub-dominant with respect to the energy
density of the inflaton but it will eventually dominate because the one of the
inflaton decreases during kination as $a^{-6}(\tau)$.

Finally, in order to understand better these results it is useful to recall,
as we have already explained, that the authors of the diagonalization method
assume that, during the whole evolution of the universe, quanta named
quasiparticles are created and annihilated due to the interaction with the
quantum field with gravity gmmbook , i.e., this is a vacuum polarization
effect where pairs are created and annihilated. And, following this
interpretation, the number density of the created quasiparticles at time
$\tau$ is given by $\langle
N(\tau)\rangle=\frac{1}{2\pi^{2}a^{3}(\tau)}\int_{0}^{\infty}k^{2}|\beta_{k}(\tau)|^{2}dk$.
However, one has to be very careful with this interpretation and especially
keep in mind that, as we have already pointed out in our toy model (19), real
particles are only created when the adiabatic regime breaks. Effectively,
before the beginning of kination, i.e., before the break of the adiabatic
evolution, the main term of the $\beta$-Bogoliubov coefficient is given by
$-\frac{\omega^{\prime}_{k}(\tau)}{4i\omega_{k}^{2}(\tau)}$, whose
contribution to the energy density is $\frac{m_{\chi}^{2}H^{2}}{96\pi}$, and
to the number density of the quasiparticles $\frac{m_{\chi}H^{2}}{512\pi}$,
and thus, at any time $\tau$ before the beginning of kination
$\langle\rho(\tau)\rangle\not=m_{\chi}\langle N(\tau)\rangle$, meaning that
the quasiparticles do not evolve as real massive particles. On the contrary,
during kination the polarization effects disappear and the leading term of
$\langle N(\tau)\rangle$ is given by
$10^{-5}\left(\frac{m}{m_{\chi}}\right)^{3}m^{3}\left(\frac{a_{kin}}{a(\tau)}\right)^{3}$,
having $\langle\rho(\tau)\rangle=m_{\chi}\langle N(\tau)\rangle$ and a decay
of $a^{-3}(\tau)$, which justifies the interpretation of massive particle
production.

###  III.2  Gravitational particle creation in Lorentzian Quintessential
Inflation

Now we are ready to calculate the particle production for our QI models.
Coming back to the LQI potential (13), first of all we have integrated
numerically the conservation equation for the inflaton field, namely

$\displaystyle\ddot{\varphi}+3H\dot{\varphi}+V_{\varphi}=0,$ (27)

where
$H=\frac{1}{\sqrt{3}M_{pl}}\sqrt{\frac{\dot{\varphi}^{2}}{2}+V(\varphi)}$,
with initial conditions at the horizon crossing, i.e., when the pivot scales
leaves the Hubble radius. Recall that in that moment the system is in the
slow-roll phase and, since this regime is an attractor, one only has to take
initial conditions in the basin of attraction of the slow-roll solution, for
example, $\varphi_{*}=-0.154M_{pl}$ and
$\dot{\varphi}_{*}=-\frac{V_{\varphi}(\varphi_{*})}{3H_{*}}$, where the “star”
denotes that the quantities are evaluated at the horizon crossing.

Once we have obtained the evolution of the background, and in particular the
evolution of the Hubble rate, we compute the evolution of the scale factor,
which is given by

$\displaystyle a(t)=a_{*}e^{\int_{t_{*}}^{t}H(s)ds},$ (28)

where we have chosen as the value of the scale factor at the horizon crossing
$a_{*}=1$.

From the evolution of the scale factor, we can see in Figure 3 that a spike
appears in the plot of the quantity $\omega_{k}^{\prime}/\omega_{k}^{2}$
during the phase transition from the end of inflation to the beginning of
kination, that is, at that moment when the adiabatic evolution is broken and
particles are gravitationally produced.

Figure 3: Plot of the adiabatic evolution for a heavy field with mass
$m_{\chi}\cong 10^{15}$ GeV, when the background is given by the Lorentzian
Quintessential Inflation potential. The value $k=a_{kin}H_{kin}$ has been used
for the quantities $a_{kin}$, $H_{kin}$ of this model. On the left, with the
same scale as the one used in Figure 1 and on the right with a wider range of
values of the scalar field in order to appreciate the extent of the non-
adiabatic region.

Then, from the knowledge acquired studying our toy model (19), we have
numerically solved the equation (8), with initial conditions
$\alpha_{k}(\tau_{*})=1$ and $\beta_{k}(\tau_{*})=0$ at the horizon crossing
(there were neither particles nor polarization effects at that moment because
during the slow-roll regime the derivatives of the Hubble rate are negligible
compared with the powers of $H$, i.e., the system is in the adiabatic regime).

In order to get rid of complex exponentials we have transformed equation (8)
into a second order differential equation, namely

$\displaystyle\left\\{\begin{array}[]{cc}&\alpha_{k}^{\prime\prime}(\tau)=\alpha_{k}^{\prime}(\tau)\left(\frac{\omega_{k}^{\prime\prime}(\tau)}{\omega_{k}^{\prime}(\tau)}-\frac{\omega_{k}^{\prime}(\tau)}{\omega_{k}(\tau)}+2i\omega_{k}(\tau)\right)+\left(\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}\right)^{2}\alpha_{k}(\tau)\\\
&\beta_{k}^{\prime\prime}(\tau)=\beta_{k}^{\prime}(\tau)\left(\frac{\omega_{k}^{\prime\prime}(\tau)}{\omega_{k}^{\prime}(\tau)}-\frac{\omega_{k}^{\prime}(\tau)}{\omega_{k}(\tau)}-2i\omega_{k}(\tau)\right)+\left(\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}\right)^{2}\beta_{k}(\tau)\end{array}\right..$
(31)

Given that $\alpha_{k}(\tau_{*})=1$ and $\beta_{k}(\tau_{*})=0$ leads to
$\alpha_{k}^{\prime}(\tau_{*})=0$, we are interested in solving the equation
for $\alpha_{k}(\tau)$, which can be split into the real and imaginary form in
the following way,

$\displaystyle\left\\{\begin{array}[]{cc}&\alpha_{k,Re}^{\prime\prime}(\tau)=\alpha_{k,Re}^{\prime}(\tau)\left(\frac{\omega_{k}^{\prime\prime}(\tau)}{\omega_{k}^{\prime}(\tau)}-\frac{\omega_{k}^{\prime}(\tau)}{\omega_{k}(\tau)}\right)-2\omega_{k}(\tau)\alpha_{k,Im}^{\prime}(\tau)+\left(\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}\right)^{2}\alpha_{k,Re}(\tau)\\\
&\alpha_{k,Im}^{\prime\prime}(\tau)=\alpha_{k,Im}^{\prime}(\tau)\left(\frac{\omega_{k}^{\prime\prime}(\tau)}{\omega_{k}^{\prime}(\tau)}-\frac{\omega_{k}^{\prime}(\tau)}{\omega_{k}(\tau)}\right)+2\omega_{k}(\tau)\alpha_{k,Re}^{\prime}(\tau)+\left(\frac{\omega_{k}^{\prime}(\tau)}{2\omega_{k}(\tau)}\right)^{2}\alpha_{k,Im}(\tau)\end{array}\right.,$
(34)

and then $|\beta_{k}(\tau)|^{2}=|\alpha_{k}(\tau)|^{2}-1$ because of the well-
known conservation property of the Wronskian. For the value
$k=a_{kin}H_{kin}$, we obtain in Figure 4 that $|\beta_{k}(\tau)|^{2}$
stabilizes soon to a non-zero value after the beginning of kination,
containing only particle production effects. We have numerically verified that
this happens for the range $0.05\lesssim\frac{k}{a_{kin}H_{kin}}\lesssim
7\times 10^{4}$, which leads to values of $|\beta_{k}|^{2}$ of the order of
$10^{-10}$ and $10^{-11}$.

Figure 4: Evolution of $|\beta_{k}(\tau)|^{2}$

Then, introducing these values of the $\beta$-Bogoliubov coefficient in energy
density (equation (9)), we have obtained a vacuum energy density of the order
of $10^{44}\mbox{ GeV}^{4}$. So, the energy density of the produced particles
evolves as

$\displaystyle\langle\rho(\tau)\rangle=\langle\bar{\rho}\rangle\left(\frac{\bar{a}}{a(\tau)}\right)^{3},$
(35)

where $\langle\bar{\rho}\rangle\cong 10^{44}\mbox{ GeV}^{4}\cong 3\times
10^{-30}M_{pl}^{4}$ and $\bar{a}$ are, respectively, the energy density of the
produced particles and the value of the scale factor at the end of the non-
adiabatic phase ($\varphi\sim M_{pl}$, which coincides with the stabilization
of the $\beta$-Bogoliubov coefficient). Finally, the energy density of the
background at this moment is given by
$\bar{\rho}_{\varphi}=3\bar{H}^{2}M_{pl}^{2}\cong 2\times 10^{57}\mbox{
GeV}^{4}\cong 7\times 10^{-17}M_{pl}^{4}$, showing that the energy density of
the produced particles is sub-leading close to the beginning of kination, but
will eventually be dominant because, during the kination regime, the energy
density of the inflaton field decreases as $a^{-6}$.

###### Remark III.1

Note that the values of the wavenumber $k$, that leads to significant values
of the $\beta$-Bogoliubov coefficient, are not in the ultraviolet regime. On
the contrary, analytic calculations only deal with the ultraviolet spectrum
leading to an insignificant value of the $\beta$-Bogoliubov coefficient for
heavy masses. For this reason, without taking into account the long
wavelengths, it is usual to assume that the heavy masses suppress the particle
production, thus leading to an inefficient reheating temperature. However, as
we will see, this is not the case when one performs the numerical calculations
considering all the spectrum of values of the wavenumber.

###### Remark III.2

In the case of the $\alpha$-attractors in the context of Quintessential
Inflation, we have numerically obtained that the corresponding energy
densities are of the order $\langle\bar{\rho}\rangle\sim 3\times
10^{-27}M_{pl}^{4}$ and $\bar{\rho}_{\varphi}\sim 6\times 10^{-17}M_{pl}^{4}$
(see Figure 5)

Figure 5: Plot of the adiabatic evolution(left) and the $\beta$-Bogoliubov
coefficient(right) for a heavy field with mass $m_{\chi}\cong 10^{15}$ GeV,
when the background is given by the exponential $\alpha$-attractor potential.
Here we have used the value $k=a_{kin}H_{kin}$, since we have observed that
particle creation takes place for modes in the range
$1\lesssim\frac{k}{a_{kin}H_{kin}}\lesssim 10^{3}$.

We would like to finish this subsection showing more details about the
function
$\frac{\omega_{k}^{\prime}}{\omega_{k}^{2}}=\frac{m_{\chi}^{2}a^{3}H}{\omega_{k}^{3}}$.
We have calculated its temporal derivative, obtaining:

$\displaystyle\frac{d}{dt}\left(\frac{\omega_{k}^{\prime}}{\omega_{k}^{2}}\right)=\frac{3H^{2}a^{3}m_{\chi}^{2}}{\omega_{k}^{3}}\left(\frac{1-w_{eff}}{2}-\frac{m_{\chi}^{2}a^{2}}{\omega_{k}^{2}}\right),$
(36)

where $w_{eff}$ denotes the effective Equation of State (EoS) paramenter, i.e,
the ratio of the pressure to the energy density. We will easily see that
during kination, i.e. when $w_{eff}=1$, this is a decreasing function, and is
only an increasing function for the modes satisfying

$\displaystyle k>\sqrt{\frac{1+w_{eff}}{1-w_{eff}}}m_{\chi}a.$ (37)

Then, for the relevant models in LQI that contribute to the energy density of
produced particles ($0.05\lesssim\frac{k}{a_{kin}H_{kin}}\lesssim 7\times
10^{4}$), the condition (37) is fulfilled during all the inflationary regime
(see Fig 3 and recall that in LQI inflation ends when
$\varphi_{END}\cong-0.078$), and thus one obtains the spike given in Fig 3.

In addition, for the relevant modes, at the horizon crossing the quantity
$\frac{\omega_{k}^{\prime}}{\omega_{k}^{2}}$ is extremely small and, since it
is an increasing function during inflation, the quantity
$\frac{\omega_{k}^{\prime}}{\omega_{k}^{2}}$ is negligible before the horizon
crossing, that is, during this period there is no particle production and the
polarization effects are negligible, for this reason we can ensure that before
the horizon crossing the relevant modes of the field $\chi$ are in the vacuum.
After the horizon crossing, polarization effects appear and, later, during the
phase transition, particles are created, as we will see in Fig 4. Finally, the
$\beta$-Bogoliubov coefficient stabilizes (at the same time that the non-
adiabatic regime finishes (right plot of Fig 3)) during kination, that is, the
polarization effects are negligible and it only encodes the particle creation.

And exactly the same happens in the $\alpha$-QI scenario.

##  IV The reheating process

After the production of the heavy massive particles with masses around
$10^{15}$ GeV, they have to decay into lighter ones which after the
thermalization process form a relativistic plasma that depicts our hot
universe. Two different situations may arise, as follows:

  1. 1.

The decay, which occurs at time $\tau_{dec}$, is before the end of the
kination regime, when the energy density of the inflaton becomes equal to the
one of the $\chi$-field.

  2. 2.

The decay is after the end of the kination regime.

###  IV.1 Decay before the end of kination

Let us begin the discussion with the LQI potential. In this case, the energy
density of the background, i.e. the one of the inflaton field, and the one of
the relativistic plasma, when the decay is finished, that is when
${\Gamma}\sim H_{dec}=\bar{H}\left(\frac{\bar{a}}{a_{dec}}\right)^{3}\cong
10^{10}\left(\frac{\bar{a}}{a_{dec}}\right)^{3}$ GeV $\sim 5\times
10^{-9}\left(\frac{\bar{a}}{a_{dec}}\right)^{3}M_{pl}$ , will be

$\displaystyle\rho_{\varphi,dec}=3{\Gamma}^{2}M_{pl}^{2}\quad\mbox{and}\quad\langle\rho_{dec}\rangle\cong
3\times 10^{-30}\left(\frac{\bar{a}}{a_{dec}}\right)^{3}M_{pl}^{4}\sim
10^{-21}\Gamma M_{pl}^{3}.$ (38)

Imposing that the end of the decay precedes the end of kination, that means,
$\langle\rho_{dec}\rangle\leq\rho_{\varphi,dec}$, one gets $\Gamma\geq
10^{-21}M_{pl}.$

In addition, since the decay is after the beginning of the kination, one has
$\Gamma\leq H_{kin}\cong 4\times 10^{-8}M_{pl}$. So, we have the following
bound for the decay rate,

$\displaystyle 10^{-21}M_{pl}\leq\Gamma\leq 4\times 10^{-8}M_{pl}.$ (39)

Finally, the reheating temperature, i.e., the temperature of the universe when
the relativistic plasma in thermal equilibrium starts to dominate, which
happens when $\rho_{\varphi,reh}=\langle\rho_{reh}\rangle$, can be calculated
as follows: Since after the decay the evolution of the respective energy
densities is given by

$\displaystyle\rho_{\varphi,reh}=\rho_{\varphi,dec}\left(a_{dec}/a_{reh}\right)^{6},\qquad\langle\rho_{reh}\rangle=\langle\rho_{dec}\rangle\left(a_{dec}/a_{reh}\right)^{4},$
(40)

we will have
$\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,dec}}=\left(a_{dec}/a_{reh}\right)^{2},$
and thus, the reheating temperature will be

$\displaystyle
T_{reh}=\left(\frac{30}{\pi^{2}g_{reh}}\right)^{1/4}\langle\rho_{reh}\rangle^{\frac{1}{4}}=\left(\frac{30}{\pi^{2}g_{reh}}\right)^{1/4}\langle\rho_{dec}\rangle^{\frac{1}{4}}\sqrt{\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,dec}}}\sim
4\times 10^{-17}\left(\frac{M_{pl}}{\Gamma}\right)^{1/4}M_{pl},$ (41)

where $g_{reh}=106.75$ is the effective number of degrees of freedom for the
Standard Model. So, taking into account the bound (39), the reheating
temperature ranges between $7\times 10^{3}$ GeV and $2\times 10^{7}$ GeV.

###  IV.2 Decay after the end of kination

In the case that the decay of the $\chi$-field is after the end of kination,
one has to impose ${\Gamma}\leq H(\tau_{end})\equiv H_{end}$, where we have
denoted by $\tau_{end}$ the time at which kination ends. Taking this into
account, one has

$\displaystyle
H^{2}_{end}=\frac{2\rho_{\varphi,end}}{3M_{pl}^{2}}\quad\mbox{and}\quad\rho_{\varphi,end}=\bar{\rho}_{\varphi}\left(\frac{\bar{a}}{a_{end}}\right)^{6}=\frac{\langle\bar{\rho}\rangle^{2}}{\bar{\rho}_{\varphi}},$
(42)

where we have used that the kination ends when
${\langle{\rho}(\tau_{end})\rangle}={{\rho}_{\varphi,end}}$, meaning
$\left(\bar{a}/a_{end}\right)^{3}=\frac{\langle\bar{\rho}\rangle}{\bar{\rho}_{\varphi}}$.
So, the condition ${\Gamma}\leq H_{end}$ leads to the bound

$\displaystyle\Gamma\leq 10^{-23}M_{pl}.$ (43)

On the other hand, assuming once again instantaneous thermalization, the
reheating temperature (i.e., the temperature of the universe when the
thermalized plasma starts to dominate) will be obtained when all the
superheavy particles decay, i.e. when $H\sim\Gamma$, obtaining

$\displaystyle
T_{reh}=\left(\frac{30}{\pi^{2}g_{reh}}\right)^{1/4}\langle\rho_{dec}\rangle^{1/4}=\left(\frac{90}{\pi^{2}g_{reh}}\right)^{1/4}\sqrt{{\Gamma}M_{pl}}~{},$
(44)

where we have used that, after the end of the kination regime, the energy
density of the produced particles dominates the one of the inflaton field.

Consequently, since the BBN epoch occurs at the $1$ MeV regime and taking once
again $g_{reh}=106.75$, one can find that, in that case, the reheating
temperature is bounded by

$\displaystyle 1\mbox{ MeV}\leq T_{reh}\leq 5\times 10^{6}\mbox{ GeV }.$ (45)

###### Remark IV.1

Since for $\alpha$-attractors in QI the number of produced particles is higher
than in LQI, then the reheating temperature also increases. In fact, following
step by step the same calculations that we have done above, we have obtained a
very high reheating temperature when the decay is before the end of kination,
between $10^{7}$ GeV and $10^{9}$ GeV. And, when the decay is after the end of
kination, the reheating temperature belongs to the range between $1$ MeV and
$10^{9}$ GeV.

##  V Gravitational Waves and the BBN success

###  V.1 The overproduction of GWs

In the seminal paper pv , the authors pointed out that in Quintessential
Inflation a reheating due to the gravitational production of light particles
is incompatible with the overproduction of Gravitational Waves. However, as we
will show in this subsection, for our LQI model, (and also in $\alpha$\- QI,
because in this scenario, as we have already explained, the production of
particles is greater than in LQI), when the reheating is due to the
gravitational creation of superheavy particles, then the overproduction of
GW’s does not disturb the BBN success.

To prove our statement, first of all we recall that the energy density of the
produced GWs during the phase transition from the end of inflation to the
beginning of kination is given by Giovannini99

$\displaystyle\langle\rho_{GW}(\tau)\rangle\cong\frac{H_{kin}^{4}}{2\pi^{3}}\left(\frac{a_{kin}}{a(\tau)}\right)^{4}\cong
10^{-2}H_{kin}^{4}\left(\frac{a_{kin}}{a(\tau)}\right)^{4}.$ (46)

Then, taking into account that the success of the BBN demands that the ratio
of the energy density of GWs to the one of the produced particles at the
reheating time satisfies dimopoulos

$\displaystyle\frac{\langle\rho_{GW,reh}\rangle}{\langle\rho_{reh}\rangle}\leq
10^{-2},$ (47)

we will see that the constraint (47) is overcome when the decay of the
superheavy particles is previous to the end of kination. Effectively, if the
decay occurs before the end of kination one has

$\displaystyle\frac{\langle\rho_{GW,reh}\rangle}{\langle\rho_{reh}\rangle}=\frac{\langle\rho_{GW,dec}\rangle}{\langle\rho_{dec}\rangle},$
(48)

because the energy density of light relativistic particles evolves as the one
of the GW’s. Next, noting that during kination the energy density of the
background scales as $a^{-6}$, it yields that

$\displaystyle\left(\frac{a_{kin}}{a_{dec}}\right)^{4}=\left(\frac{\rho_{\varphi,dec}}{\rho_{\varphi,kin}}\right)^{2/3},$
(49)

and thus, using our previous results (see formula (38)) and recalling that the
value of the Hubble rate at the beginning of kination is $H_{kin}\cong 4\times
10^{-8}M_{pl}$, we get

$\displaystyle\frac{\langle\rho_{GW,reh}\rangle}{\langle\rho_{reh}\rangle}\cong
10^{-1}\left(\frac{\Gamma}{M_{pl}}\right)^{1/3}.$ (50)

Therefore, the bound (47) is overcome when

$\displaystyle\Gamma\leq 7\times 10^{-4}M_{pl},$ (51)

which is completely compatible with the bound (39).

On the other hand, when the decay is produced after the end of kination, and
assuming once again instantaneous thermalization, the reheating time will
coincide with the decay one. Then, since
$\langle\rho_{dec}\rangle=3{\Gamma}^{2}M_{pl}^{2}$ and

$\displaystyle
H_{dec}=H_{end}\left(\frac{a_{end}}{a_{dec}}\right)^{3/2}\Longrightarrow\left(\frac{a_{end}}{a_{dec}}\right)^{3/2}=\sqrt{\frac{3}{2}}\frac{\Gamma
M_{pl}\sqrt{\bar{\rho}_{\varphi}}}{\langle\bar{\rho}\rangle},$ (52)

where we have used that
$H_{end}=\sqrt{\frac{2}{3}}\frac{\langle\bar{\rho}\rangle}{M_{pl}\sqrt{\bar{\rho}_{\varphi}}}$.

Thus, we will have

$\displaystyle\langle\rho_{GW,dec}\rangle=\langle\rho_{GW,end}\rangle\left(\frac{a_{end}}{a_{dec}}\right)^{4}=\langle\rho_{GW,end}\rangle\left(\sqrt{\frac{3}{2}}\frac{\Gamma
M_{pl}\sqrt{\bar{\rho}_{\varphi}}}{\langle\bar{\rho}\rangle}\right)^{8/3}$
(53)

and, using that

$\displaystyle\left(\frac{a_{kin}}{a_{end}}\right)^{4}=\frac{\rho_{\varphi,end}}{\rho_{\varphi,kin}}\quad\mbox{and}\quad\rho_{\varphi,end}=\frac{\langle\bar{\rho}\rangle^{2}}{\bar{\rho}_{\varphi}},$
(54)

we get

$\displaystyle\langle\rho_{GW,dec}\rangle=\left(\frac{2}{16}\right)^{1/3}\left(\frac{H_{kin}}{M_{pl}}\right)^{2}\left(\frac{\bar{\rho}_{\varphi}}{\langle\bar{\rho}\rangle^{2}}\right)^{1/3}(\Gamma
M_{pl})^{8/3}\cong 10^{-2}\Gamma^{8/3}M_{pl}^{4/3},$ (55)

and thus,

$\displaystyle\frac{\langle\rho_{GW,reh}\rangle}{\langle\rho_{reh}\rangle}=\frac{\langle\rho_{GW,dec}\rangle}{\langle\rho_{dec}\rangle}\cong
3\times 10^{-3}\left(\frac{\Gamma}{M_{pl}}\right)^{2/3}\leq 10^{-17},$ (56)

where we have used the bound (43). So, the constraint (47) is clearly
overcome.

###  V.2 BBN constraints from the logarithmic spectrum of GWs

As we have already explained, GWs are produced in the post-inflationary
period, and its logarithmic spectrum of GWs, namely $\Omega_{GW}$ defined as
$\Omega_{GW}\equiv\frac{1}{\rho_{c}}\frac{d\rho_{GW}(k)}{d\ln k}$ (where
$\rho_{GW}(k)$ is the energy density spectrum of the produced GWs;
$\rho_{c}=3H_{0}^{2}M_{pl}^{2}$, where $H_{0}$ is the present value of the
Hubble rate, is the so-called critical density) scales as $k^{2}$ during
kination rubio , 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 maggiore ),

$\displaystyle I\equiv h_{0}^{2}\int_{k_{BBN}}^{k_{end}}\Omega_{GW}(k)d\ln
k\leq 10^{-5},$ (57)

where $h_{0}\cong 0.678$ parametrizes the experimental uncertainty to
determine the current value of the Hubble constant and $k_{BBN}$, $k_{end}$
are the momenta associated to the horizon scale at the BBN and at the end of
inflation respectively. As has been shown in Giovannini1 , the main
contribution of the integral (57) comes from the modes that leave the Hubble
radius before the end of the inflationary epoch and finally re-enter during
kination, that means, for $k_{end}\leq k\leq 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 Giovannini  (see also rubio ;
Giovannini2 ; Giovannini3 ; Giovannini:2016vkr ),

$\displaystyle\Omega_{GW}(k)=\tilde{\epsilon}\Omega_{\gamma}h^{2}_{GW}\left(\frac{k}{k_{end}}\right)\ln^{2}\left(\frac{k}{k_{kin}}\right),$
(58)

where $h^{2}_{GW}=\frac{1}{8\pi}\left(\frac{H_{kin}}{M_{pl}}\right)^{2}$ is
the amplitude of the GWs; $\Omega_{\gamma}\cong 2.6\times 10^{-5}h_{0}^{-2}$
is the present density fraction of radiation, and the quantity
$\tilde{\epsilon}$, 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 rubio ;
Giovannini1  for more details). As has been derived in Giovannini1 , the
specific form of the expression above comes from the behavior of the Hankel
functions for small arguments. Now, plugging expression (58) into (57) and
disregarding the sub-leading logarithmic terms, one finds

$\displaystyle
2\tilde{\epsilon}h_{0}^{2}\Omega_{\gamma}h^{2}_{GW}\left(\frac{k_{kin}}{k_{end}}\right)\leq
10^{-5}\Longrightarrow
10^{-2}\left(\frac{H_{kin}}{M_{pl}}\right)^{2}\left(\frac{k_{kin}}{k_{end}}\right)\leq
1\Longrightarrow 10^{-17}\left(\frac{k_{kin}}{k_{end}}\right)\leq 1,$ (59)

because in our LQI model at the beginning of kination the Hubble rate is
$H_{kin}\sim 4\times 10^{-8}M_{pl}$, and the same happens for our $\alpha$-QI
model.

So, we continue with the LQI, but taking into account that all the reasoning
is also valid for the other model. Then, to calculate the ratio
$k_{kin}/k_{end}$, we will have to study the following two different
situations:

  1. 1.

When the decay occurs before the end of kination.

In this case the reheating time coincides with the end of kination, and thus,
a simple calculation leads to

$\displaystyle\frac{k_{kin}}{k_{end}}=\frac{k_{kin}}{k_{reh}}.$ (60)

Using the formulas

$\displaystyle\rho_{\varphi,reh}=\rho_{\varphi,kin}\left(\frac{a_{kin}}{a_{reh}}\right)^{6}\qquad\mbox{and}\qquad\langle\rho_{reh}\rangle=\langle\rho_{dec}\rangle\left(\frac{a_{dec}}{a_{reh}}\right)^{4},$
(61)

we get

$\displaystyle\frac{a_{kin}}{a_{reh}}=\left(\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,kin}}\right)^{1/6}\left(\frac{a_{dec}}{a_{reh}}\right)^{2/3}=\left(\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,kin}}\right)^{1/6}\left(\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,dec}}\right)^{1/3},$
(62)

where we have used the relation
$\left(\frac{a_{dec}}{a_{reh}}\right)^{2}=\frac{\langle\rho_{dec}\rangle}{\rho_{\varphi,dec}}$.
Then, taking into account that
$H_{reh}=\Gamma\left(\frac{a_{dec}}{a_{reh}}\right)^{3}$, we obtain

$\displaystyle\frac{k_{kin}}{k_{end}}=\frac{k_{kin}}{k_{reh}}=\frac{H_{kin}a_{kin}}{H_{reh}a_{reh}}=\frac{H_{kin}}{\Gamma}\frac{\rho_{\varphi,dec}}{\langle\rho_{dec}\rangle}\left(\frac{\rho_{\varphi,dec}}{\rho_{\varphi,kin}}\right)^{1/6}.$
(63)

Finally, from our previous results (38)

$\displaystyle\rho_{\varphi,dec}=3\Gamma^{2}M_{pl}^{2},\quad\langle\rho_{dec}\rangle\sim
10^{-21}\Gamma
M_{pl}^{3}\quad\mbox{and}\quad\rho_{\varphi,kin}=3H_{kin}^{2}M_{pl}^{2},$ (64)

we arrive at

$\displaystyle\frac{k_{kin}}{k_{end}}\sim
10^{15}\left(\frac{\Gamma}{M_{pl}}\right)^{1/3},$ (65)

and consequently the constraint (59) becomes

$\displaystyle\left(\frac{\Gamma}{M_{pl}}\right)^{1/3}\leq 10^{2},$ (66)

which is obviously overcome for all the viable values of $\Gamma$, i.e, for
all values between $10^{-21}M_{pl}$ and $10^{-8}M_{pl}$.

  2. 2.

When the decay occurs after the end of kination.

Now, since during kination we have
$H_{end}=H_{kin}\left(\frac{a_{kin}}{a_{end}}\right)^{3}$, one gets

$\displaystyle\frac{k_{kin}}{k_{end}}=\left(\frac{H_{kin}}{H_{end}}\right)^{2/3},$
(67)

and taking into account that

$\displaystyle
H_{end}=\sqrt{\frac{2\langle\bar{\rho}\rangle^{2}}{3M_{pl}^{2}\bar{\rho}_{\varphi}}}\sim
3\times 10^{-22}M_{pl}\qquad\mbox{and}\qquad H_{kin}\sim 4\times
10^{-8}M_{pl},$ (68)

the bound (59) is completely overcome.

##  VI Conclusions

In the present work we have numerically studied the gravitational particle
production of superheavy particles conformally coupled to gravity for two
classes of QI scenarios, namely LQI and $\alpha$-QI. To calculate the energy
density of the produced particles we have used the well-known diagonalisation
method, where the key point is the calculation and interpretation of the time-
dependent $\beta$-Bogoliubov coefficient.

In fact, this coefficient encodes all the polarization effects (creation and
annihilation of pairs named, in the Russian literature, quasiparticles) and
the produced superheavy particles -the real ones which after decaying into
lighter ones form a relativistic plasma which reheats the universe- during the
phase transition from the end of inflation to the beginning of kination.
However, and this is the main observation of the present work, the
polarization effects disappear soon after the end of the phase transition
(during the kination regime), i.e., when the evolution of the universe gets
adiabatic again. So, we have numerically checked that effectively the value of
the $\beta$-Bogoliubov coefficient stabilizes during the kination, which
allows us to compute it numerically, and thus, to calculate numerically the
energy density of the particles created during this phase transition.

Once these superheavy particles have been created, they must decay into
lighter ones to form a relativistic plasma which eventually becomes dominant
and matches with the hot Big Bang universe. Then two different situations
arise, namely when the decay occurs before the end of kination regime and when
the decay occurs after the end of the kination regime. Thus, since we have
numerically computed the energy density of the superheavy particles, for both
situations we have been able to calculate the reheating temperature of the
universe, which depends on the decay rate of these superheavy particles and
whose maximum value is quite big -more or less around $10^{7}$ GeV in LQI and
around $10^{9}$ GeV in $\alpha$-QI-, which demystifies the belief, never
checked numerically, that heavy masses suppress the particle production, thus
leading to an abnormally low reheating temperature. In fact, the main
contribution of particle production is to a long wavelength regime, which is
impossible to quantify analytically because only ultraviolet effects can be
calculated with analytic methods, and this is the reason why in many papers
the production of superheavy particles is simply disregarded.

Finally, we have checked that the overproduction of GWs does not distrub the
BBN, because all the bounds preserving its success are clearly overcome.

## Acknowledgments

The investigation of JdH has been supported by MINECO (Spain) grant
MTM2017-84214-C2-1-P, and in part by the Catalan Government 2017-SGR-247.
L.A.S thanks the School of Mathematical Sciences (Queen Mary University of
London) for the support provided.

## References

  * (1) A. Guth, Phys. Rev. D 23, 347 (1981).
  * (2) A. Linde, Phys. Lett. B 108, 389 (1982).
  * (3) G. Chibisov and V. Mukhanov, Mon. Not. Roy. Astron. Soc. 200, 535 (1982).
  * (4) A. A. Starobinsky, Phys. Lett. B 117, 175 (1982).
  * (5) A. H. Guth and S-Y. Pi, Phys. Rev. Lett. 49, 1110 (1982).
  * (6) J. M. Bardeen, P. J. Steinhardt and M. S. Turner, Phys. Rev. D 28, 679 (1983).
  * (7) A. D. Linde, Phys. Lett. 116 B, 335 (1982).
  * (8) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594, A20 (2016) [arXiv:1502.02114 [astro-ph.CO]].
  * (9) K. Bamba, S. Capozziello, S. Nojiri and S.D. Odintsov, Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests, Astrophysics and Space Science 342:155-228 (2012) [arXiv:1205.3421 [gr-qc]].
  * (10) M. Li, X-D Li, S. Wang and Y. Wang, Dark Energy, Commun.Theor.Phys. 56, 525-604 (2011) [arXiv:1103.5870 [astro-ph]].
  * (11) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006) [arXiv:hep-th/0603057].
  * (12) P. J. E. Peebles and A. Vilenkin, Phys. Rev. D 59, 063505 (1999) [arXiv:astro-ph/9810509].
  * (13) J. Haro, J. Amorós and S. Pan, Eur.Phys.J. C79 no. 6, 505 (2019) [arXiv:1901.00167 [gr-qc]].
  * (14) J. de Haro, J. Amorós and S. Pan, Phys. Rev. D 93, no. 8, 084018 (2016) [arXiv:1601.08175 [gr-qc]].
  * (15) J. de Haro and E. Elizalde, Gen. Rel. Grav. 48, no. 6, 77 (2016) [arXiv:1602.03433 [gr-qc]].
  * (16) J. de Haro, Gen. Rel. Grav. 49, no. 1, 6 (2017) [arXiv:1602.07138 [gr-qc]].
  * (17) K. Dimopoulos, Nucl. Phys. Proc. Suppl 95, 70 (2001) [arXiv:astro-ph/0012298].
  * (18) K. Dimopoulos and J. W. F. Valle, Astropart. Phys. 18, 287 (2002) [arXiv:astro-ph/0111417].
  * (19) M. Giovannini, Phys. Rev. D 67, 123512 (2003) [arXiv:hep-ph/0301264].
  * (20) Md. Wali Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, Phys. Rev. D 89, 123513 (2014) [arXiv:1404.1445 [gr-qc]]
  * (21) Md. Wali Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, Int. J. Mod. Phys. D 24, 1530014 (2015) [arXiv:1410.6100 [gr-qc]].
  * (22) J. de Haro, J. Amorós and S. Pan, Phys. Rev. D 94, 064060 (2016) [arXiv:1607.06726 [gr-qc]].
  * (23) C. Q. Geng, C. C. Lee, M. Sami, E. N. Saridakis and A. A. Starobinsky, JCAP 06, 011 (2017) [arXiv:1705.01329 [gr-qc]].
  * (24) L. Aresté Saló and J. de Haro, Eur. Phys. J. C 77, no. 11, 798 (2017) [arXiv:1707.02810 [gr-qc]].
  * (25) J. Haro and S. Pan, Int. J. Mod. Phys. D 27, 1850052 (2018) [arXiv:1512.03033 [gr-qc]].
  * (26) J. Haro, W. Yang and S. Pan, JCAP 01, 023 (2019) [arXiv:1811.07371 [gr-qc]].
  * (27) K. Dimopoulos and C. Owen, JCAP 1706, 027 (2017) [arXiv:1703.00305 [gr-qc]]
  * (28) K. Dimopoulos, Phys. Rev. D68, 123506 (2003) [arXiv:astro-ph/0212264].
  * (29) J.C. Bueno Sanchez and K. Dimopoulos, Phys. Lett. B642, 294-301 (2006); Erratum-ibid. B647, 526 (2007) [arXiv:hep-th/0605258].
  * (30) J.C. Bueno Sanchez and K. Dimopoulos, JCAP 0710, 002 (2007) [arXiv:hep-th/0606223].
  * (31) J.C. Bueno Sanchez and K. Dimopoulos, JCAP 0711, 007 (2007) [arXiv:0707.3967 [hep-ph]].
  * (32) K. Dimopoulos and T. Markkanen, JCAP 06, 021 (2018) [arXiv:1803.07399 [gr-qc]].
  * (33) K. Dimopoulos and L. Donaldson-Wood, Phys. Lett. B796, 26-31 (2019) [arXiv:1906.09648 [gr-qc]].
  * (34) K. Dimopoulos, M. Karciauskas and C. Owen, Phys. Rev. D100, 083530 (2019) [arXiv:1907.04676 [hep-ph]].
  * (35) K. Dimopoulos and S. Sánchez López, Phys. Rev. D103, 043533 (2021) [arXiv:2012.06831 [gr-qc]].
  * (36) L. Parker, Phys. Rev. Lett. 21, 562 (1968); L. Parker, Phys. Rev. 183, 1057 (1969); L. Parker, Phys. Rev. D 3, 346 (1970).
  * (37) V. M. Folov, S. G. Mamayev and V. M. Mostepanenko, Phys. Lett. A 55, 7 (1976).
  * (38) A. A. Grib, B. A. Levitskii and V. M. Mostepanenko, Theoreticheskaya i Matematicheskaya Fizika 19, 59 (1974).
  * (39) A. A. Grib, S. G. Mamayev and V. M. Mostepanenko, Gen. Rel. Grav. 7, 535 (1976); A. A. Grib, S. G. Mamayev and V. M. Mostepanenko, Soviet Physics Journal 17, 1700 (1974).
  * (40) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * (41) Ya B. Zeldovich and A. A. Starobinsky, JETP Lett. 26, 252 (1977).
  * (42) T. Damour and A. Vilenkin, Phys. Rev. D 53, 2981 (1996) [arXiv:hep-th/9503149].
  * (43) M. Giovannini, Phys. Rev. D 58, 083504 (1998) [arXiv:hep-ph/9806329].
  * (44) D. J. H. Chung, E.W. Kolb and A. Riotto, Phys. Rev. D 59, 023501 (1998) [arXiv:hep-ph/9802238].
  * (45) D. J. H. Chung, P. Crotty, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 043503 (2001) [arXiv:hep-ph/0104100].
  * (46) N. D. Birrell and P. C. W. Davies, J. Phys. A: Math. Gen. 13, 2109 (1980)
  * (47) S. Hashiba and J. Yokoyama, JCAP 01, 028 (2019) [arXiv:1809.05410 [gr-qc]].
  * (48) S. Hashiba and Y. Yamada, JCAP 05, 022 (2021) [arXiv:2101.07634 [hep-th]].
  * (49) B. Spokoiny, Phys. Lett. B 315, 40 (1993) [arXiv:gr-qc/9306008].
  * (50) Y. Akrami, R. Kallosh, A. Linde and V. Vardanyan, JCAP 1806, 041 (2018) [arXiv:1712.09693 [hep-th]].
  * (51) J. Haro, J. Amorós and S. Pan, Eur.Phys.J. C 79 no.6, 505 (2019) [arXiv:1901.00167 [gr-qc]].
  * (52) J. de Haro and L. Aresté Saló, Phys. Rev. D 95, 123501 (2017) [arXiv:1702.04212 [gr-qc]].
  * (53) S. Hashiba and J. Yokoyama, Phys. Rev. D99, 043008 (2019) [arXiv:1812.10032 [hep-ph]].
  * (54) L. Kofman, A. Linde and A. Starobinsky, Phys. Rev. D56 3258-3295 (1997) [arXiv:hep-ph/9704452].
  * (55) D. Benisty and E.I. Guendelman, Int.J.Mod.Phys. D29, no. 14, 2042002 (2020) [arXiv:2004.00339 [astro-ph.CO]].
  * (56) D. Benisty and E.I. Guendelman, Eur. Phys. J. C 80, 577 (2020) [arXiv:2006.04129 [astro-ph.CO]].
  * (57) L. Aresté Saló, D. Benisty, E. I. Guendelman and J.d. Haro, JCAP 07, 007 (2021) [arXiv:2102.09514 [astro-ph.CO]].
  * (58) L. Aresté Saló, D. Benisty, E. I. Guendelman and J.d. Haro, Phys. Rev. D103, 123535 (2021) [arXiv:2103.07892 [astro-ph.CO]].
  * (59) M. Joyce, Phys. Rev. D 55, 1875 (1997) [arXiv:hep-ph/9606223].
  * (60) A. A. Grib, S.G. Mamayev and V. M. Mostepanenko, Friedmann Laboratory Publishing for Theoretical Physics, St. Petersburg (1994).
  * (61) G. F. Giudice, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 023508 (2001) [arXiv:hep-ph/0005123].
  * (62) J. Ellis, D.V. Nanopoulos and S. Sarkar, Nuc. Phys. B 259, 175 (1985).
  * (63) J. Ellis, A. Linde and D. Nanopoulos, Phys. Lett. B 118, 59 (1982).
  * (64) J. de Haro, S. Pan, L. Aresté Saló, JCAP 06, 056 (2019) [arXiv:1903.01181 [gr-qc]].
  * (65) Ya B. Zeldodovich and A. A. Starobinsky, Sov. Phys JETP 34, 1159 (1972).
  * (66) N. D. Birrell and C. P. W. Davies, Quantum Fields in Curved Space (Cambridge: Cambridge University Press) (1982).
  * (67) T. S. Bunch, J. Phys. A 13, 1297 (1980).
  * (68) A. Linde, Contemp. Concepts Phys. 5, 1 (2005) [arXiv:hep-th/0503203].
  * (69) M. Giovannini, Phys. Rev. D 60, 123511 (1999) [arXiv:astro-ph/9903004].
  * (70) K. Dimopoulos, L.D. Wood and C. Owen, Phys. Rev. D 97, 063525 (2018) [arXiv:1712.01760 [astro-ph.CO]].
  * (71) J. Rubio and C. Wetterich, Phys. Rev. D 96, 063509 (2017) [arXiv:1705.00552 [gr-qc]].
  * (72) M. Maggiore, Phys. Rept. 331, 283 (2000) [arXiv:gr-qc/9909001].
  * (73) M. Giovannini, Phys. Rev. D 60, 123511 (1999) [arXiv:astro-ph/9903004].
  * (74) M. Giovannini, Phys. Lett. B 759, 528 (2016) [arXiv:1603.09217 [astro-ph.CO]].
  * (75) M. Giovannini, Class. Quant. Grav. 26, 045004 (2009) [arXiv:0807.4317 [astro-ph]].
  * (76) M. Giovannini, Phys. Lett. B 668, 44 (2008) [arXiv:0807.1914 [astro-ph]].