Constraining f(R) gravity in solar system, cosmology and binary pulsar systems

The $f(R)$ gravity can be cast into the form of a scalar-tensor theory, and scalar degree of freedom can be suppressed in high-density regions by the chameleon mechanism. In this article, for the general $f(R)$ gravity, using a scalar-tensor representation with the chameleon mechanism, we calculate the parameterized post-Newtonian parameters $\gamma$ and $\beta$, the effective gravitational constant $G_{\rm eff}$, and the effective cosmological constant $\Lambda_{\rm eff}$. In addition, for the general $f(R)$ gravity, we also calculate the rate of orbital period decay of the binary system due to gravitational radiation. Then we apply these results to specific $f(R)$ models (Hu-Sawicki model, Tsujikawa model and Starobinsky model) and derive the constraints on the model parameters by combining the observations in solar system, cosmological scales and the binary systems.


Introduction
Since the discovery of cosmic acceleration in 1998 [1,2], considerable efforts have been devoted in cosmology to understand the physical mechanism responsible for it. The ΛCDM model interprets the acceleration of the universe as a consequence of the cosmological constant. Although this model matches cosmological observations well [3], the cosmological constant suffers from some theoretical problems. If the cosmological constant originates from the vacuum energy in quantum field theory, extreme fine-tuning is required to explain its smallness [4]. It is also difficult to explain its closeness to the present matter density of the universe [4]. This motivates the search for alternative explanations for the cosmic acceleration.
Two types of approaches have been considered. One can either introduce a new kind of matter whose role is to trigger acceleration, or modify the behavior of gravity on cosmological scales [5,6]. In the first approach, dark energy is introduced as a new energy form, which has positive energy density but negative pressure. In the second approach, various attempts to modify gravity have been presented. For recent reviews on modified gravity, see [7,8,9,10].
Lovelock's theorem states that General Relativity (GR) represents the most general theory describing a single metric that in four dimensions has field equations with at most second-order derivatives [11]. As a result of this theorem, one way to modify Einstein's field equations is to permit the field equations to be higher than second order. In this paper, we will consider the so-called f (R) gravity which has fourth order field equations. The Ricci scalar R in the gravity action is replaced by a general function of Ricci scalar f (R). For reviews on f (R) gravity, see [12,13]. The f (R) gravity does not introduce any new type of matter and can lead to the late time acceleration of the universe [14,15]. When cast into the scalar-tensor theory, the f (R) gravity implies a strong coupling between the scalar field and matter. This would violate all experimental constraints on deviations from Newton's gravitation [16]. Certain constraints have to be imposed on the function f (R) for the model to be linearly stable [17,18] and pass local gravitational tests [19]. The first attempt f (R) = R − µ 2(n+1) /R n proposed by Carroll et al. in [20] failed these constraints right away. However, since then, models that evade them have been found [21,22]. Fortunately, the chameleon mechanism can alleviate these constraints. Imposing the chameleon mechanism, the scalar field can develop an environment dependent mass [23,24,25]. When the ambient matter density is large enough, its mass becomes large, and the corresponding fifth force range is short. Thus the scalar field can be hidden in the high density environment and the fifth force cannot be detected [16].
The parametrized post-Newtonian (PPN) formalism is useful to study different theories of gravity [26,27,28,29]. In the PPN formalism, the PN (weak field and slow motion) limit of different theories are characterized by a set of PPN parameters and the most important two parameters are γ and β. These two parameters can be directly measured by the solar system experiments. The GR prediction (γ = 1 and β = 1) is consistent with the observations [30], which provide constraints on various modified gravity models [31,32]. Meanwhile, the binary pulsar systems can emit gravitational waves and provide a good test for gravitational theories [26,27,33,34,35]. Since these systems lose energy due to gravitational radiation, the orbital period of these systems will decay 1 . Several authors have considered this effect in f (R) gravity [41,42,43] for some specific models. However, in these works, the authors have ignored the chameleon mechanism . Although some authors have applied the chameleon mechanism to f (R) gravity when they study the PN limit, they only calculate the PPN parameter γ [44,45].
In this paper, we give a comprehensive investigation on various constraints on the general f (R) gravity with chameleon mechanism. Following the method developed in our previous work [46], we first calculate the PPN parameters γ and β, the effective cosmological constant, and the effective gravitational constant in the general f (R) gravity. Considering the current observations in solar system and cosmological scales, we derive the combined constraint for the general f (R) gravity. Binary pulsar system is a good testing ground for alternative theories of gravity. In the previous work [47], we have derived the orbital period derivative for quasicircular binary systems in scalar-tensor gravity with chameleon mechanism. Here, applying the similar analysis to f (R) gravity, we obtain the orbital period derivative for quasicircular neutron star-white dwarf (NS-WD) binary systems. Using the observational data of PSR J0348 +0432 [48] and PSR J1738 +0333 [49], we also obtain the binary pulsar constraints on f (R) gravity. We find that the chameleon mechanism cannot apply to Palatini f (R) gravity. Thus, in the paper, we mainly focus on metric f (R) gravity. Applying the general results to the specific f (R) models, including Starobinsky model, Hu-Sawicki model and Tsujikawa model, we obtain the constraints on the model parameters.
The paper is organized as follows: In Sec. 2, we review f (R) gravity and chameleon mechanism. In Sec. 3, we study various observational constraints on f (R) gravity, and obtain the parameter constraints on the specific models. We conclude in Sec. 4.

