Consistency of $f(R)$ gravity models around solar polytropes

It is stated that a class of $f(R)$ gravity models seem to obtain $\Lambda$CDM behaviour for high redshifts and general relativistic behaviour locally at high curvatures. In the present paper, we numerically study polytropic configurations that resemble stars like young sun with Hu and Sawicki $f(R)$ gravity field equations and compare the spacetime at the boundary to the general relativistic counterpart. These polytropes are stationary spherically symmetric configurations and have regular metrics at the origin. Since Birkhoff's theorem does not apply for modified gravity, the solution outside may deviate from Schwarzschild-de Sitter spacetime. At the boundary, Post-Newtonian parametrization was used to determine how much the studied model deviates from the general relativistic $\Lambda$CDM model.


I. INTRODUCTION
According to observations [1,2], a dark energy (DE) component is needed to explain the apparent acceleration of the Universe. In the ΛCDM model, Einstein's general relativity (GR) is modified by a tiny constant Λ. With this deviation from the bare Ricci scalar in the Einstein-Hilbert action, the universe accelerates according to cosmic background, supernovae and large scale observations [1][2][3][4]. Various other means to introduce the accelerated expansion into the field equations include e.g. cosmological scalar fields with gravitational interactions [5], as well as vector fields [6], string theory and higher dimension inspired models [7] as well as several type of gravity modifigations [8].
f (R) modified gravity theories [9,10], where the Einstein-Hilbert action is generalized as a function of the Ricci scalar, provide natural alternatives for dark energy that needs new physics to explain the accelerated expansion of the universe. In these models some subdominant terms in the action functional become essential in very large scales.
In the present paper, the f (R) theories of modified gravity are studied. that include an action function with a suitable function of the Ricci scalar. The observed large scale defects could then be explained by the deviation from Einstein-Hilbert gravity. However, Einstein-Hilbert gravity is well tested within solar system length scales and at low energies [11]. Therefore we, focus on models considered as alternatives for the concordance ΛCDM model with GR like behaviour at solar system scale.
When the modification in the action functional is written using auxiliary scalar fields (as in quintessence and * Electronic address: kajohe@utu.fi † Electronic address: vilja@utu.fi scalar-tensor teories [5,8]), the new effects of gravity can be interpreted as scalar-field -matter interactions. Chameleon f (R) models [12] allow the mass of the scalar field to vary according to the local matter density and give rise to non-linear self-interactions. With chameleon behaviour, correction to the Newtonian gravitational potential may be small in high curvature regimes (e.g. near the sun) and the metric therefore close to Schwarzschildde Sitter solution.
There are many f (R) models that can produce correct shift from matter domination to radiation domination and to the currently accelerating epoch (see e.g. [13,14]). Also, constraints from structure formation and stability conditions are satisfied by various models [15][16][17][18][19]. Since gravity is most rigorously tested at local scales and with low energies, tangible f (R) model must be able to settle to the observed GR values within this regime. f (R) models that pass all the above mentioned tests are scarce.
Among current viable f(R) models are [20][21][22][23][24]. In the present paper we study Hu and Sawicki f (R)-model (HS) numerically to test if it locally deviates from GR. It has non-linear rational action functional and seems to follow locally the chameleon mechanism.
While the metric outside a spherically symmetric body in GR must be the Schwarzchild-de Sitter space-time according to Birkhoff's theorem [25], this does not hold for a spherically symmetric matter distributions in nontrivial f (R) theories, like in the HS-model. In these theories, the way space-time outside the mass distribution depends on the solution inside, is more complex [26,27]. Therefore, we have compared static spherically symmetric bodies with a polytropic equation of state within the HS model and in Einsteinian gravity.
We also calculate the mass of auxiliary scalar degree of freedom and find it to be high enough to evade the "fifth force" problem [12]. This paper is organized as follows: first in Section II the gravitational framework is discussed, in Section III the polytropic configuration and the surrounding space-time is constrained and in Section IV the numerical work is represented. Finally in Section V we discuss the results.

