A review of Quintessential Inflation

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.

Understanding the universe's evolution is one of the greatest mysteries in the history of humanity. It is always the primary question: "where we come from and where we are going". In particular, its early and late expansions have been studied a great deal at present time. Looking at the scientific literature, one can find two popular and well accepted theories -though they are not observationally proved-, namely the inflation (the early surprisingly fast accelerated expansion of our universe) and the quintessence as a form of dark energy (the current cosmic acceleration). The inflationary paradigm [1][2][3] is actually a very fast accelerating phase of the early universe that lasted for an extremely tiny time and became able to solve a number of shortcomings associated with the standard Big Bang cosmology, such as the horizon problem, flatness problem or the primordial monopole problem. The predictive power of inflation was soon recognized due to its ability to explain the origin of inhomogeneities in the universe as quantum fluctuations during this epoch [4][5][6][7][8], because such an explanation greatly matches with the recent observational data from Planck's team [9]. Thus, it is remarkable to note that inflation, which appeared at the beginning of the 80's, is still considered the best way to explain the recent observational data, because it is nowadays the simplest viable theory that describes almost correctly the early universe in agreement with the recent observations [9].
On the other hand, one of the most accepted explanations for the current cosmic acceleration comes through the introduction of some quintessence field [10]. In fact, soon after the discovery of the current cosmic acceleration at the end of the last century [11,12], a class of pioneering cosmological models attempting to unify the early-and late-time accelerating expansions were introduced. By construction, unlike the standard quintessence models (see [16] for a review), these models -named as Quintessential Inflation (QI) models [13][14][15]-only contain one classical scalar field, also named inflaton as in standard inflation [1][2][3]5], and it is shown that they succeed in reproducing the two accelerated epochs of the universe expansion (see also [17][18][19][20][21][22][23][24][25][26][27][28] for other interesting QI models). This idea to unify Inflation with Quintessence was indeed a novel attempt by Peebles and Vilenkin in their seminal paper [14], and the novelty of their proposal comes through the introduction of a single potential that at early times allows inflation while at late times provides quintessence. Thus, a unified picture of the universe was effectively proposed connecting the distant early phase to the present one. Thanks to this proposal, the origin of the scalar responsible for the current inflation of the universe can be determined and fine-tunning problems are reduced [29]. In addition, the majority of models we deal with only depend on two parameters, which are determined by observational data. Hence, because of the behavior of the slow-roll regime as an attractor, the dynamics of the model are obtained with the value of the scalar field and its derivative -initial conditions-at some moment during inflation. This shows the simplicity of Quintessential Inflation, which from our viewpoint is simpler than standard quintessence, where a minimum of two fields are needed to depict the evolution of the universe, the inflaton and a quintessence field. Thus, one needs two different potentials and two different initial conditions: one for the inflaton, which has to be fixed during inflation, and another one for the quintessence field, whose initial conditions normally have to be fixed at the beginning of the radiation era.
This enhanced more investigation in order to connect Quintessential Inflation with the observational data [18-23, 25, 27, 28, 30-32] and, consequently, this particular topic has become a popular area of research at the present time. However, an important difference occurs with respect to the standard inflationary paradigm, where the potential of the inflaton field has a local minimum (a deep well) and, thus, the inflaton field releases its energy while it oscillates, producing enough particles [33][34][35][36][37] to reheat our universe. 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, the mechanism of reheating is completely different: once the inflationary phase is completed, a reheating mechanism keeping "alive" the inflaton field is needed to match inflation with the Hot Big Bang universe [1] because the particles existing before the beginning of this period were completely diluted at the end of inflation resulting in a very cold universe.
On this way, the most accepted idea to reheat the universe in the context of QI comes through a phase transition of the universe from inflation to kination (a regime where all the energy density of the inflation turns into kinetic [38]) where the adiabatic regime is broken, which allows to create particles. The mechanism to produce particles is not unique in this context since a number of distinct ones are available and can be used. The first one is the well-known Gravitational Particle Production studied long time ago in [39][40][41][42][43][44], at the end of the 90's in [45,46] and more recently applied to Quintessential Inflation in [13,14,47,48] for massless particles. A second important mechanism is the so-called Instant Preheating introduced in [49] and applied for the first time to inflation in [50] and recently in [48,51] in the context of α-attractors in supergravity. Other less popular mechanisms are the Curvaton Reheating applied to Quintessence Inflation in [52,53], the production of massive particles self-interacting and coupled to gravity [54] and the reheating via production of heavy massive particles conformally coupled to gravity [32,[55][56][57][58].
The production of superheavy massive particles conformally coupled to gravity is one of the most important concerns of this review. These particles are created during the phase transition from the end of inflation to the beginning of kination and should decay into lighter ones to form a thermal relativistic plasma, which eventually dominates to match this early period with the Hot Big Bang universe. In fact, the main motivation for using a conformally-coupled scalar field is its simplicity, which allows to employ the well-known Hamiltonian diagonalization method (see [59] for a review) to calculate the energy density of the gravitationally produced particles, showing that before the beginning of kination the vacuum polarization effects, which are geometric objects associated to the creation and annihilation of the so-called quasi-particles [59], are sub-dominant and have no relevant effect in the Friedmann equation. On the contrary, after the abrupt phase transition to kination, heavy massive particles are produced and, since their energy density decreases as a −3 before decaying in lighter particles and as a −4 after that, they will eventually dominate the energy density of the inflation whose decrease is as a −6 . And, thus, the universe will become reheated.
Coming back to the original Peebles-Vilenkin model [14], the inflationary part is described by a quartic potential and, according to the recent observations, this does not suit well. To be explicit, for the quartic potential in the inflationary part of this potential the two-dimensional contour of (n s , r), where n s is the scalar spectral index and r is the ratio of tensor to scalar perturbations, does not enter into the 95% confidence-level of the Planck results [9]. However, a simple change in the inflationary piece −quartic potential to a plateau one− can solve this issue (see [28] for a detailed discussion and also see [32]). On the other hand, the reheating mechanism followed in [14] is the gravitational production of massless particles that results in a reheating temperature of the order of 1 TeV. This reheating temperature is not sufficient to solve the overproduction of the Gravitational Waves (GWs). As a result, the Big Bang Nucleosynthesis (BBN) process could be hampered.
Dealing once again 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 [60] and on Gravitational Particle Production [61]. 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 [62]. As a consequence, the reheating temperature has to be greater than 1 MeV (see also [63] 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 bound 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 [64], but this problem can be successfully removed if the reheating temperature is of the order of 10 9 GeV (see for instance [65]). This is the reason why we will restrict the reheating temperature to remain, more or less, between 1 MeV and 10 9 GeV.
Finally, one has to take into account that a viable reheating mechanism has to deal with the effect of the Gravitational Waves (GWs), which are also produced during the phase transition, in the BBN success by satisfying the observational bounds coming from the overproduction of the GWs [14] or related to the logarithmic spectrum of their energy density [66]. As we will see throughout this review, the overproduction of GWs does not disturb the BBN success in the models studied (Lorentzian Quintessential Inflation and α-attractors in the context of QI), when the reheating is due to the production of superheavy particles, which shows their viability.
In summary, due to the lack of works that explain all the outstanding contributions made in Quintessential Inflation since the appearance of the seminal Peebles-Vilenkin article until now, the main goal of this review is to show the reader the most relevant aspects obtained on this topic during this time period. To do it we have chosen the most prominent papers on the subject from our viewpoint, and we have used the results provided by them to write this review.
The review is organized as follows: In Section 1 we review, with great detail, the seminal Peebles-Vilenkin model, calculating the value of the two parameters on which the model depends, the evolution of the inflaton field throughout history, the number of e-folds as a function of the reheating temperature, and additionally introducing models that improve the original one. Section 2 is devoted to the study of a class of Exponential Quintessential Inflation models. We show that, with initial conditions when the pivot scale leaves the Hubble radius (recall that this happens in the slow-roll regime which is an attractor, so one can safely take initial conditions close to the slow-roll solution) during the radiation phase, the solution is never in the basin of attraction of the scaling solution. Therefore, in order to match with the current observational data the potential should be modified introducing a pure exponential term. In Section 3 we review the so-called Lorentzian Quintessential Inflation (LQI), which is based on the assumption that the main slow-roll parameter as a function of the number of e-folds evolves as the Lorentzian distribution. The α-attractors in the context of Quintessential Inflation are revisited in Section 4, showing that we obtain the same results as in the case of LQI. In Section 5 we review other Quintessential Inflation models and in particular some aspects of the work of K. Dimopoulos in QI. Next, in Section 6 we study some reheating mechanisms in Quintessential Inflation, namely via Gravitational Particle Production of massless and superheavy particles, via Instant Preheating and via Curvature Reheating, obtaining bounds for the reheating temperature. In Section 7 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, showing that for the LQI model, when the reheating is via gravitational production of superheavy particles, all these bounds are easily overpassed. The last section is devoted to the conclusions of the review.
The units used throughout the paper are = c = 1, and the reduced Planck's mass has been denoted by M pl ≡ 1 √ 8πG ∼ = 2.44 × 10 18 GeV.

THE PEEBLES-VILENKIN MODEL
The first Quintessential Inflation scenario was proposed by Peebles and Vilenkin in their seminal paper [14] at the end of the last century, soon after the discovery of the current cosmic acceleration [11,12]. The corresponding potential is given by where λ is a dimensionless parameter and M M pl is a very small mass compared with the Planck's one. At this point, note that the quartic potential is the responsible for inflation and the inverse power law leads to dark energy (in that case quintessence) at late times.
One can see that the model is very simple and only depends on two parameters which can be calculated as follows: First of all, we calculate the main slow-roll parameters So, the spectral index and the ratio of tensor to scalar perturbations are given by [67] where the star denotes that the quantities are evaluated when the pivot scale leaves the Hubble radius (at the horizon crossing). Thus, one has the important relation 16(1 − n s ) = 3r between both quantities.
On the other hand, using the formula of the power spectrum of scalar perturbations (see for instance [67]) for the Peebles-Vilenkin model it yields that where we have used the formulas (2) and (3). And, taking for example n s ∼ = 0.96, which is approximately its observational value [68,69], we get The calculation of the other parameter is more involved and for that one needs to know the evolution of the inflaton field up to the present time. As we will see, the value of the inflaton at the present time is of the order of 40M pl (see Fig. 2). So, taking into account that in order to reproduce the current cosmic acceleration the kinetic energy has to be negligible compared with the potential one at the present time, in order to match with the present energy density we will have where H 0 ∼ 10 −61 M pl is the present value of the Hubble rate. Finally, we get Unfortunately, the original Peebles-Vilenkin model provides a spectral index and a tensor/scalar ratio that do not enter in the two dimensional marginalized joint confidence contour at 2σ CL for the Planck TT, TE, EE + low E+ lensing + BK14 + BAO likelihood (see Fig. 1). For this reason, we need to improve the part of the potential which provides inflation, that is, the quartic potential.

Improvements
For that purpose, one would need to consider plateau potentials [27] or α-attractors [48,70], such as an Exponential SUSY Inflation type potential or a Higgs Inflation-type potential For the first model, the slow-roll parameters are given by and, taking into account that when the pivot scale leaves the Hubble radius e αϕ * /M pl 1, one gets * ∼ = and thus, the spectral index and the tensor/scalar ratio are equal to obtaining the relation (1 − n s ) 2 = α 2 2 r.
We can also calculate the number of e-folds from the horizon crossing to the end of inflation, namely N = t EN D t *

Hdt.
In the slow-roll regime, where it can be approximated by which for the first potential leads to In the same way, for the second potential, the slow-roll parameters are given by and, taking into account again that when the pivot scale leaves the Hubble radius e αϕ * /M pl 1, we get * ∼ = 2α 2 e 2αϕ * /M pl , and thus, the spectral index and the tensor/scalar ratio are equal to obtaining the same relation as for the other potential, namely (1 − n s ) 2 = α 2 2 r.
In addition, for the number of e-folds, one also obtains the same result and what is important to note is that for both potentials the spectral index and the ratio of tensor to scalar perturbations have the same relation with the number of e-folds which implies that for α ∼ O(1) and for a number of e-folds greater than 60, which is typical in quintessential inflation due to the kination phase, the ratio of tensor to scalar perturbations is less than 0.003. Thus, the spectral index and the tensor/scalar ratio enter perfectly in the two dimensional marginalized joint confidence contour at 2σ CL for the Planck TT, TE, EE + low E+ lensing + BK14 + BAO likelihood (see Fig. 1). Now we want to understand the evolution of the Quintessential Inflation models after inflation, which only depends on the tail of the potential. We will focus on the Peebles-Vilenkin one, and for that case, inflation ends when = 1, which occurs for ϕ EN D = 2 √ 2M pl . Taking into account that an alternative expression of this slow-roll parameter is = −Ḣ H 2 , we conclude that at the end of inflation the effective Equation of State (EoS) parameter, namely w ef f , satisfies where we have used the Friedmann and Raychaudhuri equations being ρ the energy density and P the corresponding pressure.
Then, from (22) one gets the following relation at the end of inflation,φ 2 EN D = V (ϕ EN D ), meaning that the energy density at the end of inflation is Next, we have to assume, as usual, that there is no drop of energy between the end of inflation and the beginning of kination, which for these models occurs when ϕ kin ∼ = 0. So, at the beginning of kination the energy density is approximately given by ρ kin ∼ = 96λM 4 pl =⇒ H kin ∼ = 4 √ 2λM pl and, since all the energy is practically kinetic, the effective EoS parameter is very close to w ef f ∼ = 1. For this value, combining the Friedmann and Raychaudhuri equations (23), one getsḢ whose solution is given by H(t) = 1 3t .
Then, coming back to the Friedmann equation and taking into account that during kination all the energy density is kinetic, we get which in terms of the Hubble rate could be written as follows: At this point, we have to take into account that during the phase transition from the end of inflation to the beginning of kination the adiabatic regime is broken and particles are created. Here we will assume that these particles, which we consider very massive and conformally coupled to gravity, are gravitationally produced. So, they must decay into lighter ones in order to form a relativistic plasma which will eventually dominate the energy density of the universe and will become reheated. And then, since the kination regime ends when the energy density of the inflaton field is of the same order than the one of the produced particles, two different scenarios appear: 1. Decay before the end of the kination period.
2. Decay after the end of the kination period.

Decay before the end of kination
In this case, since the thermalization of the decay products is nearly instantaneous and occurs before the end of kination, the reheating time coincides with the end of the kination. Hence, we will havė Taking into account that at the reheating time the energy density of the inflaton field is the same as the one of the radiation, the Friedmann equation will become H 2 reh = 2ρ reh 3M 2 pl and, using the Stefan-Boltzmann law ρ reh = π 2 30 g reh T 4 reh , where g reh denotes the effective number of degrees of freedom at the reheating time, we get As we will see, by calculating the value of Θ for several reheating mechanisms, we can safely do the approximation which only depends on the value of the parameter Θ. Effectively, using for example Instant Preheating (see subsection 7.2), one gets Θ = g 2 2π 2 , which for the viable value of the coupling g ∼ 10 −4 yields Θ ∼ 10 −9 , and one can easily check that the approximation holds. After the matter-radiation equality the dynamical equations cannot 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 ≡ − ln(1 + z) = ln a a0 . Now, using the variable N , one can recast the energy density of radiation and matter respectively as and where ρ eq = ρ m,eq + ρ r,eq and the value of the energy density of the matter or radiation ρ m,eq = ρ r,eq at the matterradiation equality can be calculated as follows: Using the present value of the ratio of the matter energy density to the critical one Ω m,0 = 0.308 and the present value of the Hubble rate H 0 ∼ = 68 Km/sec/Mpc ∼ 10 −33 eV, one can deduce that the present value of the matter energy density is ρ m,0 = 3H 2 0 M 2 pl Ω m,0 ∼ = 3 × 10 −121 M 4 pl , and at the matter-radiation equality ρ m,eq = ρ m,0 (1 + z eq ) 3 .
In order to obtain the dynamical system for this scalar field model, we introduce the following dimensionless variables, where H 0 ∼ 1.42 × 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φ + 3Hφ + V ϕ = 0, one can construct the following non-autonomous dynamical system, where the prime represents the derivative with respect to N ,H = H H0 andV = V . Moreover, the Friedmann equation now looks asH where we have introduced the dimensionless energy densitiesρ r = ρr Then, we have to integrate the dynamical system, starting at N eq = − ln(1 + z eq ), with initial conditions 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 (50) to 1, i.e., imposingH(0) = 1.
Numerical calculations show that M ∼ 10 5 GeV for a reheating temperature of 10 5 GeV, as in our heuristic reasoning. In addition, it remains of the same order for the other values of the reheating temperature within the allowed range. In Figure 2 we display the evolution of the scalar field ϕ and in Figure 3 we plot the graph of the parameter Ω i ≡ ρi 3H 2 M 2 pl for matter, radiation and dark energy, and also the evolution of the effective EoS parameter.

The number of e-folds
In this subsection we aim to find the relation between the number of e-folds from the horizon crossing to the end of inflation and the reheating temperature T reh . We start with the well-known formula [71] where, as previously, a is the scale factor and a EN D , a kin , a reh , a eq and a 0 denote respectively its value at the end of inflation, at the beginning of kination, radiation and matter domination, and finally at the present time.
Taking into account the evolution during kination and radiation we have a kin a reh 6 = ρ reh ρ kin and a reh a eq 4 = ρ eq ρ reh (52) and, noting that H 0 ∼ 2 × 10 −4 Mpc −1 ∼ 6 × 10 −61 M pl and k * = a 0 k phys , where we have chosen as usual k phys = 0.02Mpc −1 , we have obtained where we have used once again that after the matter-radiation equality the evolution is adiabatic, that is, a 0 T 0 = a eq T eq , as well as the Stefan-Boltzmann law ρ eq = π 2 30 g eq T 4 eq and ρ reh = π 2 30 g reh T 4 reh being g eq = 3.36 the degrees of freedom at the matter-radiation equality. And we have chosen as degrees of freedom at the reheating time the ones of the Standard Model, i.e., g reh = 106.75. Now, from the formula of the power spectrum we infer that H * ∼ 4π × 10 −4 * 10 M pl ∼ 4 × 10 −4 √ * M pl , obtaining and introducing the current value of the temperature of the universe T 0 ∼ 10 −31 M pl we get Taking into account that for the Peebles-Vilenkin model one has H kin ∼ 4 √ 2λM pl ∼ 4 √ 2 × 10 −7 M pl and * ∼ = (1 − n s )/3, we get and, using the value of the spectral index n s ∼ = 0.96 and assuming that ln a EN D a kin is negligible, we finally obtain Since the scale of nucleosynthesis is 1 MeV and in order to avoid the late-time decay of gravitational relic products such as moduli fields or gravitinos which could jeopardise the nucleosynthesis success, one needs temperatures in the range 1 MeV ≤ T reh ≤ 10 9 GeV, which leads to constrain the number of e-folds to 64 N 74.
A final remark is in order: Note that in Quintessential Inflation we have obtained a number of e-folds greater than in standard inflation, which normally ranges between 50 and 60 e-folds. This fact is due to the kination regime appearing in QI, which increases significantly the number of e-folds.

EXPONENTIAL QUINTESSENTIAL INFLATION
In this section we will consider more complex models as the Peebles-Vilenkin one: the same Exponential Inflationtype potentials studied for the first time in [24], where γ is a dimensionless parameter and n is an integer.
In this case the power spectrum of scalar perturbations, its spectral index and the ratio of tensor to scalar perturbations are given by (see for a detailed derivation [27]) and An important relation is obtained combining the equations (60) and (61), which leads to the formula for the power spectrum and thus, which, for the viable values of n s ∼ = 0.96 (the central value of the spectral index) and r = 0.02 (whose value, as we will immediately see, guarantees that the model provides a viable a number of e-folds), leads to Next, it is important to realize that a way to find theoretically the possible values of the parameter γ is to combine the equations (61) and (62) to get and, using the theoretical values r ≤ 0.1 and n s = 0.9649 ± 0.0042 (see for instance [68,69]), one can find the candidates of γ at 2σ C.L. These values have to be checked for the joint contour in the plane (n s , r) at 2σ C.L., when the number of e-folds is approximately between 60 and 75, which is what usually happens in Quintessential Inflation due to the kination phase [20,71] -the energy density of the scalar field is only kinetic [13,38]-after the end of the inflationary period.
For this kind of potentials, in order to compute the number of e-folds we need to calculate the main slow-roll parameter whose value at the end of inflation is EN D = 1, meaning that at the end of this epoch the field reaches the value Then, the number of e-folds is given by Thus, combining the equations (60), (61) and (68), one obtains the spectral index and the tensor/scalar ratio as a function of the number of e-folds and the parameter γ.
On the other hand, since inflation ends at ϕ EN D = parameter is equal to −1/3, meaning thatφ 2 EN D = V (ϕ EN D ), the energy density at the end of inflation is (69) and the corresponding value of the Hubble rate is given by which will constrain very much the values of the parameter γ because in all viable inflationary models at the end of the inflation the value of the Hubble rate is of the order of 10 −6 M pl [2]. In fact, when (70) is of the order of 10 −6 M pl one gets Therefore, to perform numerical calculations, we will use the values of n = 10 and r = 0.02 and, thus, for γ ∼ 10 −16 we obtain approximately 67 e-folds, which is a viable value in Quintessential Inflation.
After this, we have to find out the time when kination starts, which can be chosen at the moment when the effective EoS parameter is very close to 1. Assuming for instance that kination starts at w ef f ∼ = 0.99, we have numerically obtained ϕ kin ∼ = 47M pl andφ kin ∼ = 6 × 10 −9 M 2 pl , hence H kin ∼ = 2 × 10 −9 M pl .
Next, we want to calculate the value of the scalar field and its derivative at the reheating time (for example, in the case that the decay was before the end of kination). Using for our model the formulas (29) and assuming that the effective number of degrees of freedom is the same as in the Standard Model, i.e. g reh = 106.75, we get for a viable reheating temperature T reh ∼ 10 7 GeV that at the beginning of the radiation the initial conditions of the scalar field are

The scaling and tracker solutions
The scaling and tracker solutions appear in pure exponential potentials and approximately in more generic exponential potentials as (58) (see for instance [27,67]).

The scaling solution
Assuming as usual the flat Friedmann-Lemaître-Robertson-Walker (FLRW) geometry of our universe, we consider a pure exponential potential, V (ϕ) = V 0 e −γϕ/M pl , and a barotropic fluid with Equation of State (EoS) parameter w = 1/3, i.e., a radiation fluid, whose energy density we continue denoting by ρ r . Then, the dynamical system is given by the equations where The scaling solution is a solution with the property that the energy density of the scalar field scales as the one of radiation, meaning that the Hubble parameter evolves as H(t) = 1 2t . At this point, we look for a solution of the form where C andt are constants. The equationρ r + 4Hρ r = 0 is satisfied for any value of C 2 and the equation On the other hand, equating H(t) = 1 2t with H = 1 √ 3M pl √ ρ ϕ + ρ r , which yields the relation shows that the scaling solution only exists for γ > 2, and the corresponding energy densities evolve as follows, and thus, one could derive the corresponding density parameters, An alternative approach to get this solution goes as follows (see for instance [72]): one can introduce the dimensionless variablesx which enable us to write down the following autonomous dynamical system, where w is the Equation of State (EoS) parameter, which in our case is equal to 1/3, together with the constraint One can see that for w = 1/3 the dynamical system (80) has the attractor solutionx = γ , whose energy density scales as radiation, and Ω ϕ = 4 γ 2 . So, since during radiation H = 1 2t , thus, fromx = 2 3 2 γ , we havė Finally, in order to obtaint, we use the equationỹ = 2

The tracker solution
A tracker solution [73,74] is an attractor solution of the field equation which describes the dark energy domination at late times. In fact, for an exponential potential, it is the solution ofφ + 3Hφ + V ϕ = 0 when the universe is only filled with the scalar field. Once again, we look for a solution of the type Inserting this equation intoφ + 3Hφ + V ϕ = 0, one gets which means that 0 < γ < √ 6 and, thus, the tracker solution is given by Moreover, for this equation, it is not difficult to show that meaning that the corresponding effective EoS parameter is given by which proves that, in order to impose the late-time acceleration of the universe, we need to restrict γ as 0 < γ < √ 2.
Finally, in terms of the time N = − ln(1 + z), taking into account that the relation between N and the cosmic time is given by t = 2 γ 2 H0 e γ 2 H/2 , the tracker solution has the following expression (see Section V of [75] for a detailed derivation),

Numerical simulation during radiation
In this subsection we will show numerically that in Quintessential Inflation, for exponential types of potentials, during radiation the scalar field is never in the basin of attraction of the scaling solution and, thus, does not evolve as the scaling solution.
To prove it, first of all we heuristically calculate the value of the red-shift at the beginning of the radiation epoch a eq a rh = (1 + z eq ) ρ r,reh ρ r,eq where we have used that ρ r,eq = ρ r,reh a reh aeq 4 and the adiabatic evolution after the matter-radiation equality. Then, taking z eq ∼ 3 × 10 3 and using the current temperature T 0 ∼ 10 −31 M pl , for the reheating temperature T reh ∼ 10 9 GeV we get z reh ∼ 10 22 .
In addition, at the beginning of radiation the energy density of the matter will be ρ m,reh = ρ m,eq a eq a reh 3 = ρ m,eq ρ r,reh ρ r,eq and, for radiation, In this way, the dynamical equations after the beginning of radiation can be easily obtained using once again as a time variable N ≡ − ln(1 + z) = ln a a0 . Recasting the energy density of radiation and matter respectively as functions of N , we get and where N reh = − ln(1 + z reh ) ∼ −50 denotes the value of the time N at the beginning of radiation for a reheating temperature around 10 9 GeV.
To obtain the dynamical system for this scalar field model, we will introduce the dimensionless variables where K is a parameter that we will choose accurately in order to facilitate the numerical calculations. So, taking into account the conservation equationφ + 3Hφ + V ϕ = 0, one arrives at the following dynamical system, where the prime is the derivative with respect to N , . It is not difficult to see that one can writeH where we have defined the dimensionless energy densities asρ r = ρr Next, to integrate from the beginning of radiation up to the matter-radiation equality, i.e., from N reh ∼ −50 to N eq ∼ −8, we will choose KM pl = 10 −17 GeV 2 , yieldinḡ andV Finally, with the initial conditions for the field being x reh ∼ 70 and y rh ∼ 3 × 10 32 (these initial conditions are obtained in the equations (73)) and integrating numerically the dynamical system, we conclude that, for values of the parameter n greater than 10, the value of the EoS parameter w ϕ remains 1 during radiation domination, namely between N rh and N eq , thus proving that the inflaton field does not belong to the basin of attraction of the scaling solution for large values of n. In addition, integrating the system up to now, we can see that the model never matches with the current observational data.

A viable model
To depict the late-time acceleration, one has to modify the original potential because it does not explain the current observational data (for the potential (58), at the present time, the density parameter Ω ϕ = ρϕ 3H 2 M 2 pl is far from its observational value, namely 0.7). For this reason and following the spirit of the Peebles-Vilenkin model [14], to match with the current observational data it is mandatory to introduce a new parameter M with units of mass which must be calculated numerically. In that case, we will consider the following modification of the potential (58) by introducing a new exponential term containing the parameter M , with 0 < γ < √ 2 in order to obtain, as we have already seen, that at late times the solution is in the basin of attraction of the tracker solution, which evolves as a fluid with effective EoS parameter w ef f = γ 2 3 − 1.
Remark.-An important remark is in order. Contrary to the potential used in [24,27] where ρ ν =ρ ν e βϕ/M pl is the neutrino energy density with constant bare mass m ν and β > 0 is the non-minimal coupling of neutrinos with the inflaton field (see for details [27]), which has a minimum where the inflaton ends its evolution, thus acting as an effective cosmological constant which stands for the current cosmic acceleration, the potential considered here does not have a minimum and the scalar field continues rolling down the potential following the spirit of the quintessential inflation.
On the other hand, in this case one can also justify this choice considering that the inflaton field is non-minimally coupled with neutrinos but with a negative coupling constant, namely β = −γ < 0. Therefore, M 4 =ρ ν − 3p ν and the effective mass of neutrinos is given by m ν,ef f (ϕ) = m ν e −γϕ/M pl , which tends to zero for large values of the field, so the neutrinos become relativistic, contrary to what happens when β is positive, where the neutrinos acquire a heavy mass becoming non-relativistic particles. Finally, to match with the current observational data, M 4 =ρ ν − 3p ν must be very small compared with the Planck's energy density. In fact, for n = 10 and γ = 0.1 we have numerically obtained M ∼ 10 −2 eV. Now, one has to solve the dynamical system (96) with initial conditions at the beginning of the radiation era. Choosing now K = H 0 (the present value of the Hubble rate) we have to take into account that to obtain the parameter M one has to impose the value ofH = H(N )/H 0 to be 1 at N = 0.
The numerical results are presented in Figure 5, where one can see that at very late times the energy density of the scalar field dominates and the universe has an eternal accelerated expansion because its evolution is the same as that of the tracker solution. State parameter w ef f , from kination to future times, for n = 10 and γ = 0.1. As one can see in the picture, after kination the universe enters in a large period of time where radiation dominate. Then, after the matter-radiation equality, the universe becomes matter-dominated and, finally, near the present, it enters in a new accelerated phase where w ef f approaches 0.1 2 3 − 1 ∼ = −1 , that is, it has the same effective EoS parameter as the tracker solution, meaning that the solution is in the basin of attraction of the tracker one.

LORENTZIAN QUINTESSENTIAL INFLATION
Based on the Cauchy distribution (Lorentzian in the physics language) the authors of [76,77] considered the following ansatz, where is the main slow-roll parameter, N = ln a a0 denotes the total number of e-folds (do not confuse with N , the number of e-folds from the horizon crossing to the end of inflation), ξ is the amplitude of the Lorentzian distribution and Γ is its width. From this ansatz, we can find the exact corresponding potential of the scalar field, namely where λ is a dimensionless parameter and the parameter γ is defined by This potential can be derived by using equations (37) in [78]. However, in this work we are going to use a more simplified potential, keeping the same properties as the original potential but not coming exactly from the ansatz (102). However, from this potential we do recover the ansatz by using the suitable approximations, which are valid if we properly choose the parameters such that γ 1 and 2ξ Γπ > 1. For this reason and also because the results of the original potential happen to be exactly the same as shown in [83], it is better to work with the simplified potential, namely We can see the shape of the potential on Fig. 6, where the inflationary epoch takes place on the left-hand side of the graph, while the dark energy epoch occurs on the right-hand side.

Calculation of the value of the parameters involved in the model
We start with the main slow-roll parameter, which is given by and, since inflation ends when EN D = 1, one has to assume that 2ξ Γπ > 1 to guarantee the end of this period.
In fact, we have and we can see that, for large values of γ, one has that ϕ EN D is close to zero. Thus, we will choose γ 1 =⇒ Γξ 1, which is completely compatible with the condition 2ξ Γπ > 1.
On the other hand, the other important slow-roll parameter is given by Both slow-roll parameters have to be evaluated when the pivot scale leaves the Hubble radius, which will happen for large values of cosh (γϕ/M pl ), obtaining * = with ϕ * < 0. Then, since the spectral index is given in the first approximation by n s ∼ = 1 − 6 * + 2η * , one gets after some algebra where r = 16 * is the ratio of tensor to scalar perturbations. Now, we calculate the number of e-folds from the horizon crossing to the end of inflation, which is given by so we have that meaning that our model leads to the same spectral index and tensor/scalar ratio as the α-attractors models with α = 2 3γ 2 (see for instance [79]). Finally, it is well-known that the power spectrum of scalar perturbations is given by Now, since in our case V (ϕ * ) ∼ = λM 4 pl e ξ , meaning that H 2 * ∼ = λM 2 pl 3 e ξ , and taking into account that * ∼ = (1−ns) 2 8γ 2 , one gets the constraint where we have chosen as a value of n s its central value 0.9649.
Summing up, we will choose our parameters satisfying the condition (129), with γ 1 and ξ 1, which will always fulfill the constraints Γξ 1 and 2ξ Γπ that we have imposed. Then, to find the values of the parameters one can perform the following heuristic argument: Taking for example γ = 10 2 , the constraint (129) becomes λe ξ ∼ 7 × 10 −15 . On the other hand, at the present time we will have γϕ 0 /M pl 1 where ϕ 0 denotes the current value of the field. Thus, we will have V (ϕ 0 ) ∼ λM 4 pl e −ξ , which is the dark energy at the present time, meaning that So, for the value H 0 = 67.81 Km/sec/Mpc = 5.94 × 10 −61 M pl , we get the equations whose solution is given by ξ ∼ 122 and λ ∼ 10 −69 .
If we choose γ ∼ 10 2 , we see that the values of ξ and γ could be set equal in order to obtain the desired results from both the early and late inflation. From now on we will set ξ = γ, since it may help to find a successful combination of parameters because it reduces the number of effective parameters. As we will see later, numerical calculations show that, in order to have Ω ϕ,0 ∼ = 0.7 (observational data show that, at the present time, the ratio of the energy density of the scalar field to the critical one is approximately 0.7), one has to choose ξ = γ ∼ = 121.8.
Next, we aim to find the relation between the number of e-folds and the reheating temperature T reh . We start once again with the formula a kin a reh a reh a eq a eq a 0 (115) and, following the same steps as in the case of the Peebles-Vilenkin model, we get Now, from the formula of the power spectrum of scalar perturbations (127) we infer that H * ∼ 4 × 10 −4 √ * M pl , obtaining and introducing the current value of the temperature of the universe T 0 ∼ 9.6 × 10 −32 M pl we get So, since we have numerically checked that H kin ∼ 4 × 10 −8 M pl , we get where we have used that * = 1 2γ 2 N 2 and we have also numerically computed that ln a EN D a kin ∼ = −0.068.
Finally, taking into account that the scale of nucleosynthesis is 1 MeV and in order to avoid the late-time decay of gravitational relic products such as moduli fields or gravitinos which could jeopardise the nucleosynthesis success, one needs temperatures lower than 10 9 GeV. So, we will assume that 1 MeV ≤ T reh ≤ 10 9 GeV, which leads to constrain the number of e-folds to 58 N 67. And for this number of e-folds, 0.966 n s 0.970, which enters within its 2σ CL range.

Present and future evolution: Numerical analysis
To understand the present and future evolution of our universe, we have integrated numerically the dynamical system (96) with K = H 0 and initial conditions at the reheating time (obtained numerically) ϕ kin ∼ −0.003M pl anḋ ϕ kin ∼ 8 × 10 −8 M 2 pl . The results are presented in Figure 8, where one can see that in LQI the universe accelerates forever at late times with an effective EoS parameter equal to −1. kination to future times, for a reheating temperature around 10 9 GeV. Right: The effective Equation of State parameter w ef f , from kination to future times. As one can see in the picture, after kination the universe enters in a large period of time where radiation dominates. Then, after the matter-radiation equality, the universe becomes matter-dominated and, finally, near the present, it enters in a new accelerated phase where w ef f approaches −1.

α-ATTRACTORS IN QUINTESSENTIAL INFLATION
The concept of α-attractor, in the context of standard inflation, was introduced for the first time in [80], obtaining models which generalize the well-known Starobinsky model (see [81] for a detailed revision of the Starobinsky model).
In this section we consider α-attractors in the context of Quintessential Inflation. For this reason we deal with the following Lagrangian motivated by supergravity and corresponding to a non-trivial Kähler manifold (see for instance [47] and the references therein), combined with a standard exponential potential, where φ is a dimensionless scalar field, and κ and λ are positive dimensionless constants.
In order that the kinetic term has the canonical form, one can redefine the scalar field as follows, obtaining the following potential, where we have introduced the dimensionless parameter n = κ √ 6α. Similarly to [76,77,83], the potential satisfies the cosmological seesaw mechanism, where the left side of the potential gives a very large energy density -the inflationary side-and the right side gives a very small energy density -the dark energy side. The asymptotic values are V ± = λ exp(±n). The parameter n is the logarithm of the ratios between the energy densities, as ξ in the earlier versions [83].

Calculation of the value of the parameters involved in the model
Dealing with this potential at early times, the main slow-roll parameter is given by where we must assume that n 2 12α > 1 because inflation ends when EN D = 1, and the other slow-roll parameter is Both slow-roll parameters have to be evaluated at the horizon crossing, which will happen for large values of cosh with ϕ * < 0.
Next, we calculate the number of e-folds from the horizon crossing to the end of inflation, which for small values of α is given by so we get the standard form of the spectral index and the tensor/scalar ratio for an α-attractor [2], Finally, it is well-known that the power spectrum of scalar perturbations is given by and, since in our case V (ϕ * ) ∼ = λM 4 pl e n and, thus, pl 3 e n , taking into account that * ∼ = 3α 16 (1 − n s ) 2 one gets the constraint where we have chosen as the value of n s its central value given by the Planck's team, i.e., n s = 0.9649 [68].
2. If one takes α very small -for example of the order 10 −2 -, then, as has been shown in [48], a much simpler model than the exponential one is possible with a linear potential. In fact, the authors of [48] showed that the potential V (φ) = λ(φ + √ 6α)M 4 pl + ΛM 4 pl ,, which in terms of the canonically normalized field ϕ has the form is viable for values of α ∼ 10 −2 and Λ ∼ 10 −120 .
3. Another important important application of the α-attractors is its use to alleviate the current Hubble tension [82]. In that work, the authors, in the framework of α-attractors, include to the model an Early Dark Energy (EDE) component that adds energy to the Universe, because the success in easing the Hubble tension crucially depends on the shape of this energy injection. In fact, in [82] the authors use the following potential for EDE, obtaining for p = 2 and n = 4 the following value of the Hubble rate at the present time: H 0 = (70.9 ± 1. Once again, we have integrated numerically the dynamical system (96) with K = H 0 and initial conditions at the reheating time (obtained numerically) ϕ kin ∼ −0.5M pl andφ kin ∼ ×10 −6 M 2 pl . The results presented in Figure 9 are the same as the ones obtained in the LQI scenario. So, both models predict the same evolution of the universe.

OTHER QUINTESSENTIAL INFLATION MODELS
In the present section, some QI models, different from the ones studied in the previous sections, are presented and revisited.

Dimopoulos work in Quintessential Inflation
In this subsection we briefly review some points of the vast contribution of K. Dimopoulos in this field, starting with the class of models introduced in [84], whose potential is given by where M M pl and m are two masses and k is a positive integer.
It is interesting to note that the potential has the following asymptotic behavior, meaning that at very early times the potential approaches to a constant non-vanishing false vacuum responsible for inflation, and at late times it is approximately a quasi-exponential one which leads to quintessence.
The value of the parameter m is obtained, as in standard inflation, using the power spectrum of scalar perturbation, which for k ≤ 4 leads to the value m ∼ 10 15 GeV. And, in order to match with the current observational data, one has to choose which has only sub-Planckian values for k ≤ 4. In fact, for k = 4 one has M ∼ 10 18 GeV.
Another important contribution comes from Modular Quintessential Inflation [85], where the author introduced the following toy potential, where m and M are masses and q a positive integer. Once again, the asymptotic behavior is given by showing that the inflationary piece of the potential is like a top-hill one and the tail, which leads to quintessence, is a pure exponential one. For this modular model the minimum value of m is around 10 7 GeV and, as we have already seen, for pure exponential potentials, in order that the tracker solution, which is an attractor, reproduces an accelerated universe, one needs that q and M satisfy the constraint 0 < q In addition, it is not difficult to show that the spectral index is given by and thus, taking n s ∼ = 0.96, we can deduce that q M pl M 2 ≤ 2 × 10 −2 . Finally, using the power spectrum of scalar perturbation one gets the following important relation between the parameters involved in the model: Next, we review the treatment of α-attractor done in [47]. First of all, one has to introduce a negative Cosmological Constant (CC) (recall that Einstein's CC is positive) with the following form, Λ = −λM 2 pl e −n , and thus, adding to the Lagrangian (120) the term ΛM 2 pl leads to the following effective potential, whose main difference with the potential (122) is that it vanishes at very late times.
During inflation ϕ < 0 and the potential becomes as (122) because e −n 1. So, as we have shown for the potential (122), for small values of α inflation works also well when this CC is introduced in the model.
In the same way, for large values of the scalar field the potential will become with γ ≡ 2 3α . Since for an exponential potential a late-time eternal acceleration is achieved when γ < √ 2, one has to choose α > 1/3. However, since in this case one has to choose α > 1/3, the calculation of the spectral index and the ratio of tensor to scalar perturbations changes a little bit with respect to the case α 1, obtaining (see for details [? ]) Another important contribution of Dimopoulos is the introduction of the concept of Warm Quintessential Inflation in [86], which is applied to the original Peebles-Vilenkin model, obtaining the unification of early-and late-time acceleration, and providing a very natural mechanism to reheat our universe. In the same way, an interesting idea is to consider Quintessential Inflation in the context of Palatini R 2 -gravity [87], adapting the well-known Starobinsky model to QI.
About the reheating mechanism in Quintessential Inflation, in [54] the authors introduced a very interesting concept: an originally subdominant, non-minimally coupled and massive scalar field with quartic self-interactions, and without any interaction with the inflaton field. During inflation the field is at zero vacuum, but during kination the minimum of the potential displaces from zero because the field becomes tachyonic, and thus, it is "kicked-off" of the origin (which becomes a potential hill) due to quantum fluctuations. Finally, when it reaches the new minimum of potential it starts to oscillate coherently realising its energy and creating particles as in standard inflation. However, as it is highlighted in [88][89][90][91], there is an important limitation of the treatment considered in [54], namely the fact that reheating through this mechanism is far from being an homogeneous, coherent and perturbative process. This observation has far reaching consequences, completely absent in the original proposal, such as the production of topological defects, the generation of a detectable gravitational wave background or baryogenesis. Moreover, it allows the onset of radiation domination for arbitrary potentials, and not only for quartic interactions, as expected from the homogeneous treatment.
Another effective reheating can also be obtained from a combination of supergravity and string theory [92]. The model contains the inflaton and a Peccei-Quinn field, which will be the responsible for particle production. In fact, both fields are coupled as in the theory of Instant Preheating, which we will study in next section, and thus, when the adiabatic evolution is broken the total kinetic energy of the inflaton decays into radiation through resonant production of Peccei-Quinn particles, and the residual potential density of the inflaton field can be the responsible for dark energy.
Finally, the Curvaton rehating, which we deal with in next section, is applied to QI in [93], where bounds on the parameters of curvaton models are found. In addition, using a minimal curvaton model, the authors showed that the allowed parameter space is considerably larger than in the case of the usual oscillatory inflation models, i.e, than in standard inflation.

Quintessential Inflation with non-canonical scalar fields
In [24] (see also [23,27]) the authors study Quintessential Inflation using the following action, where S m and S r are the actions for matter and radiation and S ν is the action for neutrinos (see [94] and references therein for a discussion of the coupling of neutrinos and QI). Here, the Lagrangian density for massive neutrinos is given by [27] which depicts a non-minimal coupling between the inflaton and the neutrinos, the potential is a pure exponential one, and the coupling k 2 is given by being α,α, γ and β the parameters of the model.
For this model the slow-roll parameters as a function of the number of e-folds N are given by * =α and thus, the spectral index and the tensor/scalar ratio are At first glance the model seems more complicated than the ones we have studied previously because it contains four parameters, but it simplifies very much because at late times, i.e, for large values of ϕ, the coupling k 2 goes to 1, recovering the canonical form of the action. In addition, the effect of neutrinos, at late times, is the modification of the potential V (ϕ), becoming the following effective potential, where ϕ 0 and ρ ν,0 are the current values of the inflaton and the energy density of the massive neutrinos.
Therefore, since the potential V (ϕ) is an exponential one, we are in the same situation studied in the subsection 3.3, which, as we have already seen, leads to a viable model.

Gauss-Bonnet Quintessential Inflation
One of the problems of Quintessential Inflation is that due to the huge difference between the energy density scale of inflation and the current energy density (which is over a hundred orders of magnitude), in standard Quintessential Inflation the inflaton field typically rolls over super-Planckian distances in field space, resulting into a multitude of problems. Firstly, the flatness of the quintessential tail may be lifted by radiative corrections. Also, because the associated mass is so small, the quintessence field may give rise to a so-called 5th force problem, which can lead to a violation of the equivalence principle. To avoid these problems, it is desirable to keep the field variation sub-Planckian. In this case, however, to bridge the huge difference between the inflation and dark energy density scales, the quintessential tail must be steep. But, if the quintessential tail is too steep, when the field becomes important today, it unfreezes and rolls down the steep potential not leading to accelerated expansion at all. One way to overcome this problem is to make sure the scalar field remains frozen today even though the quintessential tail is steep. To this end, a solution (see [95] for a more detailed discussion) is the coupling of the field with the Gauss-Bonnet term, because such coupling impedes the variation of the field even if the potential is steep. Thus, a reason to use the GB scalar is that in some models the GB coupling could become important at late times making sure that the field freezes with sub-Planckian displacement, such that it becomes the dark energy today without the aforementioned problems. Fortunately, as we have already shown, the potential in Lorentzian Quintessential Inflation and the one of α-attractors is steep enough and the field also freezes (w ef f approaches to −1 at late times) at the present time leading to the current acceleration. Therefore, here we will analyse the paper [95] (see also [96] for another paper about Gauss-Bonnet QI) where the authors consider the following action, where F (ϕ) = F 0 e −qϕ/M pl (q > 0) is the coupling with the field and G = R 2 − 4R µν R µν + R ρµσν R ρµσν is the Gauss-Bonnet scalar, which for the flat FLRW geometry has the simple form G = 24H 2 (Ḣ + H 2 ).
On the other hand, the authors have chosen as a potential a mathematically convenient prototype for situations in which an early time plateau, favoured by Planck, as well as an exponential quintessential tail are present. For example, For this potential the number of e-folds is N ∼ = 1 4p 2 e 2pϕ/M pl , and the spectral index and the ratio of tensor to scalar perturbations as a function of N are given by which have the same form of the α-attractors for α = 1 6p 2 , so they clearly enter in the 2σ C.L. Finally, since our universe must be reheated by other means than inflaton decay, the authors employ the instant preheating mechanism, in which, as we will see in next section, the field is coupled to some other degrees of freedom. So, as the field is rapidly rolling down the quintessential tail of its runaway potential, it induces massive particle production, which after the decay into lighter ones produces the radiation bath of the hot Big Bang. Soon afterwards, the inflaton field freezes at some value with small residual energy density, which becomes important at present, playing the role of dark energy.

Simple models of Quintessential Inflation
In this subsection we will discuss very simple models of Quintessential inflation that either do not have enough physical motivation or do not agree with the experimental data.
We start with a model based on the U (1) Peccei-Quinn symmetry proposed in [97], where the real part of a complex field plays the role of the inflaton whereas the imaginary part is the quintessence field.
The Lagrangian is given by where arg(Φ) denotes the argument of the field Φ. Writing the field as follows, Φ = 1 √ 2 φe iϕ/f , one finds the following equations for the real and imaginary part of the field, where g(φ) = φ 2 f 2 and the corresponding potentials are given by The inflaton field φ moves in a potential like a "mexican-hat" and when it arrives at one of its minima it starts to oscillate and releases its energy producing particles as in standard inflation. On the contrary, the quintessence field ϕ rolls in the potential V (ϕ) producing enough dark energy to match with the current data provided that the parameters involved in the model satisfy Another simple model was built in [98,99] in the context of brane-world context, where the modified Friedmann equation for the flat FLRW geometry becomes where M pl,4 is the 4D reduced Planck's mass and the brane tension is given by λ = 3 In this context it has been shown that a potential composed by a sum of exponentials or a hyperbolic cosine leads to Quintessential Inflation. For this reason the authors, adopting natural units (M pl,4 = 1), choose as a potential For this potential, it is not difficult to calculate the spectral index and the tensor/scalar ratio as a function of the number of e-folds, which leads to the relation 6(n s − 1) = r, which is, as in the original Peebles-Vilenkin model, clearly incompatible with the recent Planck's data (see Figure 1). Anyway, we continue with the model. The chosen reheating mechanism is the Instant Preheating, thus obtaining, as we will see in next section, a reheating temperature around 10 9 GeV. The last simple Quintessential Inflation model that we deal with was introduced in [18], where the idea is to find an unstable non-singular solution of the equations of General Relativity (the Friedmann and Raychaudhuri equations), as in the Starobinsky model [100] or the famous Einstein's static model [101], in QI. The idea to build the model goes as follows: First of all, we look at early times for a Raychaudhuri equation of the formḢ = F (H) with F (H) = −α 2 H β , and note that it has a finite time singularity for β > 1. So, for large values of H, we choose F as a linear function of the Hubble rate. Secondly, to obtain a phase transition, one assumes that the derivative of F is discontinuous at some point H E (this discontinuity enhances the particles production of heavy massive particles obtaining a finite reheating temperature compatible with the BBN success [18]), choosing for H H E a kination regime, which corresponds to w ef f = 1 =⇒ F (H) = −3H 2 . And finally, if the model has to take into account the current cosmic acceleration, the simplest way is to assume that our dynamical systemḢ = F (H) has a fixed point H f , i.e., F (H f ) = 0, meaning that H f is a de Sitter solution, which can be modeled for An example with all the properties mentioned above iṡ where the parameters involved satisfy H e H f , and in order for F to be continuous one has to choose For that model the effective EoS parameter w ef f = −1 − 2Ḣ 3H 2 is given by which shows that for H H e one has an early quasi de Sitter period with w ef f ∼ = −1, when H e H H f the universe is in a kination phase (w ef f = −1), and finally, for H H f one has a late quasi de Sitter regime with w ef f ∼ = −1.
On the other hand, integrating the equationḢ = F (H), we can see that the non-singular solution is To end, the potential can be obtained using the Raychaudhuri equation, which leads to Then, for the model studied here, when H > H E , one gets and for H < H E a simple calculation shows that Finally, using the relation V (H) = (3H 2 +Ḣ)M 2 pl , one obtains the following potential, which depicts, at early times, a 1-dimensional Higgs potential, also called Double Well Inflationary potential, and, at late times, an exponential-quintessence potential.

REHEATING MECHANISMS
Different reheating mechanisms, such as gravitational particle production, instant preheating and curvaton reheating are revisited, with all the details, in this section.

Massless particle production
Here, we consider a massless quantum field χ (intensively considered in the early literature, see for instance [14,[43][44][45][46]), which will be the responsible for particle production. This kind of field only interacts with gravity and we only assume that the particles are nearly conformally coupled with gravity, i.e., the coupling constant is approximately 1/6 (ξ ∼ = 1 6 ) but not equal to 1/6, because free massless spinor and gauge fields are conformally-invariant, so they do not contribute to the energy density of the relativistic plasma formed by the produced particles.
Then, the modes, in Fourier space, satisfy the Klein-Gordon equation where τ is, once again, the conformal time and R(τ ) is the Ricci scalar curvature.
To define the vacuum modes before and after the phase transition, we assume that at early and late times the term a 2 R will vanish fast enough at early and late times, then its behavior at early and late times is respectively Therefore, the vacuum modes at early ("in" modes) and late times ("out" modes) (exact solutions of (170)) will be given by [43] and, since we are considering particles nearly conformally coupled to gravity, we can consider the term (ξ − 1/6)a 2 (τ )R(τ ) as a perturbation, and we can approximate the "in" and "out" modes by the first order Picard's iteration, i.e., inserting (171) in the right-hand side of (172), as which will represent, respectively, the vacuum before and after the phase transition.
Thus, after the phase transition, we could write the "in" mode as a linear combination of the "out" mode and its conjugate as follows, Imposing the continuity of χ and its first derivative at the transition time we get, up to order (ξ − 1/6) 2 , the value of these coefficients [44,57], Finally, in a simple model of gravitational production of quanta with negligible rest mass the energy density of the produced particles due to the phase transition is given by [102] where N S is the number of fermion and boson fields involved in the model, which for simplicity we will assume of the order 1, and the integral of the β-Bogoliubov coefficient (175) is convergent because at early and late times the term a 2 (τ )R(τ ) converges fast enough to zero. In addition, if at the transition time the first derivative of the Hubble parameter is continuous one has β k ∼ O(k −3 ), which means that the energy density of the produced particles is not ultraviolet divergent. Then, the energy density of the produced massless particles approximately becomes [43] ρ whereN is a dimensionless numerical factor.
Remark.-The numberN is clearly model dependent. In the scenario proposed by Ford in [43] where there is a transition from de Sitter to a matter domination modeled by a 2 (τ )R(τ ) ≡ 12 τ 2 +τ 2 0 , the numberN can be calculated analytically giving as a result 9 8 . However, note that in this case reheating is impossible because the energy density of the produced particles decreases faster than those of the background. To go beyond, we have calculated numerically this factor for some simple models that have a transition from a de Sitter regime to a kination one, and in all cases N is of the order 1 (see for instance [103]).
On the other hand, it is well-known that the thermalization of the produced particles is a very fast but not instantaneous process [104] (see also the Section III of [32]). In fact, following the work of Spokoiny [13] (see also [14]), we can consider the thermalization rate Γ th = n χ,kin σ 2→2 , where the most important process for kinetic equilibrium are 2 → 2 scatterings with gauge boson exchange, whose typical energy energy is E ∼ ρ 1/4 χ,kin , in the t-channel, which is given by σ 2→2 = α 3 E 2 (see for instance the section IV of [104]), where, as usual, α 2 ∼ 10 −3 [13].
Given that where for many models one finds [105] M ≡ 1 which is finite, because in all the considered models R = 6(Ḣ + 2H 2 ) is continuous and, as we have already explained, a 2 R vanishes fast enough at early and late times. So, we get that The relativistic fluid reaches the thermal equilibrium when H(t th ) ∼ H th ∼ Γ th [13,14] and, since during kination the Hubble rate scales as a −3 , one has meaning that the temperature when the thermalization is achieved is of the order T th = 30 π 2 g th Taking for instance H kin ∼ 10 −6 M pl (the typical value of the Hubble rate at the beginning of kination for many models), g th = 106.75 (the effective number of degrees of freedom of the Standard Model) and ξ − 1 6 ∼ 10 −2 , one gets a thermalization temperature of the order of 10 9 GeV.
Finally, the reheating occurs when the energy density of the background and the one of the created particles is of the same order (ρ χ (t reh ) ∼ ρ(t reh ) =⇒ a kin a reh ∼ ξ − 1 6 H kin M pl ), which implies, for the typical value H kin ∼ 10 −6 M pl , a reheating temperature of the order which for ξ − 1 6 ∼ 10 −2 and g reh = 106.75 leads to a low reheating temperature around 10 3 GeV.
6.1.2. Superheavy particle production: Calculations using the WKB approximation The study of the massive particle production in QI is a bit more involved than for massless particles. Therefore, to warm up, we will consider, once again, a toy model based on that of Peebles-Vilenkin, namely where m is the inflaton mass and M is another very small mass whose value can be calculated in the same way as in the Peebles-Vilenkin model: during the radiation and matter domination epochs the inflaton field is all the time of the order M pl (see [14] for a detailed discussion). Then, in the model the field will dominate at late times when where we have used that the current value of the Hubble parameter is H 0 ∼ 10 −61 M pl . Remark The mass of the inflaton field m can be calculated using the formula of the power spectrum of scalar perturbations P ζ = H 2 * 8π 2 M 2 pl ∼ 2 × 10 −9 , obtaining after a simple calculation and, since * ∼ 10 −2 , one finally gets m ∼ 10 −5 M pl .
On the other hand, regarding the creation of superheavy massive particles conformally coupled with gravity that have no interaction with the inflaton field, the time-dependent frequency of the χ-particles in the k-mode is ω k (τ ) = k 2 + m 2 χ a 2 (τ ), where m χ denotes the mass of the quantum field χ, that is, the χ modes, in Fourier space, satisfy the Klein-Gordon (KG) equation Before continuing, the following observation should be made: Remark.-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 of the universe is of the order of T ∼ 10 −3 M pl ∼ 10 15 GeV), the beginning of the classical Hot Big Bang (HBB) scenario is possible. Since for the flat FLRW universe the energy density of the universe, namely ρ, and the Hubble parameter H are related through the Friedmann equation ρ = 3H 2 M 2 pl and the temperature of the universe is related to the energy density via the Stefan-Boltzmann law ρ = (π 2 /30)g * T 4 (where g * = 106.75 is the number degrees of freedom for the energy density in the Standard Model), one can conclude that a classical picture of our universe might be possible when H ∼ 10 −5 M pl ∼ 10 13 GeV. Then, if inflation starts at this scale, i.e. taking the value of the Hubble parameter at the beginning of inflation as H B ∼ 10 −5 M pl , we will assume as a natural initial condition that the quantum χ-field is in the vacuum at the beginning of inflation. We will also choose the mass of the χ-field two orders greater than this value of the Hubble parameter (m χ ∼ 10 −3 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 [106]) because, as we will immediately see, the polarization terms will be sub-dominant and do not affect the dynamics of the inflaton field. So, we will choose 10 13 GeV ∼ m ∼ H B m χ ∼ 10 15 GeV M pl ∼ 10 18 GeV.
Thus, during the adiabatic regimes, i.e., when H(τ ) m χ =⇒ ω k (τ ) ω 2 k (τ ), one can use the WKB approximation [107] of the equation (187) where n is the order of the approximation to calculate the k-vacuum mode. When some high order derivative of the Hubble parameter is discontinuous, one has to match the k-vacuum modes before and after this moment, and it is precisely at this moment when one needs to use positive frequency modes after the breakdown of the adiabatic regime in order to perform this matching, which, following Parker's viewpoint [39], is the cause of the gravitational particle production.
For our toy model, note that the derivative of the potential is discontinuous at ϕ = −M pl , which means, due to the conservation equation, that the second derivative of the inflaton field is discontinuous at the transition time and, consequently, from the Raychaudhuri equationḢ = −φ 2 2M 2 pl one can deduce that the second derivative of the Hubble parameter is also discontinuous at this time.
In this case one only needs the first order WKB solution to approximate the k-vacuum modes before and after the phase transition where [108] Before the transition time, namely τ kin , which is very close to the beginning of kination (ϕ ∼ = −M pl ), the vacuum is depicted by χ W KB 1,k (τ ), but after the phase transition this mode becomes a mix of positive and negative frequencies of the form α k χ W KB . Then, the β-Bogoliubov coefficient could be obtained matching both expressions at τ = τ kin , obtaining where is the Wronskian of the functions f and g at the transition time, and being F (τ ± kin ) = lim τ →τ ± kin F (τ ).
The square modulus of the β-Bogoliubov coefficient will be given approximately by [28] |β where we have used that at the transition time all the energy at the end of inflation, which is approximately 1 2 m 2 M 2 pl because ϕ EN D = − √ 2M pl , was converted into kinetic energy.
Thus, for our model, the number density of the produced particles and its energy density will be Now, one has to note that there are two different situations, namely, when the superheavy massive particles decay before and after the end of the kination regime.
1. Decay of the superheavy particles into lighter ones before the end of kination.
Let Γ be the decay rate of the superheavy particles. The decay is practically finished when Γ is of the same order of the Hubble rate, i.e., when Γ ∼ H(t dec ) = H kin a kin a dec a kin a dec

3
, and thus, the corresponding energy densities will be Since the decay is before the end of kination, one has Γ < H kin and ρ χ,dec < ρ ϕ,dec , which leads to the following bound, 10 −5 m 2 m χ M pl which for m ∼ 10 13 GeV, m χ ∼ 10 15 GeV and H kin ∼ 10 12 GeV becomes 10 −6 < Γ/GeV < 10 12 . (197) Then, assuming as usual nearly instantaneous thermalization, the reheating temperature, i.e., the temperature of the universe when the relativistic plasma in thermal equilibrium starts to dominate (ρ χ,reh ∼ ρ ϕ,reh ), will be obtained taking into account that and thus, which for the values of our parameters leads to a reheating temperature of the order GeV.
Finally, using the bound of the decay rate, we deduce that the reheating temperature is bounded by 3 × 10 2 GeV < T reh < 10 6 GeV.

Decay of the superheavy particles into lighter ones after the end of kination
First of all, recall that we have already seen that at the end of kination one has Taking into account that the decay is after the end of kination we will have the bound Γ < H end ∼ 10 −6 GeV, and thus, since the thermalization is nearly instantaneous, the reheating time coincides when the decay is completed, i.e, when H ∼ Γ, and consequently, the reheating temperature will be which is bounded by where the lower bound is obtained taking into account that the reheating will be before the Big Bang Nucleosynthesis which occurs around 1 MeV.

The diagonalization method
In the last subsection we have used the WKB method to calculate the particle production when some derivative of the potential has a discontinuity, but when the potential is smooth enough it does not work because one needs the exact solution of the equation (187).
So, here we will explain the so-called diagonalization method, which is essential to calculate the massive particle creation in more realistic scenarios such as LQI and α-QI.
The idea of the method goes as follows: Given the quantum scalar field χ of superheavy particles conformally coupled to gravity satisfying the KG equation (187), the modes that define the vacuum state, at a given initial time τ i , must satisfy the condition and the energy density of the vacuum is given by [109] ρ where in order to obtain a finite energy density [59] we have subtracted the energy density of the zero-point oscillations of the vacuum Following the method developed in [110] (see also Section 9.2 of [59]), we will write the modes as follows, where α k (τ ) and β k (τ ) are the time-dependent Bogoliubov coefficients. Now, imposing that the modes satisfy the conditions one easily shows that the Bogoliubov coefficients must satisfy the fundamental system in order for the expression (207) to be a solution of the equation (187).
Finally, inserting (207) into the expression for vacuum energy density (206), and taking into account that the Bogoliubov coefficients satisfy the equation |α k (τ )| 2 − |β k (τ )| 2 = 1, one finds that At this point, it is very important to notice that |β k (τ )| 2 encodes the vacuum polarization effects and also the particle creation, which only happens when the adiabatic evolution breaks. In fact, the quantity was named in the Russian literature as the number density of quasi-particles, which is very different from the number density of the produced particles because, as we will see immediately, it also contains the vacuum polarization effects.
As an application, we will consider our toy model (184) once again. Then, in order to obtain the value of the β-Bogoliubov coefficient we come back to the equation (209) and in the first approximation we take α k (τ ) = 1, getting After integrating by parts before the beginning of kination, it yields However, after kination the β-Bogoliubov coefficient must be given by where the constant C has to be chosen in order that the β-Bogoliubov coefficient becomes continuous at τ kin because the equation (209) is a first order differential equation, and then one has to demand continuity to the solution. Thus, after some cumbersome calculation [111] one has where we have denoted a kin ≡ a(τ kin ).
Notice that the terms that do not contain C lead to sub-leading geometric quantities in the energy density. Effectively, the term − ω k (τ ) 4iω 2 k (τ ) leads to the following contribution to the energy density: 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 seen are always sub-leading, are vacuum polarization effects, and they rapidly disappear in the adiabatic regime. That is, soon after the beginning of kination |β k (τ )| 2 approaches to |C| 2 , whose value coincides with the one obtained using the WKB approach (see formula (192)), getting once again which at the beginning of kination is sub-dominant with respect to the energy density of the background but it will eventually dominate because the one of the background decreases during kination as a −6 (τ ).
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 quasi-particles are created and annihilated due to the interaction of the quantum field with gravity [59], i.e., this is a vacuum polarization effect. And, following this interpretation, the number density of the created quasi-particles at time τ is given by N χ (τ ) = 1 2π 2 a 3 (τ ) ∞ 0 k 2 |β k (τ )| 2 dk. However, one has to be very careful with this interpretation and especially keep in mind that, as we have shown in our toy model (184), real particles are only created when the adiabatic regime breaks. Effectively, before the beginning of kination the main term of the β-Bogoliubov coefficient is given by − ω k (τ ) 4iω 2 k (τ ) , whose contribution to the energy density is , and to the number density of quasi-particles is , and thus, at time τ before the beginning of kination ρ χ (τ ) = m χ N χ (τ ), meaning that the quasi-particles do not evolve as real massive particles. On the contrary, during kination the leading term of N χ (τ ) is given by 10 −5 m mχ 3 m 3 a kin a(τ ) 3 , so we have ρ χ (τ ) = m χ N χ (τ ) and the decay follows a −3 (τ ), which justifies the interpretation of massive particle production.
So, once we have understood the gravitational particle production, we review our last work [112] about the creation of particles in the LQI scenario. Coming back to our potential (104), first of all one has to integrate numerically the conservation equation for the inflaton field, namelÿ where H = 1 √ 3M pl φ 2 2 + V (ϕ), with initial conditions at the horizon crossing (when the pivot scale leaves the Hubble radius). Recalling 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, ϕ * = −0.154M pl andφ * = 0, where the "star" denotes, as usual, that the quantities are evaluated at the horizon crossing.
Once the evolution of the background is obtained, and in particular the evolution of the Hubble rate, one can compute the evolution of the scale factor, which is given by where one can choose a * = 1 as the value of the scale factor at the horizon crossing.
From the evolution of the scale factor, one can see in Figure 9 that a spike appears during the phase transition from the end of inflation to the beginning of kination, i.e., when the adiabatic evolution is broken and particles are gravitationally produced. , with initial conditions α k (τ * ) = 1 and β k (τ * ) = 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 it is useful to transform the equation (209) into a second order differential equation, namely Given that α k (τ * ) = 1 and β k (τ * ) = 0 leads to α k (τ * ) = 0, one is interested in solving the equation for α k (τ ), which can be split into the real and imaginary form in the following way, and then |β k (τ )| 2 = |α k (τ )| 2 − 1 because of the well-known conservation property of the Wronskian. For the value k = a kin H kin , one can see in Figure 11 that |β k (τ )| 2 stabilizes soon to a non-zero value after the beginning of kination, containing only particle production effects. One can numerically check that this happens for the range 0.05 whereρ χ ∼ 10 44 GeV 4 ∼ 10 −30 M 4 pl andā are, respectively, the energy density of the produced particles and the value of the scale factor at the end of the non-adiabatic phase. Finally, the energy density of the background at this moment is given byρ ϕ = 3H 2 M 2 pl ∼ 10 57 GeV 4 ∼ 10 −17 M 4 pl , showing that the energy density of the produced particles is sub-leading at 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 .
Once again, two different situations arise: 1. Decay of the superheavy particles into lighter ones before the end of kination.
In this case, the energy density of the the inflaton field and the one of the relativistic plasma, when the decay is finished, i.e., when Γ ∼ H dec =H Imposing that the end of the decay precedes the end of kination, that means, ρ χ,dec ≤ ρ ϕ,dec , one gets Γ ≥ 10 −21 M pl , and, since the decay is after the beginning of the kination and for our LQI model H kin ∼ 10 −8 M pl one gets Γ ≤ H kin ∼ = 4 × 10 −8 M pl . So, the following bound for the decay rate is obtained, Finally, the reheating temperature, i.e., the temperature of the universe when the relativistic plasma in thermal equilibrium starts to dominate, which happens when ρ ϕ,reh ∼ ρ χ,reh , can be calculated as follows: Since after the decay the evolution of the respective energy densities is given by where g reh = 106.75 is the effective number of degrees of freedom for the Standard Model. So, taking into account the bound (251) the reheating temperature ranges between 10 3 GeV and 10 7 GeV.
2. Decay of the superheavy particles into lighter ones after the end of kination.
In the case that the decay of the χ-field is after the end of kination, one has to impose Γ ≤ H end . Now, taking into account that at the end of kination where it is used that the kination ends when ρ χ,end ∼ ρ ϕ,end , meaning ā a end 3 =ρ χ ρϕ . So, the condition Γ ≤ H end leads to the bound 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 ∼ Γ, obtaining where one has to take into account 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

Instant preheating
In this subsection we review the work of Felder, Kofman and Linde [50] (see also [49]), where Instant Preheating for QI scenarios was introduced and discussed in detail.
Essentially, instant reheating is based on the interacting part of the Lagrangian density, namely − 1 2 g 2 ϕ 2 χ 2 , where g is a coupling constant, which due to the phase transition between the end of inflation and the beginning of kination ceases to be in the vacuum state to produce superheavy massive particles. These created particles will interact decaying into light ones, becoming a relativistic plasma whose energy density eventually dominates those of the background and as a consequence the universe will thermalize and will reheat.
To deal with particle production, note that, in Fourier space, the k-mode of the quantum field χ satisfies the equation of a time dependent harmonic oscillator where the derivative is with respect to the conformal time τ and the square of the frequency is given by being once again m χ the bare mass of the quantum field, R(τ ) the scalar curvature and ξ the coupling constant with gravity.
To clarify ideas about preheating (i.e., particle production) we consider the simplest toy model although the reasoning will serve for a general class of Quintessential Inflation potentials.
Near the phase transition, which actually corresponds to the beginning of kination (ϕ kin = 0), one has ϕ(τ ) ∼ = ϕ kin τ . Then, if we simplify choosing a conformal coupling ξ = 1 6 , we can approximate the frequency ω k (τ ) by k 2 + a 2 kin (m 2 χ + g 2 (ϕ kin ) 2 τ 2 ), where we have disregarded the expansion of the universe taking a(τ ) = a kin in order to obtain the well-known overbarrier problem in scattering theory [113], which can be analytically solved. Effectively, the β k -Bogoliubov coefficient is related with the reflection coefficient via the formula (see [114][115][116] and [117,118] for a mathematical explanation) where γ is a closed path in the complex plan that contains the two turning points τ ± = ±i √ k 2 +a 2 kin m 2 χ ga kin ϕ kin . A simple calculation yields that the number of particles in the k-mode is given by [116] n k ≡ |β k | 2 = e − π(k 2 +a 2 kin m 2 χ ) ga kin ϕ kin .
Then, the number density of produced particles is given by in terms of the cosmic time t [102], namely Analogously, the energy density of the produced particles is given by [102] ρ Then, at the phase transition one has and at late times ρ χ (t) ∼ = m 2 χ + g 2 ϕ 2 (t)n χ (t), which means that at late times the χ-particles acquire an effective mass equal to Now, we will see that three constraints must be imposed [50] in order to have a viable model: 1. If we do not want an exponential suppression of the energy density, one has to choose a bare mass satisfying m χ ≤ √ gφ kin . In fact, for the sake of simplicity, we will take m χ = 0, and thus, the effective mass is given by m ef f = g|ϕ|.
2. For masses -the effective mass of the field χ is g|ϕ|-greater than the Hubble parameter the vacuum polarization energy density due to the field χ, which for H g|ϕ| can be calculated using the WKB approximation, is of the order H 6 g 2 ϕ 2 [119]. This quantity is smaller than the energy density of the background (∼ H 2 M 2 pl ) when the pivot scale leaves the Hubble radius, meaning that the polarization effects will not affect the last stages of inflation. Therefore, since H ≤ Since the decay must be before thant, the condition t dec <t , where t dec ∼ = 1 3Γ is the time when the field χ decayed, i.e., H dec ∼ Γ, leads to the relation g < 10 2 Γ M pl 2/5 , which together with the condition g m M pl ∼ 10 −6 constrains the value of the decay rate to satisfy Γ 10 −20 M pl ∼ 10 −2 GeV.
Let's now calculate the temperature at the reheating time, which as we have already seen is given by Taking into account that, if t dec > 10 6 M −1 pl =⇒ Γ < 10 −6 M pl , the effective mass is approximately gM pl , one has ρ ϕ,dec = 3Γ 2 M 2 pl and ρ χ,dec ∼ gM pl n χ,dec ∼ 10 −2 g 5/2 H kin M pl M 2 pl Γ, getting that The values of the parameters g and Γ must satisfy the constraint 10 −6 g < 10 2 Γ M pl
Finally, one can conclude that the value of the parameters involved in the theory must accomplish all the following requirements:

Curvaton reheating in Quintessential Inflation
In this section we review the so-called curvature reheating mechanism in Quintessential Inflation. To do that we follow [52], but taking into account that for the authors of that paper, contrary to our convention, the reheating temperature is the temperature of the universe when the curvaton field has totally decayed in relativistic particles. Recall that the universe is really reheated when the energy density of inflaton is of the same order than the one of the decay products, provided that the curvaton field had previously decayed.
We assume that the potential of the curvaton field, namely σ, is quadratic, V (σ) = 1 2 m 2 σ σ 2 , where the mass of the curvaton is chosen to be smaller than the value of the Hubble parameter at the end of inflation m σ H end . At this moment the curvaton field is in a slow-roll regime because the condition m σ H end ⇐⇒ V σσ H 2 end means that the curvaton potential is flat enough at the end of inflation [120]. Then, in order to avoid a second inflationary stage now driven by the curvaton, one has to impose that its energy density is subdominant when the curvaton starts to oscillate, which happens when m σ ∼ = H [120] (see also the section 5.4.1 of [121] for a detailed discussion of the quadratic potential). Then, ρ σ (t osc ) ≡ ρ σ,osc < ρ ϕ,osc = 3H 2 osc M 2 pl , where t osc is the time when the curvaton starts to oscillate. Taking into account that H osc ∼ = m σ and since at the beginning of the oscillations, i.e., at the end of the slow-roll period for the curvaton, Hσ ∼ V σ =⇒σ osc ∼ m σ σ osc , one has ρ σ,osc ∼ = m 2 σ σ 2 osc obtaining the bound σ 2 osc < 3M 2 pl .
2. The curvaton decays when the curvaton field dominated the universe.
Assuming once again instantaneous thermalization, since the curvaton decays when it dominates, that is, when ρ ϕ,dec ≤ ρ σ,dec , the reheating time will occur at the decay time (H reh = H dec ∼ Γ).
Finally, from these values and the equation (259), one can conclude that only for curvaton light masses satisfying m σ ≤ 5 × 10 −16 M pl ∼ 10 3 GeV a reheating temperature compatible with the BBN success is obtained.

Overproduction of GWs
It is well-known that a reheating due to the gravitational production of light particles is incompatible with the overproduction of Gravitational Waves because both satisfy the same equation, and thus, they scale in the same way. For this reason it is impossible that at the reheating time the ratio of the energy density of the GWs to the energy density of the relativistic plasma was smaller that 10 −2 , which is essential to prevent any complication in the BBN success [14]. In fact, the minimum value of this ratio, which is obtained when the light particles are minimally coupled with gravity (ξ = 0), is of the order 1/N s , where N s is the number of light fields. As has been explained in [14] in a minimal GUT theory N s = 4, which corresponds to the electroweak Higgs doublet, and only in supersymmetric theories N s is of the order 10 2 (see the explanation given by Peebles-Vilenkin in [14]).
However, as we will see in this subsection, for our LQI model, when the reheating is due to the gravitational creation of heavy particles, the overproduction of GW's do not modify the BBN.
So, to show this, first of all we recall that the production of Gravitational Waves during the phase transition from the end of inflation to the beginning of kination is [124] ρ GW (τ ) ∼ = 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 [70] ρ GW,reh ρ χ,reh ≤ 10 −2 .
First of all we see that the constraint (262) is overpassed when the decay of the superheavy massive particles is previous to the end of kination. Effectively, if the decay occurs before the end of kination one has ρ GW,reh ρ χ,reh = ρ GW,dec ρ χ,dec because the light relativistic particles evolve as the GW's. Now, taking into account that during kination the energy density of the background scales as a −6 , it yields that a kin a dec 4 = ρ ϕ,dec ρ ϕ,kin and then, since in LQI one has H kin ∼ 4 × 10 −8 M pl and using the results obtained previously (see the formula (222), we get ρ GW,reh ρ χ,reh ∼ = 10 −1 Γ M pl Therefore, the bound (262) is overpassed when which is completely compatible within the bound (251).
On the contrary, when the decay of superheavy particles is produced after the end of kination, and assuming once again instantaneous thermalization, the reheating time will coincide with the decay one. Then, since ρ χ,dec = 3Γ 2 M 2 pl and H dec = H end a end a dec and, using that a kin a end 4 = ρ ϕ,end ρ ϕ,kin and ρ ϕ,end =ρ we get ρ GW,dec = 2 16 and thus, where we have used the bound (258), so the constraint (262) is clearly overpassed.

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 postinflationary 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 [125], 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 [66]), 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 to the horizon scale at the BBN and at the end of inflation respectively. As has been shown in [126], the main contribution of the integral (272) 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 ≤ 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 [46] (see also [125,[127][128][129] where the graviton spectra in quintessential models have been reassessed, in a model-independent way, using numerical techniques), where is the amplitude of the GWs; Ω r ∼ = 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 [125,126] for more details). As has been derived in [126], the specific form of the expression above comes from the behavior of the Hankel functions for small arguments. Now, plugging expression (273)  To calculate the ratio k kin /k end , we will have to study the following two different situations: 1. When the decay of superheavy particles occurs before the end of kination.
In this case the reheating time coincides with the end of kination, so a simple calculation leads to Using the formulas ρ ϕ,reh = ρ ϕ,kin a kin a reh 6 and ρ χ,reh = ρ χ,dec a dec a reh .