From a bounce to the dark energy era with F(R) gravity

A major challenge in modern theoretical cosmology is to reveal whether the Universe emerged from an initial singularity or the Universe started its expansion from a non-singular bounce-like phase. In other words, whether the Universe evolve according to the standard big-bang cosmology which must end in a singularity (known as big-bang singularity) when extrapolated backwards, or, the Universe passes through a bouncing stage leading to a non-singular behavior of the spacetime curvature. The inflationary description is an appealing early Universe scenario, due to the fact that it can produce an almost scale invariant power spectrum, which is consistent with the latest observations. The inflationary scenario [1–8] requires an early period of accelerated expansion in order to solve the horizon and flatness problems. However, most inflationary scenarios have inherent connection with an initial singularity, at least most of the canonical scalar models of inflation. Moreover, the observations to date do not undoubtedly confirm that inflation indeed took place. One of the attractive alternatives to inflation is the bouncing cosmology scenario [9–53]. Apart from generating an observationally compatible power spectrum in some bouncing cosmology scenarios, bouncing cosmology has the additional advantage that it leads to a singular free evolution of the Universe. However, the big-bang singularity may be a manifestation for the failure of classical gravity theory which may be unable to describe the evolution of the Universe at such small scales. Quite possibly, the quantum generalization of the gravity theory may resolve the singularity issue; just like the classical electrodynamics fails to remove the Coulomb potential singularity at the origin, which, however is resolved by the quantum electrodynamics. However, in the absence of a fully accepted quantum gravity theory, bouncing cosmology is the most promising one to deal with the initial singularity issue.

Among the various bouncing models proposed so far, the matter bounce scenario (MBS) [14, 22, 23, 54–70] has attracted a lot of attention since it produces a scale invariant power spectrum and moreover in a MBS, the Universe passes through a matter dominated (MD) epoch at the late-time. The MBS is characterized by a symmetric scale factor where the Hubble radius diverges to infinity at late times both in the contraction and expanding eras of the bounce. This indicates that the perturbation modes are generated far away from the bouncing point, deeply inside the Hubble radius. Despite the aforementioned successes that the MBS is able to produce a scale invariant power spectrum, the matter bounce model has some serious drawbacks: (1) in scalar–tensor MBS model, the scalar power spectrum becomes exactly scale invariant leading to the spectral index for curvature perturbations as unity, which is, in fact, not compatible with the Planck observations. Such inconsistency was also confirmed in [64, 67] from a different point of view, in particular in the context of F(R) gravity theory. However, it is well known that a F(R) model can be equivalently mapped to a scalar–tensor model by a suitable conformal transformation of the spacetime metric and thus the inconsistencies of the scalar power spectrum (with regard to the Planck observational data) in MBS for both the scalar–tensor and F(R) theories are well justified. (2) The running of the scalar spectral index is constrained to lie within −0.0085 ± 0.0073 according to the Planck 2018 observations. However, for the MBS model with a single scalar field, the running of the spectral index vanishes and thus it is not compatible with the observations. This is a direct consequence of the first problem because if the spectral index is independent of the perturbation momentum, then the running index will obviously become zero. (3) The Hubble radius, defined by ${r}_{\mathrm{h}}=\frac{1}{aH}$, in the MBS monotonically increases with the cosmic time in the late era and asymptotically diverges to infinity, due to the reason that the scale factor far away from the bounce behaves as ∼t^2/3. Such increasing behavior of the Hubble radius confirms a decelerating evolution of the Universe at late times. On other hand, the supernovae observations [71–73] confirm that the present Universe is dominated by some negative pressure matter component, known as DE, which generates the accelerated expansion. Thus the dark energy (DE) dominated epoch, which is responsible for the late time acceleration, is not well described by the matter bounce scenario.

It is shown in the literature, that in the so-called quasi-matter bounce scenario, where the scale factor evolves as ${t}^{\frac{2}{3\left(1+w\right)}}$ (with w ≠ 0, note for w = 0, it becomes similar to the exact MBS), it is possible to recover the consistency of the spectral index and of the running index even in a single scalar field model [14]. However, the tensor-to-scalar ratio is found to be problematic in the case of a quasi-matter bounce model; actually the tensor and the scalar perturbation amplitudes in the quasi-matter bounce scenario are found to be comparable to each other leading to the tensor-to-scalar ratio as order of unity. On other hand, modified gravity theories [74–81] may provide successful descriptions for both the primordial and the late-time era of our Universe. In most modified gravity descriptions of bouncing cosmology, the Hubble radius monotonically increases at late times and thus the problem to describe the DE epoch of the Universe still persists, because an increasing behavior of the Hubble radius reveals a decelerating phase of the Universe. Motivated by this problem, in this paper, we intend to generalize the bouncing scenario which is also compatible with the late-time acceleration as well, in the context of F(R) gravity theory. It is clear that in order to get a DE dominated phase, the Hubble radius must decrease with respect to the cosmic time at late times. This in turn indicates that the primordial curvature perturbation modes generate near the bounce, in contrast to the usual MBS where the perturbation modes generate far away from the bounce deeply in the contracting era. Keeping these issues in mind, we try to provide a cosmological scenario which unifies certain cosmological epochs of the Universe, in particular, from a non-singular bounce to a MD era and from the MD to a late time DE epoch, in the context of F(R) gravity. In this regard, we would like to mention that the merging of bounce with the DE epoch was studied earlier [82], however in a different context, where the Universe evolution was not considered to be symmetric with respect to the bounce point, unlike to our present work where we will consider a symmetric behavior of the scale factor as a function of the cosmic time and consequently the Hubble radius asymptotically goes to zero in both sides of the bounce.

The outline of this paper is as follows: after briefly describing the essential features of F(R) gravity in section 2, we will demonstrate the behavior of F(R) gravity near the bounce and during the late-era respectively in section 3 and consequently in section 4, we will perform the scalar and tensor perturbations. Some numerical treatment of several qualitative aspects of the theory at hand are presented in section 5. Finally the conclusions follow in the end of the paper.

Let us briefly recall some basic features of F(R) gravity, which are necessary for our presentation, for reviews on this topic see [74, 75, 77]. The gravitational action of F(R) gravity in vacuum is equal to,

$$\begin{equation}S=\frac{1}{2{\kappa }^{2}}\int {\mathrm{d}}^{4}x\sqrt{-g}F\left(R\right)\end{equation} \tag{ 1 }$$

where κ² stands for ${\kappa }^{2}=8\pi G=\frac{1}{{M}_{\mathrm{P}\mathrm{l}}^{2}}$ and also M_Pl is the reduced Planck mass. By using the metric formalism, we vary the action with respect to the metric tensor g_μν , and the gravitational equations read,

$$\begin{equation}{F}^{\prime }\left(R\right){R}_{\mu \nu }-\frac{1}{2}F\left(R\right){g}_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{F}^{\prime }\left(R\right)+{g}_{\mu \nu }\square {F}^{\prime }\left(R\right)=0\end{equation} \tag{ 2 }$$

where R_μν is the Ricci tensor constructed from g_μν . Since the present article is devoted to cosmological context, the background metric of the Universe will be assumed to be a flat Friedmann–Robertson–Walker (FRW) metric,

$$\begin{equation}\mathrm{d}{s}^{2}=-\mathrm{d}{t}^{2}+{a}^{2}\left(t\right)\left[\mathrm{d}{x}^{2}+\mathrm{d}{y}^{2}+\mathrm{d}{z}^{2}\right]\end{equation} \tag{ 3 }$$

with a(t) being the scale factor of the Universe. For this metric, the temporal and spatial components of equation (2) become,

$$\begin{align}\hfill 3{H}^{2}& =-\frac{f\left(R\right)}{2}+3\left({H}^{2}+\dot {H}\right){f}^{\prime }\left(R\right)-18\left(4{H}^{2}\dot {H}+H\ddot {H}\right){f}^{{\prime\prime}}\left(R\right)\hfill \\ \hfill 0& =\frac{F\left(R\right)}{2}-\left(3{H}^{2}+\dot {H}\right){F}^{\prime }\left(R\right)+6\left(8{H}^{2}\dot {H}+4{\dot {H}}^{2}+6H\ddot {H}+ \dddot {H}\right){F}^{{\prime\prime}}\left(R\right)\hfill \\ \hfill & \quad +36{\left(4H\dot {H}+\ddot {H}\right)}^{2}{F}^{\prime \prime \prime }\left(R\right)\hfill \end{align} \tag{ 4 }$$

respectively, where $H=\dot {a}/a$ is the Hubble parameter of the Universe and f(R) is the deviation of F(R) gravity from the Einstein gravity, that is F(R) = R + f(R). Comparing the above equations with the usual Friedmann equations, it is easy to understand that F(R) gravity provides a contribution in the energy–momentum tensor (EMT), with its effective energy density (ρ_eff) and pressure (p_eff) being given by,

$$\begin{equation}{\rho }_{\text{eff}}=\frac{1}{{\kappa }^{2}}\left[-\frac{f\left(R\right)}{2}+3\left({H}^{2}+\dot {H}\right){f}^{\prime }\left(R\right)-18\left(4{H}^{2}\dot {H}+H\ddot {H}\right){f}^{{\prime\prime}}\left(R\right)\right]\end{equation} \tag{ 5 }$$

$$\begin{align}\hfill {p}_{\text{eff}}& =\frac{1}{{\kappa }^{2}}\left[\frac{f\left(R\right)}{2}-\left(3{H}^{2}+\dot {H}\right){f}^{\prime }\left(R\right)+6\left(8{H}^{2}\dot {H}+4{\dot {H}}^{2}+6H\ddot {H}+ \dddot {H}\right)\right.\hfill \\ \hfill & \quad {\times}\left.{f}^{{\prime\prime}}\left(R\right)+36{\left(4H\dot {H}+\ddot {H}\right)}^{2}{f}^{{\prime\prime\prime}}\left(R\right)\right]\hfill \end{align} \tag{ 6 }$$

respectively. Thus, the effective EMT depends on the form of F(R), as expected. Therefore, different forms of F(R) will lead to different evolution of the Hubble parameter. In the present context, we will use such effective EMT of F(R) gravity to realize the cosmological evolution of the Universe.

In this section, we reconstruct the F(R) gravity from equation (4) by considering a certain form of the scale factor of the Universe. In particular, we reconstruct the F(R) gravity in the two distinct eras, namely near the bouncing point and at the late-time epoch, which are the primary subjects in the following two subsections respectively. Clearly such forms of F(R) will describe the gravity theory in the two extreme stages of the Universe and thus will not be able to reveal the Universe evolution as a whole. However later in section 5, we will numerically solve equation (4) and will determine the form of F(R) for a wide range of cosmic times, which will provide an unification of certain cosmological epochs of the Universe, particularly from a bounce to MD era, followed by a late time accelerating phase. During the numerical solution of equation (4), the form of F(R) in the two distinct eras, which will be analytically evaluated in the following two subsections, will act as boundary conditions.

3.1. Reconstruction near the bounce

In this section, we reconstruct the form of the F(R) gravity near the bounce of the Universe. The Universe's evolution in a general bouncing cosmology, consists of two eras, an era of contraction and an era of expansion. Some of the well known functional forms of the scale factor which correspond to a non-singular bounce, have the form, $a\left(t\right)={\text{e}}^{\alpha {t}^{2}}$, a(t) =  cosh t, $a\left(t\right)={\left({a}_{0}{t}^{2}+1\right)}^{n}$ and so on. In general, the non-singular bounce of a scale factor is characterized by

$$\begin{equation*}a\left({t}_{\mathrm{b}}\right)\ne 0,\qquad \dot {a}\left({t}_{\mathrm{b}}\right)=0,\qquad \ddot {a}\left({t}_{\mathrm{b}}\right){ >}0\end{equation*}$$

where t_b is the cosmic time when the bounce occurs. The above conditions lead to a finite Kretschmann scalar (K = R_μναβ R^μναβ ) and thus the Universe becomes free of the initial singularity (known as big-bang singularity). Keeping these conditions in mind, we consider the scale factor near the bounce as,

$$\begin{equation}{a}_{\mathrm{b}}\left(t\right)=1+\alpha {t}^{2}\end{equation} \tag{ 7 }$$

