A Review of Quintessential Inflation

Some of the most important quintessential inflation scenarios, such as the Peebles–Vilenkin model, are described in detail. These models are able to explain the early- and late-time accelerated expansions of our universe, and the phase transition from the end of inflation to the beginning of kination where the adiabatic evolution of the universe was broken in order to produce enough particles to reheat the universe with a viable temperature, thereby aligning with the Hot Big Bang universe. In addition, while considering the reheating to be due to the gravitational production of superheavy particles conformally coupled to gravity, we checked that the considered scenarios do not suffer problems due to the overproduction of gravitational waves at the end of inflation, and thus the validity of Big Bang nucleosynthesis is preserved.


Introduction
Understanding the universe's evolution is one of the greatest mysteries in the history of humanity. The primary questions have always been: "where do we come from and where we are going?" In particular, its early and late expansions have been studied a great deal lately. Looking at the scientific literature, one can find two popular and well accepted theories-though they are not observationally proven-namely, the inflation (the early surprisingly quickly accelerated expansion of our universe) and the quintessence as a form of dark energy (the current cosmic acceleration). The inflationary paradigm [1][2][3] describes a very fast, accelerating phase of the early universe that lasted for an extremely small amount of time and solves a number of shortcomings associated with standard Big Bang cosmology, such as the horizon problem, the flatness problem and the primordial monopole problem. The predictive power of inflation was soon recognized due to its ability to explain the origins of inhomogeneities in the universe as quantum fluctuations during this epoch [4][5][6][7][8], because such an explanation aligns well 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 80s, is still considered the best way to explain the recent observational data, because it is nowadays the simplest theory that describes the early universe almost in agreement with the recent observations [9].
On the other hand, one of the most accepted explanations for the current cosmic acceleration relies on the 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 was introduced. By construction, unlike the standard quintessence models (see [13] for a review), these models-named quintessential inflation (QI) models [14][15][16]only contain one classical scalar field, also named inflaton as in standard inflation [1][2][3]5], and it has been shown that they succeed in reproducing the two accelerated epochs of the universe's expansion (see also [17][18][19][20][21][22][23][24][25][26][27][28] for other interesting QI models). This idea to since their energy density decreased as a −3 before decaying in lighter particles and as a −4 after that, they eventually dominated the energy density of the inflation whose decrease was as a −6 . Thus, the universe was reheated.
Coming back to the original Peebles-Vilenkin model [15], 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 [15] 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, Big Bang nucleosynthesis (BBN) 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 can be obtained; recall that the radiation dominated era was prior to the Big Bang nucleosynthesis (BBN) epoch, which occurred in the 1 MeV regime [62]. As a consequence, the reheating temperature had 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 considered in 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 gravitational waves (GWs), which are also produced during the phase transition in BBN by satisfying the observational bounds coming from the overproduction of the GWs [15] or those related to the logarithmic spectrum of their energy density [66]. As we will see throughout this review, the overproduction of GWs does not disturb BBN's 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 of this topic. 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 2 we review, in great detail, the seminal Peebles-Vilenkin model, calculating the values of the two parameters on which the model depends; the evolution of the inflaton field throughout history; and the number of e-folds as a function of the reheating temperature. Additionally, we introduce models that improve the original one. Section 3 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 the current observational data, the potential should be modified, introducing a pure exponential term. In Section 4, 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 5, showing that we obtain the same results as in the case of LQI. In Section 6, we review other quintessential inflation models and in particular some aspects of the work of K. Dimopoulos in QI. Next, in Section 7 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 8, we deal with the constraints to preserve BBN's success from the logarithmic spectrum of GWs and also from BBN's 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 surpassed. The last section is devoted to the conclusions of the review.

The Peebles-Vilenkin Model
The first quintessential inflation scenario was proposed by Peebles and Vilenkin in their seminal paper [15] 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 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: Thus, the spectral index and the ratio of tensor to scalar perturbations are given by [67] n s ∼ = 1 − 6 * + 2η 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). Furthermore, taking, for example, n s ∼ = 0.96, which is approximately its observational value [68,69], we get λ ∼ π 2 36 (1 − n s ) 3 × 10 −9 ∼ 10 −14 .
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 in the order of 40M pl (see Figure 1). Thus, 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 Figure 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,51], 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 when taking into account that when the pivot scale leaves the Hubble radius e αϕ * /M pl 1, and thus, the spectral index and the tensor/scalar ratio are equal to leading to 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 = 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 leading to 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 Figure 1).