II. GRAVITY FORMALISM
We consider the f (R) gravity action [20] of the form where κ 2 = 8πG, and L m is the matter Lagrangian. The corresponding field equations, derived in the metric formalism, are (2) Here T µν is the standard minimally coupled stress-energy tensor for perfect fluid and F (R) ≡ df /dR. Contracting the field equations we get the trace equation which can be used as one independent dynamical equation.
As usual, the matter term needs to obey the equation of continuity D µ T µν = 0. The only non-trivial component for spherically symmetric system is where comma stands for radial derivative, ≡ d/dr. Like in GR, the equation of continuity is automatically satisfied here if the field equations are satisfied [28] and no additional information is gained. In this work we chose the continuity equation to be one of the solved dynamical equations instead of using the full set of field equations. Therefore, as the set of independent equations we use the 11-component of equation (2), the trace equation (3), the definition for the Ricci scalar in terms of metric components and the continuity equation (4).
As there are still several physically acceptable f (R) models we chose the Hu-Sawicki model for numerical purposes. This model is able to reproduce ΛCDM at high curvature when the f (R) function is expanded. It also reproduces correct matter domination and the current accelerated eras and avoids Ostrogradski and matter instabilities [15][16][17]20]. This particular model also vanishes for flat space-time, f (0) = 0, and is able to produce the asymptotic limit f (R) ∼ R + Λ as R → ∞.
The action function for the Hu-Sawicki model is given by [20] f (R) This contains three essential parameters, a positive real number c 1 , c 2 and a positive integer n. The mass scale m can be always included in parameters c 1 , c 2 , but it is convenient to use it as the mass scale as in the paper [20]: where h is the reduced Hubble parameter and h = 0.72 in our work. The authors of [20] study the solar system behaviour of the model by an approximated f (R) function. They conclude, that the model should correctly arrive at the solar system value outside a sun like matter distribution because the scalar degree of freedom is locked to the general relativistic value out through the solar corona [20]. The approximate solution for (5) in high curvature (i.e. when m 2 /R → 0) gives accordingly an effective ΛCDM model when the f (R) function is expanded. By approximating the f (R) function to begin with, the non-linearity of the modified equations that allow the Birkhoff's theorem to be violated, is not taken fully into account.
In the present paper the above assumption is tested by using the full non-linear, rational action function with the non-linear field equations. We, therefore, numerically study whether the approximation of the rational function (5), that makes the solution fully merge with the general relativistic space-time around a polytrope, is appropriate. This is important because the model is currently used as a viable DE model (see e.g. [29]).
One useful way to study this family of f (R) models is to write (5) as: Now one finds that only when the Ricci scalar R is near the vacuum value m 2 , can the third term contribute to ΛCDM . However, because of the non-linearity of the field equations, behaviour of the solution may be different from ΛCDM .
The system was first required to give the correct static vacuum value R vac = 12H 2 0 with the trace equation (3). This restricted the initial-parameter space considerably. Moreover only analytical and regular solutions at the origin were considered, by demanding the solutions to reach a finite value at the origin [30].

III. STATIC SPHERICALLY SYMMETRIC SOLUTIONS
We consider spherically symmetric, static configurations and adopt the metric for this configuration: SdS-metric describes an expanding vacuum solution of the Einstein field equations. In order to test whether the studied modified gravity polytrope can really reproduce the general relativity weak field limit, the metric of the HS-polytrope is compared to the SdS-metric on the boundary of the spherical matter distribution, Further, parametrized Post-Newtonian parameters are compared to observations in order to determine the relevance of the SdS-solutions in the solar system. Derivation of the γ P P N parameter by fitting the general metric (8) to the PPN SdS-solution, was done in [30]. The used parameter is here the cosmological term H is neglected because it is vanishingly small within solar system.

A. Polytropic considerations
We use standard polytropic configurations that resemble young sun-like solar mass main sequence stars. We chose a solar polytrope as the reference star with polytropic index n p = 3. This is the so called Eddington model and it is given by the general relativistic Tolman-Oppenheimer-Volkoff (TOV) equation. This model is able to reproduce the observed solar mass and radius rather well with a simple set of equations and initial conditions.
As self gravitating globes of gas, the stars are kept in hydrodynamic equilibrium with the equilibrium condition Here F α is the gravitational force [25]. The TOV equation is defined with the hydrostatic equilibrium (10) and spherically symmetric metric (8).
As static stellar objects, we consider perfect fluid solutions with polytropic equation of state and parametrization The index n p = 3 gives the Eddington polytropic solution for the Lane-Emden (LE) equations for solar mass M and radius R . According to [31] the solar polytropic model with index n p = 3 is not stabile with respect to Jacobi stability analysis. The solar polytropic index is stabile within linear stability analysis and Kosambi-Cartan-Chern theory (see [31] and references therein). The polytropic index values:  [33], Very Long Baseline Interferometry (VLBI) [34], Lunar Laser Ranging (LLR) tests of general relativity [35] and the Mercury perihelion shift measurements [36,37]. The constraints for βP P N are not discussed in this paper, as the Cassini mission results give the most stringent limits for SdS here even for βP P N = 1 would conform to also Jacobi stabile polytropes. However, when departing from the solar index the radiusmass relation changes dramatically. The behaviour of the polytropic star with the used field equations would also change for this index and new degrees of freedom would be needed in the analysis. The interconnected Lane-Emden coefficients (K, γ in (11) and scale length l in LE equations [25]) need to be chosen correctly in order to obtain a star that best describes the sun. Lane-Emden equations for a polytropic configuration are derived from the general relativistic Tolmann-Oppenheimer-Volkoff equations. We also derived general relativistic ΛCDM model, for reference, using the ansaz described in Section (II). This model is of course identical to the LE polytrope. The standardized Eddington model was chosen for this study, because it is considered to be in good resemblance with the standard solar model [32] with the central density We varied the central density of the Eddington polytrope in order to find statistically significant differences or similarities between the modified gravity polytropes and general relativistic ones. Also, other realistic main sequence stars (M HS = M ) were searched from larger polytropic parameter space with f (R) gravity. But as in the case of general relativistic n p = 3 polytropes, the only possible configurations are those with the solar mass. This occurs because (see e.g. (11.3.18) in [25]) The mass M is constant for the Eddington model n p = 3.