(the suffix 'b' stands for near bounce scale factor) with α being a free parameter having mass dimension [+2] and the bounce happens at t = 0. The above form of the scale factor can be thought as a Taylor series expansion of a(t) around t = 0 and keeping up-to quadratic order in cosmic time (t). We neglect the higher orders of t in the Taylor expansion of a(t) as, in this section, we are interested to reconstruct the F(R) gravity near the bouncing point. The linear order of t in the Taylor expansion vanishes due to the condition $\dot {a}=0$, necessary for a bounce. Moreover the condition $\ddot {a}\left(0\right){ >}0$ indicates that the parameter α must be positive and thus we take α > 0 in the present context. Here we would like to mention that in the next sections, we will also provide a complete form of the scale factor which is valid for a wide range of cosmic time and may yield the behavior (7) near the bouncing point. The scale factor of equation (7) immediately leads to the following Hubble parameter ($H=\frac{\dot {a}}{a}$) and Ricci scalar (R(t)) as,

$$\begin{align}\hfill H\left(t\right)& =\frac{2\alpha t}{1+\alpha {t}^{2}}\simeq 2\alpha t,\hfill \\ \hfill R\left(t\right)& =12{H}^{2}+6\dot {H}=\frac{12\alpha \left(1+3\alpha {t}^{2}\right)}{{\left(1+\alpha {t}^{2}\right)}^{2}}\simeq 12\alpha +12{\alpha }^{2}{t}^{2}\hfill \end{align} \tag{ 8 }$$

respectively, with the H(t) and R(t) being considered up to $\mathcal{O}\left({t}^{2}\right)$, similar to the case of the scale factor. However equation (8) clearly indicates that the Hubble parameter varies linearly with t and goes to zero at the bouncing point, while the Ricci scalar, on the other hand, becomes R(0) = 12α. Later, during the calculations of scalar and tensor perturbations, we will show that the parameter α should be at the order ∼10⁻⁸/κ² to make the model compatible with the Planck constraints and thus the Ricci scalar becomes ∼10²⁸ GeV² (with $\frac{1}{{\kappa }^{2}}={M}_{\mathrm{P}\mathrm{l}}^{2}=1{0}^{36}\enspace {\mathrm{G}\mathrm{e}\mathrm{V}}^{2}$) at the bouncing point. In order to reconstruct the F(R) gravity, we will use equation (4), for which we determine the following quantities (appearing in equation (4)) with the scale factor given in equation (7),

$$\begin{align*}\hfill {H}^{2}+\dot {H}& =\frac{2\alpha }{1+\alpha {t}^{2}}\simeq 2\alpha -2{\alpha }^{2}{t}^{2},\hfill \\ \hfill 4{H}^{2}\dot {H}+H\ddot {H}& =\frac{8{\alpha }^{3}{t}^{2}\left(1-3\alpha {t}^{2}\right)}{{\left(1+\alpha {t}^{2}\right)}^{4}}\simeq 8{\alpha }^{3}{t}^{2}\hfill \end{align*}$$

With the above expressions, equation (4) becomes,

$$\begin{equation}24\alpha \left(R-12\alpha \right){f}_{\mathrm{b}}^{\prime \prime }\left(R\right)+\left(R-24\alpha \right){f}_{\mathrm{b}}^{\prime }\left(R\right)+{f}_{\mathrm{b}}\left(R\right)+2\left(R-12\alpha \right)=0\end{equation} \tag{ 9 }$$

Solving the above equation for f(R), we get,

$$\begin{equation}{f}_{\mathrm{b}}\left(R\right)=\left(\frac{12\alpha D}{\sqrt{\mathrm{e}}}-1\right)R-D\sqrt{6\alpha \pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}{\left(R-12\alpha \right)}^{3/2}\enspace \mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]\end{equation} \tag{ 10 }$$

where Erfi[z] is the imaginary error function defined as Erfi[z] = −iErf[iz] with Erf[z] being the error function and 'i' is the imaginary unit. Moreover D is an integration constant having mass dimension [−2]. The above solution of f_b(R) immediately leads to the form of F_b(R) = R + f_b(R) as,

$$\begin{equation}{F}_{\mathrm{b}}\left(R\right)=\frac{12\alpha D}{\sqrt{\mathrm{e}}}R-D\sqrt{6\alpha \pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}{\left(R-12\alpha \right)}^{3/2}\enspace \mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]\end{equation} \tag{ 11 }$$

The presence of ${\mathrm{e}}^{-\frac{R}{24\alpha }}$ and $\mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]$ in the right-hand side of equation (11) confirm that the F_b(R) contains all the positive integer powers of R where the coefficient of the linear Ricci scalar is other than unity. This indicates a deviation from Einstein gravity. However recall, the form of F(R) in equation (11) represents the gravity theory near the bounce and in a later section 5, during the numerical reconstruction of F(R) gravity for a wide range of cosmic times, we will show that the late time F(R) indeed matches with the Einstein gravity.

3.2. Reconstruction in the late time epoch

In this section, we apply a reconstruction scheme to describe the behavior of the Universe at late times. Recall, in the previous section during the reconstruction near the bounce, we considered a scale factor (suitable for bounce) and then, by solving equation (4), we determine the corresponding form of F(R). However the late-time reconstruction method which we are going to apply in this section, is slightly different in comparison to that of the earlier one, in particular, we will start with a form of F(R) (rather than a scale factor) suitable for a DE model and then reconstruct the corresponding Hubble parameter from the gravitational equation of motion.

In order to reproduce the current acceleration of the Universe, several versions of viable modified gravity including the so-called one-step models have been proposed in references [83–87]. The simplest one is given by the so-called exponential gravity which includes an exponential function of the Ricci scalar in the action. We consider such exponential F(R) model in the present section, which is known to provide a good DE model as described in [88, 89]. In particular, we consider,

$$\begin{equation}{F}_{\mathrm{l}}\left(R\right)=R-2{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)\end{equation} \tag{ 12 }$$

with β and Λ are two model parameters having mass dimensions [0] and [+2] respectively and the suffix 'l' is for 'late' time epoch. The exponential correction over the usual Einstein–Hilbert action becomes important at cosmological scales and at late-times, providing an alternative to the DE problem. Due to the supernovae Ia (Sne-Ia) [90, 91], baryon acoustic oscillations (BAO) [92–95], cosmic microwave background (CMB) [96–99] and H(z) [100–103] data, the parameters β and Λ are well constrained, particularly the F(R) model (12) is best fitted with Sne-Ia + BAO + H(z) + CMB data for the parametric regimes given by: $\beta =3.9{8}_{-2.46}^{+\infty }$ and Λ = 1.2 × 10⁻⁸⁴ GeV² [88]. Such parametric ranges also make the model compatible with Sne-Ia + BAO + H(z) data. In effect, we stick with β = 4 and Λ = 1.2 × 10⁻⁸⁴ GeV² throughout the paper. With these viable ranges of the model parameters along with the F(R) of equation (12), now we are going to solve the gravitational equation of motion to reconstruct the evolution of the Hubble parameter during the DE epoch. Plugging the above form of F(R) into equation (4), we get,

$$\begin{align}\hfill & -3{H}_{\mathrm{l}}^{2}+{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)-3\beta {\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\left({H}_{\mathrm{l}}^{2}+{\dot {H}}_{\mathrm{l}}\right)\hfill \\ \hfill & \qquad -\frac{9{\beta }^{2}}{{\Lambda}}\enspace {\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\left(4{H}_{\mathrm{l}}^{2}{\dot {H}}_{\mathrm{l}}+{H}_{\mathrm{l}}{\ddot {H}}_{\mathrm{l}}\right)=0\hfill \end{align} \tag{ 13 }$$

where the 'dot' represents $\frac{\mathrm{d}}{\mathrm{d}t}$. However getting an analytic solution of equation (13) is troublesome, so we proceed to solve it numerically and for this purpose, we introduce a dimensionless time scale (t_r) and a dimensionless Hubble parameter (H_r) defined as follows:

$$\begin{equation*}{t}_{\mathrm{r}}=\frac{t}{{t}_{\mathrm{s}}},\quad {H}_{\mathrm{r}}=\frac{1}{{a}_{\mathrm{l}}}\frac{\mathrm{d}{a}_{\mathrm{l}}}{\mathrm{d}{t}_{\mathrm{r}}}\end{equation*}$$

with t_s = 10¹⁶ s i.e. at the order of the present age of the Universe. The actual Hubble parameter (H_l(t)) is related with the rescaled one as ${H}_{\mathrm{l}}\left(t\right)=\frac{1}{{t}_{\mathrm{s}}}{H}_{\mathrm{r}}\left({t}_{\mathrm{r}}\right)=1{0}^{-16}{H}_{\mathrm{r}}$ (s)^−1, moreover the derivative of the Hubble parameter in the two time coordinates are connected as ${\dot {H}}_{\mathrm{l}}\left(t\right)=\frac{1}{{t}_{\mathrm{s}}^{2}}\frac{\mathrm{d}{H}_{\mathrm{r}}}{\mathrm{d}{t}_{\mathrm{r}}}$ (note, $\frac{\mathrm{d}}{\mathrm{d}t}$ is represented by a 'dot' while the $\frac{\mathrm{d}}{\mathrm{d}{t}_{r}}$ is mentioned by as it is). In view of these relations, equation (13) can be expressed in terms of the rescaled quantities as follows,

$$\begin{align}\hfill & -3{H}_{\mathrm{r}}^{2}+{\Lambda}{t}_{\mathrm{s}}^{2}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)-3\beta {\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\left({H}_{\mathrm{r}}^{2}+\frac{\mathrm{d}{H}_{\mathrm{r}}}{\mathrm{d}{t}_{\mathrm{r}}}\right)\hfill \\ \hfill & \qquad -\frac{9{\beta }^{2}}{{\Lambda}{t}_{\mathrm{s}}^{2}}{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\left(4{H}_{\mathrm{r}}^{2}\frac{\mathrm{d}{H}_{\mathrm{r}}}{\mathrm{d}{t}_{\mathrm{r}}}+{H}_{\mathrm{r}}\frac{{\mathrm{d}}^{2}{H}_{\mathrm{r}}}{\mathrm{d}{t}^{2}}\right)=0\hfill \end{align} \tag{ 14 }$$

where $\frac{R}{{\Lambda}}=\frac{1}{{\Lambda}{t}_{\mathrm{s}}^{2}}\left(12{H}_{\mathrm{r}}^{2}+6\frac{\mathrm{d}{H}_{\mathrm{r}}}{\mathrm{d}{t}_{\mathrm{r}}}\right)$. With β = 4 and ${\Lambda}{t}_{\mathrm{s}}^{2}=1.5{2}^{2}{\times}1{0}^{-4}$ (the conversion 1 s = 1.52 × 10²⁴ GeV⁻¹ may be useful), equation (14) is solved numerically and is given in the left part of figure 1, where the Hubble parameter is plotted for 60 ⩽ t_r ⩽ 90 (or equivalently 6 × 10¹⁷ ⩽ t ⩽ 9 × 10¹⁷ s) during the late era. Moreover we use H_r(t_r = 60) = 0.1 and ${\left.\frac{\mathrm{d}{H}_{\mathrm{r}}}{\mathrm{d}{t}_{\mathrm{r}}}\right\vert }_{{t}_{\mathrm{r}}=60}=-0.1$ as the boundary conditions to solve the above equation numerically. The x-axis of the left part of figure 1 is the rescaled time t_r and by using the relation H_l = 10⁻¹⁶ H_r s⁻¹, we give the actual Hubble parameter H_l(t) (in the unit of (sec)⁻¹) along the y-axis. The solution of the Hubble parameter immediately leads to the evolution of the Hubble radius (${r}_{\mathrm{h}}^{\left(\mathrm{l}\right)}=\frac{1}{{H}_{\mathrm{l}}\enspace \mathrm{exp}\left[\int {H}_{\mathrm{l}}\left(t\right)\enspace \mathrm{d}t\right]}$) as shown in the right part of figure 1 where once again, the x-axis is the rescaled time coordinate and the y-axis represents the actual Hubble radius i.e. ${r}_{\mathrm{h}}^{\left(\mathrm{l}\right)}\left(t\right)={\left({a}_{\mathrm{l}}\left(t\right){H}_{\mathrm{l}}\left(t\right)\right)}^{-1}$.

Zoom In Zoom Out Reset image size

Figure 1 clearly reveals that the Hubble radius decreases with time and leading to an accelerating stage of the Universe at late time. Therefore as a whole, the F(R) (12) successfully provides a viable DE model in respect to Sne-Ia + BAO + H(z) + CMB data where the evolution of the Hubble parameter and the Hubble radius are given in figure 1. Moreover, here it may be mentioned that the exponential F(R) model has a Schwarzschild de-Sitter solution in the regime R ≫ Λ ∼ 10⁻⁸⁴ GeV². For example, in the Solar System regime R* ≃ 10⁻⁷⁸ GeV², the exponential F(R) model can be approximated as F(R*) = R* − 2Λ (note that R* is greater than the Λ) and thus in this case, the Schwarzschild de-Sitter solution can be obtained, which is indeed a stable solution as depicted by some of our authors in [86].

Therefore equations (11) and (12) represent the F(R) gravity in the two situations, one depicting the early bounce of the Universe where the scale factor behaves as a(t) = 1 + αt² (=a_b(t)) and the other providing a late time accelerating phase when the scale factor (a_l(t)) or more explicitly the Hubble parameter H_l(t) goes according to the figure 1 (the corresponding scale factor can be easily obtained by using ${a}_{\mathrm{l}}\left(t\right)={\mathrm{e}}^{\int {H}_{\mathrm{l}}\left(t\right)\mathrm{d}t}$). Till now, we have described the bounce and the DE era in a disjoint manner i.e. F_b(R) describes the Universe near the bounce but not the whole cosmological epochs, while the F_l(R) reveals the behavior only at present time. Thus the immediate question that springs to mind is, what is the form of F(R) for the entire range of cosmic time, which will further demonstrate the unification of bounce and DE era followed by a deceleration stage in the intermediate region? As a deceleration era, we may consider a MD epoch in between the bounce and the present era. With this consideration, in section 5, we will reconstruct the full F(R) gravity, however numerically, by taking the forms of F_b(R) and F_l(R) as boundary conditions. In the method of numerical reconstruction, we need an appropriate scale factor (or equivalently Hubble radius) describing the Universe evolution for a wide range of cosmic time. For this purpose, we may consider an ansatz of the scale factor as follows,

$$\begin{align}\hfill a\left(t\right)& ={\text{e}}^{\alpha {t}^{2}}\enspace \mathrm{tanh}\left[{\beta }_{1}\left(1-\frac{t}{{t}_{\mathrm{p}}}\right)\left(1-\frac{t}{{t}_{\mathrm{I}}}\right)\right]+{\left(\frac{t}{T}\right)}^{2/3}\hfill \\ \hfill & \quad {\times}\enspace \mathrm{tanh}\left[{\beta }_{2}\left(\frac{t}{{t}_{\mathrm{b}}}-1\right)\left(1-\frac{t}{{t}_{\mathrm{p}}}\right)\right]\hfill \\ \hfill & \quad +{a}_{\mathrm{l}}\left(t\right)\enspace \mathrm{tanh}\left[{\beta }_{3}\left(\frac{t}{{t}_{\mathrm{b}}}-1\right)\left(\frac{t}{{t}_{\mathrm{I}}}-1\right)\right]\hfill \end{align} \tag{ 15 }$$

where t_b is the bouncing time or even the horizon crossing time, t_p is the present time instance and t_I is an intermediate time. The above scale factor behaves as a_b(t) = 1 + αt² near t ≃ t_b ≪ t_I ≪ t_p i.e. near the bounce, however goes as ${a}_{\mathrm{I}}\left(t\right)={\left(t/T\right)}^{2/3}$ near t_b ≪ t = t_I ≪ t_p and as ${a}_{\mathrm{l}}\left(t\right)={\mathrm{e}}^{\int {H}_{\mathrm{l}}\left(t\right)\mathrm{d}t}$ near the present time i.e. t ≃ t_p. The parameters β_i can be fixed by considering the corresponding hyperbolic tangent function tending to unity at the respective time.

Though the scale factor (15) provides the correct evolution of the Universe near the bounce and at the late stage and also at a certain instance of the intermediate stage (t ≃ t_I), but it does not lead to correct evolution of the Universe for the whole intermediate epoch. Thus to address the Universe evolution smoothly for a wide range of cosmic time, we may consider a(t) as,

$$\begin{align}\hfill a\left(t\right)& =1+\alpha {t}^{2}={a}_{\mathrm{b}}\left(t\right),\hfill & \hfill \mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\enspace \mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}& \hfill \quad t\simeq {0}^{+}\\ \hfill a\left(t\right)& ={\left(\frac{t}{T}\right)}^{2/3}={a}_{\mathrm{I}}\left(t\right),\hfill & \hfill \mathrm{i}\mathrm{n}\enspace \mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{m}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}& \hfill \quad {0}^{+}\lesssim t\lesssim 5{\times}1{0}^{17}\enspace \mathrm{sec}\\ \hfill a\left(t\right)& ={\mathrm{e}}^{\int {H}_{\mathrm{l}}\left(t\right)\mathrm{d}t}={a}_{\mathrm{l}}\left(t\right),\hfill & \hfill \mathrm{d}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\enspace \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}& \hfill \quad t\enspace \gtrsim \enspace 6{\times}1{0}^{17}\enspace \mathrm{sec}\end{align} \tag{ 16 }$$

