Bayesian Analysis of a Generalized Starobinsky Model with Reheating Constraints

While the origin of the Universe is largely unknown, a lot of understanding has been gained through the hypothesis known as the "inflationary paradigm" [1, 2, 3, 4, 5]. According to this hypothesis, our Universe experienced a rapid expansion at its inception, followed by a reheating period that led to the hot Big Bang (reviews on reheating are found e.g., un Refs. [6, 7, 8]). Inflation is typically described by a scalar field referred to as the "inflaton" field. After approximately 380,000 years, the early universe became transparent to photons which then produced what we observe today as the Cosmic Microwave Background radiation (CMB). This probe of the early universe provides valuable information regarding the physics of the inflationary era [9]. In particular, the inflationary paradigm predicts the presence of primordial gravitational waves; tensor fluctuations in spacetime [10]. Numerous experiments have been conducted to detect such waves, and have managed to impose upper limits on their amplitude [11]. Understanding the dynamics of the inflaton field and its realizations is crucial for interpretation the observed CMB data.

Among the variety of inflationary models lies the scalar field potential belonging to the Starobinsky type. The Starobinsky model [12] emerged as one of the early proposals for inflation. It is characterized by an action that includes the Einstein-Hilbert term, together with a scalar curvature term of higher order. The action takes the form

In this equation, $R$ denotes the Ricci scalar, $M_{Pl}=2.44\times 10^{18}\,\mathrm{GeV}$ is the reduced Planck mass, and $M$ corresponds to the Starobinsky free parameter, expressed in terms of mass units. Similar to any $f(R)$ theory, the Starobinsky action can be transformed into a form that includes the Einstein-Hilbert term and an action for a scalar field. This is achieved by considering a conformal transformation of the metric

Upon this transformation, we obtain the following action for the Starobinsky model:

where $V(\phi)$ represents the potential of the scalar field, given by

Here $V_{0}$ coreresponds to the product $\frac{3}{4}M_{Pl}^{2}M^{2}$.

In this paper we use Bayesian analysis to explore a generalization of this model. Specifically the generalized Starobinsky model which adds a term of the from $R^{2p}$ to the Einstein-Hilbert action [13]. The study includes an examination of a range of values for the exponent, as part of the parameter space. Our analysis presents a comprehensive exploration of the parameter space, with a particular emphasis on constraining the potential amplitudes and curvature through numerical evolution of the Boltzmann equation. The statistical analysis yields posterior distributions for the model parameters, together with cosmological parameters such as the number of $e$–folds, spectral index, and the tensor-to-scalar ratio.

The article is organized as follows: Section 2 provides an overview of the $R^{2p}$ generalized Starobinsky model. Section 3 presents the results of the statistical analysis using Bayes’ theorem, including the posterior distributions and the Bayes’ factor for model selection. We provide a discussion with final remarks in Section 4.

The generalized model in question stems from a straightforward extension of the Starobinsky action

which recovers the original Starobinsky model when $p=1$. In the Einstein Frame, the model gives rise to a generalised potential as follows

where $V_{0}=\frac{2p-1}{2}\left(\frac{1}{2p}\right)^{\frac{2p}{2p-1}}\left(6M^{2}\right)^{\frac{1}{2p-1}}M_{Pl}^{2}$ (see Fig. 1).

This generalization of the Starobinsky model has been studied in various contexts. For example, Motohashi et al. [13] explored the consistency relations of the generalized model at first order, while Renzi et al. [14] investigated the stability of predictions for $r$, taking into account experimental uncertainties on $n_{s}$ and assuming the validity of $\Lambda$CDM. Here we follow a similar approach to that of Ref. [15].

Cosmological observables in the context of inflationary models can be described through expressions derived from the slow-roll (SR) approximation, as demonstrated in previous studies (see, for example, [2], [16]). These are given by

In these equations, $n_{t}$ represents the tensor spectral index, $r$ is the tensor-to-scalar ratio, $n_{s}$ corresponds to the scalar spectral index and $\alpha$ its running. The amplitude of density perturbations at a specific wavenumber $k$ is denoted by $A_{s}(k)$. These quantities are evaluated at the time when the wavenumber $k$ crosses the horizon. The SR parameters introduced above are defined as

Of course, primes on $V$ indicate derivatives with respect to the inflaton field $\phi$.