Dynamical Evolution of the Peebles-Vilenkin Model: From the Beginning of Kination to the Matter-Radiation Equality
Now we want to understand the evolution of the quintessential inflation models after inflation, which only depends on the tail of the potential. We focus on the Peebles-Vilenkin one, and for that case, inflation ends when = 1, which occurs for ϕ END = 2 √ 2M pl .
By 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 e f f , satisfies where we have used the Friedmann and Raychaudhuri equations ρ being the energy density and P being the corresponding pressure. Then, from (22) one gets the following relation at the end of inflation:φ 2 END = V(ϕ END ). Thus, 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. Thus, 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 e f f ∼ = 1. For this value, combining the Friedmann and Raychaudhuri Equation (23), one getṡ 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 must be broken and particles created. Here we assume that these particles, which we consider to have been very massive and conformally coupled to gravity, were gravitationally produced. Thus, they had to decay into lighter ones in order to form a relativistic plasma which eventually dominated the energy density of the universe and became reheated. Then, since the kination regime ended when the energy density of the inflaton field was of the same order as that of the produced particles, two different scenarios could have come about: 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 scenario, 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 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 becomes 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 The next step is to calculate the value of the field at the beginning of the matter-radiation equality. During the radiation domination the effective EoS parameter is w e f f = 1/3. Thus, combining once again the Friedmann and Raychaudhuri equations, we haveḢ = −2H 2 , whose solution is H(t) = 1 2t . Then, an approximate solution of the conservation equation can be obtained disregarding once again the potential, because during the radiation domination the potential energy is irrelevant compared with the kinetic one. Hence, the equation has the following solution,φ (t) =φ reh t reh t

Dynamical Evolution of the Peebles-Vilenkin Model: From the Matter-Radiation Equality up to Now
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 chose to be the number of e-folds up to the present epoch, namely, N ≡ − ln(1 + z) = ln a a 0 . Now, using the variable N, one can recast the energy densities of radiation and matter, respectively, as ρ r (a) = ρ r,eq a eq a 4 =⇒ ρ r (N) = ρ eq 2 e 4(N eq −N) (46) and ρ m (a) = ρ m,eq a eq a where ρ eq = ρ m,eq + ρ r,eq and the value of the energy density of the matter or radiation ρ m,eq = ρ r,eq during the matter-radiation 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 during 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, Moreover, the Friedmann equation now looks as such: 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 were obtained analytically in the previous subsection. The value of parameter M can be obtained at N = 0, by making Equation (50) equal 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 and , from kination to future times. Right: The effective equation of state parameter w e f f , from kination to future times. As one can see in the picture, after kination the universe entered into a large period of time where radiation dominated. Then, after the matter-radiation equality, the universe became matter-dominated, and finally, near the present, it entered into a new accelerated phase where w e f f approaches −1.

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 [70] k * a 0 H 0 = e −N H * H 0 a END a kin a kin a reh a reh a eq a eq a 0 , where, as previously, a is the scale factor and a END , 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, respectively 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.02 Mpc −1 , we obtain N ∼ = −5 + ln H * H 0 + ln a END a kin + 1 4 ln g eq g reh 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 , and the Stefan-Boltzmann law ρ eq = π 2 30 g eq T 4 eq , ρ reh = π 2 30 g reh T 4 reh , g eq = 3.36 being the degrees of freedom during the matter-radiation equality.
and by 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 and using the value of the spectral index n s ∼ = 0.96 and assuming that ln a END 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 constraining 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 significantly increases the number of e-folds.