(with T being an arbitrary parameter with dimension sec) and the transition from the bounce to the MD and from MD to the DE era will be obtained by the method of numerical interpolation. However, the process of interpolation requires appropriate choices for the values of the free parameters present in our model. It may be observed that there are three parameters hanging around in the theory, α, Λ and β. The parameter α appears from the near-bounce scale factor and thus can be named as 'early stage parameter', while Λ and β appear during the reconstruction in the late-time and thus the name—'late stage parameters'. The late stage parameters have been already estimated in view of the viability of the exponential F(R) as a successful DE model. However, the early stage parameter α is still remaining to be determined and in the next section, we will estimate α from the latest Planck 2018 observational results.

Being the early stage parameter, α can be determined from various primordial observable quantities like the scalar spectral index (n_s), the tensor-to-scalar ratio (r) by directly confronting their theoretical expectations with the Planck 2018 observations. In this section, we consider the spacetime fluctuations over the FRW metric and consequently calculate various observable quantities like the scalar spectral index and tensor-to-scalar ratio. In a bouncing Universe, the primordial perturbation modes (relevant to the present day observation) generate either near the bounce or at a distant past far away from the bouncing point, depending upon the late-time behavior of the Hubble radius. Thereby before moving to the explicit perturbation calculation, it is important to analyze when the perturbation modes generate in the present context of bouncing Universe. In all the bouncing models, the Hubble parameter becomes zero and thus the comoving Hubble radius, defined by ${r}_{\mathrm{h}}=\frac{1}{aH}$ (where H is the Hubble parameter), diverges at the bouncing point. However, the asymptotic behavior of the comoving Hubble radius makes a difference in various bouncing models. In this regard, (1) for some bouncing scale factors, the Hubble radius decreases monotonically at both sides of the bounce and finally shrinks to zero size asymptotically, which corresponds to an accelerating late-time Universe. Therefore, in such cases, the Hubble horizon goes to zero at large values of the cosmic time, and only for cosmic times near the bouncing point the Hubble horizon has an infinite size. So the primordial perturbation modes relevant for present time era generate for cosmic times near the bouncing point, because only at that time all the primordial modes are contained in the horizon. As the horizon shrinks, the modes exit the horizon and become relevant for present time observations. On the other hand, (2) some bouncing scale factors lead to a divergent Hubble radius asymptotically, which corresponds to a decelerating Universe at late times. In such cases, the perturbation modes generate at very large negative cosmic times, corresponding to the low curvature regime of the contracting era, unlike to the previous situations, where the perturbation modes generate near the bouncing era. More explicitly, in the latter case, the comoving wave number k begins its propagation through spacetime at large negative cosmic times, in the contracting phase on sub-Hubble scales, and exits the Hubble radius during this phase, and re-enters the Hubble radius during the low-curvature regime in expanding phase and thus being relevant for present time observations. Therefore, the physical picture in the two cases is very different with regard to when the perturbation modes generate. However as mentioned earlier, in the present work, we will deal with an unified scenario of bounce with late time acceleration connected by an intermediate deceleration epoch. Thus, the late-time Universe is characterized by an accelerating scale factor and hence the Hubble radius should shrink to zero asymptotically at both 'sides' of the bounce (as the scale factor is considered to be symmetric around the bounce). This indicates that the perturbation modes generate near the bounce, rather than far away from the bouncing point, because near the bouncing regime the Hubble radius has an infinite size and all the perturbation modes are contained inside the horizon. Hence we solve the perturbation equations near the bounce i.e. near t = 0, which is the primary subject in the remaining part of this section.

4.1. Scalar perturbation

In principle, perturbations should always be expressed in terms of gauge invariant quantities. In the present work, we shall work in the comoving gauge defined by the vanishing of the momentum density δT_0i = 0 (where T_μν is the effective EMT and the symbol 'δ' denotes the corresponding perturbation). In this gauge, the scalar metric perturbation is expressed as,

$$\begin{equation}\delta {g}_{ij}={a}^{2}\left(t\right)\left[1-2\k{e}\right]{\delta }_{ij}\end{equation} \tag{ 17 }$$

where $\k{e}\left(t, \overrightarrow {x}\right)$ denotes the scalar perturbation and known as the comoving curvature perturbation which is indeed a gauge invariant quantity. The additional metric perturbations δg₀₀ and δg_0i can be obtained in terms of  ę from perturbed gravitational equations and as a result, the second order action of $\k{e}\left(t, \overrightarrow {x}\right)$ is given by [104–106],

$$\begin{equation}\delta {S}_{\k{e}}=\int \mathrm{d}t\enspace {\mathrm{d}}^{3} \overrightarrow {x}a\left(t\right)z{\left(t\right)}^{2}\left[{\dot {\k{e}}}^{2}-\frac{1}{{a}^{2}}{\left({\partial }_{i}\enspace \k{e}\right)}^{2}\right],\end{equation} \tag{ 18 }$$

with z(t) has the following expression,

$$\begin{equation}z\left(t\right)=\frac{a\left(t\right)}{\kappa \left(H\left(t\right)+\frac{1}{2{F}^{\prime }\left(R\right)}\frac{\mathrm{d}{F}^{\prime }\left(R\right)}{\mathrm{d}t}\right)}\sqrt{\frac{3}{2{F}^{\prime }\left(R\right)}{\left(\frac{\mathrm{d}{F}^{\prime }\left(R\right)}{\mathrm{d}t}\right)}^{2}}\end{equation} \tag{ 19 }$$

Equation (18) clearly indicates that the sound speed of the scalar perturbation (${c}_{\mathrm{s}}^{2}$) is unity, which in turn guarantees the absence of superluminal modes or equivalently one may argue that the model is free from gradient instabilities. The unit sound speed for scalar perturbation is, in fact, a generic nature of F(R) theory and also of a scalar–tensor theory. This equivalence, in respect of sound speed, between F(R) and scalar tensor theory is however expected as both the theories can be mapped to one another (in the action level) by a conformal transformation of the metric. Coming back to the action (18), the scalar perturbation has positive kinetic terms if z²(t) > 0 holds, which is equivalent to the condition F'(R) > 0 as evident from equation (19). We will show that the F(R) obtained in the earlier sections indeed satisfies F'(R) > 0 and ensures the stability of the scalar perturbation. The action δS_Re leads to the equation for the perturbed variable $\k{e}\left( \overrightarrow {x},t\right)$ as,

$$\begin{equation}\frac{1}{a\left(t\right){z}^{2}\left(t\right)}\frac{\mathrm{d}}{\mathrm{d}t}\left[a\left(t\right){z}^{2}\left(t\right)\dot {\k{e}}\right]-\frac{1}{{a}^{2}}{\partial }_{i}{\partial }^{i}\enspace \k{e}=0\end{equation} \tag{ 20 }$$

In terms of the Fourier transformed scalar perturbation variable ${\mathfrak{R}}_{k}\left(t\right)=\int \mathrm{d} \overrightarrow {x}\enspace {\text{e}}^{-\text{i} \overrightarrow {k}. \overrightarrow {x}}\enspace \k{e}\left( \overrightarrow {x},t\right)$, the above equation can be written as,

$$\begin{equation}\frac{1}{a\left(t\right){z}^{2}\left(t\right)}\frac{\mathrm{d}}{\mathrm{d}t}\left[a\left(t\right){z}^{2}\left(t\right){\dot {\k{e}}}_{k}\right]+\frac{{k}^{2}}{{a}^{2}}{\mathfrak{R}}_{k}\left(t\right)=0\end{equation} \tag{ 21 }$$

where k is the wave number for the kth perturbation mode. As mentioned earlier, we solve the above equation for cosmic times near the bouncing point as the perturbation modes themselves are generated close to the bounce and thereby we can use the form of F(R) reconstructed in equation (11). With this expression of near-bounce-F(R), we will determine z(t) (an essential ingredient of perturbation equation (21)) and for this purpose, we first need to evaluate F'(R) and F''(R) which are given by,