f (R) gravity and Chameleon mechanism
The f (R) gravity comes about by a straightforward generalization of the Ricci scalar R to become a general function f (R) in the action for gravity. When varying the action, there exist two formalisms: the metric formalism and the Palatini formalism. In the Palatini formalism, the connection is not taken to be the Levi-Civita connection of the metric a priori and one varies the action assuming that the metric and the connection are independent variables. Although these two formalisms lead to the same field equations in GR [50], this is no longer true for f (R) gravity. We will investigate these two formalisms respectively.

Metric f (R) gravity
The total action for metric f (R) gravity takes the form [12] S = 1 16πG where Ψ m denotes all the matter fields. Variation with respect to the metric g µν gives the field equations [12] f where a prime denotes differentiation with respect to R and = ∇ µ ∇ µ . Since the field equations contain the second derivative of R and R includes second derivatives of the metric, the field equations are fourth order partial differential equations in the metric.
Handling fourth order equations can be troublesome, but f (R) gravity can be recast as a scalar-tensor theory via a conformal transformation and the corresponding field equations become second order. Conformal transformation of the metric can also show the scalar degree of freedom explicitly. Introducing a new field χ, we obtain a dynamical equivalent action [12] Varying this action with respect to χ, we have f (χ)(R − χ) = 0. If f (χ) = 0, we have R = χ. Substituting this into Eq. (3) leads to Eq. (1). Redefining the field by θ = f (χ) and setting U (θ) = θχ(θ) − f (χ(θ)), we have The action (4) is in the Jordan frame, which should be transformed into the Einstein frame to utilize the results of the prior studies [46,47], although the chameleon mechanism also works in the Jordan frame [51]. Defining the metric in Einstein frame asg µν = θg µν , we get the Einstein frame action as follows [12], where (∂θ) 2 =g µν ∂ µ θ∂ ν θ andR is the Ricci scalar ofg µν . To change the kinetic term into the standard form, we introduce another scalar field φ that satisfies the The scalar field φ can be directly related to the Jordan frame Ricci scalar by Therefore, the action in the Einstein frame has the form [12], where the bare potential is The conformal coupling function is [12] A(φ) = 1 with the conformal coupling parameter ξ = 1/ √ 6. Variation of S E with respect tog µν and φ gives the field equations whereT µν ≡ (−2/ √ −g)δS m /δg µν is the energy-momentum tensor of matter in the Einstein frame, and˜ ≡g µν∇ µ∇ν . The covariant derivatives∇ µ obeỹ ∇ µgαβ = 0. The scalar field equation can be rewritten as follows: with the effective potential Here the matter is assumed to be nonrelativistic, and ρ ≡ −T /A is the conserved energy density in the Einstein frame, which is independent of φ [24].

