Matter Power Spectra in Viable $f(R)$ Gravity Models with Massive Neutrinos

We investigate the matter power spectra in the power law and exponential types of viable $f(R)$ theories along with massive neutrinos. The enhancement of the matter power spectrum is found to be a generic feature in these models. In particular, we show that in the former type, such as the Starobinsky model, the spectrum is magnified much larger than the latter one, such as the exponential model. A greater scale of the total neutrino mass, $\Sigma m_{\nu}$, is allowed in the viable $f(R)$ models than that in the $\Lambda$CDM one. We obtain the constraints on the neutrino masses by using the CosmoMC package with the modified MGCAMB. Explicitly, we get $\Sigma m_{\nu}<0.451 \ (0.214)\ \mathrm{eV}$ at 95% C.L. in the Starobinsky (exponential) model, while the corresponding one for the $\Lambda$CDM model is $\Sigma m_{\nu}<0.200\ \mathrm{eV}$. Furthermore, by treating the effective number of neutrino species $N_{\mathrm{eff}}$ as a free parameter along with $\Sigma m_{\nu}$, we find that $N_{\mathrm{eff}} = 3.78^{+0.64}_{-0.84} (3.47^{+0.74}_{-0.60})$ and $\Sigma m_{\nu} = 0.533^{+0.254}_{-0.411}$ ($<0.386) \ \mathrm{eV}$ at 95% C.L. in the Starobinsky (exponential) model.


I. INTRODUCTION
The observations of the Type-Ia supernovae [1,2] in the last decade of the 20th century indicated that our universe is undergoing an accelerating expansion. Since then, the phenomenon has been further verified by several succeeding experiments [3][4][5]. To explain this interesting phenomenon, people have tried various methods. One of which is to introduce a homogeneous and isotropic energy density with negative pressure into the theory of General Relativity (GR), so-called "Dark Energy." [6] The other way is to modify Einstein's gravity theory by extending the Ricci scalar R in the Einstein-Hilbert action to an arbitrary function f (R) [7,8]. Several viable f (R) models have been proposed to satisfy the constraints from theoretical considerations as well as cosmological observations [8].
On the other hand, the oscillations between the three flavors of neutrinos in the standard model of particle physics have been detected [9,10], suggesting that either two or three of the active neutrinos have tiny masses. Clearly, from the cosmological point of view, it is necessary to consider the effect of massive neutrinos on the evolution of our universe [11][12][13][14][15][16]. For example, massive neutrinos will suppress the matter power spectrum in the small scale [11,12]. In other words, cosmology offers strong constraints on the mass scales of neutrinos. However, such cosmological constraints are highly model dependent. It is know that the simplest model in cosmology, the ΛCDM model, permits only a small range for the sum of the active neutrino masses. For example, the constraint from Planck [5] allows Σm ν < 0.23eV at 95% confidence level. In the viable f (R) models, their matter power spectra are normally larger than that of the ΛCDM [12,13]. This enhancement then can be used to compensate for the suppression due to massive neutrinos so that the neutrino mass constraint is relaxed to a broader window. Recently, Σm ν < 0.5 eV has been extracted in Refs. [14][15][16] for the chameleon type of f (R) gravity. In this paper, we will examine two typical viable f (R) models with their exact forms.
In general, the viable f (R) models can be categorized into power-law and exponential types. In this paper, we focus on the Starobinsky and exponential gravity models, which belong to these two types, respectively. Without loss of generality, we consider these f (R) models with one massive neutrino along with the other two being massless. Using the modified Code for Anisotropies in the Microwave Background (MGCAMB) [17,18] and the CosmoMC package [19], we study the constraints on the neutrino masses from the latest cosmological observational data, including those of the cosmic microwave background (CMB) from Planck [5] and WMAP [20], baryon acoustic oscillation (BAO) from Baryon Oscillation Spectroscopic Survey (BOSS) [21], Type-Ia supernova (SNIa) from Supernova Legacy Survey (SNLS) [22], and matter power spectrum from Sloan Digital Sky Survey (SDSS) [23] and WiggleZ Dark Energy Survey [24]. The constraint on the effective number of neutrino species, N eff , is also acquired in order to examine the non-standard properties of neutrinos. Since MGCAMB uses the parametrized framework to include f (R) gravity into CAMB, we only consider the linear perturbation and assume that the background evolution is the same as the ΛCDM model. This paper is organized as follows: In Sec. 2, we first give a brief review on the f (R) modification in the linear perturbation theory and then show the matter power spectra P (k) in the two types of the viable f (R) models. The effect of massive neutrinos is examined.
In Sec. 3, we show the results of the constraints on massive neutrinos using the CosmoMC package. Finally, we present our conclusions in Sec. 4.
The action of f (R) gravity is given by where κ 2 ≡ 8πG, g is the determinant of metric tensor, f (R) is an arbitrary function of Ricci scalar and L M denotes the matter Lagrangian density. By varying the action (1) with respect to the metric g µν , we obtain where the subscript "R" denotes the derivative of R, i.e. f R ≡ ∂f /∂R, = g µν ∇ µ ∇ ν is the d'Alembertian operator, and T µν is the energy-momentum tensor, defined by To deal with the dark energy problem, we must check whether f (R) gravity satisfies the following viable conditions: (i) f R > 0 for R > R 0 ≡ R(z = 0), which keeps the positivity of the modified Newton constant and avoids the attractive gravitational force; (ii) f RR > 0 for R > R 0 , which guarantees that the mass of the scalaron is real defined; (iii) f (R) → R −2Λ in the high redshift region (R ≫ R 0 ), which reproduces the ΛCDM behavior in the early universe; (iv) a late-time stable solution exists, which eliminates the appearance of singularity in the future; and (v) it should pass the local gravity tests, including those from the equivalence principle and solar system. Several viable f (R) models with these conditions have been proposed, such as Hu-Sawicki [25], Starobinsky [26], Tsujikawa [27], exponential [28][29][30][31], and Appleby-Battye [32,33] gravity models. These models can be grouped into power law and exponential types, denoted as Type-I and II, respectively [34].
In the following discussion, in order to investigate the behaviors of these two types of models, we concentrate on the Starobinsky and exponential models, given by respectively, where R c represents the characteristic curvature.