$$\begin{align}\hfill {F}^{\prime }\left(R\right)& =\frac{12\alpha D}{\sqrt{\mathrm{e}}}+\frac{D}{24}\left[-\frac{6\left(R-12\alpha \right)}{\sqrt{\mathrm{e}}}\right.\hfill \\ \hfill & \quad +\left.\frac{\sqrt{6\pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}\left(R-48\alpha \right)\sqrt{R-12\alpha }\enspace \mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]}{2\sqrt{\alpha }}\right]\hfill \end{align} \tag{ 22 }$$

and

$$\begin{align}\hfill {F}^{{\prime\prime}}\left(R\right)& =\frac{D}{576{\alpha }^{3/2}}\left[\frac{12\sqrt{\alpha }\left(R-72\alpha \right)}{\sqrt{\mathrm{e}}}\right.\hfill \\ \hfill & \quad -\left.\frac{\sqrt{6\pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}\left({R}^{2}-96\alpha R+1440{\alpha }^{2}\right)\enspace \mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]}{\sqrt{R-12\alpha }}\right]\hfill \end{align} \tag{ 23 }$$

respectively. Recall, the Ricci scalar at the bounce becomes R(0) = 12α and thus the function $\mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]$ can be expressed as Taylor series expansion in the powers of $\sqrt{R-12\alpha }$ near the bounce. However the terms with $\mathcal{O}\left({\left(R-12\alpha \right)}^{5/2}\right)$ in the Taylor expansion can be neglected as such terms will eventually lead to a higher power of t in comparison to t² and all the near-bounce quantities have been or will be evaluated up to $\mathcal{O}\left({t}^{2}\right)$, as we also did for the near-bounce scale factor in section 3.1. Using the Taylor expansion of $\mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]$ along with the expression of R(t) (see equation (8)), we evaluate the F'(R) and F''(R) up to $\mathcal{O}\left({t}^{2}\right)$ as follows,

$$\begin{equation}{F}^{\prime }\left(R\left(t\right)\right)=\frac{12\alpha D}{\sqrt{\mathrm{e}}}-\frac{12\alpha D}{\sqrt{\mathrm{e}}}\enspace \alpha {t}^{2}\end{equation} \tag{ 24 }$$

and

$$\begin{equation}{F}^{{\prime\prime}}\left(R\left(t\right)\right)=-\frac{1}{\alpha }\left[\frac{\alpha D}{\sqrt{\mathrm{e}}}-\frac{4\alpha D}{3\sqrt{\mathrm{e}}}\enspace \alpha {t}^{2}\right]\end{equation} \tag{ 25 }$$

respectively. Equation (24) indicates that F'(R) is positive near the bounce and hence ensures the stability of the scalar perturbation. On other hand, F''(R) becomes negative near the bounce, which is clear from equation (25). Actually the condition F''(R) < 0 is connected with the violation of null energy condition in the F(R) gravity theory and thus such condition should hold in order for a bounce to be generated, and the demonstration goes as follows: recall, as mentioned in section 2 that the F(R) gravity contributes an effective EMT where the effective energy density (ρ_eff) and the pressure (p_eff) are given in equations (5) and (6) respectively. Using these expressions of ρ_eff and p_eff along with the Hubble factor H(t) = 2αt, it is easy to show that at the bounce ρ_eff + p_eff becomes $=2\dot {H}\left({F}^{\prime }\left(R\left(t\right)\right)-1\right)+24{\dot {H}}^{2}{F}^{{\prime\prime}}\left(R\left(t\right)\right)$. The F'(R) and F''(R), as we obtained near t = 0 in the present context, leads to ρ_eff + p_eff < 0 i.e. the violation of energy condition which in turn ensures a bouncing phenomena at t = 0. Plugging the expressions of F'(R), F''(R) into equation (19), we get,

$$\begin{equation}a\left(t\right){z}^{2}\left(t\right)=\frac{72\alpha D}{{\kappa }^{2}\sqrt{e}}\left[1+\frac{2}{3}\alpha {t}^{2}\right]={U}_{1}+{V}_{1}\enspace \alpha {t}^{2}\end{equation} \tag{ 26 }$$

where ${U}_{1}=\frac{72\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$ and ${V}_{1}=\frac{48\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$. With the above expression of a(t)z²(t), the scalar perturbed equation (21) takes the following form,

$$\begin{equation}{U}_{1}{\ddot {\k{e}}}_{k}+2\alpha {V}_{1}{\dot {\k{e}}}_{k}\enspace t+{k}^{2}{U}_{1}{\mathfrak{R}}_{k}\left(t\right)=0\end{equation} \tag{ 27 }$$

at leading order in t. Solving equation (27) for ${\mathfrak{R}}_{k}\left(t\right)$, we get,

$$\begin{equation}{\mathfrak{R}}_{k}\left(t\right)={b}_{1}\left(k\right)\enspace {\mathrm{e}}^{-\frac{{V}_{1}}{{U}_{1}}\alpha {t}^{2}}\enspace H\left[-1+\frac{{k}^{2}}{2\alpha }\frac{{U}_{1}}{{V}_{1}},\sqrt{\frac{\alpha {V}_{1}}{{U}_{1}}}\enspace t\right]\end{equation} \tag{ 28 }$$

with H[n, x] is the nth order Hermite polynomial. b₁(k) is the integration constant which can be determined by the initial state of the canonical Mukhanov–Sasaki variable v_k (t) defined by ${v}_{k}\left(t\right)=z\left(t\right){\mathfrak{R}}_{k}\left(t\right)$, which is the adiabatic vacuum state. Thus near the bounce (i.e. near t ≃ 0), the field v_k (t) satisfies

$$\begin{equation}{\mathrm{lim}}_{t\to 0}{\enspace v}_{k}=\frac{1}{\sqrt{2k}}\enspace {\mathrm{e}}^{-\mathrm{i}k\tau }\end{equation} \tag{ 29 }$$

where τ is the conformal time defined by $\mathrm{d}\tau =\frac{\mathrm{d}t}{a\left(t\right)}$. However near t ≃ 0, the conformal and cosmic time become same as,

$$\begin{equation*}\tau =\int \frac{\mathrm{d}t}{a\left(t\right)}=\int \frac{\mathrm{d}t}{\left(1+\alpha {t}^{2}\right)}\simeq t.\end{equation*}$$

Furthermore for cosmic times that satisfy t ≃ 0, the Hubble horizon is infinitely large and thus the primordial modes are well inside the Hubble horizon i.e. the perturbation modes satisfy k ≫ aH. As a result, the equation for the field v_k (t) becomes

$$\begin{equation}\frac{{\mathrm{d}}^{2}{v}_{k}}{\mathrm{d}{\tau }^{2}}+{k}^{2}{v}_{k}=0\enspace \enspace .\end{equation} \tag{ 30 }$$

This is the equation of a simple harmonic oscillator with time-independent frequency. For this case, the solution for the minimum energy state is given by ${v}_{k}=\frac{1}{\sqrt{2k}}\enspace {\mathrm{e}}^{-\mathrm{i}k\tau }$. Hence we impose the initial condition as shown in equation (29). This justifies the choice of the adiabatic vacuum state in the present context. On other hand, the function $H\left[-1+\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}},\sqrt{\frac{\alpha {V}_{1}}{{U}_{1}}}t\right]$ has a limiting value for t → 0 and is given by,

$$\begin{equation*}{\left.H\left[-1+\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}},\sqrt{\frac{\alpha {V}_{1}}{{U}_{1}}}\enspace t\right]\right\vert }_{t\to 0}=\frac{{2}^{\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}}}\sqrt{\pi }}{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{1}}{4\alpha {V}_{1}}\right)}\end{equation*}$$

with Γ(z) symbolizes the 'gamma function'. Therefore, by considering the adiabatic vacuum as the initial condition of ${\mathfrak{R}}_{k}\left(t\right)$, the integration constant turns out to be,

$$\begin{equation}{b}_{1}\left(k\right)=\frac{1}{z\left(t\to 0\right)}\left[\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{1}}{4\alpha {V}_{1}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}}}}\right]=\frac{1}{\sqrt{{U}_{1}}}\left[\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{1}}{4\alpha {V}_{1}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}}}}\right]\end{equation} \tag{ 31 }$$

where in the second equality, we used $z\left(t\to 0\right)=\frac{1}{\sqrt{{U}_{1}}}$ from equation (26). Thus as a whole, the solution for the scalar perturbation variable becomes,

$$\begin{align}\hfill {\mathfrak{R}}_{k}\left(t\right)& =\left(\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{1}}{4\alpha {V}_{1}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}}}\sqrt{{U}_{1}}}\right)\enspace {\mathrm{e}}^{-\frac{{V}_{1}}{{U}_{1}}\alpha {t}^{2}}\enspace H\left[-1+\frac{{k}^{2}{U}_{1}}{2\alpha {V}_{1}},\sqrt{\frac{\alpha {V}_{1}}{{U}_{1}}}\enspace t\right]\hfill \\ \hfill & =\left(\frac{2\kappa \enspace {\Gamma}\left(1-\frac{3{k}^{2}}{8\alpha }\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{3{k}^{2}}{4\alpha }}\sqrt{12\alpha D}}\right)\enspace {\mathrm{e}}^{\left[\frac{1}{4}-\frac{2}{3}\alpha {t}^{2}\right]}\enspace H\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]\hfill \end{align} \tag{ 32 }$$

where we have used ${U}_{1}=\frac{72\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$ and ${V}_{1}=\frac{48\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$ in the second equality. Consequently, the above solution of ${\mathfrak{R}}_{k}\left(t\right)$ immediately leads to the scalar power spectrum for kth mode as follows,

$$\begin{align}\hfill {P}_{\k{e}}\left(k,t\right)& =\frac{{k}^{3}}{2{\pi }^{2}}\enspace {\left\vert {\mathfrak{R}}_{k}\left(t\right)\right\vert }^{2}\hfill \\ \hfill & =\frac{{k}^{2}}{12\alpha D{\pi }^{3}}\enspace \frac{{\left(\kappa \enspace {\Gamma}\left(1-\frac{3{k}^{2}}{8\alpha }\right)\right)}^{2}}{\enspace {2}^{\frac{3{k}^{2}}{2\alpha }}\enspace }\enspace {\mathrm{e}}^{\left[\frac{1}{2}-\frac{4}{3}\alpha {t}^{2}\right]}\enspace \hfill \\ \hfill & \quad {\times}{\left\{H\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]\right\}}^{2}\hfill \end{align} \tag{ 33 }$$

Using the horizon crossing relation k = aH, the scalar power spectrum can be expressed in terms of the quantities at horizon crossing as,

$$\begin{align}\hfill {\left.{P}_{\k{e}}\left(k,t\right)\right\vert }_{\mathrm{h}.\mathrm{c}}& =\frac{\alpha {t}_{\mathrm{h}}^{2}}{3D{\pi }^{3}}\enspace \frac{{\left(\kappa \enspace {\Gamma}\left(1-\frac{3\alpha {t}_{\mathrm{h}}^{2}}{2}\right)\right)}^{2}}{\enspace {2}^{6\alpha {t}_{\mathrm{h}}^{2}}\enspace }\enspace {\mathrm{e}}^{\left[\frac{1}{2}-\frac{4}{3}\alpha {t}_{\mathrm{h}}^{2}\right]}\enspace \hfill \\ \hfill & \quad {\times}{\left\{H\left[-1+3\alpha {t}_{\mathrm{h}}^{2},\sqrt{\frac{2\alpha }{3}}\enspace {t}_{\mathrm{h}}\right]\right\}}^{2}\hfill \end{align} \tag{ 34 }$$

where t_h is the horizon crossing time. With equation (33), we can determine the spectral index of the primordial curvature perturbations. However before proceeding to calculate the spectral index, we will perform first the tensor perturbation, which is necessary for evaluating the tensor-to-scalar ratio.

4.2. Tensor perturbation

Let us now focus on the tensor perturbations, and the tensor perturbation on the FRW metric background is defined as follows,

$$\begin{equation}\mathrm{d}{s}^{2}=-\mathrm{d}{t}^{2}+a{\left(t\right)}^{2}\left({\delta }_{ij}+{h}_{ij}\right)\mathrm{d}{x}^{i}\enspace \mathrm{d}{x}^{j},\end{equation} \tag{ 35 }$$

where ${h}_{ij}\left(t, \overrightarrow {x}\right)$ is the tensor perturbation. The variable ${h}_{ij}\left(t, \overrightarrow {x}\right)$ is itself a gauge invariant quantity, and the tensor perturbed action up to quadratic order is given by,

$$\begin{equation}\delta {S}_{\mathrm{h}}=\int \mathrm{d}t\enspace {\mathrm{d}}^{3} \overrightarrow {x}a\left(t\right){z}_{T}{\left(t\right)}^{2}\left[{\dot {\mathrm{h}}}_{ij}{\dot {\mathrm{h}}}^{ij}-\frac{1}{{a}^{2}}{\left({\partial }_{l}{h}_{ij}\right)}^{2}\right],\end{equation} \tag{ 36 }$$

where z_T(t) has the following form [104],

$$\begin{equation}{z}_{\mathrm{T}}\left(t\right)=\frac{a\left(t\right)}{\kappa }\sqrt{{F}^{\prime }\left(R\right)},\end{equation} \tag{ 37 }$$

Therefore the speed of the tensor perturbation is ${c}_{\mathrm{T}}^{2}=1$ i.e. the gravitational waves propagate with the speed of light which is unity in the natural units.

Coming back to equation (37), the stability of the tensor perturbation will be guaranteed by the condition F'(R) > 0—the same condition by which the scalar perturbation is ensured to be stable. Recall, in the previous section, we showed the positivity of F'(R) in equation (24) and thus the tensor perturbation is indeed stable in the present context. The action (36) leads to the equation of the tensor perturbed variable as,

$$\begin{equation}\frac{1}{a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right)}\frac{\mathrm{d}}{\mathrm{d}t}\left[a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right){\dot {\mathrm{h}}}_{ij}\right]-\frac{1}{{a}^{2}}{\partial }_{\mathrm{l}}{\partial }^{\mathrm{l}}{h}_{ij}=0\end{equation} \tag{ 38 }$$