Chameleon mechanism
An important consequence of the conformal coupling function A(φ) is that matter will generally feel a fifth force mediated by the scalar field. Since the conformal coupling parameter ξ is of order unity, the fifth force will have a significant impact on the motion of particles [16]. In order to evade the fifth force constraints, the mass of the field should be sufficiently large in high density environment [52]. Since scalar field needs to have cosmological effects to accelerate the expansion of the Universe, on cosmological scales, the magnitude of the scalar mass can be Hubble scale to cause the acceleration of the universe. Thus a mechanism is needed to screen the scalar field in local environment while let the scalar field accelerate the Universe on large scale [16,45].
The behavior of the scalar field is governed by the effective potential V eff (φ). An essential element of the model is the fact that V eff (φ) depends explicitly on the matter density, as seen in Eq. (14). The shape of the effective potential is determined by the function f (R). For a suitably chosen function f (R), the effective potential can have a minimum. We denote by φ min the value at the minimum, that is [46], dV eff dφ φmin = 0.
Whilst the mass of small fluctuations around the minimum is [16], It can be observed that the scalar field has a density dependent mass. When the density of the environment is large enough, the mass becomes large, and the corresponding fifth force range is so small that it cannot be detected by gravitational experiments [16]. Laboratory constraints can be greatly alleviated if the mass develops a strong dependence on the ambient density of matter.
Theories in which such a dependence is realized are called to have a chameleon mechanism. Therefore, if the following three conditions can be satisfied in some regions of φ space, the f (R) model can have a chameleon mechanism [6]: The mass can increase with density. Using Eq. (6), these conditions can be translated into [16]

Palatini f (R) gravity
Previous discussions have focused on the metric formalism. We now discuss the Palatini formalism. The action in the Palatini formalism is formally the same as in the metric formalism. However, the Ricci tensor is constructed from the independent connection and is not related to the metric tensor. The Palatini action takes the form [12] Here R ≡ g µν R µν and the Ricci tensor R µν is determined by the independent connection Γ µ αβ . Variations with respect to the metric and the connection can yield the following formulae respectively [12], and Transforming the action into the Einstein frame, we obtain [12] which follows the scalar field equation, Note that, the scalar field θ is algebraically related toT , i.e., θ = θ(T ), which is non-dynamical and cannot propagate in spacetime. Therefore, we cannot define a mass of the scalar field θ, as discussion above. As a result of the non-dynamical nature of the scalar field, the chameleon mechanism does not apply to Palatini f (R) gravity. There exists another significant difference between Palatini f (R) gravity and the chameleon theory: Since the fifth force is produced by the gradient of the scalar field, and in chameleon theories a compact object in a homogeneous background can generate a scalar field with Yukawa profile. A test particle in the homogeneous background can feel the fifth force. While in Palatini f (R) gravity, the scalar field does not have gradient in a homogeneous background and does not mediate a fifth force. In addition, there are other serious shortcomings of Palatini f (R) gravity [12]. So, in the rest of this paper, we will only focus on metric f (R) gravity.