B. Space-time constraints
If the studied f (R) modification (5) is to produce general relativistic gravitation within the weak field of solar system, the Schwarzschild-de Sitter space-time is anticipated as the resulting space-time around the calculated matter distributions. This metric is also a fair approximation around slowly rotating bodies, like the sun, so we calculate the Post-Newtonian γ P P N parameter at the boundary to see how well the general relativistic value, γ P P N = 1, is reproduced.
The Post-Newtonian parametrization (PPN) gives the strongest bound for a modified gravity theory within the solar system. Observational constraints are presented in the Table I. The Cassini mission [33,37] measurement constrains the space-time parameter γ P P N most stringently within the solar system: Therefore, this Cassini constraint is used as reference for choosing the acceptable solutions that mimic general relativistic ΛCDM polytropes.

IV. NUMERICAL STUDY
We study the local effects of the f (R) model by integrating the matter distribution, ρ(r) from the center of the polytrope to the boundary. In the numerical integration the center is defined to be under one meter. SI units with c = 1 were used throughout this study. The regularity at the origin was separately tested for the considered set of equations.
The calculations were highly dependent on the equation ansatz, thus we used one of the modified field equations (2), equation for the Ricci scalar and the continuity equation (4) to overcome numerical difficulties.
We solved the dynamical equations for the metric components A(r), B(r), the energy density ρ(r) and the scalar curvature R(r), with respect to the radial coordinate r (as was done in [26,30]).
The space-time around the polytrope is found by integrating from non-singular central region with the most general spherically symmetric metric (8) outwards until a boundary is found at ρ(r) = 10 −5 kg m −3 , as sun's photosphere lies at ρ ps = 2 × 10 −4 kg m −3 .
This definition gives a discontinuous surface, since we assume the space-time to settle to the vacuum SdS value at the defined surface of a polytropic star. The Solar system potential or galactic potential is not discussed in this paper, because all the observations here are made either inside the polytrope or at the surface. For reference, the same set of equations in II were also used with f (R) = R + Λ with the same numerical parameters. This is the Eddington polytrope and is also denoted as ΛCDM model.
We parameterized the polytropes (11) as is done when deriving the Lane-Emden equations: with (14), and could therefore variate around the solar polytrope parameters. i.e. ρ i and γ were determined by the Eddington model that fits the standard solar model well and θ 0 was varied. This was done in order to find out how well the f (R) polytropes satisfy the stringent observational condition for γ P P N (16) when compared to the Eddington polytrope. By variating θ 0 we were able to easily find statistically significant differences between these models, because γ P P N parameter depends on the central density.
The model depends effectively only on the central density θ 0 in (17) as is shown in Chapter IV A, since the scalar curvature effectively follows the energy density R(r) ∼ κρ(r).
Similar parameterization was also done with the ΛCDM model and the resulting polytropes, although very similar, differ systematically.