In this Section we perform the corresponding statistical analysis of the generalized Starobinsky potential (2.2). Whereas the main parameters to be constrained with MCMC method are $\log V_{0}\,,p\,,N_{k}\,,$ and $n_{s}$, we will use the consistency relations (cr) from Eqs. (2.12) and (2.13) to obtain the posterior probability distribution for the tensor–to–scalar ratio $r^{cr}$ and the scalar running $\alpha^{cr}$, respectively. For that matter, we have used the Boltzmann code class [22], implementing the generalized Starobinsky potential the primordial module, and the software Monte Python [23] to perform the Bayesian inference.

Particularly we will be interested in comparing the theoretical bounds from reheating with the observational constraints from the latest Planck Satellite Collaboration data [24], including the likelihoods for temperature ($TT$) power spectra, the polarisation spectra ($TE$ and $EE$), and CMB lensing measurements. Besides, we combine such Planck data with Baryon Acoustic Oscillation (BAO) observations¹¹1For reference, qe will refer to the combination of all these sets of experiments as Planck+BAO.: the final galaxy clustering data set of the Baryon Oscillation Spectroscopic Survey [25], the BAO signal from large-scale correlation function of the 6dF Galaxy Survey (6dFGS) [26], and a sample of galaxies with low redshift ($z<0.2$) from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) [27].

Table 2 shows the priors on the parameters. The MCMC was performed for those parameters dubbed with role “Cosmological” (see last column), whereas the scalar spectral index $n_{s}$ is obtained as a “Derived” parameter from the previous ones. On the other hand, the tensor–to–scalar ratio $r^{cr}$ and the scalar running $\alpha^{cr}$ are obtained from the consistency relations (2.12) and (2.13) respectively.

In order to compare with results from previous analysis in the literature, we use the same priors on $p$ and $N_{k}$ as those used in [14]. For the rest of the standard cosmological parameters ($\Omega_{b}h^{2}\,,\Omega_{cdm}h^{2}\,,100\theta_{s}\,,\tau_{reio}$) we have considered fixed values from Planck 2018²²2From [28] we take $\Omega_{b}h^{2}=0.0224\,,\Omega_{cdm}h^{2}=0.120\,,100\theta_{s}=1.0411\,,\tau_{reio}=0.054\,.$, since such parameters do not change considerably from the observed values even when the priors are not specified. That is, the dynamics of the inflaton with the generalized Starobinsky potential does not alter the observed values for the matter content (baryonic $\Omega_{b}h^{2}$ and cold dark matter $\Omega_{cdm}h^{2}$), angular diameter distance $100\theta_{s}$, or the reionization optical depth $\tau_{reio}$.

The statistical constraints on the parameters are shown in Table 3. The only parameter with no restrictions within the broad prior we have considered is (the logarithm of) the amplitude of the inflationary potential (first row). For the rest of the parameters we found their corresponding mean values and confidence levels (CL). It is important to notice that, while the consistency relation (2.13) does not show explicit dependence on $p$, by expressing $r$ as function of $p$ and $n_{s}$ in equation (2.11), and substituting in equation (2.13), we obtain the function $\alpha^{cr}=\alpha^{cr}(p,n_{s})$.

In Figure 2 the 1D and 2D posteriors are shown. Besides, the theoretical bounds from reheating presented in Table 1 are shown in red stripes for each parameter. It can be seen that observational constraints for $\log V_{0}$ are broader than those imposed by reheating. On the other hand, for the potential parameter $p$, and the spectral index $n_{s}$ it can be observed that the theoretical bounds from reheating are thoroughly consistent with the constraints imposed by Planck+BAO. In the case of the tensor–to–scalar ratio $r^{cr}$, the lower theoretical bound of $0.0024$ is within the $1\sigma$ region. However combined data Planck+BAO indicate that lower values than that are allowed.

It can be seen that the $e$-folds number $N_{k}$ is constrained within a region that excludes part of the theoretical bounds from reheating for such parameter. As already seen in Ref. [29], $N_{k}$ depends monotonically on $\omega_{\rm re}$ such that lower values of $N_{k}$ are preferred for small $\omega_{\rm re}$. Similarly occurs for the running of the spectral index $\alpha_{s}^{cr}$, where the bounds from reheating are in a low probability region of the posterior.

In order to extend our analysis for the parameters obtained from the consistency relations, we have performed MCMC for $r$ and $\alpha_{s}$. The consistency relation (2.13) as well as the MCMC method give a negative mean value for the running. Specifically, we have obtained that the mean value for $\alpha_{s}$ from the MCMC is given by $\alpha_{s}=-5.26\times 10^{-4}\,.$ On the other hand, the tensor–to–scalar ratio as a cosmological parameter gets constrained with the value $r=2.42\times 10^{-3}\,,$ which is in agreement with the value obtained from the consistency relation (2.12) (see Table 3). For the sake of comparison, we show in Figure 3 the 1D posteriors for $r$ and $\alpha_{s}$ obtained from the MCMC, as well as $r^{cr}$ and $\alpha_{s}^{cr}$ obtained from the consistency relations (2.12) and (2.13).