The Fourier transformed tensor perturbation variable is defined as ${h}_{ij}\left(t, \overrightarrow {x}\right)=\int \mathrm{d} \overrightarrow {k}\enspace {\sum }_{\gamma }{{\epsilon}}_{ij}^{\left(\gamma \right)}\enspace {h}_{\left(\gamma \right)}\left( \overrightarrow {k},t\right){\text{e}}^{\text{i} \overrightarrow {k}. \overrightarrow {x}}$, where γ = '+' and γ = '×' represent two polarization modes. Moreover ${{\epsilon}}_{ij}^{\left(\gamma \right)}$ are the polarization tensors satisfying ${{\epsilon}}_{ii}^{\left(\gamma \right)}={k}^{i}{{\epsilon}}_{ij}^{\left(\gamma \right)}=0$. In terms of the Fourier transformed tensor variable h_k (t), equation (38) can be expressed as,

$$\begin{equation}\frac{1}{a\left(t\right){z}_{T}^{2}\left(t\right)}\frac{\mathrm{d}}{\mathrm{d}t}\left[a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right){\dot {\mathrm{h}}}_{k}\right]+\frac{{k}^{2}}{{a}^{2}}{h}_{k}\left(t\right)=0\end{equation} \tag{ 39 }$$

The two polarization modes obey the same equation (39) and thus we omit the polarization index. Moreover, both the polarization modes even follow the same initial condition and that is why they have the same solution. So in the following, what we will do in the context of tensor perturbation is applicable for both the polarization modes. Finally, in the expression of the tensor power spectrum, we will introduce a multiplicative factor 2 due to the contribution from both the polarization modes. Using the expression of F'(R) from equation (24), we determine $a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right)$ as,

$$\begin{equation}a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right)=\frac{12\alpha D}{{\kappa }^{2}\sqrt{e}}\left[1+2\alpha {t}^{2}\right]={U}_{2}+{V}_{2}\enspace \alpha {t}^{2}\end{equation} \tag{ 40 }$$

with ${U}_{2}=\frac{12\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$ and ${V}_{2}=\frac{24\alpha D}{{\kappa }^{2}\sqrt{\mathrm{e}}}$. Plugging back this expression of $a\left(t\right){z}_{\mathrm{T}}^{2}\left(t\right)$ into equation (39) and after some algebra, we get the following equation,

$$\begin{equation}{U}_{2}{\ddot {\mathrm{h}}}_{k}+2\alpha {V}_{2}{\dot {\mathrm{h}}}_{k}\enspace t+{k}^{2}{U}_{2}{h}_{k}\left(t\right)=0\end{equation} \tag{ 41 }$$

at leading order in t, which is a viable consideration because the perturbation modes generate near the bouncing phase in the present scenario. Solving equation (41) for h_k (t), we get,

$$\begin{equation}{h}_{k}\left(t\right)={b}_{2}\left(k\right)\enspace {\mathrm{e}}^{-\frac{{V}_{2}}{{U}_{2}}\alpha {t}^{2}}\enspace H\left[-1+\frac{{k}^{2}}{2\alpha }\frac{{U}_{2}}{{V}_{2}},\sqrt{\frac{{V}_{2}}{{U}_{2}}\alpha }\enspace t\right]\end{equation} \tag{ 42 }$$

where b₂(k) is an integration constant and can be determined from an initial condition. As an initial condition, we consider that the tensor perturbation field starts from the adiabatic vacuum, more precisely the initial configuration is given by, ${\mathrm{lim}}_{t\to 0}\left[{z}_{\mathrm{T}}\left(t\right){h}_{k}\left(t\right)\right]=\frac{1}{\sqrt{2k}}$. This immediately leads to the expression of b₂(k) as,

$$\begin{equation}{b}_{2}\left(k\right)=\frac{1}{{z}_{\mathrm{T}}\left(t\to 0\right)}\left[\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{2}}{4\alpha {V}_{2}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{2}}{2\alpha {V}_{2}}}}\right]=\frac{1}{\sqrt{{U}_{2}}}\left[\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{2}}{4\alpha {V}_{2}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{2}}{2\alpha {V}_{2}}}}\right]\enspace \enspace \enspace .\end{equation} \tag{ 43 }$$

In the second equality of the above equation, we use ${z}_{\mathrm{T}}\left(t\to 0\right)=1/\sqrt{{U}_{2}}$ from equation (40). Plugging this expression of b₂(k) into equation (42) yields the final solution of h_k (t) as follows,

$$\begin{align}\hfill {h}_{k}\left(t\right)& =\left(\frac{2{\Gamma}\left(1-\frac{{k}^{2}{U}_{2}}{4\alpha {V}_{2}}\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}{U}_{2}}{2\alpha {V}_{2}}}\sqrt{{U}_{2}}}\right)\enspace {\mathrm{e}}^{-\frac{{V}_{2}}{{U}_{2}}\alpha {t}^{2}}\enspace H\left[-1+\frac{{k}^{2}{U}_{2}}{2\alpha {V}_{2}},\sqrt{\frac{{V}_{2}}{{U}_{2}}\alpha }\enspace t\right]\hfill \\ \hfill & =\left(\frac{2\kappa \enspace {\Gamma}\left(1-\frac{{k}^{2}}{8\alpha }\right)}{\sqrt{2\pi k}\enspace {2}^{\frac{{k}^{2}}{4\alpha }}\sqrt{12\alpha D}}\right)\enspace {\mathrm{e}}^{\left[\frac{1}{4}-2\alpha {t}^{2}\right]}\enspace H\left[-1+\frac{{k}^{2}}{4\alpha },\sqrt{2\alpha }\enspace t\right]\hfill \end{align} \tag{ 44 }$$

where we have used ${U}_{2}=\frac{12\alpha D}{{\kappa }^{2}\sqrt{e}}$ and ${V}_{2}=\frac{24\alpha D}{{\kappa }^{2}\sqrt{e}}$. Equation (44) represents the solution of the tensor perturbation for both the polarization modes. The solution of h_k (t) immediately leads to the tensor power spectrum as,

$$\begin{align}\hfill {P}_{\mathrm{h}}\left(k,t\right)& =\frac{{k}^{3}}{2{\pi }^{2}}\enspace \sum _{\gamma }{\left\vert {h}_{k}^{\left(\gamma \right)}\left(t\right)\right\vert }^{2}\hfill \\ \hfill & =\frac{{k}^{2}}{6\alpha D{\pi }^{3}}\enspace \frac{{\left(\kappa \enspace {\Gamma}\left(1-\frac{{k}^{2}}{8\alpha }\right)\right)}^{2}}{\enspace {2}^{\frac{{k}^{2}}{2\alpha }}}{\mathrm{e}}^{\left[\frac{1}{2}-4\alpha {t}^{2}\right]}\enspace {\left\{H\left[-1+\frac{{k}^{2}}{4\alpha },\sqrt{2\alpha }\enspace t\right]\right\}}^{2}\hfill \end{align} \tag{ 45 }$$

It may be noticed that γ = '+' and γ = '×' modes contribute equally to the power spectrum, as expected because their solutions behave similarly. At the horizon crossing, the tensor power spectrum turns out to be,

$$\begin{align}\hfill {\left.{P}_{\mathrm{h}}\left(k,t\right)\right\vert }_{\mathrm{h}.\mathrm{c}}& =\frac{2\alpha {t}_{\mathrm{h}}^{2}}{3D{\pi }^{3}}\enspace \frac{{\left(\kappa \enspace {\Gamma}\left(1-\frac{\alpha {t}_{\mathrm{h}}^{2}}{2}\right)\right)}^{2}}{\enspace {2}^{2\alpha {t}_{\mathrm{h}}^{2}}}\enspace {\mathrm{e}}^{\left[\frac{1}{2}-4\alpha {t}_{\mathrm{h}}^{2}\right]}\enspace \hfill \\ \hfill & \quad {\times}{\left\{H\left[-1+\alpha {t}_{\mathrm{h}}^{2},\sqrt{2\alpha }\enspace {t}_{\mathrm{h}}\right]\right\}}^{2}\hfill \end{align} \tag{ 46 }$$

Now we can explicitly confront the model at hand with the latest Planck observational data [107], so we shall calculate the spectral index of the primordial curvature perturbations n_s and the tensor-to-scalar ratio r, which are defined as follows,

$$\begin{equation}{n}_{\mathrm{s}}-1={\left.\frac{\partial \enspace \mathrm{ln}\enspace {P}_{\k{e}}}{\partial \enspace \mathrm{ln}\enspace k}\right\vert }_{\mathrm{H}.\mathrm{C}},\qquad r={\left.\frac{{P}_{\mathrm{h}}\left(k,t\right)}{{P}_{\k{e}}\left(k,t\right)}\right\vert }_{\mathrm{H}.\mathrm{C}}\end{equation} \tag{ 47 }$$

As evident from these expressions, n_s and r are evaluated at the time of the first horizon exit near the bouncing point, for positive times (symbolized by 'H.C' in the above equations), when k = aH i.e. when the mode k crosses the Hubble horizon. Using equation (33), we determine $\frac{\partial \enspace \mathrm{ln}\enspace {P}_{\k{e}}}{\partial \enspace \mathrm{ln}\enspace k}$ as follows,

$$\begin{align}\hfill \frac{\partial \enspace \mathrm{ln}\enspace {P}_{\k{e}}}{\partial \enspace \mathrm{ln}\enspace k}& =2-\frac{3{k}^{2}}{2\alpha }\enspace {\psi }^{\left(0\right)}\left[1-\frac{3{k}^{2}}{8\alpha }\right]-\frac{3{k}^{2}}{\alpha }\left(\mathrm{ln}2\right)\hfill \\ \hfill & \quad +\frac{3{k}^{2}}{\alpha }\left\{\frac{{H}^{\left(1,0\right)}\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]}{H\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]}\right\}\hfill \end{align} \tag{ 48 }$$

where ψ⁽⁰⁾[z] is the zeroth order Polygamma function or equivalently the digamma function and H^(1,0)[z₁, z₂] is the derivative of H[z₁, z₂] with respect to its first argument. Moreover the following expressions were used to calculate equation (48),

$$\begin{align}\hfill \frac{\partial }{\partial k}\left\{{\Gamma}\left(1-\frac{3{k}^{2}}{8\alpha }\right)\right\}& =-\frac{3k}{4\alpha }\enspace {\Gamma}\left(1-\frac{3{k}^{2}}{8\alpha }\right)\enspace {\psi }^{\left(0\right)}\left(0,1-\frac{3{k}^{2}}{8\alpha }\right)\hfill \\ \hfill \frac{\partial }{\partial k}\left\{H\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]\right\}& =\frac{3k}{2\alpha }\enspace {H}^{\left(1,0\right)}\left[-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right]\hfill \end{align} \tag{ 49 }$$

