Unifying inflation with early and late dark energy epochs in axion F ð R Þ gravity

We provide a theoretical model of F ð R Þ gravity in which it is possible to describe in a unified way inflation, an early and a late dark energy era, in the presence of a light axion particle which plays the role of the dark matter component of the Universe. Particularly, the early time phenomenology is dominated by an R 2 term, while the presence of the other terms f ð R Þ ensure the occurrence of the early and late-time dark energy eras. The inflationary phenomenology is compatible with the Planck 2018 data for inflation, while the late-time dark energy era is compatible with the Planck 2018 constraints on the cosmological parameters. Also, the model exhibits an early dark energy era, at z ∼ 2 . 5 approximately, followed by a deceleration era, which starts at approximately z ∼ 1 . 5 , which in turn is followed by a late-time dark energy era for redshifts z ∼ 0 . 5 , which lasts for approximately 5 billion years up to present time. A notable feature of our model is that the dark energy era is free from dark energy oscillations, at least in the redshift interval z ¼ ½ 0 ; 10 (cid:2) . In addition, we also discuss several features related to observational data at z ∼ 2 . 34 , at which redshift intricate observational data exist in the literature. Moreover, the numerical code for the dark energy phenomenology, written in PYTHON 3 , is presented in the end of the article. Finally, the model has another interesting characteristic, a sudden jump of the value of the Hubble rate in the redshift interval z ∼ ½ 2 ; 2 . 6 (cid:2) where its value suddenly increases and then decreases until z ∼ 0 .