Exponential Quintessential Inflation
In this section we will consider more complex models than the Peebles-Vilenkin one: the same exponential inflation-type 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 by 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 theoretically find the possible values of the parameter γ is to combine the Equations (61) and (62) to get n−1 r 8n 2 n/2 (66) 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,70]-the energy density of the scalar field is only kinetic [14,38]-after the end of the inflationary period.
For potentials of this kind, 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 END = 1, meaning that at the end of this epoch the field reaches the value ϕ END = Then, the number of e-folds is given by Thus, by 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 ϕ END = 2 n 2 γ 2 1 2n−2 M pl , when the effective equation of state (EoS) parameter is equal to −1/3, meaning thatφ 2 END = V(ϕ END ), the energy density at the end of inflation is (69) and the corresponding value of the Hubble rate is given by which greatly constrains the values of the parameter γ because in all viable inflationary models at the end of the inflation the value of the Hubble rate is in the order of 10 −6 M pl [2]. In fact, when (70) is in the order of 10 −6 M pl , one gets Therefore, to perform numerical calculations, we used the values of n = 10 and r = 0.02, and thus, for γ ∼ 10 −16 we obtained approximately 67 e-folds, which is a viable value in quintessential inflation.
Next, we have to find out the time when kination started, which can be chosen by the moment when the effective EoS parameter is very close to 1. Assuming, for instance, that kination starts at w e f f ∼ = 0.99, we have numerically obtained ϕ kin ∼ = 47M pl anḋ ϕ 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 H = 1 √ 3M pl ρ ϕ + ρ r . 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φ + 3Hφ On the other hand, equating H(t) = 1 2t with H = 1 √ 3M pl ρ ϕ + ρ r , which yields the relation which 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, [10]): one can introduce the dimensionless variables 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 constraintx One can see that for w = 1/3 the dynamical system (80) has the attractor solutioñ γ , whose energy density scales as radiation, and Ω ϕ = 4 γ 2 .

The Tracker Solution
A tracker solution [71,72] 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. Note also that the system (80) when w = 0 (in the matter domination era) has the fixed point (x,ỹ) = γ/ √ 6, 1 − γ 2 /6 , which of course corresponds to the solution (86). 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 H 0 e γ 2 H/2 , the tracker solution has the following expression (see Section V of [73] 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 a eq 4 and the adiabatic evolution after the matterradiation 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 3/4 = π 2 30 g eq g reh g 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 a 0 . 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 introduce the dimensionless variables where K is a parameter that we will choose accurately in order to facilitate the numerical calculations. Thus, 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 writē 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 use KM pl = 10 −17 GeV 2 , yieldinḡ 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 by numerically integrating 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 , as seen in Figure 4, thereby proving that the inflaton field does not belong to the basin of attraction of the scaling solution for large values of n. In addition, considering 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 [15], 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 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 e f f = γ 2 3 − 1.

Remark 1.
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, thereby 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 ν,e f 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 of H = 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. and (green curve), from kination to future times, for a reheating temperature around 10 9 GeV. To perform numerical calculations we have taken n = 10 and γ = 0.1. Right: The effective equation of state parameter w e f f , from kination to future times, for n = 10 and γ = 0.1. As one can see in the picture, after kination the universe theoretically entered into a long period of time where radiation dominated. Then, after the matter-radiation equality, the universe became matter-dominated; and finally, near the present, it entered into a new accelerated phase where w e f f approaches 0.1 2 3 − 1 ∼ = −1that 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 [74,75] considered the following ansatz: where is the main slow-roll parameter; N = ln a a 0 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 Equation (37) in [76]. However, in this work we 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 those shown in [77], it is better to work with the simplified potential, namely, We can see the shape of the potential in Figure 6, where the inflationary epoch is shown on the left-hand side of the graph, and the dark energy epoch is on the right-hand side. The left-hand side shows the inflationary energy density and the right-hand side shows the late dark energy density.

Calculating the Values of the Parameters Involved in the Model
We start with the main slow-roll parameter, which is given by and since inflation ends when END = 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 ϕ END is close to zero. Thus, we chose γ 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, [2]). 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−n s ) 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 chose our parameters to satisfy 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 have γϕ 0 /M pl 1 where ϕ 0 denotes the current value of the field. Thus, we have V(ϕ 0 ) ∼ λM 4 pl e −ξ , which is the dark energy at the present time, meaning that Thus, 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 set ξ = γ, since it could 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 by 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 by introducing the current value of the temperature of the universe T 0 ∼ 9.6 × 10 −32 M pl we get Thus, 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 END 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 jeopardize the nucleosynthesis, one needs temperatures lower than 10 9 GeV. Thus, we assume that 1 MeV ≤ T reh ≤ 10 9 GeV, which leads to constraining the number of e-folds to 58 N 67. Furthermore, for this number of e-folds, 0.966 n s 0.970, which is within its 2σ CL range, as can be seen in Figure 7.

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 andφ 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. and , from kination to future times, for a reheating temperature around 10 9 GeV. Right: The effective equation of state parameter w e f f , from kination to future times. As one can see in the picture, after kination the universe entered in a large period of time where radiation dominated. Then, after the matter-radiation equality, the universe became matter-dominated; and finally, near the present, it entered a new accelerated phase where w e f f approaches −1.