B. Matter Power Spectrum
We review the perturbation equations of f (R) gravity in the Newtonian gauge and study the effect of f (R) gravity on the matter power spectrum with the parametrization used in MGCAMB [17,18].
The perturbed FLRW metric in the Newtonian gauge is given by where a(η) is the scale factor, η is the conformal time, and Ψ and Φ are the scalar perturbations.
It is worth noting that the viable f (R) gravity theory is indistinguishable from the ΛCDM model at the high-redshift stage. Although the specific distinction between f (R) and ΛCDM depends on the forms of f (R), it always happened later than z = 10, at which the nonrelativistic matter dominated our universe. Therefore, we only consider the perturbation equations in the matter dominated region. The perturbed energy-momentum tensor is where v m is the velocity field. By following the similar procedure in Ref. [35] with the energy-momentum conservation, δT µ ν;µ = 0, we can derive the perturbation equations in f (R) theory: where the "dot" denotes the time derivative and k is a comoving wavenumber under the Fourier transform. Besides, it is convenient to use the gauge-invariant matter density perturbation δ m , defined by From Eqs. (9) -(14) with the sub-horizon limit (k 2 ≫ a 2 H 2 ), the relations between the metric potential and the perturbation of the matter density are given by where Eqs. (15) and (16) are the parametrizations used in MGCAMB to incorporate f (R) gravity.
In other words, MGCAMB modifies the linear perturbation to include the effect of f (R) that has the ΛCDM limit.
In MGCAMB, the f (R) gravity effect is introduced by Bertschinger and Zukin (BZ) in Ref. [36] through the form where λ 1 is defined by the Compton wavelength in units of the Hubble scale [18,36], where However, comparing Eq. (17) to (19), the BZ approach describes the f (R) effect as a function of time, indicating f RR ∝ a 6 in the high redshift regime (f R ∼ 1). This approach can only be used in the Starobinsky model with n = 1 in the matter dominated era, but not in the Starobinsky (n = 1) and exponential models since they deviate from the ΛCDM model as f RR ∝ (R c /R) 2n and e −R/Rc , respectively. Instead of the BZ approach used in Refs. [15,18], we further modify MGCAMB and take the exact f (R) form in Eq. (17) to represent the matter power spectrum in f (R) gravity. Note that in Refs. [14,16] the authors have done the analysis by using the full linear perturbation equations of f (R) gravity [37]. Through the viable condition (iii), the characteristic curvature R c is defined by the dark energy density ρ Λ and the parameter λ (β), in the Starobinsky (exponential) model, given by Within this framework, we conduct the program of the CosmoMC package with MGCAMB to perform our calculations.
Instead of investigating the numerical result from MGCAMB directly, we can figure out the effect of f (R) gravity on the matter power spectrum, P (k) ∼ δ 2 m , based on the equation, where µ 1 = f −1 R and µ 2 = (1 + 4k 2 f RR /a 2 f R ) / (1 + 3k 2 f RR /a 2 f R ). For the Starobinsky and exponential models, their first-order derivatives, f (s) are both smaller than unity, so that µ 1,2 are greater than unity. The first two viable conditions, f R > 0 and f RR > 0, guarantee that the numerator of µ 2 is larger than its denominator.
Furthermore, µ 2 is enhanced by a large wavenumber k, corresponding to the matter power spectrum in the small scale. In all, the matter power spectra of the viable f (R) gravity theories are always larger than that of the ΛCDM model as explicitly shown in Fig. 1, where the density of dark energy is fixed to be Ω DE ∼ 73% and neutrino masses are taken to be zero.
The deviation between the Starobinsky (exponential) and ΛCDM models is proportional to (R c /R) 2n (e −R/Rc ). If the dark energy density ρ Λ is fixed at the present time (z = 0), the characteristic curvature R c is inverse proportional to the model parameter, denoting that a smaller model parameter corresponds to a larger deviation for the matter power spectrum. As a result, we consider the model parameter space close to the lower bound of the viable condition (λ, β ≥ 1 and n ≥ 2) in the following calculations. In Fig. 2, we show the matter power spectra in the Starobinsky (n = 2) and exponential models, where the ΛCDM result is also given. As indicated above, a large enhancement of f (R) gravity occurs at the large wavenumber k. There is an interesting phenomenon: the magnification of P (k) in the Starobinsky model is more greater than that in the exponential one within the same allowed viable model parameter (e.g., λ = β and R star c = R exp c ). Since there is an exponential decay in the Type-II models, it converges toward the ΛCDM result much faster than that in the Type-I ones as the curvature increases. Accordingly, it is possible to exist a much larger enhancement the Type-I model than that in the Type-II one. To study the effect of massive neutrinos on the matter power spectrum, we first use the relation between the contributions of massive neutrinos to the total energy density, Ω ν , and the total mass, Σm ν , in unit of eV, given by where h is the reduced Hubble constant. The upper value at 95% C.L. for the total neutrino mass constrained by the Planck data in the ΛCDM model [5], Σm ν < 0.23 eV, leads to Ω ν ≤ 5 × 10 −3 . Although it is a rather small ratio to the total energy density, its effect on the matter power spectrum is still detectable, as mentioned in Sec. 1. As demonstrated in Fig 3, we see that the free streaming massive neutrino suppresses the growth of P (k) in the subhorizon scale [11]. This effect of massive neutrinos on the matter power spectrum is opposite to that of f (R) models [12]. As a result, the Type-I f (R) models would allow a higher scale for the total neutrino mass than Type-II.