A. Constraining the parameter space
For studying the possible space-times around the f (R) polytropes we scanned a large parameter space of (n, c 1 , c 2 , R 0 , θ 0 ). The initial values A 0 , B 0 for the spherically symmetric metric coefficients (8) were chosen by using the regularity condition at the center when r → 0. The metric parameters θ and φ (8) do not enter the equation ansatz.
We first constrained the parameter space of the f (R) function by considering only those parameters c 1 and c 2 , that give the correct static vacuum value for the Ricci scalar R = 12H 2 0 . The allowed values include only solutions that reproduce the correct vacuum with the trace equation (3) for The Hu-Sawicki parameters c 1 and c 2 , due to this requirement, have a fixed ratio. The parameters n, c 1 , c 2 were scanned over values n ∈ [1,12], c 1 and c 2 ∈ [1, 10 10 ] and the results seem highly degenerate over these parameter values. θ 0 was varied, in order to gain statistical significance, through the values that produce physical central densities of red giants ρ 0 /ρ 0,Edd 0.01 to a very dense star with central density ρ 0 /ρ 0,Edd 1000.
The Hu-Sawicki model has also been compared to the Type Ia supernova (SneIa) and gamma ray burst (GRB) data, by [38]. To note, the parameter values favoured by SneIa and GRB data with 95% CL in that work, do not conform to the correct static vacuum solutions for this model.

B. Behaviour of HS polytropes
For the scanned range of initial scalar curvature values the space-time outside the sphere is near to, but not exactly, either SdS or Einsteinian (corresponding to f (R) = R) space-time This is easy to see if the function f (R) is written as in (7) (17), such that the mass of the polytrope is always the solar mass M . The density parameter gives the solar value with θ0 = 1. There also exists f (R)HS-polytropes almost identical to GR LE solution with γP P N within observationally acceptable limits (16). The f (R)-polytropes shown in this figure, however, have the γP P N value out of this observable range and are not observationally favorable.
Hu-Sawicki polytropic solutions are, indeed, very similar to general relativistic polytropes. The energy density always tracks to the general relativistic ΛCDM solution throughout most of the stellar radius and, therefore, the radius and mass of the HS-polytrope are similar to the GR solution.
A key issue that distinguishes GR and f (R) HS polytropes, is the space-time outside the sphere. The γ P P N parameter is not sufficiently close to 1 around the f (R) polytropes for the majority of the solutions.
The solutions show unstabile numerical behaviour for the scalar curvature R in the central regions and near the surface. The solution, however tracks to fulfill R(r) κρ(r) for almost all the radius. Because of the tracking behaviour, initial scalar curvature can be almost anything and still produce a GR-like solution. Generally with a non-SdS space-time, the only effective initial parameter for distinguishing Hu-Sawicki model from the ΛCDM model is, therefore, the central density parameter θ 0 in (17). For the central density of the Eddington polytrope (ρ = 76.5 g cm −3 ), the space-time parameter outside the polytrope varies around the SdS-value with a mean value close to the ΛCDM value (see FIG.2).
Because the equations in II are of higher order than the field equations of general relativity, the conditions at the surface of the star need to be more stringent. For example, the metric parametersB andȦ in (8), need not be continuous at the surface of the star in general relativity but in the case of fourth order gravity, are required to be continuous. The reference value for a general relativistic Eddington polytrope at the defined surface approaching from inside is γ P P N,ΛCDM = 0.9999991717. It is exactly one once the vacuum is reached, since the polytropes surface is discontinuous at the border by the definition of our polytropes.