α-Attractors in Quintessential Inflation
The concept of the α-attractor, in the context of standard inflation, was introduced for the first time in [78], providing models which generalize the well-known Starobinsky model (see [79] 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.
providing the following potential: where we have introduced the dimensionless parameter n = κ √ 6α. Similarly to [74,75,77], 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 [77].

Calculating the Values 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 END = 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, H 2 * ∼ = λM 2 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].
To end, some remarks are in order:

2.
If one takes α to be very small-for example, in the order of 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 with values of α ∼ 10 −2 and Λ ∼ 10 −120 .

3.
Another important important application of the α-attractors is its use to alleviate the current Hubble tension [80]. In that work, the authors, in the framework of αattractors, included in 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 [80] the authors used the following potential for EDE:

Present and Future Evolution: Numerical Analysis
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 anḋ ϕ kin ∼ ×10 −6 M 2 pl . The results presented in Figure 9 are the same as the ones obtained in the LQI scenario. Thus, both models predict the same evolution of the universe.
, from kination to future times for a reheating temperature around 10 9 GeV. Right: The effective equation of state parameter w e f f , from kination to future times. As one can see in the picture, after kination the universe entered a large period of time where radiation dominated. Then, after the matter-radiation equality, the universe became matter-dominated; and finally, near the present, it entered into a new accelerated phase where w e f f approaches −1.

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 contributions of K. Dimopoulos in this field, starting with the class of models introduced in [30], 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 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. Furthermore, in order to align 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 [81], 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 requires that q and M satisfy the constraint 0 < q M pl M < √ 2. 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) in the following form: Λ = −λM 2 pl e −n . 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 (122) because e −n 1. Thus, as we have shown for the potential (122), for small values of α inflation also work 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 [47]) Another important contribution of Dimopoulos is the introduction of the concept of warm quintessential inflation in [82], which is applied to the original Peebles-Vilenkin model, justifies the unification of early-and late-time acceleration and providing a very natural mechanism for reheating our universe. In the same way, an interesting idea is to consider quintessential inflation in the context of Palatini R 2 -gravity [83], 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 realizing its energy and creating particles as in standard inflation. However, as is highlighted in [84][85][86][87], there is an important limitation of the treatment considered in [54], namely, the fact that reheating through this mechanism is far from being a 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 mechanism can also be obtained from a combination of supergravity and string theory [88]. Said model contains the inflaton and a Peccei-Quinn field, which is responsible for particle production. In fact, both fields are coupled as in the theory of instant preheating, which we 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 considered responsible for dark energy.
Finally, the curvaton rehating, which we deal with in next section, is applied to QI in [89], 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 studied 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 [90] 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 α,α, γ and β being the parameters of the model. For this model, the slow-roll parameters, as functions 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 in a multitude of problems. Firstly, the flatness of the quintessential tail may be lifted by radiative corrections. Furthermore, 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. However, if the quintessential tail is too steep, regarding the present, 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 [91] 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 e f f approaches to −1 at late times) for the present time, leading to the current acceleration.
Therefore, here we analyze the paper [91]. See also [92] 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, favored by Planck, and 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 employed the instant preheating mechanism, in which, as we will see in next section, the field is coupled to some other degrees of freedom. Thus, 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 is important for the present, playing the role of dark energy.