Thereby equation (48) immediately leads to the spectral index as,

$$\begin{align}\hfill {n}_{\mathrm{s}}& =\left[3-\frac{3{k}^{2}}{2\alpha }\enspace {\psi }^{\left(0\right)}\left(0,1-\frac{3{k}^{2}}{8\alpha }\right)-\frac{3{k}^{2}}{\alpha }\left(\mathrm{ln}2\right)\right.\hfill \\ \hfill & \quad +{\left.\frac{3{k}^{2}}{2\alpha }\left\{\frac{{H}^{\left(1,0\right)}\left(-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right)}{H\left(-1+\frac{3{k}^{2}}{4\alpha },\sqrt{\frac{2\alpha }{3}}\enspace t\right)}\right\}\right]}_{\mathrm{H}.\mathrm{C}}\hfill \end{align} \tag{ 50 }$$

As mentioned earlier, the perturbation modes are generated and also cross the horizon near the bounce. Thus we can safely use the near-bounce scale factor in the horizon crossing condition to determine k = aH = 2αt_h (where t_h is the horizon crossing time). Using this relation, equation (50) turns out to be,

$$\begin{align}\hfill {n}_{\mathrm{s}}& =3-6\alpha {t}_{\mathrm{h}}^{2}\enspace {\psi }^{\left(0\right)}\left(0,1-\frac{3}{2}\enspace \alpha {t}_{\mathrm{h}}^{2}\right)-12\alpha {t}_{\mathrm{h}}^{2}\left(\mathrm{ln}\enspace 2\right)\hfill \\ \hfill & \quad +12\alpha {t}_{\mathrm{h}}^{2}\left\{\frac{{H}^{\left(1,0\right)}\left(-1+3\alpha {t}_{\mathrm{h}}^{2},\sqrt{\frac{2\alpha }{3}}\enspace {t}_{\mathrm{h}}\right)}{H\left(-1+3\alpha {t}_{\mathrm{h}}^{2},\sqrt{\frac{2\alpha }{3}}\enspace {t}_{\mathrm{h}}\right)}\right\}\hfill \end{align} \tag{ 51 }$$

Furthermore, the tensor-to-scalar ratio is given by,

$$\begin{equation}r={\left\vert \frac{{h}_{k}\left(t\right)}{{\mathfrak{R}}_{k}\left(t\right)}\right\vert }_{k=a\left({t}_{\mathrm{h}}\right)H\left({t}_{\mathrm{h}}\right)}^{2}\end{equation} \tag{ 52 }$$

where the solutions of h_k (t) and ${\mathfrak{R}}_{k}\left(t\right)$ are shown in equations (44) and (32) respectively. Equations (51) and (52) clearly indicate that n_s and r depend on the dimensionless parameter $\alpha {t}_{\mathrm{h}}^{2}$ which is further connected to the Ricci scalar at horizon crossing by $\alpha {t}_{\mathrm{h}}^{2}=\left(\frac{{R}_{\mathrm{h}}}{12\alpha }-1\right)$. Thereby, we can argue that the observable quantities n_s and r depend on R_h/α. With this information, we now directly confront the theoretical expectations of spectral index and tensor-to-scalar ratio with the Planck 2018 constraints [107]. In figure 2 we present the estimated spectral index and tensor-to-scalar ratio of the present scenario for three choices of $\frac{{R}_{\mathrm{h}}}{\alpha }$, on top of the 1σ and 2σ contours of the Planck 2018 results [107]. As we observe, the agreement with observations is efficient, and in particular well inside the 1σ region. At this stage it is worth mentioning that in an F(R) gravity theory, the MBS, in which case the perturbations are generated far away from the bouncing point deeply in the contracting regime, is not consistent with the Planck results as shown in [64]. However, here we demonstrate that a F(R) gravity model indeed leads to a viable bouncing model where the primordial perturbations are generated near the bounce. Thereby, we can argue that the viability of a F(R) model (with respect to the Planck constraints) is interrelated with the generation era of the perturbations modes. Furthermore the scalar perturbation amplitude (A_s) is constrained by $\mathrm{ln}\left[1{0}^{10}{A}_{\mathrm{s}}\right]=3.044{\pm}0.014$ due to the Planck results [107]. Equation (33) along with the consideration of the integration constant $D=\frac{1}{\alpha }$ (recall D has mass dimension [−2]) depict that the scalar perturbation amplitude depends on the dimensionless parameters $\frac{{R}_{\mathrm{h}}}{\alpha }$ and ακ². We may choose $\frac{{R}_{\mathrm{h}}}{\alpha }=16$ (which is within the range that makes the n_s and r compatible with the Planck results) and consequently the scalar perturbation amplitude in the present context becomes ${A}_{\mathrm{s}}=\frac{1}{20{\pi }^{3}}\alpha {\kappa }^{2}$. Therefore the theoretical expectation of A_s will match with the Planck observations provided ακ² lies within ακ² = [1.283 × 10⁻⁷, 1.320 × 10⁻⁷]. Thus as a whole, the observable quantities n_s, r and A_s are simultaneously compatible with the Planck constraints for the parameter ranges: $14{\leqslant}\frac{{R}_{\mathrm{h}}}{\alpha }{\leqslant}20$ and ακ² = [1.283 × 10⁻⁷, 1.320 × 10⁻⁷] respectively. However the viable range of ακ² depends on the consideration D = 1/α, i.e. a different D will eventually lead to a different range of viability of the parameter ακ². As an example, for $D=\frac{1}{4\alpha }$, the scalar perturbation amplitude becomes ${A}_{\mathrm{s}}=\frac{1}{5{\pi }^{3}}\alpha {\kappa }^{2}$ while n_s and r have the same form as given in equations (51) and (52) respectively and therefore the parameters should lie within $14{\leqslant}\frac{{R}_{\mathrm{h}}}{\alpha }{\leqslant}20$ and ακ² = [3.208 × 10⁻⁸, 3.300 × 10⁻⁸] in order to make the theoretical values of the observable quantities compatible with the latest Planck 2018 results. Such parametric ranges make the horizon crossing Ricci scalar of the order R_h ∼ 10⁻⁸/κ² = 10²⁸ GeV².

Zoom In Zoom Out Reset image size

Thereby, the viable ranges of the early and late stage parameters are given by $14\lesssim \frac{{R}_{\mathrm{h}}}{\alpha }\lesssim 20$, ακ² = [3.208 × 10⁻⁸, 3.300 × 10⁻⁸], $\beta =3.9{8}_{-2.46}^{+\infty }$ and Λ = 1.2 × 10⁻⁸⁴ GeV² respectively. As a result, we can interpolate the scale factor or the Hubble radius in the two transition eras (i.e. from the bounce to MD and from the MD to DE epoch), leading to a smooth evolution of the Hubble radius for a large range of cosmic time.

Before moving to the next section we would like to mention that in the present context, the non-singular bounce Universe is time symmetric and the initial adiabatic vacuum condition is set at the bouncing point, which seems that the model is more similar to a scenario of the creation of the Universe from nothing. The spontaneous birth of the Universe from 'nothing' has been discussed earlier in the context of quantum cosmology where the Universe is described by a wave function satisfying the well known Wheeler–DeWitt equation [108]. With the development of quantum cosmology theory, it has been suggested that the Universe can be created spontaneously from nothing, where 'nothing' means there is neither matter nor space or time, and the problem of singularity can be avoided naturally [109, 110]. In particular, the author of [109] showed that a quantum Universe of zero size, or, indeed, the cosmological 'nothing' may quantum mechanically tunnel to a Universe of a finite size. However the argument of [109] was established for a closed type Universe, unlike to the present bounce scenario which has been presented for a flat FRW Universe.

Given the structure of the scale factor in the early and late stages of the Universe (as a_b(t) = 1 + αt² in equation (7) and ${a}_{\mathrm{l}}\left(t\right)={\mathrm{e}}^{\int {H}_{\mathrm{l}}\left(t\right)\mathrm{d}t}$ respectively, where the H_l(t) is given in figure 1), we would like to provide a complete picture by considering the intermediate region to be a MD epoch, and moreover the transitions from the bounce to MD and from the MD to the DE era will be obtained by the method of numerical interpolation. Due to complicated nature of the equations governing the evolution of the scale factor in a general context, we will determine the interpolating function using numerical techniques and will illustrate the same. Let us briefly point out the methods one may use in order to generate such interpolating solutions. In the transition regions (i.e. from the bounce-MD and from the MD–DE), one approximates the behavior of the physical quantity of interest (e.g., the scale factor a(t) or the Hubble radius ${r}_{\mathrm{h}}=\frac{1}{aH}$) by a polynomial function of time, with degree of the polynomial kept arbitrary. Then, in the early bounce epoch one uses the given behavior of the desired physical quantity (for example for the scale factor, the behavior near the bounce is a_b(t) = 1 + αt²) to generate numerical estimates of the respective quantity at various time instants till the description is reliable. Similar numerical estimations are being made at the intermediate MD and at the late stage as well. With these sets of data and the polynomial function one can use any standard interpolation software package to end up getting the desired plots. The structure of the plot, of course, depends on the degree of the polynomial and desired accuracy level. All the plots in this paper are for an accuracy level of $\mathcal{O}\left(1{0}^{-8}\right)$. Since our main aim is to merge certain cosmological epochs of the Universe, a good physical quantity to start with is the Hubble radius (r_h(t)) rather than the scale factor, because the accelerating or decelerating stage of the Universe is easily realized by the decreasing or increasing behavior of the Hubble radius (with respect to cosmic time) respectively. Following equation (16), we may construct the Hubble radius as,

$$\begin{align}\hfill {r}_{\mathrm{h}}\left(t\right)& ={\left[{\dot {a}}_{\mathrm{b}}\left(t\right)\right]}^{-1}=\frac{1}{2\alpha t},\hfill & \hfill & \mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\enspace \mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}\hfill & \hfill & t\simeq {0}^{+}\hfill \\ \hfill {r}_{\mathrm{h}}\left(t\right)& ={\left[{\dot {a}}_{I}\left(t\right)\right]}^{-1}=\frac{3t}{2}{\left(\frac{T}{t}\right)}^{2/3},\hfill & \hfill & \mathrm{i}\mathrm{n}\enspace \mathrm{t}\mathrm{h}\mathrm{e}\enspace \mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{m}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}\hfill & \hfill & {0}^{+}\lesssim t\lesssim 5{\times}1{0}^{17}\enspace \mathrm{sec} .\hfill \\ \hfill {r}_{\mathrm{h}}\left(t\right)& ={\left[{\dot {a}}_{l}\left(t\right)\right]}^{-1}\quad \left(\mathrm{s}\mathrm{e}\mathrm{e}\enspace \mathrm{fi}\mathrm{g}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{1}\right),\hfill & \hfill & \mathrm{d}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\enspace \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\enspace \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e},\enspace \mathrm{f}\mathrm{o}\mathrm{r}\hfill & \hfill & t\enspace \gtrsim \enspace 6{\times}1{0}^{17}\enspace \mathrm{sec} .\hfill \end{align} \tag{ 53 }$$

and then with the procedure described above, we numerically interpolate the Hubble radius in the two transitions i.e. from the bounce to MD and from the MD to DE era respectively. Moreover the Hubble radius in the contracting regime can also be obtained by assuming a symmetric behavior of $\left\vert {r}_{\mathrm{h}}\left(t\right)\right\vert$ around the bounce. However, the details of the interpolation of the curve connecting the early bounce to the late accelerating stage is an artifact of the procedure followed and admits possible variations depending on the process of interpolation by numerical techniques. However such indeterminacy in determining the interpolating function would not affect our main conclusions in the present context. Thus, having explained the details of the interpolating procedure, we now turn to the corresponding implications and present the variations of all the relevant parameters with time.