C. Numerical results
For a wide range of Hu-Sawicki parameters and initial values for the energy density parameter θ 0 in (17), the solution is similar to the general relativistic polytropes. The notable difference is the space-time around the configurations, that turns out to be non-SdS for roughly half the f (R) stars per scan.
All the polytropic modeling was done with the full non-linear field equations and non-linear rational function (5). The considered numerical accuracy did not suffice for finding any solutions for the used set of equations (2,3,4) with n ≥ 9.
Overall, for n = 1 . . . 8 in (5) with the chosen ansatz, higher than 13 decimal-accuracy was not feasible for a large set of solutions. This accuracy, however is sufficient for finding an adequate γ P P N value. With higher accuracy most of the equations were not solvable because of the high degree of non-linearity. When the accuracy is improved, the amount of successful polytropes drops dramatically because of numerical difficulties. Even with increasing accuracy the mean γ P P N value does not converge to the ΛCDM value, as is shown in FIG.3.
The solutions also have numerical instabilities in the Ricci scalar, R(r), near the center and surface. The solutions, however, always track to fulfill R(r) = κρ(r) effectively as the resulting configuration.
Outside the correct parameter area confronting the static vacuum solution for the trace equation, f (R) polytropes do not yield γ P P N = 1/2 outside, as in the case of the simpler polynomial f (R) models. Rather, the γ P P N value lies near 1 everywhere except if the initial value for the scalar curvature is that of the vacuum.
All the solutions were required to be regular at the origin. The regularity condition was not, however, sufficient to force the solution to be the Schwarzchild-de Sitter space-time This is because around 60% of the polytropes do not fulfill the observational requirement of SdS given by (16).
The Hu-Sawicki parameters are degenerate when considering the mass, radius and γ P P N of the polytrope. This value is, however, dependent on the central density of the polytrope. To get statistically relevant comparison to the ΛCDM model with constant γ P P N , polytropic solutions with modified gravity were scanned over a wide range of central densities.
As a result, the PPN parameter values typically lie within |γ P P N − 1| < 1 × 10 4 , e.g. the solutions usually tend to a SdS-like space-time. SdS space-time is, Only the observationally acceptable values are plotted. In FIG. 3 the whole set of polytropes is shown. In both figures, the ΛCDM value is also plotted for reference. It is always closer to the SdS value, γ P P N = 1, than the Hu-Sawicki mean value is just inside the surface. Mean γ P P N value for the Hu-Sawicki polytropes converges to a value, that lies within the observational bounds, but is always different than the ΛCDM value.
We found that the mass of a f (R) polytrope is always closer to the observed solar mass than the mass of the Eddingon polytrope. The difference in the masses M HS /M 1.002 and M GR /M 1.048 is systematic for all the n p = 3 polytropes. Although the solutions are similar with respect to the initial values and obtained radii and masses, the solutions never identical.
The f (R)-model is stabile with respect to shorttimescale instabilities [40] if f (R) is positive. For the present model, this stability condition is satisfied with ex-   (17). f (R) polytropes with θ0 = 0.5 and 2.0 are also shown in Fig.(1). ΛCDM values were computed for reference with the same code. The Lane-Emden solution is very similar to θ0 = 1.0 ΛCDM solution in the table. Also, the spacetime parameter is not dependent on the initial density parameter, bacause of the instabilities in the curvature near the center. Note, the space-time parameter for ΛCDM model is (γP P N -1)=−8 × 10 −7 in our work.
tremely small f (R) everywhere inside and near the polytrope. The mass term within the corresponding Brans-Dicke framework, is dependent on a 1/f (R) term. The mass of the scalar, m 2 σ , is allways big enough for the "fifth force" problem to be evaded for these polytropic configurations.

V. CONCLUSIONS AND DISCUSSION
Although, Hu-Sawicki polytropes are very similar to the young sun like stars in GR, only one third of the HS polytropes in a typical 29 000 polytrope run have γ P P N parameter within the experimental Cassini limits [33].
Many f (R) models result in γ P P N = 1/2 [26,39], far from the GR value γ P P N = 1. The studied HS-model is different in this respect although not fully consistent with the observations. For the Hu-Sawicki -model, γ P P N value is nearly one everywhere in our model parameter space. However, the f (R) HS space-time around a spherically symmetric configuration with a polytropic equation of state nearly always converges to a different space-time than is allowed by GR. The non-linearity of the system poses difficulties, and very high accuracy was not feasible. With a sufficient numerical accuracy all the n < 9 models were solvable. Furthermore, even as the numerical accuracy was increased γ P P N,HS value does not converge to γ P P N,ΛCDM .
With only 40% of the γ P P N values reaching the observational limits with the Hu-Sawicki gravity, it is obvious that HS-model does not easily reach the Schwarzchild-de Sitter solution outside a polytropic stellar configuration.
This model, therefore, seems to require some finetuning of the initial scalar curvature in order to reproduce general relativistic SdS space-time outside the surface of a sun-like HS polytrope. The scalar curvature is always unstabile near the center, so no specific value can be attached to the HS polytrope, that would exactly correspond to the general relativistic polytropic sun.
To get away from this fine-tuning a physical super selection rule might exist. We do not suggest any mechanism for this, however, we would like to note that the sun is the only star for which the SdS space-time is mea-sured to be correct with good accuracy. Could it be, that other sun like young stars that can be modeled with a polytropic equation of state, might arrive at a slightly different vacuum than SdS? In this case HS -model would, indeed, describe the space-time around young main sequence stars more accurately.
The short-timescale instabilities are avoided with extremely small and positive f (R). Also, the scalar field mass within the Brans-Dicke scenario becomes large, as m σ contains a term proportional to 1/f (R).
The solutions do not change notably in the mass, radius and γ P P N when the Hu-Sawicki parameters are varied. This indicates that only the ratio c 1 /c 2 , which is fixed, is effective for R m 2 . However, the γ P P N value is naturally dependent on the initial central density ρ 0 , that was varied throughout the numerical studies.
With a different polytropic equation of state, some white dwarfs and neutron stars can be modeled. Work is under way to test the space-time specifically around f (R) low density white dwarfs and neutron star -polytropes. Moreover it would be interesting to study numerically also other proposed well-behaving action functions.