Simple Models of Quintessential Inflation
In this subsection we 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 [93], 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 Φ. By writing the field as Φ = 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 [94,95] 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 mass and the brane tension is given by λ = 3 , M pl,5 being the 5D reduced Planck mass. 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), chose 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 Collaboration data (see Figure 1). Anyway, we continue with the model. The chosen reheating mechanism is the instant preheating, providing, as we will see in next section, a reheating temperature around 10 9 GeV. Finally, the model exhibits transient acceleration at late times for 0.96 ≤ Aα 2 ≤ 1.26 and 271 ≤ αϕ 0 ≤ 273, whereas eternal acceleration is necessarily obtained for 2.3 × 10 −8 ≤ Aα 2 ≤ 0.98 and 255 ≤ αϕ 0 ≤ 273.
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 [96] or the famous Einsteinian static model [97], 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. Thus, 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's success [18]). Then for H H E , one chooses a kination regime, which corresponds to w e f f = 1 =⇒ F(H) = −3H 2 . 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 For that model the effective EoS parameter w e f 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 e f f ∼ = −1, when H e H H f the universe is in a kination phase (w e f f = −1), and finally, for H H f one has a late quasi de Sitter regime with w e f f ∼ = −1.
On the other hand, integrating the equationḢ = F(H), we can see that the nonsingular solution is which can be approximated by when t > 0. (165) 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 pl , one obtains the following potential: which depicts, at early times, a 1-dimensional Higgs potential, also called the 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, [15,[43][44][45][46]), which is considered 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 in early and late times; thus, 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)) are given by [43] and since we are considering particles to be 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 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: By imposing the continuity of χ and its first derivative at the transition time, we get, up to order (ξ − 1/6) 2 , the values 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 [98] where N S is the number of fermion and boson fields involved in the model, which for simplicity we assume to be in 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 2.
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 , 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 that of the background. To go beyond, we have calculated this factor numerically for some simple models that include a transition from a de Sitter regime to a kination one, and in all casesN is of the order 1 (see, for instance, [99]).
On the other hand, it is well-known that the thermalization of the produced particles is a very fast but not instantaneous process [100] (see also the Section III of [32]). In fact, following the work of Spokoiny [14] (see also [15]), we can consider the thermalization rate Γ th = n χ,kin σ 2→2 , where the most important process for kinetic equilibrium is 2 → 2 scattering 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 [100]), where, as usual, α 2 ∼ 10 −3 [14].
Given that where for many models, one finds [18] 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. Thus, we get that The relativistic fluid reaches the thermal equilibrium when H(t th ) ∼ H th ∼ Γ th [14,15]; and since during kination the Hubble rate scales as a −3 , one has meaning that the temperature when the thermalization is achieved is roughly T th = 30 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 in 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 in the same order (ρ which implies, for the typical value H kin ∼ 10 −6 M pl , a reheating temperature of roughly which for ξ − 1 6 ∼ 10 −2 and g reh = 106.75 leads to a low reheating temperature of around 10 3 GeV.

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 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 [15] for a detailed discussion). Then, in the model, the field will dominate at late times when where we have assumed that the current value of the Hubble parameter is H 0 ∼ 10 −61 M pl .

Remark 3.
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 4.
It is well-known that at temperatures in the order of 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 in 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 assume as a natural initial condition that the quantum χ-field is in the vacuum at the beginning of inflation. We also chose a 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 in the same order as those of the vector mesons responsible for transforming quarks into leptons in simple theories with SU(5) symmetry [101]), because as we will immediately see, the polarization terms will be sub-dominant and do not affect the dynamics of the inflaton field. Thus, we chose 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 [102] of the Equation (187) where n is the order of the approximation used 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 [103] Before the transition time, namely, τ kin , which is very close to the beginning of kination (ϕ ∼ = −M pl ), the vacuum is depicted by χ WKB 1,k (τ), but after the phase transition, this mode becomes a mix of positive and negative frequencies of the form α k χ WKB 1,k (τ) + β k (χ WKB 1,k ) * (τ). Then, the β-Bogoliubov coefficient could be obtained by matching both expressions at τ = τ kin , obtaining where is the Wronskian of the functions f and g at the transition time, and F(τ ± kin ) = lim τ→τ ± kin F(τ). The square modulus of the β-Bogoliubov coefficient will be given approximately by [28] becausë where we have assumed that at the transition time all the energy at the end of inflation, which is approximately 1 2 m 2 M 2 pl because ϕ END = − √ 2M pl , was converted into kinetic energy.
Thus, for our model, the number density of the produced particles and their energy density will be n χ (t) ∼ , 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, which for m ∼ 10 13 GeV, m χ ∼ 10 15 GeV and H kin ∼ 10 12 GeV becomes 10 −6 < Γ/GeV < 10 12 .
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. 2.
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 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 must have been before Big Bang nucleosynthesis, which theoretically occurred at around 1 MeV.