I. INTRODUCTION
The observation that utterly changed the way of thinking in modern theoretical cosmology was made in the late 90's [1] and indicated that the Universe is expanding in an accelerated way, contrary to the standard model of cosmology. This late-time acceleration era is known as the dark energy era and to date still remains a mystery on its driving force. In the context of general relativity (GR), the dark energy era can be generated by using a cosmological constant, and the model that is successful in fitting the observational data is the Λ-cold-dark-matter (ΛCDM) model. However, the ΛCDM model strongly relies on the presence of a tiny cosmological constant, which is added by hand in the model, and of course in the presence of the other mysterious component of the dark sector, dark matter [2][3][4][5][6][7]. The fact that the cosmological constant has such a small value, is in conflict with the vacuum energy predictions of quantum field theory, thus such a shortcoming has to be alleviated in a theoretical way. Moreover, in the context of GR, the early and late-time era are usually described by two distinct theories; in the case of inflation, the usual description is given by a canonical scalar field with a sufficiently flat potential, while the late-time era is generated by the cosmological constant. Thus, a unified description is lacking in the context of GR. Modified gravity bridges the gap between the early and late-time phenomenology, since a unified description can be achieved; see, for example, the pioneer article [8] and Refs. [9][10][11][12][13][14][15][16] and also Refs. [17][18][19][20][21][22] for reviews. Apart from the unified description of the early and the late-times acceleration eras, modified gravity has the attribute of generating a dark energy era with nonconstant equation of state (EOS) parameter, since for the ΛCDM model, the dark energy era has the EOS P ¼ −ρ.
In this paper, we shall present an FðRÞ gravity model for which the unification of the early and the late-time acceleration eras can be achieved. In addition, the model also generates an early dark energy era around z ∼ 5, followed by a deceleration era, and the latter is followed by the present day acceleration. In addition to the FðRÞ gravity model, we assume the presence of a light misalignment axion field [15, which has a broken Peccei-Quinn symmetry during inflation. During inflation, the axion remains frozen at its vacuum expectation value, but after the inflationary era it starts oscillating and behaves as cold dark matter. We analyze the inflationary phenomenology of the model at hand in detail and we demonstrate that it is compatible with the Planck constraints on inflation [79]. Also, we provide a thorough numerical analysis of the dark energy era, utilizing a numerical code written in PYTHON 3 which is freely available and can be found in [80] (see the Appendix for details on the code). For the dark energy era, we quantify our analysis by expressing all the physical quantities as functions of the redshift and as functions of an appropriate statefinder quantity. As we demonstrate, the results are compatible with the latest Planck constraints on the cosmological parameters [81]. Furthermore, an early dark energy era is produced by our theoretical model, starting approximately at z ∼ 2.5 and ending at z ∼ 1.5. Accordingly, the model decelerates until z ∼ 0.5 and after that it accelerates again until the present time era. In the literature, several models of early dark energy models have been worked out; see, for example, Refs. [82][83][84][85][86]. Furthermore, we analyze an interesting feature of the model, a sudden jump in the Hubble rate at z ∼ 2. Interestingly enough, the Hubble rate at the aforementioned redshift has a local minimum followed by a local maximum, and after that it decreases until the present day. Finally, we briefly discuss our results in view of the observational data [87] for redshifts z ∼ 2.
The motivation for using a combined scalar field and fðRÞ gravity framework is twofold. First, the scalar field in our case is basically a misalignment axion field, which is a primordial scalar field coming as a remnant of the quantum theory preceding the inflationary era. Second, the inflationary era emerged directly from a quantum epoch of our Universe, where all interactions were unified and possibly spacetime could be higher dimensional. However, during the inflationary era, the Universe was four-dimensional and classical; nevertheless, there inflationary effective Lagrangian may contain higher order curvature terms and even string corrections as remnants of the quantum era of our Universe. Such terms could be fðRÞ gravity terms, such as the case we shall study in this paper, or even Gauss-Bonnet terms coupled to the scalar field, thus effectively having an Einstein-Gauss-Bonnet theory . Einstein-Gauss-Bonnet theories in the presence of R 2 terms were studied in Ref. [128], and it would be interesting to study these in the case that the scalar field is actually the axion; however, this is out of the scope of this paper. Also let us notice that if higher order curvature terms are present in the inflationary effective Lagrangian, these terms would apparently dominate the evolution at early times compared to the Einstein-Hilbert term or even the frozen axion, since even for the low scale inflationary scenario, the Hubble rate during inflation is of the order of H I ∼ Oð10 13 Þ GeV. We shall explicitly show this in the following sections.

II. ESSENTIAL FEATURES OF THE f ðRÞ-AXION MODEL
We shall assume that the gravitational action of our model has the following form: with where G is Newton's gravitational constant and M p is the reduced Planck mass. In addition, L m describes the Lagrangian density of the perfect matter fluids present, which in our case is radiation. The light axion field is quantified by the scalar field ϕ and its dynamics will be explained shortly. The FðRÞ gravity model has the following form: where m s is m 2 m is the present day energy density of cold dark matter, and 0 < δ < 1, while ζ and γ are dimensionless constants to be specified latter. Also, M is chosen to be M ¼ 1.5 × 10 −5 ð N 50 Þ −1 M p for inflationary phenomenological reasoning [129], where N is the e-foldings number. Moreover, the parameter Λ takes values of the order of the cosmological constant and has dimensions eV 2 , while λ is a dimensionless parameter. For a flat Friedmann-Robertson-Walker (FRW) metric, the field equations for the fðRÞ gravity in the presence of axion dark matter and radiation, and the scalar field reads with F R ¼ ∂F ∂R , and the "dot" and the "prime" denote differentiation with respect to the cosmic time and the scalar field, respectively.
Let us turn our focus to the axion scalar field, and let us discuss how it evolves during the evolution of the Universe. Here we shall briefly outline the axion dynamics, which is described in detail in Refs. [15,26]. We shall consider an axion scalar field with a broken Uð1Þ Peccei-Quinn symmetry during the inflationary era, with a potential during inflation of the form [26] VðϕÞ where ϕ i is the large vacuum expectation value of the axion field during inflation. The axion field is frozen during the inflationary era to its vacuum expectation value and does not evolve dynamically, so the following initial conditions describe it as [26] ϕðt i Þ ≃ 0; where t i is the cosmic time during the inflationary era, f a is the axion decay constant, θ a denotes the initial misalignment angle, and m a is the axion mass. The above frozen dynamical evolution continues to occur as long as H ≫ m a . When m a ∼ H, the axion starts to oscillate and thus for the eras that m a ≽ H; the axion energy density evolves as [15,26] where ρ ð0Þ m ¼ 1 2 m 2 a ϕ 2 i ; hence, it evolves as a cold dark matter component for the eras for which the condition m a ≽ H holds true. In the following, we shall study the dynamics of the model (1) for both the early and late-time eras.

III. INFLATIONARY EVOLUTION OF THE MODEL
During the inflationary era, the Hubble rate is of the order of H I ¼ 10 13 GeV for most grand unified theories with low-scale inflation. Thus, the Ricci scalar R ∼ H 2 I takes quite large value; hence, at leading order, the effective FðRÞ of the model (2) becomes and thus at leading order, the effective FðRÞ gravity during inflation is equal to Now, in order to analyze the inflationary dynamics, we need to understand which terms in the equations of motion (4) dominate the evolution. First, m s which was defined below Eq. (2) is m 2 s ≃ 1.87101 × 10 −67 eV 2 and the param- [129]; hence, for N ∼ 60, M is equal to M ≃ 3.04375 × 10 22 eV. Moreover, since the inflationary era is a slow-roll era, meaning that _ H ≪ H 2 , we have R ≃ 12H 2 , and for H ¼ H I ∼ 10 13 GeV, we have R ∼ 1.2 × 10 45 eV 2 . Also the reduced Planck mass is M p ≃ 2.435 × 10 27 eV, and the parameter Λ will be assumed to take the value Λ ≃ 11.895 × 10 −67 eV 2 , which is close to the value of the cosmological constant at present day. Finally, for phenomenological reasoning, the values of ϕ i and m a are chosen to be ϕ i ¼ Oð10 15 Þ GeV and m a ≃ Oð10 −14 Þ eV. Let us now compare the terms appearing in the equations of motion, having in mind that the leading order terms in the FðRÞ gravity are ∼R and R 2 . First, the radiation density term during inflation is highly subdominant since κ 2 ρ r ∼ e −N during inflation, and also the kinetic term of the axion scalar ∼ _ ϕ 2 i is also subdominant, since the axion during inflation is frozen and obeys the initial conditions (7). Now let us proceed with the potential term which is of the order κ 2 Vðϕ i Þ ∼ Oð8.41897 × 10 −36 Þ eV 2 , while the terms R and R 2 are R ∼ 1.2 × Oð10 45 Þ eV 2 and also R 2 =M 2 ∼ Oð1.55 × 10 45 Þ eV 2 . Finally, the power-law curvature terms, for δ ¼ 0.1 and ζ ¼ 0.2 (which are the values for ζ and δ we shall also assume for the late-time analysis), take values of the order Þ eV 2 and the rest of the terms are highly subdominant. Thus during the inflationary era, the dynamical evolution of the cosmological system is controlled by the fðRÞ gravity, By looking the effective form of the FðRÞ gravity appearing in Eq. (11), one would expect that the term (1 þ λ) will affect the phenomenology of this deformed R 2 model. We shall work out in detail the inflationary phenomenology of the model (11) and to our surprise, for the R 2 model there is no difference between the λ ≠ 0 and λ ¼ 0 theories. Let us explicitly show this, so by substituting Eq. (11) in the Friedmann equation in Eq. (4), we get without any approximation, the following differential equation: Since the inflationary era is governed by the slow-roll approximation, quantified by the relations, in view of the relations (13), the Friedmann equation (12) becomes which is solved easily and yields a quasi-de Sitter evolution, where H I is an integration constant, which is basically the inflationary scale. Now, we can obtain easily the phenomenology of the model and the slow-roll parameters are [17,88,130] and the spectral index of the primordial scalar curvature perturbations and the tensor-to-scalar ratio are [130] n s ≃ 1 − 6ϵ 1 − 2ϵ 4 ; ð17Þ Now, the exact expression for the slow-roll index can easily be found [130], and it is equal to Hence, for the deformed R 2 model of Eq. (11), the slow-roll index ϵ 4 is ϵ 4 ≃ −ϵ 1 . Hence, the spectral index of the primordial curvature perturbations and the tensor-to-scalar ratio are The slow-roll index ϵ 1 can easily be found by using the analytic quasi-de Sitter solution for the Hubble rate (14), and it is equal to and from it we can easily find the time instance at which inflation ends, namely, t f , by solving the equation ϵ 1 ðt f Þ ¼ 1, and the solution is We can now use the definition of the e-foldings number N, in order to find t i as a function of the e-foldings number and the rest of the parameters, so we have Now if we calculate the slow-roll index ϵ 1 at the first horizon crossing time instance t i , we get remarkably a λ-independent result, so the spectral index and the tensor-to-scalar ratio are at leading order n s ∼ 1 − 2 N and r ∼ 12 N 2 , which are identical to the R 2 model with λ ¼ 0. The model is compatible with the latest Planck data [81]. Thus, the inflationary phenomenology of the model (2) is viable and compatible with the Planck data. As a final comment, the fact that the parameter λ does not affect at all the phenomenology of the model is an artifact of the R 2 model and of the fact that the axion scalar is not dynamically evolving during the inflationary era. In fact, it can be shown that if a dynamically evolving canonical scalar field is present, this rescaled Einstein-Hilbert R term affects the inflationary phenomenology significantly in some cases [131].

IV. EARLY AND LATE DARK ENERGY ERAS
For the FðRÞ gravity theory in the presence of radiation and the axion dark matter fluid, the field equations for the flat FRW metric can be cast in the Einstein-Hilbert form as follows: where ρ tot ¼ ρ a þ ρ DE þ ρ r stands for the total energy density, ρ a is the axion field energy density, which recalls that it scales as ∼a −3 , while ρ r is the radiation energy density and ρ DE is the dark energy density, a purely geometric term since it is equal to Also, P tot ¼ P r þ P a þ P DE denotes the total pressure of the cosmological fluid, and the dark energy pressure is In this section, we aim to study the late-time behavior of the model (2), by solving numerically the Friedmann equation (4). To this end, we shall use the redshift as a dynamical parameter quantifying the evolution instead of the cosmic time, and also we shall introduce an appropriate statefinder quantity, which will make the dark energy effects more transparent. The redshift z is defined as and the present time scale factor, that is, at z ¼ 0, is assumed to be equal to unity. The statefinder we shall use is the function y H ðzÞ [15,17,[132][133][134][135], with ρ ð0Þ m being the present time energy density of cold dark matter. Obviously, y H ðzÞ is a dark energy-dependent quantity, and it is different from zero, when geometric terms appear in the gravitational action. We can write the first Friedmann equation in terms of the statefinder quantity y H ðzÞ and by recalling that for the axion field where we introduced the parameter χ ¼ ρ ð0Þ and ρ ð0Þ r is the radiation energy density, and recall that Apparently, the statefinder y H ðzÞ clearly shows deviations from the standard model of cosmology, and it is constant for the ΛCDM model. Hence, it is indeed a statefinder for dark energy since it shows deviations from the Einstein gravity and also shows if the dark energy is dynamical or not, and the latter issue is still a mystery too, along with dark energy itself. In terms of the statefinder y H ðzÞ, the Friedmann equation reads [133] with the dimensionless functions J 1 , J 2 , and J 3 being defined in the following way: with F RR ¼ ∂ 2 F ∂R 2 . Our aim is to solve numerically the differential equation (33) in the redshift interval z ¼ ½0; 10, by using appropriate physical motivated initial conditions which correspond to the late stages of the matter domination era. These are We developed a numerical code appropriately constructed to integrate the differential equation (33) backward from z ¼ 10 to z ¼ 0, using PYTHON and the total EOS parameter is Finally, the deceleration parameter reads Also for ΛCDM model, the Hubble rate equals to with H 0 being the value of the Hubble rate at present time, which according to the latest Planck data is H 0 ≃ 1.37187 × 10 −33 eV according to the latest Planck data [81]. Also Ω Λ ≃ 0.681369 and Ω M ∼ 0.3153 [81], while Ω r =Ω M ≃ χ, with χ being presented below Eq. (32). Also, for the numerical analysis, we shall choose . Now, let us discuss the results of our numerical analysis and we gather the most characteristic examples of the latetime behavior for the model (2) in Fig. 1, where we plot the statefinder y H ðzÞ (left plot), the total EOS parameter (right plot) and in Fig. 2 the deceleration parameter, as functions of the redshift. Also, in Fig. 3, we plot the Hubble rate as a function of the redshift. The behavior of the statefinder y H ðzÞ clearly shows that y H ðzÞ is essentially negative for the whole interval z ¼ ½0; 10, and this has interesting consequences in view of the observational data reported in Ref. [87], so we comment on the y H ðzÞ behavior in the next section. The behavior of the deceleration parameter can clearly indicate when the Universe is accelerating and when the Universe is decelerating. In our case, the Universe decelerates until z ∼ 2.5, then it enters an acceleration phase, until z ∼ 1.5, followed by a decelerating epoch which lasts for nearly 4 billion years until z ∼ 0.5. After z ∼ 0.5, the Universe enters the final acceleration era which lasts until present day, nearly 5 billion years. In Fig. 2, we also quote the ΛCDM behavior of the deceleration parameter (red), and as it can be seen, the FðRÞ model is indistinguishable from the ΛCDM model only for the last 2 billion years approximately, from z ∼ 0.2 to z ¼ 0. The same conclusion applies in the behavior of the total EOS parameter in the right plot of Fig. 1, where the red curve represents the behavior of the ΛCDM model. Thus, the FðRÞ gravity model clearly has two dark energy epochs, one early and one late dark energy epoch, and in between these two dark energy eras, the Universe is decelerating. The same behavior can be verified by looking the right plot of Fig. 1, that is, the total EOS parameter for the FðRÞ gravity model (blue curve), while the red curve corresponds to the ΛCDM model. An interesting behavior in the Hubble rate which occurs for z ∼ 2.6 until z ∼ 2 is that the Hubble rate at z ∼ 2.2 has a local minimum, followed by a local maximum at z ∼ 2 and then the Hubble rate evolves normally to its present day value for the FðRÞ gravity model. This intriguing behavior is presented in Fig. 3, where we plot H=Hð0Þ as a function of the redshift. The curious behavior is highlighted on the right plot, and it seems that it is an inherent characteristic of the model; we could not attribute this behavior to any other qualitative reasoning. Also we should note that the fraction of the present day value of the Hubble rate for the FðRÞ gravity model, over the Planck value for the Hubble rate at present day is Hð0Þ H 0 ¼ 0.563636. Finally, let us discuss the viability of the FðRÞ gravity model in a more quantitative way, so let us give the values of the dark energy density parameter Ω DE and of the dark energy EOS parameter at present day, and we compare these with the Planck 2018 constraints on these quantities. The PYTHON code using the "LSODA" integration technique yields ω DE ð0Þ ≃ −0.9961 which is within the viability limits of the Planck constraint ω DE ¼ −1.018 AE 0.031, and accordingly the dark energy density parameter is found equal to Ω DE ð0Þ ≃ 0.6872, which is an acceptable value when compared to the Planck constraint Ω DE ¼ 0.6847 AE 0.0073. Also the deceleration parameter value at present day for the FðRÞ gravity model is qð0Þ ≃ −0.5267 while for the ΛCDM model is q ¼ −0.527; hence, the two values are nearly identical.
A feature of the present FðRÞ gravity model that is worth mentioning is the fact that the model is free from dark energy oscillations which are known to plague FðRÞ gravity models for the redshift interval z ¼ ½5; 10. This in an inherent characteristic of the model as it proves, and thus we may conclude that the occurrence of dark energy oscillations in FðRÞ gravity models is a model-dependent feature.
In conclusion, the FðRÞ gravity model (2) in the presence of a light axion field produces a unified phenomenological picture for the inflationary era and the dark energy eras, with the intriguing characteristic of producing actually two distinct dark energy era, one early and one late dark energy eras, with a brief deceleration epoch occurring between these two deceleration eras. Finally, let us finally comment that the inverse integration PYTHON code we used and the Mathematica 11 results agree to a great extent.

A. The behavior at z ∼ 2.34
Before closing, let us briefly discuss another interesting feature of the model, which is related to the measurement of the Hubble rate at z ¼ 2.34 reported by [87]. Assuming that the observation is correct, we shall discuss how such an observation is supported by the FðRÞ gravity model (2). Basically, the observation of [87] indicates that the Hubble rate at redshift z ∼ 2.34 is Hðz ¼ 2.34Þ ¼ 222 km=Mpc= sec, a result which is also discussed in [136][137][138][139]. An interesting way to explain the result of [87] is to make the assumption that ρ DE < 0 in Eq. (32), since otherwise one would get Ω m h 2 ¼ 0.142 which obviously does not agree with the cosmic microwave background related value of the Planck Collaboration data Ω c h 2 ¼ 0.12 AE 0.001 [81]. Although that a negative dark energy density sounds strange, this concept has appeared in the literature [140,141]. In our case, negative dark energy density would mean negative values for the statefinder y H ðzÞ, and obviously by looking the upper left plot of Fig. 1, we can see that y H ðzÞ is negative for nearly the whole interval z ¼ ½0; 10. In fact, it only becomes positive at z ∼ 1.4, where it has the value y H ð1.4Þ ¼ 0.227363 and has the present day value y H ð0Þ ¼ 2.19527. Let us also comment that negative values of y H ðzÞ in FðRÞ gravity models frequently occur, since most fðRÞ gravity models are plagued with dark energy oscillations for the interval z ¼ ½5; 10 [15,17]. The present model is free from oscillations though; however, y H is negative for nearly the whole interval z ¼ ½0; 10.

B. The ghost issue
Before closing, we need to discuss the important issue of ghosts in the essentially fðR; ϕÞ theory we discussed in this paper. The ghost instabilities may be developed in a modified gravity theory if the wave speed of the cosmological perturbations, which we will denote as c A , is larger than unity. As was shown in Ref. [88], however, for fðR; ϕÞ theories, the wave speed of the scalar perturbations is unity; thus, no ghost instabilities occur in the curvature perturbations. In fact, for a flat spacetime, which is our case, the sound speed is equal to the wave speed of the cosmological perturbations (see the last table of Ref. [88]). Let us elaborate on this issue a bit further in order to clarify this important issue. We shall obtain the cosmological perturbations by perturbing the flat FRW metric g ð3Þ αβ in the following way [88]: Hð0Þ as a function of the redshift. In the right plot, the intriguing behavior of the Hubble rate for redshifts z ∼ 2.6 until z ∼ 2 is shown in a more transparent way. The green line in the right plot corresponds to the H=Hð0Þ value 2.56196. g ð3Þ αβ dx α dx β ¼ dr 2 þ r 2 ðdθ 2 þ sin 2 θdϕ 2 Þ: ð43Þ In addition, η and a in Eq. (42) denote the conformal time and the scale factor. Moreover α, β, γ, and φ in Eq. (42) denote the scalar type order variables of the cosmological perturbations, while the tensor C αβ stands for the traceless and transverse tensor perturbation. The dynamical evolution of the scalar type perturbation quantified by the variable Φ ¼ φ δϕ , and which is directly related to the scalar type gauge invariant quantity δϕ φ ¼ − _ ϕ H φ δϕ , is determined by the differential equation that follows: with Δ being the Laplacian of the spatial section of the FRW metric. It is apparent from the right-hand side of the above evolution equation that the wave speed of the scalar perturbations is equal to unity; thus, no ghost instabilities are expected in the cosmological perturbations. But why does the wave speed determine actually whether or not ghost instabilities occur in the first place? Basically, by looking at Eq. (44), the wave speed of the perturbation Φ ¼ φ δϕ is equal to unity, but apart from the perturbation Φ for a general fðR; ϕÞ theory, there is another scalar perturbation, denoted as Ψ, which is defined as follows: the evolution of which is determined by the differential equation that follows: with the quantity S being defined below, Hence, by looking at both Eqs. (44) and (46), which are essentially wave equations, the propagation wave speed is equal to unity, so the wave speed of the propagation of both the fluctuating field and of the perturbed metric is equal to unity. In fact, since the FRW metric we use has a flat spatial part, the sound wave speed defined as c 2 ρ is identical to the wave speed, which is equal to unity. In order to make this more apparent, let us bring both Eqs. (44) and (46) to the Mukhanov-Sasaki form. We definez ¼ c A z, with c A ¼ 1 in our case, and z is while E is In addition, we introducev ¼ zΦ and also u ¼ a κ 2 H 1 z Ψ; hence, the wave equations appearing in Eqs. (44) and (46), respectively, can be brought in the well-known Mukhanov-Sasaki formsv where in our case c A ¼ 1. Thus, the wave speed of the fluctuating fluid and of the perturbed metric is equal to unity; hence, no ghost instabilities occur in the theoretical framework we used.

V. CONCLUDING REMARKS
In this work, we investigated the phenomenology of an FðRÞ gravity model in the presence of a primordial light axion scalar field with broken Uð1Þ Peccei-Quinn symmetry during the inflationary era. Due to the fact that the axion scalar is frozen to its primordial vacuum expectation values, it has a highly subdominant role during the inflationary era, which is predominantly controlled by the R 2 term of the FðRÞ gravity. The resulting effective Lagrangian during inflation has also a rescaled Einstein-Hilbert R term at leading order during inflation. As we demonstrated by performing an explicit calculation, the resulting inflationary phenomenology is identical with that corresponding to the R 2 model; thus, surprisingly, the rescaled R term does not affect the dynamical evolution during inflation, at least at leading order. However, if the axion was dynamically evolving during inflation, this result would not be valid. Actually, the rescaled Einstein-Hilbert gravity can significantly affect the inflationary phenomenology of a canonical scalar field as we will show in [131]. After the inflationary era, and specifically when the axion mass m a becomes of the same order as the Hubble rate, and for all cosmic times for which m a ≽ H, the axion field starts to oscillate and its energy density evolves as ρ a ∼ a −3 , thus behaves as a dark matter perfect fluid. Accordingly, we examined the dark energy phenomenology of the FðRÞ gravity model in the presence of the axion dark matter perfect fluid and of the radiation perfect fluid. After expressing the Friedmann equation as a function of the redshift and of a suitable statefinder function, we solved numerically the resulting equation, using suitable physically motivated initial conditions corresponding to the late stages of the matter domination era. For the numerical analysis, we developed a PYTHON 3 numerical code which is freely available here [80], along with a pedagogical description of the FðRÞ gravity dark energy phenomenology. The resulting picture is quite interesting phenomenologically, since three novel features for FðRÞ gravity appear in the theory. Firstly, the Universe experiences an early dark energy era, starting at z ∼ 2.5 and ending at z ∼ 1.5, followed by a decelerating epoch which lasts until z ∼ 0.5, and for the approximately remaining 5 billion years, from z ∼ 0.5 to z ¼ 0, the Universe accelerates again. We found that the model is also viable and compatible with the 2018 Planck constraints on the cosmological parameters, at least when the dark energy density parameter and the dark energy EOS parameters are considered. Also the model mimics the ΛCDM model significantly at the last stages of the evolution prior to present day redshift z ¼ 0 and at least qualitatively. Second, the model is free from dark energy oscillations, which are known to plague the FðRÞ gravity dark energy phenomenology. Thus, our result indicates that the dark energy oscillations might be a model-dependent feature of FðRÞ gravity dark energy phenomenology. Third, the model also complies with the z ∼ 2.34 observation of [87]. Finally, the model has also another interesting characteristic occurring at z ∼ 2, where the Hubble rate has a local minimum, followed by a local maximum, after which evolves normally up to z ¼ 0. This however is another model-dependent feature. An interesting phenomenological concept which emerged from the effective Lagrangian of the present FðRÞ gravity model at early times is the appearance of a rescaled Einstein-Hilbert term αR. Although for the present work, this term played a subdominant term during inflation, due to the presence of the R 2 term, in conjunction with the fact that the axion scalar was frozen in its vacuum expectation value, it would be interesting to consider the same FðRÞ gravity model without the R 2 term, and in the presence of a dynamically evolving canonical scalar field. As we shall show in [131], the results can be quite intriguing in some cases of interest.