Stability issues
More recent attention has focused on the stability issues about metric f (R) gravity. These include Ostrogradski instability [53], Frolov instability [54], Dolgov-Kawasaki instability [55] and instability of de Sitter space [18]. A scrutiny of these issues is needed to make sure that f (R) gravity is viable. The first two stability issues can be bypassed in the specific models discussed below [53,56].
Dolgov and Kawasaki [55] found that the Ricci scalar is instable in the f (R) model proposed by Carroll et al. [20]. Their analysis is generalized to a general function by Faraoni [17]. The origin of this issue is that the mass squared of the scalar degree of freedom is negative. Since the mass squared has the same sign as f (R), the stability condition can be written as [12] f where R 0 is the Ricci scalar today. This condition is satisfied for all the models studied in the following section.
In order to investigate the stability of de Sitter space, we consider a spatial flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe. The vacuum field equations take the form [12] where an overdot denotes differentiation with respect to t. The stationary points of the dynamical system (26) are de Sitter space with Hubble constant H. The condition for the existence of de Sitter space is [12] Rf − 2f = 0.
The stability condition of de Sitter space with respect to inhomogeneous linear perturbations reads [18] If the solution to Eq. (27) meets the stability condition (28), the Universe will enter into a stable de Sitter phase in the future [14]. We now impose the stability condition of de Sitter space on the specific f (R) models. We will investigate the following well studied models The models (A), (B) and (C) are proposed by Hu and Sawicki [45], Tsujikawa [19] and Starobinsky [15], respectively. In the model (A), the mass scale is chosen to be [45] whereρ 0 is the average matter density in the universe today. In the models (B) and (C), R c roughly corresponds to the order of observed cosmological constant for µ = O(1). During the whole expansion history of the Universe, the Ricci scalar is in the high curvature region, i.e., R m 2 or R R c [45]. Thus, the model (A) can be approximated by and the model (C) can be approximated by It can be observed that the free parameters of the model (A) are in one-to-one correspondence with that of the model (C) through the relations m 2 c 1 /c 2 → µR c , m 2(n+1) c 1 /c 2 2 → µR 2k+1 c and n → 2k. So we only study the models (A) and (B) in the following.
The model (A) can be expressed as another useful form [19] f where α = c 1 c 1/n−1 2 and R c = m 2 c −1/n 2 . The following relation holds at the de Sitter point: [13] where x ≡ R/R c . The stability condition (28) implies the relation [13], Thus for each specific n, the above inequality gives a bound on x and this bound can be transformed into a bound on α through Eq. (36). For instance, when n = 2, one has x ≥ √ 3 and α ≥ 8 √ 3/9. In the following section, we will come back to discuss this inequality.

Constraints on f (R) gravity
In this section, we consider the observational constraints on metric f (R) gravity in cosmological scale, solar system and binary pulsar systems, respectively.

