Gravitational production of superheavy baryonic and dark matter in quintessential inflation: nonconformally coupled case

TL;DR

This paper investigates the gravitational creation of superheavy baryonic and dark matter during quintessential inflation, highlighting the importance of coupling differences and establishing a mass threshold for dark matter viability.

Contribution

It introduces a model where baryonic matter and dark matter have different gravitational couplings during inflation, providing new insights into superheavy dark matter production.

Findings

Superheavy dark matter mass must exceed 8 x 10^{15} GeV for viability.

Nonconformal coupling of baryonic matter influences reheating processes.

Dark matter conformally coupled with gravity affects production mechanisms.

Abstract

The gravitational production of superheavy dark matter is studied in the context of quintessential inflation. The superheavy particles, whose decay products are baryonic matter and are the responsible for the reheating of the universe after the end of the inflationary period, are not conformally coupled with gravity. On the contrary, dark matter is assumed to be conformally coupled with gravity. We show that the viability of these scenarios requires the mass of the superheavy dark matter to be greater than 8 x 10^{15} GeV.

Equations

$$\displaystyle V(\varphi)=\left\{\begin{array}[]{ccc}\lambda M_{pl}^{4}\left(1-e^{\sqrt{\frac{2}{3}}\frac{\varphi}{M_{pl}}}+\frac{M^{2}}{M_{pl}^{2}}\right)^{2}&\mbox{for}&\varphi\leq 0\\ \lambda\frac{M^{8}}{\varphi^{4}+M^{4}}&\mbox{for}&\varphi\geq 0,\end{array}\right.$$

$${\mathcal{L}}_{JF}=\frac{M_{pl}^2}{2}\left(R+\alpha R^2\right)a^3,$$

$$\tilde{a}=a\sqrt{1+2\alpha R},d\tilde{t}=dt\sqrt{1+2\alpha R}.$$

$${\mathcal{L}}_{EF}=\left(\frac{M_{pl}^2}{2}\tilde{R}+\frac{1}{2}({\tilde{\phi }}^{\prime }{)}^2-V(\tilde{\phi })\right){\tilde{a}}^3,$$

$$\tilde{\phi }=-\sqrt{\frac{3}{2}}{M}_{pl}\ln (1+2\alpha R),V(\tilde{\phi })=\frac{\alpha R^2{M}_{pl}^2}{2(1+2\alpha R{)}^2}.$$

$$V(\tilde{\phi })=\frac{M_{pl}^2}{8\alpha }{\left(1-e^{\sqrt{\frac{2}{3}}\frac{\tilde{\phi }}{M_{pl}}}\right)}^2.$$

$$1-n_s\cong 6{ϵ}_{\ast }-2{\eta }_{\ast }\cong \frac{8}{3}{e}^{\sqrt{\frac{2}{3}}\frac{{\phi }_{\ast }}{M_{pl}}},$$

$${\phi }_{\ast }\cong \sqrt{\frac{3}{2}}{M}_{pl}\ln \left(\frac{3}{8}(1-n_s)\right).$$

$$\lambda \sim 9{\pi }^2(1-n_s{)}^2\times 10^{-9}.$$

$$\lambda =9\times 10^{-11}\text{mbox}and{\phi }_{\ast }\cong -5.42M_{pl}.$$

$$\ddot{\phi }+3\sqrt{\frac{\frac{{\dot{\phi }}^2}{2}+V(\phi )}{3M_{pl}^2}}\dot{\phi }+V_{\phi }=0$$

$$H_{kin}=\frac{{\dot{\phi }}_{kin}}{\sqrt{6}{M}_{pl}}\cong 1.44\times 10^{-6}{M}_{pl},$$

$${\rho }_{\phi ,kin}\cong 6.26\times 10^{-12}{M}_{pl}^4.$$

$$N={\int }_{t_{\ast }}^{t_{END}}Hdt=\frac{1}{M_{pl}}{\int }_{{\phi }_{\ast }}^{{\phi }_{END}}\frac{1}{\sqrt{2ϵ}}d\phi$$

$$N\cong \frac{3}{4}\left(\frac{8}{3(1-n_s)}+\frac{1}{3-2\sqrt{3}}+\ln \left(\frac{3}{8}\frac{n_s-1}{3-2\sqrt{3}}\right)\right),$$

$$|{\beta }_k(\tau ){|}^2\cong \frac{m_Y^4{a}_{kin}^{12}(\stackrel{...}{H}({\tau }_{kin}^{-})-\stackrel{...}{H}({\tau }_{kin}^{+}){)}^2}{1024{\omega }_k^{12}({\tau }_{kin})},$$

$${\chi }_k^{\prime \prime }+{\Omega }_k^2{\chi }_k=0,$$

$$|{\beta }_k(\tau ){|}^2\cong \frac{9(\xi -1/6{)}^2{a}_{kin}^8(\stackrel{...}{H}({\tau }_{kin}^{-})-\stackrel{...}{H}({\tau }_{kin}^{+}{)}^2}{32{\omega }_k^8({\tau }_{kin})}.$$

$$(\stackrel{...}{H}({\tau }_{kin}^{-})-\stackrel{...}{H}({\tau }_{kin}^{+}){)}^2={\dot{\phi }}_{kin}^4{\left(\frac{V_{\phi \phi }(0^{-})}{M_{pl}^2}\right)}^2$$

$$=\frac{16{\lambda }^2}{9}{\dot{\phi }}_{kin}^4$$

$${\rho }_A(\tau )\cong \frac{m_A}{2{\pi }^2{a}^3(\tau )}{\int }_0^{\infty }{k}^2|{\beta }_k(\tau ){|}^{2}dk$$

$${\rho }_X(\tau )\cong \frac{{\lambda }^2}{128\pi }{\left(\xi -\frac{1}{6}\right)}^2{\left(\frac{{\dot{\phi }}_{kin}}{m_X}\right)}^4{\left(\frac{a_{kin}}{a(\tau )}\right)}^3,$$

$${\rho }_Y(\tau )\cong \frac{7{\lambda }^2}{589824\pi }{\left(\frac{{\dot{\phi }}_{kin}}{m_Y}\right)}^4{\left(\frac{a_{kin}}{a(\tau )}\right)}^3.$$

$${\rho }_Y(\tau )=\frac{7}{4608{\left(\xi -\frac{1}{6}\right)}^2}{\left(\frac{m_X}{m_Y}\right)}^4{\rho }_X(\tau ),$$

$${\rho }_Y(\tau )\cong 1.5\times 10^{-3}{\left(\frac{m_X}{m_Y}\right)}^4{\rho }_X(\tau ).$$

$${\rho }_{Y,rh}=1.5\times 10^{-3}{\left(\frac{m_X}{m_Y}\right)}^4{\rho }_{X,rh},$$

$${\rho }_X(\tau )={\rho }_{X,rh}{\left(\frac{a_{rh}}{a(\tau )}\right)}^4,{\rho }_Y(\tau )={\rho }_{Y,rh}{\left(\frac{a_{rh}}{a(\tau )}\right)}^3,$$

$$\frac{a_{rh}}{a_{eq}}=\frac{{\rho }_{Y,rh}}{{\rho }_{X,rh}}\cong 1.5\times 10^{-3}{\left(\frac{m_X}{m_Y}\right)}^4,$$

$${\rho }_{Y,eq}\cong 5\times 10^{-12}{\rho }_{X,rh}{\left(\frac{m_X}{m_Y}\right)}^{16}$$

$$=\frac{{\pi }^2{g}_{\ast }}{6}\times 10^{-12}{T}_{rh}^4{\left(\frac{m_X}{m_Y}\right)}^{16},$$

$$T_{rh}\cong 2.2\times 10^{-7}{\left(\frac{m_Y}{m_X}\right)}^4\text{mbox}GeV.$$

Full text

Superheavy particles; Dark matter; Reheating; Quintessential inflation.

I Introduction

Quintessential inflation, which was addressed for the first time by Peebles and Vilenkin (PV) in pv , is an attempt to unify inflation and quintessence via a single scalar field whose potential allows inflation at early times while at late time provides quintessence (see for instance dimopoulos and references therein). A remarkable property of the PV model is that it contains an abrupt phase transition from inflation to kination (a regime where all the energy density of the inflation turns into kinetic), where the adiabatic regime is broken and, thus, particles could be gravitationally created ford ; Damour . This leads to the possibility to explain the abundance of dark matter through the gravitational production of superheavy particles during the phase transition in quintessential inflation hashiba ; hashiba1 , or during the oscillations of the inflaton field in standard inflation kolb1 ; ema ; kolb2 .

The potential of the model presented here depends on two parameters which are determined using observational data: one with the observational value of the power spectrum of scalar perturbations and the other one taking into account that the ratio of the energy density of the scalar field to the critical energy density at the present time is approximately $0.7$. Moreover, this potential is obtained matching a Starobinsky inflationary-type potential with the inverse power law potential used in pv . The former one leads to theoretical values of the spectral index and the ratio of tensor to scalar perturbations agreeing with the recent observational data provided by the Planck’s team planck18 , and the second one is responsible for the current cosmic acceleration.

Since the potential has an abrupt phase transition at the end of the inflationary phase, we will consider the gravitational production of two kinds of superheavy particles: $X$-particles, nonconformally coupled with gravity, whose energy density after their decay into baryonic light particles and later thermalization of decay products will dominate the energy density of the scalar field in order to match with the Hot Big Bang (HBB), and dark $Y$-particles, conformally coupled with gravity, which are only gravitationally interacting massive particles (GIMP). We will show that the quintessential inflation model presented in this work preserves the Big Bang Nucleosynthesis (BBN) success, in the sense that the overproduction of Gravitational Waves (GWs) does not disturb the BBN for $X$-particles and $Y$-particles with masses in the range of $10^{15}-10^{17}$ GeV and $10^{16}-10^{18}$ GeV respectively, leading to a maximum reheating temperature in the TeV regime.

The paper is organized as follows: In Section II we present our quintessential inflation model based on a Starobinsky Inflation-type potential matched with a quartic inverse power law potential. Section III is devoted to the calculation of the energy density of the superheavy produced particles and to give viable bounds for the reheating temperature and for the masses of $X$ and $Y$ particles. In Section IV numeric calculation has been performed in order to show the viability of the model at the present time and its future evolution. Finally, we present the conclusion of the work in Section V.

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

II The quintessential inflation model

It is well-known that in quintessential inflation the number of e-folds from the pivot scale exiting the Hubble radius to the end of inflation is greater than $60$. For this reason, in order that the theoretical values of the spectral index and the ratio of tensor to scalar perturbations enter in their marginalized joint confidence contour in the plane $(n_{s},r)$ at $2\sigma$ C.L. for the Planck2018 TT, TE, EE + low E+ lensing + BK14 + BAO likelihoods planck18 , we have changed the quartic inflationary potential of the original PV quintessential inflation model pv by a Starobinsky-type potential in the Einstein Frame (EF) riotto ; aho (also named Higgs Inflation martin ), obtaining:

[TABLE]

where $\lambda$ is a dimensionless parameter which we will calculate right now and, as we will show in Section IV, $M\cong 2.6\times 10^{5}$ GeV is a small mass.

Remark II.1

To show the equivalence of $R^{2}$-gravity in the Jordan Frame (JF) and the first piece of the potential (3) in the EF, we consider, in the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, the Lagrangian of $R^{2}$-gravity in the JF (see for instance odintsov )

[TABLE]

where $\alpha$ is a positive parameter with dimension of $M_{pl}^{-2}$.

To work in the EF, we perform the change of variable aho

[TABLE]

Then, the Lagrangian in the EF becomes

[TABLE]

where ′ denotes the derivative with respect to $\tilde{t}$, the Ricci scalar in the EF is $\tilde{R}=6(\tilde{H}+2\tilde{H}^{2})$, with $\tilde{H}=\tilde{a}^{\prime}/\tilde{a}$, and the relation between both frames is given by

[TABLE]

Therefore, since $1+2\alpha R=e^{-\sqrt{\frac{2}{3}}\frac{\tilde{\varphi}}{M_{pl}}}$, we conclude that $R^{2}$-gravity in the JF is equivalent to General Relativity (GR) in the EF when the potential is given by

[TABLE]

On the other hand, the tail of the potential $(\varphi>0)$, which is motivated by SUSY QCD SUSY , is the same used by Peebles and Vilenkin in his seminal paper [1], and has been studied in many papers dealing with quintessence such as rp ; albrecht .

In this model, the kination phase starts at $\varphi_{kin}\cong 0$. Thus, to obtain the value of the Hubble parameter at that time, namely $H_{kin}$, first of all we calculate the slow roll parameters: Denoting by $\epsilon_{*}=\frac{M_{pl}^{2}}{2}\left(\frac{V_{\varphi}(\varphi_{*})}{V(\varphi_{*})}\right)^{2}$ and $\eta_{*}={M_{pl}^{2}}\frac{V_{\varphi\varphi}(\varphi_{*})}{V(\varphi_{*})}$ the values of the slow roll parameters and by $\varphi_{*}$ the value of the scalar field when the pivot scale exits the Hubble radius, since the mass $M$ satisfes $M\ll M_{pl}$, one has $\epsilon_{*}\cong\frac{4}{3}e^{2\sqrt{\frac{2}{3}}\frac{\varphi_{*}}{M_{pl}}}$ $\eta_{*}=-\frac{4}{3}e^{\sqrt{\frac{2}{3}}\frac{\varphi_{*}}{M_{pl}}},$ and thus, the spectral index is given by btw

[TABLE]

meaning that

[TABLE]

On the other hand, the observational estimation of the power spectrum of the scalar perturbations when the pivot scale leaves the Hubble radius is ${\mathcal{P}}_{\zeta}\cong\frac{H_{*}^{2}}{8\pi^{2}M_{pl}^{2}\epsilon_{*}}\sim 2\times 10^{-9}$ btw . Since during the slow roll regime the kinetic energy density is negligible compared with the potential one, we will have $H_{*}^{2}\cong\frac{\lambda}{3}M_{pl}^{2}$, and using the relation $\epsilon_{*}=\frac{3}{16}(1-n_{s})^{2}$ one gets

[TABLE]

Taking into account that the observational value of the spectral index is $n_{s}=0.968\pm 0.006$ Planck , if one chooses its central value one gets

[TABLE]

Then, once we have these quantities we can solve numerically the conservation equation

[TABLE]

with initial conditions $\varphi_{*}=-5.42M_{pl}$ and $\dot{\varphi}_{*}=0$ (obviously, one can choose other similar initial conditions and the result has to be practically the same because the inflationary dynamics are that of an attractor).

Using event-driven integration with an ode RK78 integrator one gets $\dot{\varphi}_{kin}=3.54\times 10^{-6}M_{pl}^{2}$, and thus

[TABLE]

and

[TABLE]

To end this section, let’s calculate the number of e-folds between the time when $\varphi=\varphi_{*}$ and $\varphi=\varphi_{END}$ (i.e. the end of inflation) provided by our model

[TABLE]

So, using the value of $\varphi_{*}$ above, that $\epsilon\cong\frac{4}{3}\left(\frac{s}{1-s}\right)^{2}$, where $s=e^{\sqrt{\frac{2}{3}}\frac{\varphi}{M_{pl}}}$, and that $s_{END}\cong-3+2\sqrt{3}$ (which corresponds to $\epsilon_{END}=1$), one gets that

[TABLE]

which leads to $41.34\leq N\leq 95.29$ for the values of $0.956\leq n_{s}\leq 0.98$ within its $2\sigma$ C.L. In particular, at $1\sigma$ C.L., i.e., for the values $0.969\leq n_{s}\leq 0.975$, the expected number of e-folds in quintessential inflation, is between $60$ and $75$.

III Reheating via gravitational particle production

Since the second derivative of the potential (3) is discontinuous at $\varphi=0$, from the conservation equation one can see that the third temporal derivative of the inflation field is discontinuous at the beginning of kination, and using the Raychaudhuri equation $\dot{H}=-\frac{\dot{\varphi}^{2}}{2M_{pl}^{2}}$ one can deduce that at the beginning of kination the third derivative of the Hubble parameter is discontinuous, enhancing the particle production as discussed in kolb . Then, in order that vacuum polarization effects do not disturb the dynamics of the $\varphi$-field, the mass of the superheavy particles, produced gravitationally, must be greater than $10^{15}$ GeV, where we have assumed that the beginning of inflation occurs at GUT scales, that is, when the Hubble parameter is of the order of $10^{14}$ GeV (see for instance hyp ). Therefore, for the $Y$-particles, which we assume to be conformally coupled with gravity, since $m_{Y}\gg H$ one can safely use the WKB approximation (see section $2$ of bunch for a detailed explanation) to calculate the $\beta$-Bogoliubov coefficient of the $k$-mode hap1 , leading for our model to

[TABLE]

where $\tau_{kin}$ denotes the beginning of the kination in conformal time, $\omega_{k}(\tau)=\sqrt{k^{2}+a^{2}(\tau)m_{Y}^{2}}$ is the time dependent frequency of the $k$-mode and the third derivative of the Hubble parameter is evaluated on the right $(+)$ and on the left $(-)$ of $\tau_{kin}$.

Remark III.1

In hps the calculation of the $\beta$-Bogoliubov coefficient was done using the well-known diagonalization method Grib ; Zeldovich , and the importance of the discontinuity of some derivative (in our case the second one) of the potential at the phase transition is pointed out. In fact, the greater the order of the discontinuous derivative is, the less the number density of superheavy gravitationally produced particles hap is, which is in agreement with kolb . So, for a smooth phase transition the production of superheavy particles would be suppressed and its energy density would be abnormally small, meaning that in such a model the reheating is impossible via gravitational production of superheavy particles and, thus, other mechanisms of reheating, such as ”instant preheating” fkl ; fkl1 , must be invoked.

On the contrary, for the $X$-particles, which are nonconformally coupled with gravity, we have that the $k$-mode satisfies the equation bunch

[TABLE]

where $\Omega_{k}^{2}=\omega_{k}^{2}+(\xi-\frac{1}{6})a^{2}R$, being $\xi$ the coupling constant, $\omega_{k}(\tau)=\sqrt{k^{2}+a^{2}(\tau)m_{X}^{2}}$ and $R$ the Ricci scalar. At this point, one has to note that the WKB is a perturbative approximation which holds when $m_{X}\gg|\xi-\frac{1}{6}|R$, and thus, since at the GUT scales one has $R\sim 10^{29}\mbox{ GeV}^{2}$ so that the mass $m_{X}$ is far from the Planck’s mass, one has to choose $|\xi-\frac{1}{6}|\leq 1$, and the square of the $\beta$-Bogoliubov is given by

[TABLE]

Therefore, taking into account that

[TABLE]

and the fact that the energy density of $A$-particles, with $A=X,Y$, is given by

[TABLE]

before the decay of the $X$-particles, its energy density evolves as

[TABLE]

and the one of the $Y$-particles evolves as

[TABLE]

Thus, before the decay of the $X$-particles, one will have

[TABLE]

and, assuming that $|\xi-\frac{1}{6}|\cong 1$ so that the energy density of the $X$-particles is the maximum possible, we will have

[TABLE]

Now, it is important to take into account that, when reheating is due to the gravitational production of superheavy particles, in order that the overproduction of GWs does not alter the BBN success, the decay of these particles has to take place after the end of kination hyp . Then, assuming as usual instantaneous thermalization, the reheating is produced immediately after the decay of the $X$-particles, obtaining

[TABLE]

where the subindex “rh” means that the quantities are evaluated at the reheating time. After reheating, the evolution of the corresponding energy densities will be

[TABLE]

meaning that at the matter-radiation equality

[TABLE]

and consequently

[TABLE]

where $T_{rh}$ denotes the reheating temperature and $g_{*}=106.75$ are the degrees of freedom for the Standard Model.

On the other hand, considering the central values obtained in planck of the red-shift at the matter-radiation equality $z_{eq}=3365$, the present value of the ratio of the matter energy density to the critical one $\Omega_{m,0}=0.308$, and $H_{0}=67.81\;\mbox{Km/sec/Mpc}\cong 1.42\times 10^{-33}$ eV, one can deduce that the present value of the matter energy density is $\rho_{m,0}=3H_{0}^{2}M_{pl}^{2}\Omega_{m,0}=3.26\times 10^{-121}M_{pl}^{4}$, and at the matter-radiation equality one will have $\rho_{m,eq}=\rho_{m,0}(1+z_{eq})^{3}=4.4\times 10^{-1}\mbox{eV}^{4}$. Since practically all the matter has a non-baryonic origin, one can conclude that $\rho_{Y,eq}\cong\rho_{m,eq}$, meaning that the reheating temperature is given by a function of $m_{Y}/m_{X}$ as follows:

[TABLE]

III.1 Decay after the end of the kination regime

As we have already explained in the previous section, in order that the overproduction of GWs does not alter the BBN success, the decay of the $X$-particles has to be produced after the end of kination, which occurs when the energy density of the inflaton field is equal to the one of the $X$-particles. Then, the decaying rate, namely $\Gamma$, has to satisfy ${\Gamma}\leq H(\tau_{end})\equiv H_{end}$, where we have denoted by $\tau_{end}$ the time at which kination ends. Therefore, one has

[TABLE]

and

[TABLE]

in which, taking into account that during kination the energy density of the inflaton field decays as $a^{-6}$ and the one of the produced particles as $a^{-3}$, we have introduced the so-called heating efficiency defined in rubio as

[TABLE]

Consequently, (32) leads to $H_{end}=\sqrt{2}H_{kin}\Theta$, and from the constraint $\Gamma\leq H_{end}$ one obtains the bound

[TABLE]

On the other hand, assuming once again instantaneous thermalization, the energy density of the $X$-particles at the reheating time will be $\rho_{X,rh}=3{\Gamma}^{2}M_{pl}^{2}$, and thus, the reheating temperature will be given by

[TABLE]

As a consequence, from the two expressions of the reheating temperature (31) and (36) one can write the mass of the dark matter as a function of $\Gamma$ and $m_{X}$ as follows:

[TABLE]

III.2 Overproduction of GWs

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 hossain3

[TABLE]

where the energy density of the GWs is given by $\rho_{GW}(\tau)\cong 10^{-2}H^{4}_{kin}\left(\frac{a_{kin}}{a(\tau)}\right)^{4}$ (see for instance ford ).

Thus, since

[TABLE]

we will have,

[TABLE]

meaning that the bound (38) leads to the constraint

[TABLE]

Here, it is important to recall that, in order to apply the WKB approximation, we have assumed that the mass of the particles is greater than the Hubble parameter at the beginning of inflation, which is of the order of $10^{14}$ GeV if inflation starts at GUT scales. Therefore, choosing $m_{X}\geq 10^{15}$ GeV one can easily show that the constraint (41) automatically implies (35), and thus, taking into account that $T_{rh}>1$ MeV because the BBN occurs at the MeV regime gkr , one gets that $\Gamma$ must satisfy

[TABLE]

which always holds when

[TABLE]

Consequently, from (36) and (42), for our model the reheating temperature is bounded by

[TABLE]

and from (37) and (42) the mass of the $Y$-particles by

[TABLE]

Then, choosing for example $m_{X}=10^{15}$ GeV, one gets the following bound for the reheating temperature

[TABLE]

and from (31), if one assumes that the universe reheats when the temperature is around $1$ GeV, the mass of the $Y$-particles has to be $m_{Y}\cong 4.6\times 10^{16}$ GeV. In general, for $m_{X}=10^{15}$ GeV the relation between the mass of the particles that generate dark matter and the reheating temperature is presented in Figure 1.

To end this section, a final remark is in order: When one considers that the $X$-particles are conformally coupled with gravity, the relation (31) becomes

[TABLE]

Then, for $m_{X}=10^{15}$ GeV and a reheating temperature of $1$ GeV, one gets $m_{Y}\cong 2.3\times 10^{17}$ GeV, which means that the mass of the $Y$- particles is increased in one order with respect to the nonconformally coupled case. This shows that, if one wants a model with elementary superheavy $X$ and $Y$ far from the Planck scale, one has to consider that the $X$-particles -the ones which decay into light baryonic matter- do not have to be conformally coupled with gravity.

IV Numerical calculations

In this section we want to calculate the value of the parameter $M$ as a function of the reheating temperature and the late time evolution of our model.

IV.1 Analytic results

To perform this calculation, first of all, as we have already shown at the end of Section $2$, we take as initial conditions at the beginning of kination

[TABLE]

During kination one can safely disregard the potential, so during this phase one has $a\propto t^{1/3}\Longrightarrow H=\frac{1}{3t}$, and using the Friedmann equation, the dynamics in this regime will be

[TABLE]

Then, at the end of kination, one has

[TABLE]

and using once again that $H_{end}=\sqrt{2}H_{kin}\Theta$, one gets

[TABLE]

During the period between $t_{end}$ and $t_{rh}$ the universe is matter dominated and, thus, the Hubble parameter becomes $H=\frac{2}{3t}$. Since the gradient of the potential could also be disregarded at this epoch, hence, the equation of the scalar field becomes $\ddot{\varphi}+\frac{2}{t}\dot{\varphi}=0$, and thus, at the reheating time

[TABLE]

and

[TABLE]

Note that for the allowed reheating temperatures, i.e., for temperatures satisfying $T_{rh}\leq 66$ TeV, one has ${\frac{\pi}{6}\sqrt{\frac{g_{*}}{10}}\frac{T_{rh}^{2}}{M_{pl}H_{kin}\Theta}}\ll 1$, so we can safely make the approximation

[TABLE]

During the radiation period one can continue disregarding the potential, obtaining

[TABLE]

and thus, at the matter-radiation equality one has

[TABLE]

where $g_{eq}\cong 3.36$ are the degrees of freedom at this scale gr and $T_{eq}$ is the temperature of the radiation at the matter-radiation equilibrium, which is related with the energy density via the relation $\rho_{eq}=\frac{\pi^{2}}{15}g_{eq}T^{4}_{eq}$, and thus, given by $T_{eq}\cong 7.8\times 10^{-10}$ GeV.

In the same way,

[TABLE]

After the matter-radiation equality the dynamical equations can not be solved analytically and, thus, one needs to use numerics to compute them. In order to do that, we need to use a “time” variable that we choose to be the number of $e$-folds up to the present epoch, namely, $N\equiv-\ln(1+z)=\ln\left(\frac{a}{a_{0}}\right)$. Now, using the variable $N$, one can recast the energy density of radiation and matter respectively as

[TABLE]

and

[TABLE]

where the value of the energy density at the matter-radiation equality $\rho_{eq}\cong 8.8\times 10^{-1}$ ${\mbox{e}V}^{4}$ has been obtained in the previous Section III and one can also understand that $N_{eq}$ is the value of $N$ at the matter-radiation equality.

IV.2 The dynamical system

In order to obtain the dynamical system for this scalar field model, we introduce the following dimensionless variables

[TABLE]

where $H_{0}\sim 1.42\times 10^{-33}$ eV denotes once again the current value of the Hubble parameter. Now, using the variable $N=-\ln(1+z)$ defined above and also using the conservation equation $\ddot{\varphi}+3H\dot{\varphi}+V_{\varphi}=0$, one can construct the following non-autonomous dynamical system:

[TABLE]

where the prime represents the derivative with respect to $N$, $\bar{H}=\frac{H}{H_{0}}$ and $\bar{V}=\frac{V}{H_{0}^{2}M_{pl}^{2}}$. Moreover, the Friedmann equation now looks as

[TABLE]

where we have introduced the following dimensionless energy densities $\bar{\rho}_{r}=\frac{\rho_{r}}{H_{0}^{2}M_{pl}^{2}}$ and $\bar{\rho}_{m}=\frac{\rho_{m}}{H_{0}^{2}M_{pl}^{2}}$.

Then, we have to integrate the dynamical system, starting at $N_{eq}=-8.121$, with initial condition $x_{eq}$ and $y_{eq}$ which are obtained analytically in the previous subsection. The value of the parameter $M$ is obtained equaling at $N=0$ the equation (64) to $1$, i.e., imposing $\bar{H}(0)=1$.

Numerical calculations show that $M\cong 2.6\times 10^{5}$ GeV, independently of the reheating temperature, which is a value of the same order as the one obtained in pv . On the other hand, in Figure 2 we have drawn the evolution of the different dimensionless energy densities, obtaining a frozen quintessence, that is, the energy density of the scalar field is frozen and starts to dominate very close to the present time. We have also considered the evolution of the ratio of the energy density to the critical one for the different constituents, i.e., $\Omega_{B}(t)=\frac{\rho_{B}(t)}{3H^{2}(t)M_{pl}^{2}}$ where $B=r,m,\varphi$. And in Figure 3 we have depicted the evolution of the effective EoS parameter $\omega_{eff}(t)=-1-\frac{2\dot{H}(t)}{3H^{2}(t)}$. Finally, we can see that for $N$ greater than $1.5$ one has $\Omega_{\varphi}=1$ and $\omega_{eff}=-1$ meaning that, at late times, the universe enters in a de Sitter phase and, thus, exhibits an eternal acceleration.

V Conclusions

In the present work we have considered a quintessential model whose potential, which only depends on two parameters, is composed by a Starobinsky Inflationary type-potential matched with an inverse power law potential, which is responsible for quintessence. Since the phase transition from the end of inflation to the beginning of kination is very abrupt, the adiabatic regime is broken and particles are produced. We have assumed that during this period two kind of superheavy particles are gravitationally produced: $X$-particles, which are nonconformally coupled with gravity and whose decay products form the baryonic matter, and $Y$-particles, which are conformally coupled with gravity but are only GIMP, and thus, they are responsible for the dark matter abundance. For this model we have shown that, for reasonable masses of the $X$-particles around $10^{15}$ GeV, a viable model with a reheating temperature in the GeV regime is obtained when the mass of the dark matter particles is of the order of $5\times 10^{16}$ GeV. Finally, we have shown numerically that the model leads, at late times, to a frozen quintessence, and thus, to an eternal inflation.

Acknowledgments. We want to thank Prof. Salvatore Capozziello for telling us, during the workshop Modified Gravity and Cosmology, the possibility to consider the production of superheavy particles nonconformally coupled with gravity in order to reduce the masses of the particles involved in the theory. This investigation has been supported by MINECO (Spain) grant MTM2017-84214-C2-1-P, and in part by the Catalan Government 2017-SGR-247.

Appendix: The number of e-folds in quintessential inflation

In this Appendix we will perform an accurate calculation of the number of e-folds for our model, i.e., for a quintessential inflation model where reheating is produced after the end of kination, and we will see that, due to the kination era, the number of e-folds is greater than in standard inflation, where reheating is produced due the oscillations of the inflaton field (see for a detailed calculation Cook ).

Let $k_{*}=a_{*}H_{*}$ be the value of the pivot scale in co-moving coordinates when it exits the Hubble radius and $N$ the number of e-folds from the exiting of the pivot scale to the end of inflation, i.e., $a_{END}=a_{*}e^{N}$, where, once again, $a_{END}$ denotes the value at the end of inflation.

Now we write

[TABLE]

where, as in previous sections, $a_{kin}$, $a_{end}$, $a_{rh}$, $a_{eq}$ and $a_{0}$ denote the value of the scale factor at the beginning of kination, at the end of kination, at the reheating time, at the matter-radiation equality and at present time, respectively.

Choosing, as usual, $k_{phys}\equiv\frac{k_{*}}{a_{0}}=0.02\mbox{Mpc}^{-1}$ and taking into account that $H_{0}\cong 2\times 10^{-4}\mbox{Mpc}^{-1}$ one gets

[TABLE]

where we have used that

[TABLE]

Now, using that $\frac{\rho_{X,eq}}{\rho_{X,rh}}=\frac{g_{eq}T^{4}_{eq}}{g_{*}T^{4}_{rh}}$, where, as we have already seen in the previous Section, the number of degrees of freedom at the matter-radiation equality is $g_{eq}=3.36$, and taking into account that after reheating the evolution is adiabatic, i.e., $a_{eq}T_{eq}=a_{0}T_{0}$, one gets

[TABLE]

and, from equations (15) and (33) and using that $\rho_{X,rh}=\frac{\pi^{2}}{30}g_{*}T_{rh}^{4}$, we obtain

[TABLE]

At this point, we use the observational data $T_{0}\cong 2.33\times 10^{-13}$ GeV and $H_{0}\cong 5.95\times 10^{-61}M_{pl}$ to get

[TABLE]

and finally, from equation (34) and the value $H_{*}\cong\sqrt{0.3}\pi(1-n_{s})M_{pl}\times 10^{-4}$, we conclude that

[TABLE]

Thus, since for many solvable models haro1 ; haro2 one has $\ln\left(\frac{a_{END}}{a_{kin}}\right)\cong-1$, we obtain

[TABLE]

In particular, for $m_{X}=10^{15}$ GeV and $n_{s}=0.968$, one gets

[TABLE]

which for the allowed temperatures $1\mbox{ MeV}\leq T_{rh}\leq 66\mbox{ TeV}$ (see formula (44)), leads to

[TABLE]