As mentioned earlier, we start with the Hubble radius (${r}_{\mathrm{h}}=\frac{1}{\dot {a}}$), in particular, taking the reduced Planck mass to be M_Pl = 10¹⁸ GeV and using the aforementioned construction of the Hubble radius, we obtain r_h(t) as a function of time in the expanding phase of the Universe. Further the Hubble radius in the contracting regime is obtained by considering a symmetric character of $\left\vert {r}_{\mathrm{h}}\left(t\right)\right\vert$ in the both sides of the bounce. As a result, we get the Hubble radius for −9 × 10¹⁷ ⩽ t ⩽ 9 × 10¹⁷ s, which is presented in the left part of figure 3, while the right part is a zoomed-in version of the left one near t = 0. The x-axis of the figure 3 corresponds to a 'rescaled' time coordinate obtained as $\frac{t}{{t}_{\mathrm{s}}}$ (with t_s = 10¹⁶ s) which is dimensionless, while the y axis represents the rescaled Hubble radius ${\bar{r}}_{\mathrm{h}}=1{0}^{-5}{r}_{\mathrm{h}}$ which is in the unit of 's'. Figure 3 clearly demonstrates that the Hubble radius depicts a bounce at t = 0 (as expected, because it is constructed from ${r}_{\mathrm{h}}^{\left(\mathrm{b}\right)}\left(t\right)={\left({a}_{\mathrm{b}}{H}_{\mathrm{b}}\right)}^{-1}$ near t = 0), then it increases monotonically with time up to $\frac{t}{{t}_{\mathrm{s}}}\simeq 50$ (or t = 5 × 10¹⁷ s) leading to a decelerating phase of the Universe and finally after t ≳ 5 × 10¹⁷ s, the Hubble radius starts to decrease, which in turn indicates the DE epoch of the Universe. This may provide a merging of a non-singular bounce to MD epoch, followed by a late time DE era. It may be noticed that the Hubble radius goes to zero asymptotically at both 'sides' of the bounce, which in turn confirms the fact that the relevant primordial perturbation modes generate near the bouncing point as we have considered in section 4 during the calculation of scalar and tensor perturbations. The generation era of perturbation modes in the present context makes the model different than the usual MBS where the perturbation modes generate in the distant past far away from the bounce. Having this Hubble radius, we are now going to solve the scale factor and Hubble parameter, however numerically and then by using the numerical solution of H(t) we reconstruct the form of F(R) from the gravitational equation (4). The scale factor can be obtained from $\dot {a}\left(t\right)=\frac{1}{{r}_{\mathrm{h}}\left(t\right)}$ which is a first order differential equation and thus the solution of the same requires one boundary condition. We choose a(0) = 1, because we want the scale factor to behave like a_b(t) near t = 0 and recall a_b(0) = 1. With such initial condition along with the r_h(t) given in figure 3, we solve a(t) numerically, which is presented in figure 4 where once again, the x axis is rescaled by t_r = t/t_s.

Zoom In Zoom Out Reset image size

It is evident that the scale factor is smooth everywhere and the graph in the inset shows that the scale factor acquires a minimum at t = 0. Thus we can argue that the Universe transits from a contracting phase to an expanding one through a non-singular bounce at t = 0 and thus the Universe evolution becomes free of the initial singularity. The above numerical solutions of the scale factor can be immediately differentiated providing the Hubble parameter as a function of time, shown in figure 5. The Hubble parameter becomes zero and shows an increasing behavior with time (i.e. $\dot {H}{ >}0$) at t = 0, however this is expected as the scale factor itself acquires a minimum at t = 0 i.e. a bounce. Moreover, the Hubble parameter acquires a maximum very near to bounce and then decreases rapidly, and continued to decrease until the late-time era. A zoomed-in version of the late time behavior of H(t) is shown in the right plot of figure 5, which clearly demonstrates that at the late epoch i.e. for t ≳ 7 × 10¹⁷ s, the Hubble parameter becomes of the order ∼10⁻¹⁸ s⁻¹. The occurrence of the maximum of H(t) near the bounce in the present context is in agreement with [82] where the authors explained the merging of bounce with late-time acceleration from a different viewpoint, namely from the holonomy generalizations of deformed MBS. In this regard, it may be mentioned that a deformed MBS is represented by an asymmetric scale factor where the Hubble radius diverges at t → −∞ and asymptotically goes to zero at t → +∞. Thus the primordial perturbation modes in a deformed matter bounce model generate deeply in the contracting regime far away from the bounce, unlike to our present case where the Hubble radius goes to zero asymptotically at both 'sides' of the bounce (see figure 3) leading to the generation of the perturbation modes near the bouncing regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the Hubble parameter at hand, our next task is to solve the gravitational equation (4) to reconstruct the form of F(R) which can realize such evolution of the Hubble parameter. For this purpose, we need two boundary conditions as equation (4) is a second order differential equation with respect to the Ricci scalar. As mentioned earlier, we use the forms of F_b(R) and F_l(R) obtained in section 3 for such boundary conditions and thus the boundary conditions are given by: F(R = 0) = F_l(0) and F(R_b = 12α) = F_b(12α) respectively (recall, R_b is the Ricci scalar at the bounce). As a result, we solve equation (4) numerically and the form of F(R) we get is depicted in the left part of figure 6. Actually the F(R) is demonstrated by the red curve, while the yellow one represents the Einstein gravity. Before going into the discussion regarding the behavior of the F(R) gravity, let us briefly comment about the range of the 'rescaled Ricci scalar' that we take along the x-axis in figure 6. Recall, the Ricci scalar near the bounce gets the value R(t → 0) ≃ 12α where the parameter α lies within [3.208 × 10⁻⁸/κ², 3.300 × 10⁻⁸/κ²] to make the primordial observable quantities compatible with the Planck constraints. Thereby, the Ricci scalar near the bounce becomes of the order ∼12 × 10²⁸ GeV², while in the present epoch, the scalar curvature is generally considered to take 10⁻⁸⁴ GeV². Thus in order to correctly describe the epochs from bounce to late time acceleration, we have to cover the range of R(t) from ∼12 × 10²⁸ GeV² to 10⁻⁸⁴ GeV². However in figure 6, we consider a 'rescaled Ricci scalar' given by $\bar{R}=\frac{R}{\alpha }$ along the x-axis, which is dimensionless as α has the mass dimension [+2]. Due to such rescaling, the x-axis of the plot runs from the value x_i = 10⁻¹¹² (or we may take x_i = 0) to x_f = 12. Moreover, the y-axis of figure 6 is also rescaled by $\bar{F}\left(\bar{R}\right)=F\left(R\right)/\alpha$ so that the Einstein gravity can be described by a straight line having slope of unity even in these rescaled coordinates. The left part of the figure 6 clearly demonstrates that the solution of F(R) matches with the Einstein gravity as the Ricci scalar approaches to the present value, while the F(R) seems to deviate from the usual Einstein gravity, when the scalar curvature takes larger and larger values. It is evident that near the bounce, F'(R) is positive, which in turn indicates the stability for the scalar and tensor perturbations. This finding from the numerical solution, in fact, resembles with the analytic results obtained in section 4. The comparison of F(R) gravity with the Einstein one in the low-curvature regime is explicitly demonstrated in the right part of figure 5, which reveals that the form of F(R) reconstructed in the present context is indeed comparable with Einstein gravity during the present epoch. Thus, we can safely argue that the numerical solution of F(R) obtained above indeed matches with the analytic results determined near the bounce and at the DE era. Using the form of F(R) and the Hubble parameter from figures 5 and 6 respectively, we determine the remaining bit i.e. the effective equation of state (EoS) parameter defined by,

$$\begin{equation*}{w}_{\text{eff}}=\frac{\left[\frac{f\left(R\right)}{2}-\left(3{H}^{2}+\dot {H}\right){f}^{\prime }\left(R\right)+6\left(8{H}^{2}\dot {H}+4{\dot {H}}^{2}+6H\ddot {H}+ \dddot {H}\right){f}^{{\prime\prime}}\left(R\right)+36{\left(4H\dot {H}+\ddot {H}\right)}^{2}{f}^{\prime \prime \prime }\left(R\right)\right]}{\left[-\frac{f\left(R\right)}{2}+3\left({H}^{2}+\dot {H}\right){f}^{\prime }\left(R\right)-18\left(4{H}^{2}\dot {H}+H\ddot {H}\right){f}^{{\prime\prime}}\left(R\right)\right]}\enspace \enspace ,\end{equation*}$$

as shown in figure 7, where the x-axis is the rescaled time coordinate i.e. $\frac{t}{{t}_{\mathrm{s}}}$ with t_s = 10¹⁶ s. Actually the red curve of the plot represents the w_eff for the present model while the yellow one is for the constant value $-\frac{1}{3}$ (we will keep the yellow graph to investigate the accelerating or decelerating era of the Universe). Figure 7 clearly demonstrates that at the bounce i.e. at t = 0, the EoS parameter diverges from the negative side, however this is expected because at the bounce the Hubble parameter itself becomes zero and in turn makes the ${w}_{\text{eff}}=-1-\frac{2\dot {H}}{3{H}^{2}}\to -\infty$. Then immediately after the bounce, w_eff crosses the value $-\frac{1}{3}$ leading to a transition from a bounce to a decelerating phase of the Universe, and during the deceleration epoch, the EoS parameter remains almost constant and near to zero indicating a quasi-matter dominated Universe which continues till t ≃ 5 × 10¹⁷ s, when the EoS parameter again crosses the value $-\frac{1}{3}$ and thus the Universe transits from a stage of deceleration to a stage of acceleration continued until the late-time era. Therefore, the present model may provide an unified scenario of certain cosmological epochs from bounce to late-time acceleration followed by a quasi-matter dominated epoch in the intermediate regime. Moreover, the EoS parameter during the DE era approaches to the value −0.995, as evident from the figure 7. This late time value of the w_eff is, however, in agreement with [86] where the authors considered the exponential F(R) as a DE model, similar to the present case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we consider a different F(R) DE model in comparison to the one we considered in the previous sections, in particular, we incorporate an additional logarithmic correction to the exponential F(R) gravity. Thus the form of F(R) is given by [88],

$$\begin{equation}F\left(R\right)=R-2{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)\left[1-\frac{\gamma R}{2{\Lambda}}\enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right)\right]\end{equation} \tag{ 54 }$$

where the correction over Einstein gravity is given by

$$\begin{equation*}F\left(R\right)-R=-2{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)\left[1-\frac{\gamma R}{2{\Lambda}}\enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right)\right].\end{equation*}$$

It may be observed that such correction goes to zero at R → 0 i.e. in the low curvature regime. Now in order to avoid including an implicit cosmological constant in a F(R) model or equivalently to avoid the fine tuning problem associated with the cosmological constant, the correction over the Einstein–Hilbert term of the F(R) gravity requires to vanish in the low curvature limit i.e.

$$\begin{equation}{\mathrm{lim}}_{R\to 0}\left(F\left(R\right)-R\right)=0.\end{equation} \tag{ 55 }$$

As just mentioned, this condition is satisfied for the considered F(R) model in equation (54) and thus this F(R) model can produce the vacuum solution as the Minkowski spacetime. Thereby we may argue that the model (54) is free from the problem of including an implicit cosmological constant, unlike to the scenario with the action given by $S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R-2{\Lambda}\right]$. The condition (55) for avoiding a cosmological constant is also satisfied for the exponential F(R) gravity as considered earlier in equation (12) i.e. without the logarithmic corrections. However as described in [88], the presence of the logarithmic correction, modelled by the free parameter γ by equation (54), leads to a better fit (in respect to a viable DE model) than in absence of the logarithmic correction and also better than the ΛCDM model. Hence the inclusion of the logarithmic correction provides a test for this type of F(R) theory and motivates to check how a deviation is allowed in comparison to the ΛCDM model. Note that the logarithmic correction introduces a new parameter γ i.e. for γ = 0, the F(R) model (54) becomes similar to the exponential model considered earlier in equation (12). This new F(R) model with the logarithmic term still satisfies the viability conditions (under some conditions of the free parameter) and provides an extra term in the action that evolutes smoothly along the cosmological evolution (far from the pole obviously), as addressed in the reference [88]. In particular, the logarithmic corrected exponential F(R) gravity comes as a viable DE model with respect to the Sne-Ia + BAO + H(z) + CMB data for the parametric choices: $\gamma =0.001{4}_{-0.0014}^{+0.0025}$, Λ = 1.2 × 10⁻⁸⁴ GeV² and $\beta =4.7{1}_{-2.87}^{+\infty }$ respectively. Such parametric regimes makes the F(R) model free from antigravity effects due to the condition $\gamma \enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right){< }1$ holds true. The authors of [88] clearly investigated the possible effects of the logarithmic term on the DE epoch in view of the aforementioned observational data and as a result, they found that the logarithmic correction modelled by the parameter γ leads to better fit than in the absence of logarithmic correction (i.e. the pure exponential F(R) model), obviously at the price of introducing a new parameter γ. Moreover the model (54) avoids the presence of large corrections on the Newton's law as well as the appearance of large instabilities at local systems, leading to a suitable model that recovers the well known results of general relativity (GR) at the appropriate scales. In view of this new DE model, we propose the following F(R) gravity,

