Palatini $R^2$ Quintessential Inflation

We construct a model of quintessential inflation in Palatini $R^2$ gravity employing a scalar field with a simple exponential potential and coupled to gravity with a running non-minimal coupling. At early times, the field acts as the inflaton, while later on it becomes the current dark energy. Combining the scalar sector with an ideal fluid, we study the cosmological evolution of the model from inflation all the way to dark energy domination. We interpret the results in the Einstein frame, where a coupling emerges between the fluid and the field, feeding energy from the former to the latter during the matter-dominated era. We perform a numerical scan over the parameter space and find points that align with observations for both the inflationary CMB data and the late-time behaviour. The final dark energy density emerges from an interplay between the model parameters, without requiring the extreme fine-tuning of the cosmological constant in $\Lambda$CDM.


Introduction
The content and evolution of the Universe are very well encoded and described by the concordance cosmological model, also called Λ Cold Dark Matter (ΛCDM). Its three major components, namely, ordinary matter (including photons and neutrinos), non-relativistic or cold dark matter (CDM), and the cosmological constant Λ are captured by 6 independent parameters which completely specify the model. In addition, the model assumes two phases of accelerated expansion which are taking place in the very early and late Universe: Inflation and Dark Energy.
The inflationary epoch [1][2][3][4][5][6][7] was originally proposed as a solution to the shortcomings of the hot Big Bang model. In a single stroke, inflation is able to generate a flat, homogeneous and isotropic Universe without any topological defects. Moreover, it provides a natural mechanism where vacuum quantum fluctuations of the gravitational and matter fields get amplified to cosmological perturbations [8][9][10][11][12][13], which later became the seeds for the cosmic microwave background (CMB) primordial anisotropy and the large-scale structure of the Universe. In its simplest realization, inflation is described by a scalar field which is minimally coupled to gravity, has a canonical kinetic term and is governed by a potential whose energy density drove the (quasi-)exponential expansion. However, recent observations by the Planck [14] and BICEP/Keck [15] collaborations have essentially ruled out many simple models like monomial chaotic inflation or natural inflation and have prompted the exploration of more complicated models. The constraint on the tensor-to-scalar ratio (r < 0.036) [15] suggests that the inflaton potential has to be sufficiently flat at large field values. An easy way to achieve the flattening of the potential is to either couple the inflaton field non-minimally to gravity or to add a quadratic curvature term in the action. Both of these terms are generated from quantum corrections so it is natural to include them.
Modeling inflation in modified gravity is further facilitated in the Palatini formalism. The Palatini formulation of gravity [16,17] has recently gained considerable popularity as an alternative to the usual metric formulation. It treats the metric and the connection as independent variables, which means that one has to vary the action with respect to both of them. For a minimally coupled scalar field and an action linear in R the two formulations result in the same equations of motion and the connection turns out to be the Levi-Civita one. However, when the field is non-minimally coupled to gravity  and/or quadratic or higher curvature terms are included , significant differences arise. In the case of the non-minimal coupling, the difference can be readily seen when one transforms the Jordan frame action to the Einstein frame one. Because the Riemann tensor only depends on the connection in the Palatini formalism, this means that the Ricci scalar (which is a contraction of the metric with the Riemann tensor) transforms differently under a Weyl transformation in the two formalisms. As a result, the scalar picks up an extra coefficient in its kinetic term which is absent in the Palatini version of the theory. Therefore, the field redefinition which renders the scalar field canonical is different and the resulting Einstein frame potential is usually flatter in the Palatini formulation. Similarly, when an αR 2 term is added to the action, the auxiliary field which is usually introduced in order to eliminate this term turns out to be non-dynamical in the Palatini formulation, in contrast to the metric version. Consequently, while the metric theory becomes two-field and therefore complicated to analyze, in the Palatini version the auxiliary field can be eliminated through its equation of motion and the resulting action is single-field, albeit modified. The main modification concerns the inflaton potential which is divided by a factor that again renders it asymptotically flat.
ΛCDM is very simple but extremely fine-tuned. This is why alternatives have been put forward, one of the most prominent of which is quintessence; the fifth element after baryons, CDM, photons and neutrinos. Similar to the inflaton, quintessence is also a scalar field [88] (for a review see [89]). Therefore, unlike the cosmological constant, quintessence corresponds to a dynamic degree of freedom. Its evolving equation of state is different from that of the other constituents of the Universe content at present (baryons, CDM, neutrinos, and photons). Depending on the ratio of its kinetic and potential energy it can be either attractive or repulsive. As in slow-roll inflation, a slowly varying quintessence field can lead to the current accelerated expansion of the Universe. The energy scale of the potential energy density of quintessence must be of the order of (10 −12 GeV) 4 today, which is more than a hundred orders of magnitude lower than the typical scale of inflation. Moreover, being dynamical, quintessence requires determination of its initial conditions, such that it attains the desired energy density at present. This is called the coincidence requirement.
One way to account for coincidence, is to connect quintessence with inflation so that the initial conditions of quintessence are determined by the inflationary attractor. Indeed, it is economic to combine the inflationary and quintessential epochs and describe them by a single scalar field in the context of a common theoretical framework. This idea has been dubbed quintessential inflation [90]. 1 Of course, for such a model to be consistent, the scalar field should not interfere with the thermal history of the Universe. It should be "invisible" for much of its evolution after inflation and only start dominating around the present epoch and precipitate the late-time acceleration. Naturally, the construction of such a model is a very challenging endeavour, because it has to explain simultaneously the observations of both inflation and dark energy, but many successful models do exist [79, (for recent reviews see Refs. [154,155]).
First of all, one needs to construct an inflationary phase with a successful exit. The model must bridge the enormous gap of energy density between inflation and dark energy, in a way that introduces as little fine-tuning as possible, for this is the main motivation for quintessence over ΛCDM. Indeed, a successful quintessential inflation model needs to rely on realistic theoretical foundations. Additionally, in quintessential inflation an alternative reheating method is needed since the scalar field must survive until late times to become quintessence. Therefore, conventional reheating through inflaton decay is not applicable. Fortunately, a number of successful reheating mechanisms exist, such as instant preheating [156,157], curvaton reheating [158,159] or Ricci reheating [160][161][162].
In the post-inflationary era till the present epoch, the scalar potential should be steep, allowing the radiation domination to commence, followed by the thermal history as envisaged by hot big bang. The steep potential is necessary for sending the field into hiding after the end of inflation. In particular, the post-inflationary dynamics is characterized by a field that evolves in the kinetic regime for some time [101,102], but it then overshoots the background and gets frozen due to Hubble damping. As the background energy density redshifts and becomes comparable to the field energy density, the field resumes its evolution. Today, in most scenarios, the scalar field slow-rolls down its flat potential while dominating the Universe again and driving the current accelerated expansion. However, if we consider interaction of the scalar field with matter at present, the dark energy period is more complicated.
In this paper, we study a model of quintessential inflation in the context of R 2 Palatini gravity where the scalar field has a running non-minimal coupling to gravity. Employing Palatini gravity to study quintessential inflation was first considered in Ref. [79], considering a variation of the original quintessential inflation model in Ref. [90]. This toy-model investigation demonstrated that modeling quintessential inflation with Palatini gravity is promising. In this, much more elaborated and realistic work, we consider a simple negative exponential potential in the Jordan frame. When we transform the theory to the Einstein frame, the potential becomes flat for both negative and positive field values with a steep transition region in-between, resembling a step function. The two flat regions are suitable for inflation and quintessence. Working in the Palatini formulation allows us to modify the inflationary plateau, in particular, through the R 2 term. The running non-minimal coupling allows us to obtain the correct quintessence behaviour. To study the full time evolution of the system throughout its cosmic history, we provide the equations of motion of the scalar field and an ideal fluid component representing other matter sources in the universe. We solve these equations numerically and scan over the parameter space, finding working scenarios matching both the CMB and late-time observations for parameter values that are free of fine-tuning. A preliminary study of the model can be seen in Ref. [87]; here, our treatment and findings are more complete and comprehensive.
The paper is structured as follows. In the next section, we describe our model and perform the Jordan to Einstein frame transformation. Then, in Sec. 3, we describe the model's time evolution in a cosmological setup. We employ the slow-roll approximation and discuss the inflationary behaviour of the model, adopt Ricci reheating as the mechanism responsible for reheating the Universe and describe its details, and outline the post-inflationary expansion history, namely, kination, radiation/matter domination, and quintessence. Numerical results for inflationary and late-universe observables are presented in Sec. 4, and we conclude in Sec. 5. Further computational details are relegated to the appendices.

Setup
In this section, we first present the action of the model in the Jordan frame. After a frame transformation we bring the action to its Einstein frame form. Then, we compute the equations of motion in both Jordan and Einstein frames and show how one can easily transition between them.

The model
We start by considering the action in the Palatini formalism where m P is the reduced Plank mass, ψ collectively represents the matter fields other than the inflaton ϕ, and we take them to behave as an ideal fluid 2 . The function F (ϕ, R) takes the form We let the non-minimal coupling ξ run as with ξ * > 0 and β < 0 constants, and µ an arbitrary reference scale. In the Palatini formalism, the connection Γ is independent of the metric g µν . The connection features in the Ricci tensor, which is a function of the connection Γ only, with The form of the connection is determined by constraint equations obtained by varying the action with respect to Γ, and, in the presence of the non-minimal gravitational physics introduced by the non-zero ξ and α, it will differ from the standard Levi-Civita form. The real scalar field ϕ, which plays the role of the inflaton and quintessence in quintessential inflation, is governed by an exponential potential The exponential form is well-motivated in particle physics (it usually appears in string theory and supergravity models, e.g. in gaugino condensation [163][164][165]). It can produce quintessence in agreement with observations in its flat tail at ϕ > 0, and it is also suitable for quintessence from a theoretical point of view: we do not introduce a fine-tuned cosmological constant by hand, but instead V → 0 for large ϕ, and the late time dark energy density arises dynamically from the equations of motion. The action (2.1) is dynamically equivalent (as long as ∂ 2 as can be seen by obtaining the equation of motion for the auxiliary field χ and plugging it back in Eq. (2.6). Using this, the action can be cast in the form (2.7) As is standard, we employ a conformal transformation (note that, in the Palatini formalism, this does not change Γ) to express the action in the Einstein frame where the gravitational part takes the standard Einstein-Hilbert form: (2.9) Note that, essentially, Ω 2 = ∂ R F (ϕ, R). We have introduced the short-hand notation (∂ϕ) 2 ≡ g µν∂ µ ϕ∂ ν ϕ, where∂ denotes a derivative with respect to the Einstein frame coordinates. Throughout the paper, we will use an overbar to denote Einstein frame quantities. Due to the standard form of the gravity sector, we will interpret all the usual cosmological observations in the Einstein frame.
To make the calculations that follow less cluttered, we define h(ϕ) ≡ m 2 P + ξϕ 2 . (2.10) We then get rid of the auxiliary field by obtaining its equation of motion. Let us, for a moment, ignore all matter except for the inflaton; then, we have giving (2.12) Plugging both expressions back into the action gives [63,64] Note that, because we were able to get rid of the non-dynamical auxiliary field through its equation of motion, the above action contains only one scalar field. This is in contrast to the metric version of the theory, where the auxiliary field is dynamical and the Einstein frame action contains two fields. The field can be made canonical via the redefinition (2.14) Note that for large negative ϕ, this gives dφ/dϕ ∝ e κϕ/(2m P ) , which, after integration, shows that φ approaches a constant as ϕ → −∞. We choose this constant to be equal to zero, so that φ is restricted to take positive values. The field redefinition leads finally to Note the appearance of the higher-order kinetic terms for the scalar. As we will see below, they are negligible for most of cosmological evolution. Note also the form of the Einstein frame potential,V which chiefly determines the cosmological evolution of the model. An example case is depicted in Fig. 1. The appealing features of the model are evident in the potential. For ϕ < 0, the potential decreases with increasing ϕ, but only slowly: the α term makes the potential flat and suitable for slow-roll inflation. For ϕ > 0, the α term is subleading, and the potential decreases quasi-exponentially, modified by the change of variables (2.14). The ξ contribution modifies the potential; its running enables it to fix both the inflationary CMB observables and the late-time dark energy to values that match observations. For large enough ϕ, ξ runs to negative values, causingV to first flatten and then start growing, forming a local minimum and a nearby peak when 1 + ξ(ϕ)ϕ 2 /m 2 P becomes zero. For the parameters in Fig. 1, the zero occurs at ϕ = 890.99m P , and at this point the height of the Einstein frame potential isV (890.99) = 1.14 × 10 −94 m 4 P (notice the second term in the denominator in the potential regularizes the peak). Beyond the peak, the kinetic term in (2.13) changes sign. In practice, as we will see below, dynamics never probe this region.

Equations of motion in the Jordan frame
While the Einstein frame discussed in the previous section is useful for physical interpretation of the results, the equations of motion are easier to formulate in the Jordan frame, especially when we wish to include the non-inflaton matter contribution from (2.1) and thus go beyond the simplified Einstein frame action (2.15). To obtain the equivalent of Einstein equations for our system (not a trivial task in the Palatini formulation with non-minimal gravity), we vary the action (2.1) with respect to the metric g µν and the connection Γ. The variation of the F (ϕ, R) term reads where F R (ϕ, R) ≡ ∂ R F (ϕ, R) and parentheses around indices indicate the symmetric part of a tensor. Meanwhile, the variation of the matter part gives the energy-momentum tensor, as usual: In the Jordan frame, the contributions from the field ϕ and other matter components are independent from each other (indeed, this is the reason we work in the Jordan frame to begin with). We take the matter energy-momentum tensor to be of the ideal fluid form with energy density ρ and pressure p, T (m) µν = (ρ + p)u µ u ν + pg µν . We define the fluid's barotropic parameter as w ≡ p/ρ. The energy-momentum tensor of the field takes the standard form T (ϕ) µν = ∂ µ ϕ∂ ν ϕ − g µν 1 2 (∂ϕ) 2 + V . All in all, the variation of the action with respect to the metric gives To make progress, we would like to express the left-hand side in terms of familiar geometric quantities. We start by taking the trace of (2.19), and the R-derivative of the F function reads (2.22) The trace of the energy-momentum tensor is Plugging these into (2.19), we now have a relationship between R (µν) and the matter sources. However, as explained around (2.4), R µν is defined through the Palatini connection Γ, which may differ from the Levi-Civita one, so we do not know how it relates to metric quantities such as the scale factor and the Hubble parameter. Varying the action with respect to the connection, we see (assuming symmetricity in the lower indices of the connection) that where∇ is the covariant derivative defined with Γ. After integrating by parts in the action and some manipulations (relabelling indices and tracing over λ and ν), one obtains This means that Γ is actually the Levi-Civita connection ofḡ µν = (∂ R F )g µν , that is, Using this together with the Einstein equation (2.19), we can derive an expression for the 'metric' Einstein tensor defined by the Levi-Civita connection of the Jordan frame metric [166]: Here the covariant derivatives ∇ are taken in terms of the Levi-Civita connection of the metric g µν . In other words, we have eliminated the a priori independent connection Γ, which turns out to be just an auxiliary field. We see that the energy-momentum tensor is modified by rescalings and new effective matter sources. Adopting now the flat FLRW metric, the metric Einstein tensor G µν can be written in terms of the scale factor a and its time derivatives such as the Hubble parameter H ≡ȧ/a in a standard way. Dot refers to a derivative with respect to the cosmic time. The zeroth-zeroth component of Eq. (2.27) then reads while the ij components reaḋ where R and F R are given in terms of the matter content by Eq. (2.21) and (2.22) and Eq. (2.28) is a second-order algebraic equation for H, which can be solved in terms of the field and fluid variables ϕ,φ, ρ, and w. A complication arises from the time derivatives of F R : these contain also factors such asφ andρ, which must be eliminated using the field and fluid equations introduced below. The procedure is explained in detail in Appendix A.
For the fluid, it is shown in Ref. [167] that its energy-momentum tensor is conserved in this setup, so that Finally, varying the action with respect to ϕ, we havë

Between the Jordan and Einstein frames
To give a physical interpretation for the dynamics, we want to relate the Jordan frame quantities to the Einstein frame ones. In both frames, we use a flat FLRW coordinate system with the metrics g µν = diag(−1, a 2 , a 2 , a 2 ) andḡ µν = diag(−1,ā 2 ,ā 2 ,ā 2 ). We remind the reader that we use an overbar to denote the Einstein frame quantities. The spatial coordinates of these two flat FLRW coordinate systems match, but the time coordinates are rescaled. As per our convention, we call the Jordan frame coordinates x µ = (t, x i ) and the Einstein frame coordinatesx = (t, x i ), and the conformal transformation gives for spacetime intervals where Ω 2 = F R depends on time only. We obtain the relationships as the master equations for moving between the two frames. With these, we can express various Einstein frame quantities in terms of the Jordan frame ones. In particular, where a dot still denotes a derivative with respect to the Jordan frame time, and we introduced a circle over a symbol to indicate a derivative with respect to the Einstein frame time, so that The relation between the fluid energy-momentum tensors in the Jordan and Einstein frames is [79] where we used ∂g αβ /∂ḡ µν = Ω 2 δ α µ δ β ν and √ −ḡ = Ω 4 √ −g. The Jordan frame ideal fluid is still ideal fluid in the Einstein frame; following Refs. [168,169], we write its energy-momentum tensor as where the last equations relate the Jordan and Einstein frame quantities. It follows that the barotropic parameter has the same expression in both frames: Below, we will always refer to the Einstein frame when talking of the barotropic parameter; we will omit the bar for simplicity of notation.

Equations of motion in the Einstein frame
We are now ready to examine the Einstein frame equations of motion. Their full form is complicated-in the Einstein frame action (2.9), the field and fluid components are coupled through the conformal factor Ω −2 inside S m . In a general case with α = 0, the fluid may even modify the χ constraint equation (2.11) and, as a consequence, the field transformation (2.14). We only present here approximate forms of the equations, free of some of these complications and valid during specific cosmological eras. Exact expressions can always be obtained by starting from the Jordan frame equations of Sec. 2.2 and applying the transformations of Sec. 2.3. During inflation and right after it, the fluid is subdominant and can be ignored in the field equations. Varying the action (2.15) then gives [63] with h defined in (2.10). The energy density and pressure of the field read [81] 42) and the Hubble parameter can be written as 3m 2 PH 2 =ρ φ . The higher-order kinetic terms are the only complication compared to a standard canonical scalar field.
At later times, the fluid becomes important, but the α terms turn out to be negligible. In this limit, the field transformation (2.14) can be solved explicitly to yield 3 ξϕ = m P sinh ξφ/m P , (2.43) with the Einstein frame potential The field is coupled to the fluid; action (2.15) with the fluid contribution added in gives in the α → 0 limit: Throughout the cosmic history, the fluid continuity equation in the Einstein frame can be obtained from the Jordan frame version (2.31), using the transformations of Sec. 2.3. The result is

Multiplying the field equation (2.45) by
• φ gives the continuity equation for the field energy density. The inflaton-fluid coupling terms there and in (2.47) are identical but have opposite signs: the coupling simply transfers energy from one component to the other. The coupling vanishes in the early universe when the fluid behaves like radiation, w = 1/3, but it can be non-negligible during matter domination. We will discuss the effects of this coupling in more detail in Sec. 4.
Note that these expressions are still written partly in terms of the Jordan frame field ϕ, hidden in quantities like h, V , and F R . In a general case, it is not possible to solve the field φ from ϕ analytically. This is why, in our practical numerical computations, we work in the Jordan frame. The Einstein frame expressions of this section are for the benefit of developing a physical intuition of the system.

Cosmic history with quintessential inflation
Let us now turn to the time evolution of our model in a cosmological setup. In this section, we explore the cosmic history qualitatively through its many stages, starting from inflation and ending with quintessence domination. To make contact with the standard formalism discussed in the literature, we mostly work in the Einstein frame.

Inflation
We start with the field at the plateau with ϕ < 0 and highV , with other matter components being negligible. We assume the high potential energy density dominates over the scalar's kinetic energy, giving rise to cosmic inflation, where the expansion of space accelerates. The plateau inV is suitable for slow-roll inflation, where the field slowly moves towards positive values so that the potential gradient is balanced by Hubble friction. In this limit, the Einstein frame equations of motion (2.41) take the standard form where we neglected higher-order kinetic terms as subleading slow-roll corrections [63]. The evolution is characterized by the slow-roll parameters: For slow roll to be possible, we must have V < 1 and |η V | < 1 at the corresponding field values. We can compute the slow-roll parameters for our potential (2.16) in the limit of constant ξ, that is, with β = 0. To make the computation simpler, we use a result from [63] that relates V and η V to their counterparts in the α = 0 limit (here˜ andη, respectively). The results, by using Eq. (2.14) and the chain rule, can be expressed terms of the Jordan frame field ϕ as 3) andṼ is defined in (2.44). The expression for V reveals possible extrema withṼ = 0 at κφ = −2 ± 4 − κ 2 /ξ. We demand that the potential is monotonic, i.e.Ṽ < 0 everywhere; this sets the restriction κ 2 > 4ξ on the allowed parameter space.
Asymptotically,˜ ∼ ϕ 2 , diverging for both positive and negative ϕ. However, V is suppressed by the exponential αṼ contribution so that V 1 for ϕ −m P /κ. This allows the system to undergo inflation even for large negative ϕ. Indeed, this was the motivation for us to introduce the αR 2 term to our model (2.1) in the first place. The asymptotic behaviour η V ∼η ∼ ϕ also reveals divergences for |ϕ| → ∞, this time not removed by the α terms, making slow roll impossible for ϕ −m P /(κξ). This leaves us with a range of field values near ϕ = 0 that are compatible with slow-roll inflation. We start our inflationary evolution in slow-roll in this field range. As we will see in Sec. 4, typical values of the model parameters support the 60 or so e-folds of inflation needed for a successful inflationary scenario.
The motivation for slow-roll inflation is that it produces a nearly scale-invariant spectrum of perturbations, compatible with the CMB observations [14,15] A s = 2.1 × 10 −9 , n s = 0.9649 ± 0.0042 , α s = −0.0045 ± 0.0067 , r < 0.036 . (3.4) Here A s is the scalar power spectrum amplitude, n s is the scalar spectral index, α s its running, and r is the tensor-to-scalar ratio at the CMB pivot scale k * = 0.05 Mpc −1 . In the slow-roll limit, the perturbations can be computed using the standard formalism (see e.g. [170]), giving where we also gave the forms based on the Hubble slow-roll parameters, where the approximations apply during slow roll. The expression for α s depends on higherorder slow-roll parameters, which we omit for brevity; these can be found in e.g. [170]. Using the results (3.3), we can also write down the full expression In our numerical results, we have β = 0, so the results (3.3), (3.7) will be modified slightly. We will use these expressions as guidance when scanning over the parameter space, but we will compute the CMB observables from the Hubble slow-roll parameters as laid out in (3.5). The modifications of inflation due to a non-zero β turn out to be minor; β is more important for the later evolution of the system, in particular, for fixing the final dark energy density.

Kination
Inflation ends when the field rolls down from the inflationary plateau to positive ϕ values. As the field drops off the potential 'cliff', see Fig. 1, its velocity increases and the kinetic terms in the action (2.15) start to dominate over the potential. During this stage, the extra kinetic terms proportional to α • φ 4 may play a role in the evolution. However, as the field velocity decreases due to Hubble friction, these terms die out quicker than the canonical • φ 2 kinetic term, which soon dominates. Analogously, in the Jordan frame, the αR 2 term becomes subdominant compared to the linear R term as the energy density of the universe, and thus its curvature, decreases, and it stays subdominant until today. Thus, the α term is only important during and right after inflation.
After a transition period (lasting less than 10 e-folds according to the numerics of Sec. 4), the scalar field follows standard kination [101,102,105,106,112,171] with the equations of motion Note that the exponentially suppressed potential does not play a role during this stage. The evolution (3.9) corresponds to a barotropic parameter w = 1, which is quite distinct from the standard radiation or cold matter domination (w = 1/3 and 0, respectively). The period of kination leads to a non-standard expansion history of the universe, which, in particular, shifts the number of e-folds of inflation left at the Hubble exit of the CMB pivot scale k = 0.05Mpc −1 from the standard 50-60 to 60-70. We will return to this point in Sec. 4, where we match the CMB scale based on the full expansion history.

Reheating
In many conventional models of inflation, reheating occurs through the inflaton decaying into matter particles, which then take over the energy density and start the standard Hot Big Bang era. In quintessential inflation 4 , the field condensate must be preserved and serve as dark energy later on. Therefore, radiation has to be created in some other way. There are many mechanisms which can facilitate this.
As an example, we consider one such mechanism, called Ricci reheating. Ricci reheating was first considered by Ref. [160]. Then, it was refined first by Ref. [161], which also coined the name, and further by Ref. [162]. In a nutshell, the idea behind Ricci reheating is as follows. The mechanism is based on the fact that, for a flat FRW Universe, the Ricci scalar (in the Einstein frame) isR = 3(1 − 3w tot )H 2 , where w tot is the barotropic parameter of the whole Universe. During slow-roll inflation we expect w tot = −1, while after the end of inflation during kination we have w tot = 1. This implies that the sign ofR changes in the transition from inflation to kination. If one considers also a spectator scalar field ψ with non-minimal coupling to gravity ∝Rψ 2 , then this change of sign inR would correspond to a change of sign in the effective mass-squared of ψ generated due to the non-minimal coupling. Assuming that this effective mass-squared is positive during inflation, we can safely consider that the expectation value of ψ is zero by the end of inflation. However, as we switch to kination, the effective mass of ψ becomes tachyonic and the field is displaced from zero (which corresponds to a potential hilltop, after inflation) and begins oscillations in its effective potential. The oscillating ψ has a particle interpretation and can decay into radiation, which eventually reheates the Universe, because its density is diluted less efficiently by the expansion than that of the free-falling inflaton during kination.
The mechanism has a number of advantages compared to other reheating mechanisms considered in quintessential inflation. It can be very efficient, in contrast to gravitational reheating [172,173], which means it would not challenge Big Bang Nucleosynthesis (BBN); it does not require a coupling between the spectator field and the quintessential inflaton in an enhanced symmetry point, as would be the case of instant preheating [156,157]; it does not need tuning of initial conditions for the spectator field, as does the curvaton reheating mechanism [158,159] and finally it does not presuppose a quintessential inflaton with dissipating properties as in warm quintessential inflation [142] or the generation of primordial black holes [174]. It only employs the fact that renormalisation in curved spacetime results generically in a non-minimal coupling of scalar fields to gravity.
The additional Lagrangian density of the scalar field is where V (ψ) is the part of the scalar potential which involves ψ andξ is a non-perturbative coupling, which should not be confused with ξ, the non-minimal coupling of the quintessential inflaton field.
Technically, the addition of the above in the Lagrangian density of the theory is yet another modification of gravity, which must be taken into account when switching between the Jordan and Einstein frames. However, we consider that ξ |ψ| m P always, which means that the influence of ψ on gravity remains always negligible. Thus, in effect, we can consider that the only effect of the above non-minimal coupling is to provide a contribution to the effective mass-squared of the spectator field. Additionally, the condition ξ |ψ| m P allows us to consider a perturbative scalar potential, which around the expectation value of the field during inflation, can be written as where the ellipsis denotes higher-order non-renormalisable terms, presumed negligible. We will consider at first that the non-minimal coupling overwhelms the bare effective masssquared |m 2 | |ξR| so we can ignore the first term on the right-hand-side above. This sets a limit on the mass which we discuss in Appendix B. We will also consider a positive perturbative self-coupling 0 < λ < 1, so that the potential is stabilised by the quartic term and not by non-renormalisable terms, although a modification of our results in the latter case is straightforward.
In Ref. [162] it was shown that after the end of inflation, the field ψ oscillates as determined by the terms in (3.11) that stabilize V (ψ), while the effect of the central potential hill (generated by the non-minimal coupling) is diminishing (and negligible) becauseR ∼H(t) is decreasing after inflation. If the stabilising potential is quartic, as is the case of Eq. (3.11), then the density of the oscillating condensate decays as radiation,ρ ψ ∝ā −4 [175]. However, if the potential were stabilised by a non-renormalisable term, this would not have been so. Fortunately, Ref. [162] demonstrated that it is largely irrelevant which term stabilises the potential V (ψ). This is because in Ref. [162] it was shown that the primary reheating effect is not the perturbative decay of the coherently oscillating ψ condensate, but the non-perturbative particle production on the hilltop, right after the end of inflation. At this moment, the field finds itself on top of a potential hill, leading to ample production of radiation due to a tachyon instability. In Ref. [162], it was claimed that the produced radiation dominates over the one corresponding to the oscillating condensate. Because the latter is diluted (at least) as fast as radiation, it never becomes important, at least as long as the quadratic term in Eq. (3.11) remains negligible. These considerations simplify our treatment, because they suggest that radiation is immediately produced at the end of inflation, and the further evolution of the oscillating ψ condensate is irrelevant. The only question is how much radiation is produced.
An estimate of the size of the spread of a scalar field condensate on top of a potential hill is given by ψ 2 |m eff | 2 [176], where the effective mass squared in our case is m 2 eff = −6ξH 2 during kination, which takes place near the end of inflation. Therefore, the density of radiation at the end of inflation isρ . (3.13) During kination, the total density of the Universe decreases asρ tot ∝ā −6 , while for radiation we haveρ r ∝ā −4 , which means thatρ r /ρ tot ∝ā 2 . Therefore, where 'reh' denotes reheating, which is the moment that radiation takes over and we havē ρ r ρ tot . The density of the Universe at reheating is straightforward to find, by considering thatρ tot ∝ā −6 . Indeed, we get where we used Eq. (3.14) andρ end tot = 3H 2 end m 2 P . Therefore, using that at reheatingρ tot ρ r = π 2 30 g * T 4 , the reheating temperature is where g * is the number of effective relativistic degrees of freedom at reheating. The allowed reheating efficiency for successful reheating is The lower bound in the range of the reheating efficiency in Eq. (3.17) is obtained from gravitational reheating [172,173], which challenges the process of BBN due to an over-enhancement of primordial gravitational waves during kination. Indeed, gravitational reheating suggests ρ gr r ∼ 10 −2H 4 end , whereH end ∼ 10 −8 m P and we used thatV This range includes values ofξ ∼ 1, which means that no fine-tuning is required for our mechanism to work. A similar lower bound onξ is obtained when considering the density of the primordial gravitational waves generated by inflation. Indeed, the density of the gravitational radiation at the end of inflation isρ end GW 1 4π 2H 4 end (see appendix C). We require thatρ r /ρ GW > 20 at BBN, but because bothρ r andρ GW decrease with the expansion asā −4 , we have the same requirement at the end of inflation. In view of Eq. (3.12), this requirement becomesξ > 1/6 √ 2 π 0.038, which agrees with the range in Eq. (3.18). Equations (3.16) and (3.18) suggest that the reheating temperature ranges as 10 GeV T reh < 10 14 GeV . (3.19)

Radiation and matter domination
After reheating, the universe is dominated by hot radiation, and the barotropic parameter settles to w = 1/3. As the univserse cools, particles in the thermal bath start to become non-relativistic, and this cold matter eventually takes over. We approximate this to happen instantaneously whenρ r 10 −110 m 4 P , corresponding to a temperature of ∼ 0.8 eV [170]. At the same time, the field follows the equation of motion (2.45), veering away from kination once radiation starts to take over. While the fluid is relativistic, w = 1/3, the field and fluid don't mix directly. However, in the presence of radiation the Hubble parameter is larger than it would be if induced by φ alone, and this increases the importance of the friction term. The field velocity • φ starts to decrease dramatically, until the field essentially freezes to a near-constant value φ fr . Using the known scalings of the scalar and radiation energy densities, and assuming a negligible scalar potential, we can writē where 'kin' refers to a moment at the beginning of standard kination withρ kin r /ρ kin φ ≡ Ω kin r 1. With this and 3H 2 m 2 P =ρ tot , we can solve the frozen field value as ln Ω kin r . (3.21) As the kinetic energy of the field drops, the potential again starts to play an important role in field evolution, complicating the dynamics. Two basic behaviours emerge: the field may completely freeze, so that its potential energy comes to dominate over the kinetic one and the field's barotropic parameter becomes −1, or the field may start to follow a scaling attractor with slow time evolution [177][178][179]. To estimate which fate is more likely, we can approximate the potential locally around φ = φ 0 with the exponential (3.22) If κ eff is approximately a constant, then κ eff < √ 2 leads to freezing, and κ eff > √ 2 gives the scaling solution. In our model in the examples below, we find κ eff to be small and slowly changing, leading indeed to a freezing behaviour.
After matter becomes non-relativistic with w = 1/3, time evolution is further complicated by the direct coupling between the fluid and the field in (2.45), (2.47). In practice, the dynamics have to be solved numerically; we do this in Sec. 4.

Quintessence domination
As the field rolls, ξ from (2.3) runs to smaller and smaller values, and the Einstein frame potential (2.16) becomes flatter and flatter, becoming more suitable for quintessence with a slowly rolling field. Indeed, as mentioned in Sec. 2.1, eventually ξ runs to negative values; around this point, the Einstein frame potential develops a local minimum and then starts to grow again, with a high positive peak near ξ(ϕ)ϕ 2 = −m 2 P . The coupling to matter can cause the field to overshoot the minimum and oscillate around it a few times, but eventually, as the fluid energy density dilutes away, the field will settle into the potential minimum at ϕ ≡ ϕ fin . Its barotropic parameter w φ = −1 and its energy density, given by the height of the potential, become constant. The quintessence field then behaves as dark energy. To match observations, we needV (ϕ fin ) = 7.23 × 10 −121 m 4 P , computed assuming thatH = 67.66km/s/Mpc and that roughly 70% of the energy density of the universe today is in dark energy. To be more precise, the dark energy fraction today is [180] Ω φ = Ω DE = 0.6889 ± 0.0056 . (3.23) Since we live in the transition period where both dark energy and matter have non-negligible roles, the quintessence field is not necessarily completely frozen yet. In our numerical results, we demand that the barotropic parameter of the field today respects the observational bounds of the CPL parametrisation [181], where '0' refers to today, and the limits are

Numerical results
In this section, we explain the details concerning the numerical side of our work. As it is explained above, in order to numerically solve the dynamics of the system, we work in the Jordan frame. It is then straightforward to obtain the corresponding quantities in the Einstein frame, where our intuition applies, by following the discussion in Sec. 2.3. To be more explicit, we need to solve for the scale factor a(t), the inflaton field ϕ(t) and the fluid density ρ(t) (remember that, at a classical level, homogeneity and isotropy impose that the fields depend on time only), since every other quantity depends on these. In principle this could be done by solving the system of ordinary differential equations given by Eqs. (2.29), (2.31), and (2.33). However, the Hubble factor can be algebraically solved to be where the specific forms of A and B, as well as the details of the calculation, can be found in Appendix A. For our current discussion it suffices to know that A and B depend on ϕ(t) (and its first derivative) and ρ(t) only. This means that the initial system of ordinary differential equations given by Eqs. It is also worth commenting on our choice of the Jordan frame over the Einstein frame. One obvious advantage of working in the Einstein frame is that the gravitational sector of the action is simply the Einstein-Hilbert term. However, were we to work in the Einstein frame, Eq. (2.14) would need to be solved and inverted in order to obtain ϕ(φ), to then be plugged back in the action (2.13) in order to express all the quantities in terms of the canonical Einstein frame field φ. Furthermore, the action in the Einstein frame features a quartic kinetic term and a coupling between the inflaton and the background matter fields through a conformal factor in the matter action. Although during inflation the matter action is zero, during the subsequent cosmological eras this extra coupling is present, complicating the setup. Likewise, the quartic kinetic term, which complicates the equations of motion even further, cannot be a priori discarded (although after solving the dynamics it is found to be in general negligible, see Fig. 8). All of these considerations outweigh the only hurdle in the Jordan frame: gravity is non-linear. As a matter of fact, due to working in the Palatini formalism, we can profit from further simplifications as the one explained above, where the Hubble factor can be algebraically solved in terms of the inflaton and the background fields. In this way, we find the solution of the system to be much more approachable in the Jordan frame than in the Einstein frame. Finally, it is important to keep in mind that, as we have mentioned, once the dynamics is solved in the Jordan frame it is straightforward to obtain the analogous quantities in the Einstein frame by following the discussion in Sec. 2.3.

Initial conditions
During the inflationary era the only existing field is the inflaton (even if some matter fields existed they would be inflated away), so that ρ(t) = 0. Therefore, the only equation to solve for is Eq. (2.33) (with H given by Eq. (4.1)), which is a second order ordinary differential equation. Thus, only two initial conditions are needed, ϕ(t i ) andφ(t i ). We choose ϕ(t i ) sufficiently negative to capture all the possible evolution histories when scanning over the parameter space, while respecting the bound that imposes that the field should not be much smaller than −m P /(κξ) (c.f. Sec. 3.1), for which slow-roll is not possible. This usually amounts to having ϕ(t i ) ∼ −30 m P and as it can be seen from Fig. 3 (see also the discussion in Sec. 4.3), using a smaller value would be of no help, since the region of the parameter space compatible with observations restricts ϕ(t i ) > −30 m P . Furthermore, for simplicity, since the field will eventually reach the slow-roll attractor, we chooseφ(t i ) such that slow-roll is satisfied. Effectively this means neglecting the second order derivative in Eq. (2.33). With ϕ(t i ) fixed, this equation only depends onφ(t i ), for which we can (numerically) solve to obtain the initial value.
The end of inflation gives way to kination. During this era, some reheating mechanism transfers the energy density of the inflaton to the particles of the SM, which are modelled in our setup by a perfect fluid with energy density ρ(t). In this way, the last needed initial condition is the initial energy density of radiation, at the end of inflation, ρ(t end ). It can be found by the following simple calculation where * corresponds to the time at which the CMB pivot scale exits the horizon during inflation (with k * =ā * H * = 0.05 Mpc −1 ), "end" corresponds to the end of inflation, and '0' corresponds to the present time. We have also setā 0 = 1 and made the approximation a ∝ T −1 from the end of inflation until today, where T is the temperature of radiation.
Using Eq. (2.36), we can relate the number of e-folds in the Einstein and Jordan frames as [182,183] Thus, where T 0 ≈ 2.7K. The initial energy density of radiation at the end of inflation can be written asρ where g * = 106.75 is the number of relativistic degrees of freedom. Relatingρ to ρ via (2.39) and gathering the above results together, we get written in terms of quantities that are either known or fixed by inflation. Note the cancellation of (F end R ) 2 due to the extra factor of F 2 R coming from expressing the energy density in the Jordan frame. It is important to mention that when scanning over the parameter space we require thatρ(t end ) satisfies two bounds, the upper one such that the inclusion of the radiation fluid at the end of inflation is a small perturbation to the overall dynamics, i.e., Ω end r < 0.1, and the lower one corresponding to the gravitational reheating limit, which is the least efficient reheating mechanism. Thus, we imposeρ(t end ) >ρ grav = qg * (H end ) 4 /(480π 2 ) 2.25×10 −2 (H end ) 4 [184], where we have introduced q ∼ 1 because the spectrum is not exactly thermal.

The parameter space
The model has six parameters, namely κ, ξ * , β, µ, α, and M 4 . It would be computationally costly to perform a scan over such a six-dimensional space. However, there are some simplifications that allow us to reduce the dimensionality of the parameter space.
The first thing to notice is the scaling law the model obeys. Indeed, let us rescale the coordinates, background density and parameters in the Jordan frame as 5 x µ → λx µ , ρ → λ −2 ρ, α → λ 2 α and M 4 → λ −2 M 4 .  Of course, the equations of motion are invariant under such a rescaling of the action. Furthermore, the quantity αM 4 is also invariant. Looking at the expressions for the inflationary observables in Eq. (3.6), one can see that the parameters α and M 4 only enter the expressions for n s and r through the combination αM 4 , i.e., they are invariant under the rescaling (4.7). It is not so for A s , where M 4 enters its expression alone. From this discussion we conclude that it is enough to scan over the quantity αM 4 . For each value of αM 4 we can fix M 4 such that A s satisfies the observational requirements from (3.4). In this way we have reduced the dimensionality of the parameter space to five.
There is one extra simplification that can be made by taking into account that µ in Eq. (2.3) is an arbitrary scale that can be changed by reabsorbing it into ξ * . Therefore it can be chosen to take the most convenient value, which, for us, is the field value at which the cosmological scales leave the horizon, ϕ * . This way, around this scale the effect of the running is minimal and the non-minimal coupling is roughly just ξ * . The dimensionality of the parameter space is now four.
Having defined the degrees of freedom of the system, i.e., ϕ(t), a(t) and ρ(t), the initial conditions, i.e., ϕ(t i ),φ(t i ), and ρ(t end ), and the parameters over which to scan, i.e., κ, αM 4 , ξ * and β, we first focus on the inflationary regime of the theory. In this way, we start with the initial conditions discussed above and numerically solve the system until the end of inflation, defined by the condition 6 H ≡ • φ 2 /(2H 2 m 2 P ) = 1. We take discrete slices in αM 4 , ranging from 0.0143 to 1.43 × 10 6 in steps of factor 10 and a region in β around the central value of −0.1 with a resolution of 10 −3 and scan over the parameters κ and log 10 ξ * with values in the intervals [0.2, 0.7] and [−2.5, −0.9], respectively, with resolutions of 5 × 10 −3 . The reason behind choosing such a central value for β is that we have found that a correct behaviour for quintessence is strongly peaked around it.
As the values for the field and its velocity at the end of inflation will serve as the initial conditions for the beginning of the next cosmological era, we impose a set of conditions on the points obtained from the scan through which we obtain the valid region in the parameter space. In addition to fixing A s = 2.1 × 10 −9 as discussed above, we require that: • The value of the scalar spectral index is equal to the central value obtained by Planck [185], i.e., n s = 0.9649.
• The value of the tensor-to-scalar ratio is within the latest observational bounds [15], i.e., r < 0.036.
• The value of the running of the scalar spectral index is within the 2σ bounds obtained by Planck [185], i.e., −0.0179 < α s < 0.0089.
• The initial energy density of radiation at the end of inflation, obtained via Eq. (4.6), amounts to a small perturbation of the system, i.e., Ω end r < 0.1.
• The initial energy density of radiation at the end of inflation is larger than the energy density corresponding to gravitational reheating, i.e.,ρ(t end ) > 2.25 × 10 −2 (H end ) 4 .
The last two conditions translate to the available range in the number of e-folds from the time at which the cosmological scales exit the horizon until the end of inflation (see the right panel in Fig. 2). It is usually between 60 and 75. Also note that we have not imposed a correct value for the amplitude of the power spectrum A s as a condition since every single point in the parameter space already satisfies this, by exploiting the scaling property of the model explained above. When inflation ends, and after imposing the above set of conditions to obtain the valid region of the parameter space, we use the final values of the field and its velocity as the initial conditions for the next cosmological era, as well as Eq. (4.6) for the radiation energy density, in order to solve Eqs. (2.31) and (2.33), with H given by Eq. (4.1). The barotropic parameter of the fluid is of course 1/3 up until the transition to the matter domination era, when it becomes w = 0. We model this transition by a jump from 1/3 to 0 in the barotropic parameter of the background at the time when the energy density of radiation is equal to its value at matter-radiation equality,ρ eq = 1.27 × 10 −110 m 4 P [185]. The simulation is finished when the energy density ratio of the field, corresponding now to dark energy, is equal to the central value obtained by Planck [185] of its value today, i.e., Ω 0 φ = 0.6889. At this point we impose another set of conditions, which we list here.
• The temperature of the universe at the onset of radiation domination is above T BBN 0.1MeV.
• The barotropic parameter of the field is within the latest observational bounds [185], i.e., w 0 φ < −0.95. • The running of the barotropic parameter of the field in the CPL parametrization is within the latest observational bounds [185], i.e., −0.55 < w 0 a < 0.03. Importantly, we also consider the bound on the density parameter of gravitational waves coming from BBN constraints, 20 Ω end GW < Ω end r , as discussed in Sec. 3.3. By using Eq. (3.13) this bound translates to an allowed range of values for the non-minimal coupling between the reheaton and gravityξ. The successful values ofξ for each point in the parameter space are presented in Table 1. The two points, for whichξ is the largest (ξ ∼ O(1)), are highlighted in black in Fig. 11 (while the rest are in red).
The points that satisfy this extra set of conditions are the successful points of our model. For them we have successful inflation, with correct inflationary predictions, as well as a correct evolution during the expansion history of the universe, with successful dark energy at the present time.

Numerical results for inflation
In this section, we present and analyze the obtained results for inflation. We remind the reader that the power spectrum strength at the pivot scale, A s , is fixed to its observed value in all our results. In the left panel in Fig. 2 we show an example slice of the parameter space in the (log 10 ξ * , κ) plane with fixed β = −0.1 and αM 4 /m 4 P = 1.43. The blue points have a correct value for n s while the orange points satisfy the full set of conditions for inflation stated in Sec. 4.2. In order to understand the shape of the parameter space let us consider the β = 0 case, for simplicity. First, we remember we have imposed the potential to be monotonic, i.e., κ 2 > 4ξ = 4ξ * . The lower boundary of the parameter space region corresponds to this requirement. The other consideration to take into account is Eq. (3.7), which, since in the β = 0 case the expressions for the slow-roll parameters in Eq. (3.3) are exact, is an exact expression for the scalar spectral index (in the slow-roll approximation). Since this is a quadratic equation in κ, it can be algebraically solved for, giving an expression depending on ξ * and ϕ * (the field at horizon exit), κ = κ(ξ * , ϕ * ). In Fig. 3 we plot in green the curve κ(ξ * )| ns=0.9649 , for many values of the field at horizon exit, ranging from −30 m P to 0 (in steps of 0.5 m P ). We can see that the upper boundary of the parameter space coincides with the asymptotic upper bound that the top curves form. In other words, above the upper boundary of the blue region, the value of the scalar spectral index is incorrect, for any ϕ * .
Increasing the range for ϕ * does not change the shape of the upper boundary of the parameter space in the (log 10 ξ * , κ) plane. Indeed, we also plot more κ(ξ * )| ns=0.9649 curves, now in purple, with ϕ * ranging from −200 m P to −30 m P . We find that in the range −30 m P to 0 we cover almost the entirety of the shown parameter space, while approaching more negative values simply covers a region of the parameter space discarded by observations, located at smaller and smaller values of κ.
It could also be that changing αM 4 would change the shape of the parameter space. However, we find that the main effect of this is on r. Indeed, there exists a bound, given by αM 4 /m 4 P 0.143 below which the size of the orange region in the parameter space is reduced in size (although it never fully disappears, see the left panel in Fig. 5), and above which its position shifts towards larger values of κ and ξ. This can be seen by comparing the middle and right panels in Fig. 5. It is also straightforward to see from Eqs. (3.3) and (3.5) that r can be made arbitrarily small by making αM 4 larger, as we have obtained in our numerical study (see Fig. 4). However, it is important to note that the shift in the orange region of the parameter towards larger κ and ξ can change the subsequent cosmological evolution after inflation ends, since these points serve as the initial conditions for the later evolution.
We conclude that the shape of the blue region shown in Fig. 2 is an universal feature of the model, with the caveat that the analysis concerning Fig. 3 is for the β = 0 case. We expect only minor modifications to this figure when studying the general non-zero β case, since during slow-roll inflation the value of the field barely changes and we choose the scale β=0 Figure 3: Slice of the parameter space in the (log 10 ξ * , κ) plane with β = 0 and αM 4 /m 4 P = 1.43, where we plot many curves κ(ξ * )| ns=0.9649 with ϕ * ranging from −30 m P to 0 (green) and from −200 m P to −30 m P (purple), as well as the curve κ 2 = 4ξ(= 4ξ * ) (red), so that the condition for a monotonic potential κ 2 > 4ξ(= 4ξ * ) is satisfied above it. The upper boundary of the parameter space coincides with the asymptotic upper bound from the green curves. Increasing ϕ * to more negative values explores a region of the parameter space that is not in agreement with observations, towards smaller and smaller κ, as can be seen from the purple curves. The parameter space of the theory lies between the asymptotic upper bound from the κ(ξ * )| ns=0.9649 curves and the condition κ 2 = 4ξ(= 4ξ * ), as it should. µ to be approximately equal to the field value at horizon exit ϕ * , making the running in ξ negligible. In the same spirit, it is obvious that the shapes of the blue and orange regions in Fig. 3, for which β = 0, are very similar to the analogous regions in Fig. 2, for which β = −0.1.
To conclude this section, in Fig. 6, we show an example plot of the scalar spectral index as a function of the number of e-folds before the end of inflation in the Einstein frame. The shape of n s (N ) in Fig. 6 is general and for most of the valid points of the parameter space,

Numerical results for post-inflationary evolution
In order to gain some understanding about the model, we start this section by studying one specific benchmark point of the parameter space which leads to correct dark energy predictions. After this we show the full parameter space of our quintessential inflation model.
Let which satisfies all the conditions listed above required for correct inflation and dark energy. This can be immediately confirmed by looking at Figs. 7 and 8. In the left panel in Fig. 7 we show the barotropic parameter of the inflaton and of the whole universe, which are given by where w r,m is equal to either 1/3 or 0 for a radiation (r) or a pressureless dust (m) background with energy densityρ r,m , respectively, andρ φ andp φ are given by Eq. (2.42). At the present time, which corresponds toN = 0 in both figures, the energy fraction of the field is Ω 0 φ = 0.6889 (see the right panel in Fig. 7) and its barotropic parameter and running are w 0 φ = −0.95895 and w 0 a = −0.17034, in agreement with dark energy observations. As for the energy densities at present it can be confirmed by looking at the right panel in Fig.  8, that the energy density of the field isρ φ = 1.7 × 10 −120 m 4 P while that of the fluid is ρ m = 7.5 × 10 −121 m 4 P , which are within an order of magnitude of observations. Finally, the temperature of the universe at the onset of radiation domination, i.e., when w tot = 0.36 and Ω φ = 0.05, is T = (30ρ/(π 2 g)) 1/4 2.49 × 10 −23 m P = 0.15 MeV, which is slightly above T BBN . As a far as inflationary observables and dark energy predictions go, the point given by Eq. (4.10) is fine. However, as the careful reader might have noticed, there are two issues with the matter dominated era. As it can be seen in the left panel in Fig. 7, its durationN mat = 7.5 is below what would be expected in a standard cosmology, whereN mat ∼ 8. Furthermore, the barotropic parameter of the universe is not exactly zero (although it stays below 0.1). We can explain this behaviour by taking a closer look at our model. We remind the reader that, as shown in Eqs. (2.45) and (2.47), there is a coupling between the inflaton and the fluid (the last term in both equations) coming from the conformal factor that appears in the matter action after the conformal transformation to the Einstein frame. During inflation we havē ρ r,m = 0 and during kination and the radiation dominated era we have that the barotropic parameter of the fluid is w = 1/3, so that the coupling is not present until matter-radiation equality. However, as soon as we have a pressureless dust-dominated universe, with w = 0, the coupling is turned on. In order to better understand this, after some simple algebra, one can rewrite Eq. (2.47) as where where the last equality follows from working in the matter dominated era and a prime denotes a derivative with respect to the Jordan frame number of e-folds. Thus, w eff will only be close to zero when the rate of change of F R satisfies (4.14) However, looking at the expression for F R in Eq. (2.22), and remembering that the terms coming from the α contribution are negligible at late times, the rate of change from Eq. (4.14) is approximately F R /F R ∼ ϕ /ϕ. By noticing that the field is in free fall, and, thus, has a non-negligible rate of change, during the matter dominated era (its barotropic parameter is one 7 as can be seen in Fig. 7) it immediately follows that F R /F R cannot be very small and w eff will be generally larger than zero, as we find.
As for the number of e-folds of the matter dominated eraN mat , noting that from Eq. (4.12) follows thatρ ∝ā −3(1+w eff ) , a simple calculation reveals where we have taken into account that w eff 0.1, as is the case for most of the valid parameter space. Thus,N mat will generally be smaller than its canonical value in standard Einstein-Hilbert gravity, where there is no coupling between the fluid and the inflaton so that w eff = 0.
Introducing the values of the energy density of the fluid at equality,ρ eq = 1.27 × 10 −110 m 4 P , and at the present time,ρ 0 = 3.28 × 10 −121 m 4 P , we find thatN mat could be decreased by as much as about one e-fold. We take this into account in the parameter space scans, not neglecting points that a priori would have been considered to have a too short matter dominated era. In this way, we choose six as the smallest valueN mat can take when scanning over the parameter space, although, as we will see below, for all valid pointsN mat will always be larger than seven, in agreement with the approximation w eff 0.1 that we have taken above.
In conclusion, we have obtained that the barotropic parameter of the universe during the matter dominated era will generally be larger than zero and that the length of this era will generally be shorter than in Einstein-Hilbert gravity. These effects are an inevitable consequence of working in our modified gravity setup. However, we find that for most of the parameter space w eff 0.1 (and discard the points which do not satisfy this), and in fact, around redshifts corresponding to galaxy formation, i.e.,z ∼ 4 (wherez ≡ā −1 − 1) [186], the barotropic parameter is very close to zero, thereby not significantly impeding structure formation (see Fig. 9).
Having discussed the effect of the inflaton-fluid coupling, modified gravity manifests itself in the Einstein frame through one other effect: the existence of a quartic kinetic term in the action (see the third term in Eq. (2.15)), which a priori cannot be discarded. However, as it can be seen from the left panel in Fig. 8, it remains subdominant throughout the expansion history of the universe. This is a general behaviour in all the valid parameter space. In what follows we neglect this term.
Let us next examine the evolution of the system more carefully, stage by stage. As the field approaches the end of the inflationary plateau and its velocity starts increasing, the condition H = 1 is satisfied and inflation ends. After the end of inflation there is a transition period where the field is gaining kinetic energy although its total energy density is still not dominated by it. This can be seen from Fig. 10, where we show the energy density ratios (4.16) 7 It could be that the higher order kinetic terms that appear in the Einstein frame modify the barotropic parameter of the field from its usual expression w φ = ( 1 . This is not the case, as can be seen in Fig. 8, where it is clear that the quartic kinetic term plays a subdominant role throughout the expansion history of the universe. density of the inflaton is kinetically dominated, while the energy density of the universe is still dominated by that of the field (Ω φ is still equal to one as can be seen from the right panel in Fig. 7), giving way to the kination era. This can also be seen from the left panel in Fig.  7, where the barotropic parameter of the field does not become equal to one untilN ∼ −40.
Of course, at the moment when the field becomes kinetically dominated, remembering the quartic kinetic terms are negligible, we have During kination, the radiation energy density fraction approaches that of the field, until it takes over and approaches one aroundN = −19.6, see Fig. 7. This moment corresponds to reheating. It is important to note that the scaling of the energy density of the universe betweenN = −50.6 andN = −40 is notρ ∝ā −6 (but slower), and thus the exponent 6 on the right-hand-side of Eq. (3.15) is only an approximation. Effectively this means that we are able to have a reheating temperature close to T BBN without violating the bound on Ω end GW discussed in Sec. 3.3. After reheating, the universe is dominated by the background radiation, while the field is still in free-fall, with its energy density being kinetically dominated. This can be seen from the left panel in Fig. 7, where the barotropic parameter of the field is still equal to one, as well as from the left panel in Fig. 8 and from Fig. 10. This behaviour continues until briefly before matter-radiation equality, when the field runs out of kinetic energy and starts to freeze (see Eq. (3.21)). Indeed, its barotropic parameter approaches minus one (this can also be seen from Fig. 10, where the kinetic density ratio goes from one to zero and vice versa for the potential density ratio). However, the field never fully freezes. This is due to the change in the barotropic parameter of the background from 1/3 to 0 at matter-radiation equality. As explained above, at this point the coupling between the field and the fluid is turned on and there is an energy transfer between the components. One way to understand this is by noting that w eff is larger than zero, meaning that the background dilutes faster than in the canonical case, feeding its energy into the kinetic energy of the field. Indeed, the barotropic parameter of the field jumps back to unity and the inflaton goes back into free-fall during the entirety of the matter-dominated era, only to run out of kinetic energy and freeze again (its barotropic parameter going back to minus one) at the end of it. Finally, the field does not simply slow down and freeze. If it did, we would not find the small bump in its barotropic parameter afterN = 0 in Fig. 7. The same bump can be found in Fig. 10. This is due to the local minimum of the potential in the Einstein frame, located slightly before the local maximum around 1 + ξ(ϕ max )ϕ 2 max /m 2 P = 0 (see discussion in Sec. 2.1). Indeed, the field overshoots the minimum and gains some kinetic energy, only to fall back to the minimum at ϕ min = 884.03m P and finally freeze. The present timeN = 0 corresponds to some time briefly after overshooting the minimum but before turning back.
Having characterised the dynamics of a typical valid parameter space point, including the effects of the modified gravity terms, let us now turn our attention to the location and shape of the full valid parameter space. We show some example slices in the (log 10 ξ * , κ) plane for different values of αM 4 in Fig. 11. We also scan over the parameter β, which in the orange and blue regions in the figure takes values in the interval [−0.11, −0.098] in steps of 10 −3 . We find that points giving rise to correct dark energy (shown in red and black), which satisfy the whole set of constraints given above, are only found for β = −0.099 and β = −0.105. We also show in Tab. 1 the actual parameter values all of the successful points take. We find they form no specific shape in the (log 10 ξ * , κ) plane, but expect a higherresolution scan to reveal more working points. Lowering the required minimum temperature  of the universe at the onset of the radiation dominated era, such that it is no longer larger than T BBN , makes the valid parameter space follow a curved area inside the orange region. However, imposing the appropriate bound spoils this behaviour. It is worth mentioning that although our selection criteria regarding the length of the matter dominated era is for it to be longer than 6 e-folds, allowing for a non-zero w eff to decreaseN mat , all valid points actually have at least 7 e-folds, although they are always below 8 e-folds. It is possible that the rest of constraints regarding the energy density and the barotropic parameter make the parameter space to lie in this interval.
To conclude, in this section we have characterised the behaviour, both for the field dynamics and for the modified gravity effects, of a typical successful point in the parameter space. We have also found the location of the valid points in the (log 10 ξ * , κ) plane, having scanned over β in the [−0.11, −0.098] interval in steps of 10 −3 and over αM 4 /m 4 P in the  4 , which just needs to be larger than a given lower bound αM 4 ∼ 0.1 8 . Indeed, the most successful points have κ = 0.27, log 10 ξ * = −1.9 and β = −0.099.

Conclusions
In this paper we studied a relatively simple model of quintessential inflation where a single scalar field can unify the two epochs of accelerated expansion in the history of the Universe: inflation and dark energy domination. We worked in the framework of Palatini gravity where the metric and the connection are treated as independent variables. The three main ingredients in our action are: • An exponential potential of the form M 4 e −κϕ/m P which for large positive values of the scalar field produces the quintessential tail.
• An αR 2 term which asymptotically flattens the potential for large negative values and produces inflation in agreement with observations.
• A non-minimal coupling ξϕ 2 R between the quintessence/inflaton field and gravity, where ξ ≈ ξ * is approximately constant and positive during inflation but then runs to negative values with a slope β in order to reproduce the correct late-time dark energy. Note that the region where ξ(ϕ) is negative is never probed since the field freezes before that.
The main advantage of employing the Palatini formalism is that the auxiliary field introduced in order to parametrise the R 2 term turns out to be non-dynamical and can therefore be eliminated through its equation of motion. The resulting action is then single field, but contains a quartic kinetic term and a modified effective potential. For sufficiently large values of α, the effective potential is always asymptotically flat and can therefore accommodate slow-roll inflation. In addition to the quintessence/inflaton field, we considered an ideal fluid representing the matter and radiation content of the universe. We began our analysis by examining the equations of motion for the field and the fluid in both the Jordan and Einstein frames, while at the same time relating the quantities of interest in the two frames. We determined the Jordan frame equations to be easier to solve numerically. We then studied separately all the phases arising during the time evolution of our model in a cosmological setup, namely, inflation, kination, reheating, radiation and matter domination, and finally quintessence. To produce the radiation component after inflation, we considered as an example Ricci reheating [160][161][162], where an additional scalar field with a non-minimal coupling to gravity reheats the universe during a period of kination. For quintessence, we showed that the Einstein frame scalar field potential develops a local minimum where the field eventually gets stuck, behaving like dark energy afterwards. The minimum is generated by the non-minimal coupling of the scalar field running to negative values. The dark energy density there is generated through the interplay of the different parameters, all taking natural values, avoiding the extreme fine-tuning of the cosmological constant in the standard ΛCDM scenario.
In the end, we presented a thorough analysis of our numerical procedure and results. We scanned over the inflationary parameter space and showed that, for correct choices of the parameter values, the inflationary predictions of the model match the Planck observations [14]. For late-time evolution, we noted the emergence of a coupling between the fluid and the scalar field, present in the Einstein frame during matter domination. This coupling turned out to be the biggest obstacle for our model building, threatening to disrupt the standard cosmic evolution by transferring energy from the matter fluid to the rolling field. Nevertheless, we found example points that satisfy all the criteria we set for a successful cosmological scenario, in particular for the present-time energy density and barotropic parameter of the quintessential dark energy component. We obtain definite predictions for all of the parameters of our model. The preferred parameter values which give rise to successful results are κ = 0.27, log 10 ξ * = −1.9 and β = −0.099. We did not find a preference for any specific value for the combination αM 4 , as long as it is above the threshold αM 4 /m 4 P ∼ 0.1, below which the tensor-to-scalar ratio is too large to be compatible with observations. In addition to satisfying all the available observational constraints, our model also offers testable predictions, to be probed in the future by experiments such as EUCLID [187]. Indeed, a non-zero derivative of the barotropic parameter of dark energy with respect to the scale factor (w a in the CPL parametrization), as is the case in our model, would favor dynamical dark energy models over a cosmological constant (as in ΛCDM). Our model offers specific predictions for w a , which will be useful to discern between dynamical dark energy models as measurements become more precise. It also features a non-zero barotropic parameter of the universe, probing redshifts between galaxy formation and equality, i.e.,z ∼500-1500.
To conclude, our model produces successful inflation and quintessential dark energy from the above-listed simple set of ingredients alone, without the extreme fine-tuning of ΛCDM. Our model is the first one (barring the toy-model in Ref. [79]) to produce successful quintessential inflation using modified gravity as the main ingredient.

Acknowledgments
KD is supported (in part) by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC grant: ST/T001038/1. AK was supported by the Estonian Research Council grants PSG761, MOBTT5, MOBTT86, and by the EU through the European Regional Development Fund CoE program TK133 "The Dark Side of the Universe". SSL is supported by the FST of Lancaster University. ET was supported by the Estonian Research Council grants PRG803, MOBTT5, and PRG1055 and by the EU through the European Regional Development Fund CoE program TK133 "The Dark Side of the Universe".

A Solving for the Hubble parameter
In this appendix, we solve the Jordan frame Hubble parameter H in (2.28) explicitly in terms of ϕ,φ, and the fluid energy density ρ. We begin by using (2.22) During inflation, the background matter energy density ρ = 0, and the condition is always satisfied when the potential dominates the kinetic term, 4V (ϕ) >φ 2 , in particular during slow-roll. It is also easy to satisfy later on, when ρ > 0 becomes important and α contributions become irrelevant.
B A bound on the bare mass-squared of the spectator field Let us estimate the upper bound of |m 2 |, the mass squared of the spectator field ψ from (3.11), such that it remains negligible at least until reheating. Firstly, let us obtain an upper bound of the value of ψ 2 at the end of inflation. Imposing the requirement thatρ end ψ <ρ end r , as found by Ref. [162], and using thatρ end ψ = 1 4 λ ψ 2 2 , we find where we considered Eq. (3.12). Considering that the typical value of the amplitude of the oscillating condensate at a given location is of the order |ψ| ∼ ψ 2 , we can estimate how it evolves after the end of inflation. Indeed, because 1 4 λ|ψ| 4 =ρ ψ ∝ā −4 , we have |ψ| ∝ā −1 . Thus, we obtain The above is too strict because, ifρ ψ ρ r after the end of inflation, thenρ ψ can remain subdominant until reheating even if the quadratic term in V (ψ) takes over before reheating. So the above bound is sufficient but not, strictly speaking, necessary. Its numerical value may be estimated using the range obtained in Eq. (3.18). Taking λ ∼ 1, we find 10 2 GeV m max < 10 15 GeV. (B.5) The lower bound in the above might be unrealistic because such a particle could have been already observed in the LHC. But we see that the mass range extends well above the TeV scale so there is no real conflict with the observational data.