The Diagonalization Method
In the last subsection we 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).
Hence, here we 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 with 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 [104] where in order to obtain a finite energy density [59] we have subtracted the energy density of the zero-point oscillations of the vacuum 2 ω k (τ). Following the method developed in [105] (see also Section 9.2 of [59]), we 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, by 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 consider our toy model (184) here, 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 [106] 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: 3M 2 pl H 2 . The same happens with 1 8ω k (τ) ω k (τ)/ω 2 k (τ) , which leads to a term of order H 4 , meaning that H 4 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 and vacuum polarization effects, and they rapidly disappear in the adiabatic regime-that is, soon after the beginning of kination |β k (τ)| 2 approaches |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 better understand these results it is useful to recall, as we have already explained, that the authors of the diagonalization method assumed that during the whole evolution of the universe, quanta named quasi-particles were created and annihilated due to the interactions of the quantum field with gravity [59]-i.e., a vacuum polarization effect. Furthermore, following this interpretation, the number density of the created quasiparticles 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 , so we have ρ χ (τ) = m χ N χ (τ) and the decay follows a −3 (τ), which justifies the interpretation of massive particle production.
Thus, once we have explained the gravitational particle production, we review our last work [107] 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, namely,φ 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 10 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.  (209), with initial conditions α k (τ * ) = 1 and β k (τ * ) = 0 at the horizon crossing (there were neither particles nor polarization effects at that moment because in the slow-roll regime the derivatives of the Hubble rate are negligible compared with the powers of H; i.e., the system was in the adiabatic regime).
In order to get rid of complex exponentials it is useful to transform Equation (209) into a second order differential equation, namely, Given that α k (τ * ) = 1 and β k (τ * ) = 0 lead to α k (τ * ) = 0, one should be interested in solving the equation for α k (τ), which can be split into the real and imaginary forms 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 Then, by introducing these values of the β-Bogoliubov coefficient to the Equation (210), one obtains a vacuum energy density in the order of 10 44 GeV 4 . Thus, the energy density of the produced particles evolves as 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 subleading 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 that of the relativistic plasma, when the decay is finished, i.e., when Γ ∼ H dec =H ā By imposing that the end of the decay precedes the end of kination, meaning ρ χ,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 . Thus, the following bound for the decay rate must be obtained.
Finally, the reheating temperature, i.e., the temperature of the universe when the relativistic plasma in thermal equilibrium started to dominate, which happened when ρ ϕ,reh ∼ ρ χ,reh , can be calculated as follows: Since after the decay, the evolution of the respective energy densities is given by , and thus the reheating temperature will be T reh = 30ρ χ,reh where g reh = 106.75 is the effective number of degrees of freedom for the standard model. Thus, 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 and ρ ϕ,end =ρ ϕ ā a end where it is used that the kination ends when ρ χ,end ∼ ρ ϕ,end , meaning ā a end 3 =ρ χ ρ ϕ . Thus, 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 started to dominate) was obtained when all the superheavy particles decayed-i.e., when H ∼ Γ, obtaining T reh = 30 where one has to take into account that, after the end of the kination regime, the energy density of the produced particles dominated the one of the inflaton field. Consequently, since the BBN epoch occurred in 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 then interact, decaying into light ones, becoming a relativistic plasma whose energy density eventually dominates that of the background and as a consequence the universe thermalizes and reheats.
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 with m χ once again being 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 [108], which can be analytically solved. Effectively, the β k -Bogoliubov coefficient is related to the reflection coefficient via the formula (see [109][110][111][112][113] 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 [111] Then, the number density of produced particles is given by in terms of the cosmic time t [98], namely, Analogously, the energy density of the produced particles is given by [98] ρ 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 take m χ = 0, and thus, the effective mass is given by m e f f = g|ϕ|.

2.
For masses-the effective mass of the field χ is g|ϕ|-greater than the Hubble parameterthe 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 [114]. 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 does not affect the last stages of inflation. Therefore, since H ≤

3.
The energy density of the produced χ-particles cannot dominate before its decay into light particles, which forms the relativistic plasma, because, if so, the force driving the inflation back to ϕ = 0 cannot disappear and the inflaton field would not continue its movement forward up to ∞. Effectively, the interaction term 1 2 g 2 ϕ 2 χ 2 entails that after the phase transition, the inflaton field satisfies the equation When the energy density of the χ-particles is sub-dominant, the right-hand side of (239) is negligible and the field rolls towards ∞, but when it is dominant, the right-hand side ceases to be negligible, meaning that the inflaton is under the action of the quadratic potential V(ϕ) = 1 2 g 2 ϕ 2 χ 2 , so the inflaton field will roll down to zero, which may produce a new inflationary phase. Therefore, to avoid this situation, we have to calculate when the energy density of the background and the field χ are of the same order, that is, To obtain these quantities, we use that for the model presented here, after the phase transition the universe enters in a kination regime and one haṡ and taking into account the Raychaudhuri equationḢ = −φ 2 2M 2 pl , we get Now, using that ρ χ (t) ∼ = gϕ(t)n χ (t), one gets Then, both quantities are of the same order when t ∼t ≡ 10 2 Furthermore, we obtain an important constraint for this theory: the decaying time, i.e., when the particles have decayed into a relativistic plasma, must be smaller than t, in order for the back-reaction to be subdominant so that the inflaton field rolls monotonically towards ∞.
On the other hand, by assuming as usual that there is no substantial drop of energy between the end of inflation and the phase transition time, and using that the value of the power spectrum of the curvature perturbation when the pivot scale leaves the Hubble radius is given by [67] We have used that for our model one has * = 2M 2 pl ϕ 2 * ∼ = 1−n s 4 , n s being the spectral index. Then, since the recent observational data constrain the value of the spectral index to be n s = 0.968 ± 0.006 [9], taking its central value, one gets m ∼ 5 × 10 −6 M pl , and as a consequence, H kin ∼ H end ∼ mϕ end √ 6M pl ∼ 3 × 10 −6 M pl . Thus, for values of g ≤ 10 −2 , one getst ≥ 10 10 M −1 pl and ϕ(t) ∼ M pl , meaning that for times t ∈ [10 6 M −1 pl ,t] the value of the inflaton field remains close to M pl and the effective mass of the χ-field will approximately be gM pl .
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 us 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 Γ, giving that The values of the parameters g and Γ must satisfy the constraint 10 −6 g < , and by choosing, for instance, g ∼ 10 −4 and Γ ∼ 10 −12 M pl one obtains a reheating temperature around 10 9 GeV. 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 take into account that for the authors of that paper, contrary to our convention, the reheating temperature was considered the temperature of the universe when the curvaton field had totally decayed in relativistic particles. Recall that the universe was really reheated when the energy density of inflaton was of the same order as that 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 slowroll regime because the condition m σ H end ⇐⇒ V σσ H 2 end means that the curvaton potential is flat enough at the end of inflation [115]. 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 [115] (see also Section 5.4.1 of [116] for a detailed discussion of the quadratic potential). Then, 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 .
For these values, the equation (253) becomes which means that in order to get the nucleosynthesis bounds, one has to choose light masses of the curvaton, namely, m σ ≤ 10 −14 M pl ∼ 2 × 10 4 GeV.

2.
The curvaton decays when the curvaton field dominates the universe. Assuming once again instantaneous thermalization, since the curvaton decays when it dominates, that is, when ρ ϕ,dec ≤ ρ σ,dec , the reheating time will occur during the decay (H reh = H dec ∼ Γ).
In this case the condition together with H reh < H osc , leads to the constraints Next, in this case, the reheating temperature is T reh ∼ g −1/4 reh ρ 1/4 σ,reh ∼ M pl Γ and the constraint leads to the bound On the other hand, when the curvaton decays after its domination, the power spectrum of the curvature perturbation is given by [115] P ζ ∼ = 1 9π 2 H 2 * σ 2 * ∼ 2 × 10 −9 =⇒ H * |σ * | ∼ 4 × 10 −4 M pl .

Overproduction of GWs
It is well-known that 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 than 10 −2 , which is essential to prevent any complication in the BBN's success [15]. In fact, the minimum value of this ratio, which is seen 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 [15] in a minimal GUT theory N s = 4, which corresponds to the electroweak Higgs doublet, and only in supersymmetric theories is N s in the order of 10 2 (see the explanation given by Peebles-Vilenkin in [15]).
However, as we show in this subsection, for our LQI model, when the reheating is due to the gravitational creation of heavy particles, the overproduction of GWs does not modify the BBN.
Thus, 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 [119] ρ 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 [51] ρ GW,reh ρ χ,reh ≤ 10 −2 .
First of all we see that the constraint (262) is surpassed when the decay of the superheavy massive particles is prior to the end of kination. Effectively, if the decay occurs before the end of kination, one has ρ GW,reh ρ χ,reh because the light relativistic particles evolve as the GWs. 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 surpassed 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 3/2 =⇒ a end a dec where we have used the bound (258), so the constraint (262) is clearly surpassed.

BBN Constraints from the Logarithmic Spectrum of GWs
It is well-known that during inflation GWs are produced (known as primordial GWs, in short PGWs) and in the post-inflationary period, i.e., during kination, the logarithmic spectrum of GWs, namely, Ω GW defined as Ω GW ≡ 1 ρ c dρ GW (k) d ln k (where ρ GW (k) is the energy density spectrum of the produced GWs; ρ c = 3H 2 0 M 2 pl , where H 0 is the present value of the Hubble parameter, the so-called critical density) scales as k 2 [120], 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 and k end are the momenta associated with the horizon scale at the BBN and at the end of inflation, respectively. As has been shown in [119], 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 [120][121][122][123] where the graviton spectra in quintessential models have been reassessed, in a model-independent way, using numerical techniques) .
where h 2 GW = 1 8π H kin M pl 2 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 [119,120] for more details). As has been derived in [119], the specific form of the expression above comes from the behavior of the Hankel functions for small arguments. Now, by 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.
the bound is completely surpassed.