$$\begin{align}\hfill F\left(R\right)& =R-2{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)\left[1-\frac{\gamma R}{2{\Lambda}}\enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right)\right]\hfill \\ \hfill & \quad +\mathrm{exp}\left\{1-\mathrm{cosh}\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)\right\}\left[\left(\frac{12\alpha D}{\sqrt{e}}-1-\gamma \enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right)\right)R\right.\hfill \\ \hfill & \quad -\left.D\sqrt{6\alpha \pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}{\left(R-12\alpha \right)}^{3/2}\enspace \mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]\right]\hfill \end{align} \tag{ 56 }$$

which can merge the bounce and late-time acceleration as we will show in the following, where the late time behavior is described by the F(R) of equation (54) and the bouncing scenario is still depicted by the F_b(R) determined in equation (11). In the above expression, the R_b denotes the scalar curvature at the bounce i.e. R_b = 12α ∼ 10²⁸ GeV² (see section 4 for the estimation of α). We introduce an additional cosine hyperbolic factor (let us call it 'fixing factor') with the second term of the right-hand side of equation (56) in order to avoid the effects of the second term during the late-time era. With the above F(R), the demonstration of Universe's evolution at different curvature regimes goes as follows: (1) for R ≫ Λ when the Ricci scalar is given by R(t) = 12α + 12α² t² (i.e. near the bounce), the cosine hyperbolic term behaves as,

$$\begin{align*}\hfill \mathrm{cosh}\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)& =\mathrm{cosh}\left(\frac{{\alpha }^{2}{t}^{4}}{\left(1+\alpha {t}^{2}\right)}\right)=\frac{1}{2}\left[{\mathrm{e}}^{\frac{{\alpha }^{2}{t}^{4}}{\left(1+\alpha {t}^{2}\right)}}+{\mathrm{e}}^{-\frac{{\alpha }^{2}{t}^{4}}{\left(1+\alpha {t}^{2}\right)}}\right]\hfill \\ \hfill & =1+\frac{{\alpha }^{4}{t}^{8}}{2{\left(1+\alpha {t}^{2}\right)}^{2}}+\cdots \cdots \cdots =1+\mathcal{O}\left({t}^{8}\right)\hfill \end{align*}$$

and therefore the fixing factor can be expressed by ${\mathrm{e}}^{-\mathcal{O}\left({t}^{8}\right)}=1-\mathcal{O}\left({t}^{8}\right)$. However, recall that the near bounce quantities are determined up to $\mathcal{O}\left({t}^{2}\right)$ and thus the factor ${\mathrm{e}}^{-\mathcal{O}\left({t}^{8}\right)}$ can be approximated to unity near R ≃ R_b. As a result along with the fact R_b ≫ Λ, the F(R) of equation (56) turns out to be,

$$\begin{equation}F\left(R\right)=\frac{12\alpha D}{\sqrt{e}}R-D\sqrt{6\alpha \pi }\enspace {\mathrm{e}}^{-\frac{R}{24\alpha }}{\left(R-12\alpha \right)}^{3/2}\mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]={F}_{\mathrm{b}}\left(R\right)\end{equation} \tag{ 57 }$$

in the regime R ≃ R_b, which realizes a bouncing Universe as discussed in section 3.1. However in order to ensure a bounce, it is also necessary to check F'(R) and F''(R) as they are connected to the issues of energy conditions. Being the evolution of the fixing factor as ${\mathrm{e}}^{-\mathcal{O}\left({t}^{8}\right)}=1-\mathcal{O}\left({t}^{8}\right)$ and since $\frac{\mathrm{d}R}{\mathrm{d}t}$ is proportional to t around the bounce, one can immediately write F'(R) = F_b'(R) and ${F}^{{\prime\prime}}\left(R\right)={F}_{\mathrm{b}}^{{\prime\prime}}\left(R\right)$ up to $\mathcal{O}\left({t}^{2}\right)$ near R ≃ R_b. Thereby, the F(R) we propose in equation (56) as well as its first and second derivatives match with that of F_b(R) up to $\mathcal{O}\left({t}^{2}\right)$, which clearly indicates a bouncing Universe for R ≃ R_b. Here it may be mentioned that a different fixing factor like ${\mathrm{e}}^{-\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)}$ behaves as $1+\mathcal{O}\left({t}^{4}\right)$ near the regime R ≃ R_b, unlike to our considered fixing factor where the leading order is proportional to t⁸ due to the additional cosine hyperbolic term. However in effect of fourth order term (i.e. t⁴) being the leading order one, the fixing factor ${\mathrm{e}}^{-\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)}$ makes the F'(R) different in comparison to F_b'(R) even up to $\mathcal{O}\left({t}^{2}\right)$, which in turn may violate the bounce at R ≃ R_b. In view of these arguments, we stick to the 'fixing factor' as proposed in equation (56) rather than ${\mathrm{e}}^{-\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)}$. On other hand, (2) for R ≪ R_b i.e. in the low curvature regime (or equivalently R ∼ Λ), $\frac{{R}_{\mathrm{b}}}{R}$ becomes very much larger in comparison to $\frac{R}{{R}_{\mathrm{b}}}$ and as a consequence the factor $\mathrm{exp}\left\{1-\mathrm{cosh}\left(\frac{R}{{R}_{\mathrm{b}}}+\frac{{R}_{\mathrm{b}}}{R}-2\right)\right\}$ can be approximated to $\mathrm{exp}\left[-{\mathrm{e}}^{{R}_{\mathrm{b}}/R}\right]$ which further tends to zero in the regime $\frac{R}{{R}_{\mathrm{b}}}\ll 1$. Thereby, in the low curvature regime, the form of F(R) can be written as $F\left(R\right)=R-2{\Lambda}\left(1-{\mathrm{e}}^{-\frac{\beta R}{2{\Lambda}}}\right)\left[1-\frac{\gamma R}{2{\Lambda}}\enspace \mathrm{log}\left(\frac{R}{4{\Lambda}}\right)\right]$ leading to a viable DE dominated era in respect to Sne-Ia + BAO + H(z) + CMB observations for a suitable parametric spaces, as mentioned earlier. Moreover the F(R) model (56) satisfies ${\mathrm{lim}}_{R\to 0}\left(F\left(R\right)-R\right)=0$ i.e. the correction over the Einstein–Hilbert term vanishes in the limit R → 0. This indicates that the F(R) of equation (56) is able to provide the Minkowski solution in vacuum case and thus is free from the problem of including an implicit cosmological constant. Thus, the F(R) proposed in equation (56) appropriately provides a non-singular bounce and a late acceleration near R ≃ R_b and during R ∼ Λ respectively, and hence is able to merge a bounce with DE epoch. The deceleration stage in between the bounce and the DE epochs can be obtained by the numerical interpolation, similarly as we did for the earlier case in section 5.

At this stage it deserves mentioning that we have taken the bottom-up approach to reconstruct the F(R) gravity in the present context, in which we attempt to figure out the functional form of F(R) by demanding that it should explain the observational results. However, various terms in the F(R) of equation (56) may be thought to have a fundamental origin, like—near the bouncing regime, the F(R) (denoted by F_b(R)) is given by equation (57). Thereby if we expand the Taylor series of the function $\mathrm{E}\mathrm{r}\mathrm{fi}\left[\frac{\sqrt{R-12\alpha }}{2\sqrt{6\alpha }}\right]$ (present in equation (57)) around the bouncing point i.e. around R = 12α and keeping the leading order term, then the near-bounce expression of F(R) will contain linear and quadratic power in curvature R. Such quadratic correction in the Ricci scalar plays a role for renormalizability of GR in quantum gravity. On other hand, during the late time, the F(R) behaves as of equation (54). The logarithmic correction present in the late time form of F(R) may be thought as one-loop effects of quantum gravity. Despite these arguments, it is true that due to the complicated nature of the full F(R) presented in equation (56), it is hard to realize some fundamental origin of this full form of F(R). Of course, the complete gravitational action should be defined by a fundamental theory, which, however, remains to be the open problem of modern high-energy physics. In the absence of fundamental quantum gravity, the modified gravity approach in this work is a phenomenological one that is constructed by complying with observational data.

Before concluding, we would like to mention that the present scenario merges certain cosmological epochs of the Universe, in particular, from a non-singular bounce to a MD epoch and from the MD to a late time accelerating epoch; i.e. the model is similar to a matter bounce model which is also compatible to a late DE phase of the cosmic evolution as well. However this is not the full evolution history of the Universe mainly due to the absence of the radiation era. Thereby in order to unify the entire evolutionary epochs of the Universe in the context of bouncing cosmology, one should show the unification of an early bounce with the radiation domination, followed by the photon decoupling, followed by matter domination, all the way to the late-time dominance of DE stage. The unification for the full evolution of the Universe is a longstanding problem in cosmology and we hope that our present paper may enlighten some part(s) of this bigger problem.

In the present work, we proposed a cosmological model in the context of F(R) gravity, which merges a non-singular bounce to a MD epoch and from the MD to a late time accelerating epoch; i.e. the model is similar to a generalized matter bounce model which is also compatible to a late DE dominant phase of the cosmic evolution. Using the reconstruction schemes, we reconstructed the form of F(R) near the bounce and at the late-time era respectively. However, the reconstruction techniques applied in these two eras are slightly different, in particular: near the bouncing regime, we first considered a scale factor (suitable for bounce) and then determined the form of F(R) which can realize such bouncing scale factor by using the gravitational equation of motion, while on other hand in the case of late times, we have used a reverse reconstruction procedure in comparison to that of the earlier one i.e. we started with a form of F(R) (rather than a scale factor) viable for DE model, namely the exponential F(R) gravity, and then reconstructed the corresponding Hubble parameter from the gravitational equations of motion. During such reconstructions, we got early and late stage model parameters which have been estimated from various observational constraints. The early stage parameters have been obtained from the primordial perturbations and we confronted the theoretical results of the observable quantities like the scalar spectral index, the tensor-to-scalar ratio with the latest Planck 2018 constraints. On other hand, the late stage parameters are estimated by investigating the viability of the exponential F(R) gravity as a DE model with respect to the Sne-Ia + BAO + H(z) + CMB data. Due to a late accelerating phase, the Hubble radius decreases monotonically at late times and asymptotically goes to zero, which in turn leads to the generation of the primordial perturbation near the bounce (because at that time the relevant perturbation modes are within the horizon), unlike to the usual MBS where the perturbation modes generate deeply in the contracting regime far away from the bouncing point time. Thus we have performed the perturbations near the bounce and as a result the primordial observable quantities like the spectral index for curvature perturbation, and the tensor-to-scalar ratio, are found to be simultaneously compatible with the latest Planck 2018 observations. Moreover, the scalar and tensor perturbations are stable as the condition F'(R) > 0 holds in the present context. With the early bounce and the late accelerating phase, we provided a complete evolution of the Hubble radius by considering the intermediate region to be a MD epoch and the transitions from the bounce to MD and from the MD to the DE era have been obtained by the method of numerical interpolation, where the estimated values of the model parameters are incorporated. The evolution of the Hubble radius provides the Universe's evolution for a considerable range of cosmic time, in particular, it merges a non-singular bounce to a MD epoch, followed by a late time accelerating stage. With this interpolating Hubble radius, we have numerically solved the F(R) gravitational equation to determine the form of F(R) for a wide range of cosmic time, which clearly depicts that the F(R) matches with the Einstein gravity in the low curvature regime, while it deviates from the usual Einstein gravity as the scalar curvature acquires larger and larger values. Correspondingly, the Hubble parameter and the effective EoS parameter of the Universe have been determined. As a result the EoS parameter is found to successively cross the value $-\frac{1}{3}$ twice indicating the transition of the Universe from bounce to a quasi-matter dominated phase and from the quasi-matter dominated to a late-time accelerating stage respectively. Moreover, the EoS parameter during the DE era approaches to the value −0.995. Finally, in section 6, we have proposed a different F(R) gravity form which is able to merge a non-singular bounce with a DE epoch, however the DE epoch is described by a logarithmic corrected exponential F(R) gravity, unlike to the earlier case where the DE model is depicted by the usual exponential F(R) gravity i.e. without the logarithmic corrections. Similar to the exponential F(R) model, the logarithmic generalized exponential F(R) gravity is also known to provide a viable DE model with respect to Sne-Ia + BAO + H(z) + CMB observations, in fact, the presence of logarithmic corrections leads to a better fitted model than the usual exponential F(R) model.

This work is supported by MINECO (Spain), FIS2016-76363-P (S.D.O).