III. CONSTRAINTS FROM COSMOLOGICAL OBSERVATIONS
We now perform the numerical simulations for the viable f (R) models by using the CosmoMC package with some massive neutrinos and the latest cosmological data. The dataset is as follows: the CMB data from Planck with both low-l (l < 50) and high-l (l ≥ 50) parts and WMAP with only the low-l one; the BAO data from BOSS DR11; the matter power spectral data from SDSS DR4 and WiggleZ Dark Energy Survey; and the SNIa data from SNLS. With this dataset, we explore the constraints on the Cold Dark Matter density, Ω c h 2 , and the sum of the active neutrino masses, Σm ν , in the two viable f (R) models. In Table I, we list the cosmological parameters in our analysis. Spectral index 0.9 < n s < 1.1 Reionization optical depth 0.01 < τ < 0.8 In Fig. 4, we illustrate the constraining contours of Σm ν and Ω c h 2 with the assumption of only one massive neutrino. Our results are summarized in Table II. From Fig. 4 and Table II, we see that Σm ν < 0.451 eV (95% C.L.) in the Starobinsky model, which is close to the upper bound given in Ref. [15], while the exponential model leaves the massive neutrino a rather small space, Σm ν < 0.214 eV (95% C.L.), which is very close to that of the ΛCDM model, Σm ν < 0.200 eV (95% C.L.). Clearly, the behavior the Type-II viable f (R) models is barely distinguishable from the ΛCDM one.
In addition, we can examine the effective number of neutrino species, N eff , to account for the neutrino-like relativistic degrees of freedom, defined by ρ radiation ≡ ρ γ + ρ ν = 1 + 7 8 4 11 where ρ γ is the energy density of photon, 7/8 comes from the Fermi-Dirac distribution since neutrinos are fermions, and 4/11 is due to the ratio of the neutrino to photon temperature.
The effect of N eff is mostly on the epoch of the matter-radiation equality and the expansion rate as well as the CMB power spectrum. Our results are shown in Fig. 5 and Table III. The best-fit value of N eff is 3.78 in the Starobinsky model, which is higher than the corresponding one of 3.47 in both ΛCDM and exponential models. This infers that the Starobinsky model allows more relativistic species than the other two models. However, at 95% C.L., the three models admit approximately the same range for N eff . On the other hand, the total neutrino mass in the two viable f (R) models as well and the ΛCDM one increases when N eff is treated as a free parameter.