Conclusions
In the present review we have dealt with some important quintessential inflation scenarios, namely, the original Peebles-Vilenkin one-the first quintessential inflation model introduced at the end of the 90s-and some of its improved versions; exponential scenarios; the Lorentzian quintessential inflation; and also quintessential inflation in the context of α-attractors. Studying these models, first of all we checked that for early times they provide results-spectral indices of scalar perturbations, ratios of tensor to scalar perturbations and numbers of e-folds from the horizon crossing to the end of inflationagreeing with those of the standard inflation and matching with the current observational data. Next, the models should contain a phase transition from the end of inflation to the beginning of kination where the adiabatic evolution of the universe must be broken in order to have enough particle production to reheat the universe after its cooling produced during inflation. Once we checked that the models contain this phase transition, we chose the reheating mechanism, in our case the gravitational particle production of superheavy particles, although in this review we also studied other mechanisms such as the instant preheating and curvaton reheating, and we have calculated numerically the energy density of the particles produced and the corresponding reheating temperature, whose maximum value is for the Lorentzian quintessential inflation model around 10 7 GeV.
We have also dealt with the overproduction of gravitational waves, produced as well during the phase transition from the end of inflation to the beginning of kination, showing that in the case of Lorentzian quintessential inflation-and also for the α-attractors in the context of QI-all the bounds imposed are surpassed, so they do not disturb Big Bang nucleosynthesis's success.
Finally, we found the numerical value of the other parameter (these models normally only depend on two parameters, the first one is determined using the power spectrum of scalar perturbations when the pivot scale leaves the Hubble radius), while imposing that the model matched with the current observational data at the present time. To do it, we solved numerically the corresponding dynamical system while imposing initial conditions, which were determined from the values of the inflaton field and its first derivative during the slow-roll period, from the beginning of the radiation era. In this way, we checked numerically that, for late times, the models enter in a dark energy regime (via quintessence) able to explain the current observational data.