Cosmological constraints
In order to satisfy the tests on cosmological scales, the f (R) models should mimic the ΛCDM model at the late time and provide an effective cosmological constant. Similar to the previous work [46], in this paper we do not consider the cosmological perturbations of f (R) gravity [13]. We leave this issue for the general f (R) gravity as a future work. The bare potential V (φ) in action (7) can provide the effective cosmological constant to accelerate the universe expansion, which is given by [46] where R ∞ is the background value of Ricci scalar. In order to mimic the ΛCDM model, we need that the value of Λ eff is equal to the observed cosmological constant Λ, which accelerates the cosmic expansion. Now, we can apply the cosmological constraint (38) to specific f (R) models. For model (A), we substitute Eq. (33) into Eq. (38), and obtain which, in turn, implies that Note that, we adopted the density parameters Ω m = 0.308 and Ω Λ = 0.692 [3]. This expression of c 1 /c 2 is consistent with Eq. (26) in [45]. Now it can be seen from Eq. (33) that there are two remaining parameters n and c 1 /c 2 2 in this model.
Using the relation α = c 1 c Using Eq. (33), we have For a spatial flat FLRW universe, the scalar curvature at the present epoch is R 0 = m 2 (12/Ω m − 9) [45]. Consequently, for different n, we can obtain different upper bounds on |f (R 0 ) − 1|. The results are presented in Fig. 1 with dotted line.
Similarly, for the model (B), the cosmological constraint is and the stability condition (28) implies that [13] µ > 0.905.
3.2. Solar system constraints In the solar system, the gravitational field is weak and the velocity of planets is slow compared with the speed of light. Thus we can apply the PPN formalism to solar system tests. In the PN limit, the spacetime metric predicted by different metric theory of gravity has the same structure and can be characterized by ten PPN parameters [26]. Among them, the most important parameters are γ and β.
Here, we derive the PPN parameters γ and β and the effective gravitational constant G eff in the general metric f (R) gravity with chameleon mechanism. For a scalar-tensor theory with action (7) , the solution to the scalar field equation (11) is given by [47] where the screened parameter is defined as, The parameters M E and Φ E ≡ GM E /r are the mass and the Newtonian potential at the surface of the source object in the Einstein frame, respectively. The quantity φ 0 is the field in side the source object and φ ∞ is the field in the background environment. m ∞ is the effective mass of scalar field at φ = φ ∞ . In order to solve the metric field equations, we make use of the PPN formalism introduced in [26,27]. In this formalism, the gravitational field of the source is weak GM/r 1, and the typical velocity v of the source is small, i.e. v 2 ∼ GM/r 1. Thus, we can use the perturbative expansion method to solve the field equations, and all dynamical quantities can be expanded to O(n) ∝ v 2n . The metric field g µν can be expanded around the Minkowski background as follows: We solve the field equations (10) and (11) using the PPN method [27], and transform the metric to the Jordan frame. Making use of the definitions of γ and β as follows [46], (1) where M J and χ are the mass and radial coordinate in the Jordan frame, respectively. We obtain the PPN parameters where A VEV ,A 1 and A 2 are the expansion coefficients of A(φ), i.e., Note that, here we have taken the limit m ∞ r 1, since in the solar system, the distance r is always much less than the Compton wavelength m −1 ∞ [46]. Applying to the general metric f (R) gravity, using Eqs. (6) and (9) we obtain the expansion coefficients .
Following the discussion of [45], R ∞ = 8πGρ g and ρ g = 10 −24 g cm −3 is the average galactic density in the solar vicinity. In the solar system, the source object of the scalar field is the Sun and the background is the Milky Way. Since the density of the Sun is much higher compared with the galactic background, we have φ ∞ φ 0 . Then the screened parameter can be approximated by Substituting the above parameters into Eq. (50), we obtain We find that the expression of parameter γ is consistent with Eq. (64) in [45]. The parameter β is unity no matter whatever the functional form f (R) is. As can be seen from Eq. (50), when the conformal coupling function A(φ) has the exponential form, the two terms in the bracket of the expression of β cancel each other out. And Eq. (9) shows that the conformal coupling function A(φ) always has the exponential form in metric f (R) gravity. This suggests that the experimental tests of parameter β cannot distinguish between GR and metric f (R) gravity. The relation between γ and G eff is Using the Cassini constraint |γ − 1| < 2.3 × 10 −5 [30] and the Newtonian potential at the surface of the Sun Φ E = 2.12 × 10 −6 , we obtain the constraint on general f (R) gravity as follows, Note that, this is a general constraint for any metric f (R) gravity with chameleon mechanism, which is independent of the form of f (R). We can apply this solar system constraint to the models (A) and (B). In the model (A), using Eq. (33), we have [45] 1 Here, R ∞ = 8πGρ g and ρ g = 10 −24 g cm −3 is the average galactic density in the solar vicinity. We adopted the physical matter density Ω m h 2 = 0.1415 [3].
Applying inequality (57) to the above equation, we have The equivalence principle places a bound on the parameter of model (A) [57] n > 1.8. As shown in Fig. 1, in the region 1.8 < n < 3, the solar system constraint (solid line) is fairly weak when compared with the stability condition (dotted line) and is sensitive to the value of n.
Similarly, in the model (B) we have where Eq. (44) was used to eliminate R c in terms of µ. The solar system constraint (57) yields µ > 9.5 × 10 −5 .
Compared with the stability condition (45), this shows that the solar system constraint on f (R) gravity is weaker for this model . Assuming that the cosmological constraint (38) and solar system constraint (57) are both satisfied, we have checked that in the model (A) and (B), the conditions for chameleon mechanism (17)- (19) can all be satisfied.