IV. CONCLUSIONS
We have studied the effect of massive neutrinos in the two types of viable f (R) gravity theories, the Starobinsky and exponential models, by using the CosmoMC package with the modified MGCAMB. We have considered the linear perturbations in these models by assuming that the background evolutions are the same as the ΛCDM. The enhancement of the matter power spectrum has been found to be a generic feature in the viable f (R) gravity theory. However, the biggest magnifications of the matter power spectra in the Type-I f (R) models are more significant than those in the Type-II ones. With an increasing curvature, the results in the Type-I models approach to the ΛCDM as an inverse power law, while those in the Type-II models as an exponential decay. Clearly, the Type-I viable f (R) models allow larger mass scales for massive neutrinos than those in Type-II since massive neutrinos suppress the matter power spectra. As a result, the modified f (R) gravity theory would be used to compensate the suppression from the effect of massive neutrinos. If the three neutrinos are in a large mass scale, the Type-I viable f (R) theory, such as the Starobinsky model, is favored.
Our investigation has shown that the allowed neutrino mass scale is further released when N eff is considered as a free parameter. Moreover, the best-fit value of N eff has been found to be 3.78 in the Starobinsky model, which is greater than 3.47 in the ΛCDM and exponential models. Clearly, the Starobinsky model leaves more room for a dark sector or a sterile neutrino in the universe.