In both cases we find a statistical agreement: the bounds imposed by reheating are contained within the posterior distributions obtained from any of the two methods (MCMC and consistency relations). In the case of the tensor–to–scalar ratio, the result from the MCMC $r$ as well as from the consistency relation (2.12) $r^{cr}$ are statistically consistent, and both of them contain at $1\sigma$ the lower value from the theoretical bound shown in Table 1. For the running $\alpha_{s}$ we can see that the bounds from reheating are contained within the posterior, although such bounds lie at lower probability of the posterior distribution for such parameter.

The dynamics of the early universe is ruled by an inflationary era followed by a reheating process prior to the radiation domination. Whereas the physics of inflation has been constrained by observations such as CMB anisotropies, models of reheating are not always considered as the transitional phase at the end of inflation. In this work we have considered theoretical bounds imposed by the physics of reheating in order to analyze a generalized Starobinsky potential for primordial inflation.

We foud that the information inferred from reheating leads to more restrictive values of the cosmological parameters than those imposed by Planck+BAO data. It is then important to highlight that even when analytical and semi–analytical treatment of both inflation and reheating provide insights on the most likely values of the cosmological parameters, the implementation of statistical methods using up-to-date data is mandatory if one intends to address the viability of inflationary potentials including physics from reheating. Is in this way in which we were able not only to infer the most likely value of the parameter $p$ generalizing the Starobinsky potential, but to quantify the level of agreement of the theoretical predictions from reheating with such observational constraints as well.

Particularly, in light of data from CMB and BAO we found tighter constraints for $p$ and $N_{k}$ than those presented in Ref. [14]³³3It is worth mentioning that in [14] the likelihood function for CMB data they used was Planck 2015, whereas we have used the latest data release from Planck Satellite Collaboration 2018.. The value of the parameter generalizing the standard Starobinsky potential $p$ is statistically consistent with $p=1$, and then, the original form of this inflationary potential is preferred by the combination of data Planck+BAO. In fact, we do not perform any model comparison technique due to the small departure from the Starobinsky potential. In the case of the number of $e$-folds during inflation $N_{k}$, the generalized Starobinsky potential is consistent with an approximate value of $N_{k}\sim 60$. This value lies near the maximum limit imposed by Planck observations for the standard $R^{2}$ inflationary model, where the range $49<N_{k}<59$ at $95\%$ CL is reported [9].

Concerning inflationary observational parameters such as the scalar spectral index $n_{s}$ and the scalar amplitude of perturbations $A_{s}$ we obtain the following: our inferred value for $n_{s}$ (its mean value and standard deviation) lies within the $1\sigma$ CL reported by [28] ($n_{s}^{Planck}=0.965\pm 0.004$). On the other hand, since $A_{s}$ is proportional to the amplitude of the inflaton potential $V_{0}$ (see Eq. (2.6)), we note that the values obtained from the posterior of $\log V_{0}$ are consistent with $10^{9}A_{s}=2.1^{+0.0107}_{-0.0112}$, which is also in agreement with the value reported by Planck 2018 of $A_{s}^{Planck}=2.105\pm 0.030$ (see last column of Table 2 from [28] where the $68\%$ limits is shown considering the same Planck+BAO set of data we have used in the present work).

The use of consistency relations is useful to reduce the number of degrees of freedom in the parameter space. The model–dependent nature of these relations allows us to test the generalized Starobinsky potential in a self–contained way, and to put constraints on parameters such as $r^{cr}$ and $\alpha_{s}^{cr}$ at low computational cost. Even when the considerations from reheating have been tested indirectly via theoretical constraints from the semi–analytical approach, it was possible to quantify and compare the results coming from such epoch with CMB observations, which constitutes the most robust and precise probe of the early universe. Both inflation and reheating leave characteristic signatures in the spectrum of primordial gravitational waves that could be detected in the future [30, 31], and hence the importance in constraining theoretical scenarios in order to learn more about the physics ruling the dynamics of the earliest times of our universe.

F.X.L.C. acknowledges Beca CONACHyT. J. C. H. acknowledges support from program UNAM-PAPIIT, grant IG102123 “Laboratorio de Modelos y Datos (LAMOD) para proyectos de Investigación Científica: Censos Astrofísicos". as well as the sponsorship from CONAHCyT Network Project No. 304001 “Estudio de campos escalares con aplicaciones en cosmología y astrofísica”, and through grant CB-2016-282569.