Binary pulsar constraints
It is well known that the compact binary systems can lose the orbital energy due to gravitational radiation, and the orbital period will decay. In different theories of gravity, the decay rates are different [27,26], which provides another independent opportunities to test the metric f (R) gravity. In a binary system, when the difference between the screened parameters of the two compact stars is significant, the dipole radiation dominates the orbital decay rate. Since the screened parameter is inversely proportional to surface gravitational potential, the neutron star-white dwarf (NS-WD) systems are the best testbeds to constrain the parameters of f (R) gravity. In the previous work [47], we have studied this effect in the most general scalar-tensor gravity with screening mechanism. For a quasicircular (e 1) NS-WD binary system, the orbital period derivative is given by [47]Ṗ Here, P denotes the orbital period, m = m NS + m WD is the total mass, µ = m NS m WD /m is the reduced mass, WD = φ ∞ /M Pl Φ WD and NS = φ ∞ /M Pl Φ NS are the screened parameter of the white dwarf and the neutron star respectively andṖ represents the GR prediction of the orbital period derivative. The second term in Eq. (63) corresponds to the scalar dipole radiation correction. We apply this result to the general metric f (R) gravity with chameleon mechanism. Using Eq. (6), the orbital period derivative translates intȯ It can be seen that in the special case f (R) = R−2Λ, the above result is reduced toṖ =Ṗ GR . Because Φ NS /Φ WD ∼ 10 4 , the orbital period derivative can be approximated byṖ Since all the pulsar observation agrees well with the GR prediction within the errors [35,48,49], the observation value of the period derivative can be expressed asṖ where δ is the fractional deviation of the observedṖ obs from the GR prediction, σ is the observational uncertainty. Thus the background field value f (R ∞ ) cannot deviate from unity too much, that is, (68) at 95% confidence level.
Up to now, more than 2500 pulsars have been observed [33]. However, most of them are isolated and their mass cannot be determined. Table 2 in [33] lists fifteen NS-WD systems with low-eccentricity orbits which have accurate measurement of mass. Among these fifteen NS-WD systems only PSR J0348 +0432 and PSR J1738 +0333 have accurate observation value of the radius of the white dwarf companion. Thus we use these two NS-WD systems to constrain f (R) gravity and list here the relevant parameters in Table 1.
In the PSR J0348 +0432 case (see Table 1), δ = 0.05 and σ = 0.18. Substituting the parameters into inequality (68), we obtain the upper bound at 95% confidence level. Similarly, using the observation data of PSR J1738 +0333, we obtain at 95% confidence level. Compared with the solar system constraint (57), the pulsar constraint is three orders of magnitude weaker. Applying the pulsar constraint to the model (A), we obtain The above result is also shown in Fig. 1 with dashed line. Similarly, applying the pulsar constraint to the model (B), we obtain µ > 5.4 × 10 −5 .
Consistently, we find both of them are relatively weaker than the corresponding constraints of solar system.

Conclusions
The f (R) gravity has been extensively studied to explain the accelerating expansion of the universe. In this paper, we have studied the general f (R) gravity through the scalar-tensor representation. In this theory, the chameleon mechanism is crucial for f (R) gravity to escape the fifth force constraints. However, due to the non-dynamical nature of the scalar field in Palatini f (R) gravity, this mechanism does not apply to the theory. Therefore, we focused on the metric f (R) gravity with chameleon mechanism.
We calculated the PPN parameters γ and β for the general f (R) gravity, and found that β = 1 in the limit m ∞ r 1. As a result, the observed value of β cannot constrain the parameters of f (R) models. Applying the Cassini spacecraft measurement of γ, we obtained the constraint |f (R ∞ ) − 1| < 4.9 × 10 −11 on the metric f (R) gravity, which is consistent with the previous works. To pass the cosmological test, the metric f (R) gravity should provide the effective cosmological constant. We also calculated the effective cosmological constant in f (R) gravity. In general, the cosmological constraint can reduce one free model parameter in a given specific f (R) model.
In addition, we calculated the orbital period derivativeṖ of binary pulsar systems in the metric f (R) gravity. Since GR has survived the binary pulsar test, theṖ in the metric f (R) gravity cannot deviate from that in GR too much. We found that the pulsar constraint from the observations of PSR J0348 +0432 and PSR J1738 +0333 is |f (R ∞ ) − 1| < 3.6 × 10 −8 . This is relatively weaker than the current constraints derived from the solar system observations. We also studied the stability condition of de Sitter space. Compared with the observational constraints (binary pulsar and solar system), this theoretical constraint is more stringent in Hu-Sawicki model and Tsujikawa model. With the chameleon mechanism, the metric f (R) gravity with suitable parameters can pass the cosmological test, the solar system test and the binary pulsar test at the same time.