Modified Gravity at Astrophysical Scales

Using a perturbative approach we solve stellar structure equations for low-density (solar-type) stars whose interior is described with a polytropic equation of state in scenarios involving a subset of modified gravity theories. Rather than focusing on particular theories, we consider a model-independent approach in which deviations from General Relativity are effectively described by a single parameter $\xi$. We find that for length scales below those set by stellar General Relativistic radii the modifications introduced by modified gravity can affect the computed values of masses and radii. As a consequence, the stellar luminosity is also affected. We discuss possible further implications for higher density stars and observability of the effects before described.


Introduction
It is well established that General Relativity (GR) provides an accurate description of the gravitational interaction from the submillimeter scales probed by torsion balance experiments (Adelberger et al. 2003(Adelberger et al. , 2009 to kiloparsec distances as confirmed by recent observations of strong gravitational lensing of extragalactic objects (Collett et al. 2018). These tests probe the weak-field regime of the theory, and so far no discrepancy with respect to the predictions of GR has been found. For strong fields, the theory is still poorly tested, but the recent discovery of gravitational waves produced in the merger processes of black holes (Abbott et al. 2016) or neutron stars (Abbott et al. 2017) has opened a new avenue to explore this regime of the theory. Despite the success of GR on small (subgalactic) scales, the difficulties of accommodating the observed accelerated expansion of the universe within the theory has led to the suggestion that the universal attractive character of gravity could break down on cosmological scales. Several modified gravity (MG) theories have been proposed in the last few years that introduce additional degrees of freedom, typically scalars, that mediate the gravitational interaction, thus changing its behavior on very large scales.
Despite the fact that some of the theories beyond that of Horndeski have been practically ruled out by the observation of the GW170817 event, some of them are still viable modifications of gravity (Ezquiaga & Zumalacárregui 2017. For a discussion see Baker et al. (2017), Creminelli & Vernizzi (2017), and Sakstein & Jain (2017). Thus, in such theories the equation for hydrostatic equilibrium in the nonrelativistic regime is modified in such a way that (Koyama & Sakstein 2015) where P is the pressure, m(r) is the mass inside a radius r, ρ is the energy density, G is the gravitational constant (Patrignani & Particle Data Group 2016), and ϒ is a constant parameterizing the deviation with respect to GR. Existing bounds −0.22< ϒ<1.6 are set by the Chandrasekhar limit on WD stars (Babichev et al. 2016;Jain et al. 2016) and the minimum mass of main-sequence stars (Sakstein 2015a;Babichev et al. 2016;Saltas et al. 2018). On more general grounds, effective descriptions of alternative gravity theories have been developed in recent years aiming to encapsulate in a few parameters all the relevant modifications at a given scale (Clifton et al. 2012;Silvestri et al. 2013). These parameterizations have been widely employed in the analysis of structure and lensing data of galaxy surveys (Amendola et al. 2018).
In this work, we will explore the implications for stellar structure of one of the simplest and most widely used effective parameterizations considered in the literature (Bertschinger & Zukin 2008;Silvestri et al. 2013). For a wide class of theories of gravity with one additional scalar degree of freedom, assuming the quasi-static approximation and not higher than second-order equations of motion, it can be seen that all the relevant modifications can be encoded in an effective Newton's constant parameterized in Fourier space by m , with G as the ordinary Newton's constant, and a gravitational slip parameter g f y = ( ) ( ) ( ) k k k . The μ(k) parameter changes the hydrostatic equilibrium equation introducing, generically, a new length scale in the dynamics and leading to expressions that can deviate from Equation (1). This effective approach will allow us to analyze the potential modifications in the stellar structure, including radius, mass, luminosity, and temperature of stars, in a largely model-independent way.
The manuscript is structured as follows. In Section 2, we briefly review and expose the main features of the MG theories we consider in our analysis of the changes induced in the stellar structure equations. We also present the main scenario we explore, i.e., low-density stars, and discuss the validity of our approach. The results are then detailed in Section 3, and final conclusions are drawn in Section 4.

Stellar Structure Equations
In this section, we start by presenting the standard setup of the calculation for the structure of a stellar object of mass M and radius R in GR. We consider for the sake of simplicity a nonrelativistic and nonrotating object. Following Silvestri et al. (2013) we consider a spherical Minkowski metric with linear scalar perturbations through the introduction of two radial potentials, ψ(r) and f(r), both fulfilling the condition of weakfield approximation, i.e., ψ(r)=1, f(r)=1.
is written as where the four coordinates of x μ are t, r, θ, and f. From the selected metric, g μν , and using the approximation of a perfect fluid, the energy momentum tensor can be written as r = + + the Einstein's tensor given by the Ricci tensor, mn R , and the Ricci scalar, R, for this metric. Explicitly, these expressions read to first order in metric perturbations , 5 and . 6 In addition, for the spatial diagonal component of the energy momentum tensor we can approximate to first order in metric , being i=1, 2, 3 a spatial index.
As obtained, it is clear that Equation ( It is important to note that Equations (10) and (11) along with the mass equation Equation (8) are the structure equations for the potential, pressure, and mass in the star. Besides, using Equations (5) and (6) we get y f -= + C C r 1 2 2 , with C 1 and C 2 as constants. If we impose the existence of a finite solution when  ¥ r , it follows that C 2 =0. On the other hand, for r>R we have to recover the Schwarzschild metric, which means ψ=f. In this way, C 1 =0 and ψ−f=0, making P=0. This result is due to the fact that we obtain a first-order perturbation solution, being y f~-( ) Gm r r . Then, whereas ρ and ψ are first-order functions in perturbation theory, P is a second-order function in perturbations. We can thus rewrite Equations (10) and (11)  retaining only the leading contributions, i.e., first and the second order, respectively. From Equation (12) the mass of the star can be solved as

Input from MG Theories
In the context of MG theories and for static configurations following Silvestri et al. (2013), we introduce two functions, μ(k) and γ(k), in the Fourier k-space whose effect is modifying the equations governing the solution of the potentials ψ, f. We can write y p mr  = ( ) G 4 1 5 2 and f gy = ( ) .
1 6 The physical meaning of the μ(k) function is that of providing an effective value of the gravitational constant, G. Instead, γ(k) establishes a relationship between the two potentials, ψ and f. When μ=γ=1, we recover GR equations. The most general expression for μ(k) in theories with one extra scalar degree of freedom and modified Einstein's equations involving up to second-order derivatives can be cast into the rational form (Silvestri et al. 2013) The previous step stems from the fact that it is indeed equivalent to rewrite Equation (15) as which is a local expression. Keeping this in mind we can now obtain the expression of μ in position space as 5 as constant parameters. Thus, for example, by defining b = p p 1 3 5 it has been shown (Bertschinger & Zukin 2008;Giannantonio et al. 2010) that several models such as f (R) or certain Chameleon theories correspond to β 1 =4/3. The parameter values corresponding to other models like Yukawa-type theories belong to the interval 0.75<β 1 <1.25. Furthermore, as we will show below, γ(k) plays no role in our study, so we do not provide any particular parameterization for it. In any case, we should keep in mind that more general parameterizations exist (Hojjati et al. 2014) that could also include additional vector degrees of freedom (Resco & Maroto 2018).
As can be readily seen, under the conventions used, p 3 and p 5 have units of squared length, while p 4 is dimensionless. Using Equation (19) we can rewrite Equation (15) as As we want to consider small perturbations from GR, in what follows we demand p 4 does not depart largely from unity and p 3 ∇ 2 =1 and p 5 ∇ 2 =1. In this scenario, for spherically symmetric configurations, Equation (8) now takes the form where the differential operator μ −1 is given to first order in the perturbative expansion by From integration of Equation (21) we get to first order In addition, we can also write which is the modified version of Equation (12). Finally, At this point it is important to notice that the subset of MG models we are considering only hold in the linear regime, which we assume to be the one valid inside solar-type stars (and other high-density stars such as WDs). In this sense, we must emphasize that some proposed modifications of gravity on cosmological scales that rely on screening mechanisms (Koyama & Sakstein 2015;Saito et al. 2015) through nonlinear effects on astrophysical scales are not covered by this approach.
Notice that because of the form of the effective operator μ, the structure of the pressure equation is different from that shown in Equation (1). A few comments are in order. For simplicity we have limited ourselves to the case in which the parameters p 3 , p 4 , and p 5 are just object-independent constants. In the most general case, for the spherically symmetric configurations considered in the work and unlike the cosmological case discussed in Silvestri et al. (2013), these parameters could be general functions of the stellar radius. Thus, for example, the p 4 parameter related to the effective Newton's constant could interpolate between different values in the interior and exterior of the star. However, in our case with constant p 4 , matching with the laboratory-measured value of Newton constant G imposes p 4 =1. Thus, the combination G/p 4 appearing in the obtained structure equations inside the stellar body will effectively take the value G. The possibility of changing the effective Newton's constant in astrophysical objects as compared to the laboratory value has been considered phenomenologically in previous works in nonweak-field limits; see, for instance, Velten et al. (2016) and references therein, where an α parameter is used analogously to our p 4 or in cosmological contexts, and see Planck Collaboration et al. (2016) for explicit measurements of μ 0 .

Perturbative Solution for a Polytropic Star
We have seen that the potential ψ, the density ρ, and the pressure P are related through Equation (13). In this way, given an equation of state (EoS) determining a relation between P and ρ we could, in principle, obtain an expression for ψ(ρ).
In this work, we use a polytropic EoS that can be considered a reasonable description for main-sequence (solar-type) stars (Hansen et al. 1994). In addition, we will also consider the higher-density case of WDs. No dark matter presence is assumed where n is the so-called polytropic index, which is related to the internal constituents of the star, and K is a constant with appropriate units. In this work, we use CGS units. Therefore, using Equations (13) and (27), for n>0 we have In order to solve the differential equation in Equation (29) using a perturbative approach, one can propose a solution under the form ψ=ψ 0 +ψ 1 , with y y  | | | | 1 0 . Thus, a set of two differential equations is obtained, one for each component, ψ 0 , ψ 1 : Let us mention that this same result would follow from the alternative calculation using a perturbative expansion of ψ and ρ in Equation (18). Notice that the perturbative approach allows us to solve separately the gravitational potential equation from the conservation equation shown in Equation (13), which remains unchanged with respect to the ordinary GR case.
In order to calculate the radius R and mass M we must first obtain the solution ψ(r), which will provide us the density ρ through Equation (28). By imposing ψ(R)=0, the radius R will be calculated from the usual consideration that both pressure and density vanish for r=R.
Accordingly, the total mass of the star, M(R), is obtained using Equation (28) and integrating Equation (8) It is important to note that if we impose ξ=0 we will recover the GR case. In order to be more definite we now particularize to the perturbative solution for a polytropic n=3 star focusing on the case of low-mass solar-type stars MM e . Most of the internal description is reasonably well described by a polytrope with an n=3 index. Refinements to this description would involve, at least, the use of several polytropes or, more appropriately, an improved version (Hansen et al. 1994) of the standard solar model (SSM; Turck-Chieze & Couvidat 2011). Note, however, that the goal of our calculation is not a detailed modeling of the stellar object but rather showing the effects of selected MG models, thus we will restrict to a single polytropic EoS. Therefore, in this case, it can be written as In addition, Equations (30) and ( The density profile and the mass can be obtained as from which the radius R is derived through the relation it must fulfill, ρ(R)=0, which is needed to calculate the stellar mass For reference, solutions for the radii and masses of solar-type configurations in the GR case we consider in this work have been obtained using = --K 3.8 cm g s 3 2 1 3 so that we can recover values close to those of the Sun, M 0 =1.95674×10 33 g and R 0 =6.8141×10 10 cm. A central density value r r = = -80 g cm c 0 3 has been used. In addition, stellar luminosity, L, can be derived by considering the energetics taking place inside the stellar volume through the differential law where ò is the nuclear energy generation rate in units of -erg g s 1 1 . In general, ò can display a rather complicated expression, but its main contribution can be parameterized for a low-mass solar-type star with an active pp chain under the form (Maciel 2016)  where X is the proton fraction and T is the temperature of the star. Thus, the luminosity (in erg s −1 ) can be derived as For most stars (with the exception of very low mass stars and stellar remnants) the ions and electrons can be treated as an ideal gas and quantum effects do not critically affect their behavior. In our treatment and in order to keep our modelization simplified, we will consider that the radiation pressure is much smaller than that of the gas of ions and electrons in the stellar plasma, i.e.,  P P r gas . Note that for the particular example of the Sun core~-P P 10 r 4 gas . As mentioned, a more general treatment (Koyama & Sakstein 2015) would involve considering a mixture of both pressure components.
Therefore, in our scenario, the stellar conditions are dominated by the gas pressure » =

Results
In this section, we analyze the results obtained for several magnitudes of interest we have calculated in the context of the MG modelizations under study. In order to solve Equations (34) and (35), we take = --K 3.8 cm g s 3 2 1 3 so that we can recover typical values of radius and mass of solar-type stars, R 0 and M 0 , used as normalization. Although we do not attempt to accurately model stellar values of relevant magnitudes, it is worth mentioning at this point that in our polytropic description most of the stars, i.e., in the region <  r R 0.7 , are well approximated by a polytrope of index n=3, while for the outer region, with convective behavior, it would be better described by a polytrope of index n=1.5. The latter corresponds, however, only to 0.6% of the mass.
In Figure 1, we consider the case of a low-density solar-type model with a central density r r = c 0 (r = -80 g cm 0 3 ). We show the solution potentials, i.e., the potential before introducing the perturbation, ψ 0 , and the perturbation, ψ 1 , fulfilling y y  | | | | 1 0 , as well as their sum y y y = + 0 1 , as functions of the radial coordinate r for γ=1 and ξ=−2.67×10 18 cm 2 . More generally, we obtain that only for values of  x | | 1.225 10 cm 19 2 does a perturbative correction  y y | | | | 0.1 1 0 indeed justify our framework. As we have mentioned before, the radius of the star in this approach can be obtained as the first zero of the full potential solution, ψ(R)=0. Once we know this value, the mass of the star is obtained as In Figure 2, we plot the M-R relationship for x =| | 4 10 cm 18 2 . We normalize to values obtained in the solar-type case R 0 and M 0 . The solid line corresponds to the case in which ξ>0 and the dashed line to ξ<0. We take r r Î [ ] 0.05, 5 c 0 to generate our data points. We can see that in all cases corresponding to MG the relation gets distorted from the n=3 GR case in which M is constant when varying ρ c . Objects along the M-R curve with a positive slope > dM dR 0 denote metastable configurations so that in our analysis they are discarded.
In Figure 3, the M-R diagram is shown. It has been obtained by varying the central density value in the interval r r Î [ ] 0.05, 5 c 0 . We use x =2.5 10 cm 17 2 (solid line) and x = 10 cm 18 2 (dashed line). We can observe how the shape of the M-R diagram slightly changes (compared to the flat = M M 1 0 result from GR) for different values of ξ, being more significant for objects with lower radii. The star will achieve a more massive configuration as ξ increases.
Stellar mass and radius can be described with an approximate fit as functions of (ξ, ρ c ). In this way, the stellar radius R(ξ, ρ c ) is given by For masses x r ( ) M , c , instead, the corresponding fit is In Figure 4, we plot stellar density and temperature (normalized to the central values) as a function of r/R 0 , fixing ρ c =ρ 0 for x =9 10 cm 18 2 (solid lines) and ξ=0 (dashed lines), the latter corresponding to the GR case. Central temperature is obtained using Equation (41) with m = 0.61, corresponding to hydrogen, helium fractions of solar-type stars X∼0.7, Y∼0.3 yielding a central value~T 1.21 10 K c 7 . As can be seen, the value of the density (temperature) profile for the MG case of x =9 10 cm 18 2 is only slightly changed when compared to the GR case. The temperature radial profile,   T(r), also governs the change in the luminosity of (low-mass) solar-type stars.
In Figure 5, we plot stellar luminosity as a function of ξ. We consider ρ c =ρ 0 (solid line) and ρ c =100 g cm 3 (dashed line). From considerations relative to solar emission uncertainties, a∼±1% flux accuracy in the SSM could accommodate variations of GR taking place below length scales x~10 cm 9 . Although variations expected for some objects may be indeed larger, it may be extremely difficult to disentangle the presence of such a component from an ordinary variation due to the complex dynamics of the stellar (solar) interior.
In order to size the impact of the models under study in the MG scenario on the luminosity and temperature, a prescribed correspondence with a blackbody spectrum p = L R T 4 2 eff 4 has been used in Figure 6. We plot stellar luminosity as a function of effective temperature, T eff . We consider values ξ=−4×10 18 , 0, and 4×10 18 cm 2 (blue, orange, and green lines, respectively). Values ξ<0 would correspond to metastable stellar configurations, but they are shown for the sake of completeness. We set values for the central density in the interval r r Î [ ] 0.1, 3 c 0 . We can see there is a correlation of luminosity-temperature, as expected, and thus a nontrivial dependence on the effective parameter of our MG modelization, ξ, which is clear for T eff >6000 K. Spectral types O, B, A, and F for main-sequence stars can display such high temperatures. Although we cannot expect to recover the meaningful (GR) Hertzsprung-Russell diagram with our simplified approach, the analogous logarithmic correlation we find under the form Log eff yields a weak variation in the slope, α(ξ), for this case at the ∼0.6%. See some values listed for reference in Table 1. One important difference obtained in our calculation with respect to the usual approach where ( ) T R , eff are noncorrelated variables is that for the GR case we do not recover the familiar value α=4, as in our treatment T eff is obtained from the R-dependent value of the luminosity. In other words, the MG model used determines the variation in the radius R, and thus the luminosity. Therefore, values shown in Table 1 are by no means a predictive output of the model, but rather an indication of the weak dependence of the ξ parameter.
In order to explore the variations introduced by the subset of the MG model we study in this work let us consider a higherdensity stellar object such as a WD. To model this kind of object we use a relativistic Fermi gas simplified EoS dominated by an electron component using a polytrope with n=3 and =´--K 4.8 10 cm g s 14 3 2 1 3 . In this way, we recover solutions to the structure equations Equations (34) and (35) yielding masses up to the maximum theoretical limit =  M M 1.44 (Chandrasekhar mass). Note that when using GR with an n=3 polytropic EoS, M values (but not radii) are constant when varying ρ c .
In Figure 7, we plot the WD mass-radius relationship for x = | | 10 cm 3 is a typical WD central density. As before, we keep for reference M 0 and R 0 for normalization in the axis labels using those given by solar-type stars. If we further restrict to models in which ξ0, which are those that provide    stable stellar configurations of M(R), it is clear that, given a ξ value, the maximum mass (M Ch ) will be obtained for the highest possible value of ρ c provided the perturbative solution still holds. Because the heaviest WD reported in the literature has a mass M/M e ;1.37 (Hachisu & Kato 2000), and the GR M Ch value is even higher than this, ∼1.44 M e , all values obtained in this perturbative framework are allowed. Thus, in the high-density setting presented the ξ parameter remains unconstrained.

Conclusions
Using a model-independent approach, we have studied the effects of a subset of MG theories involving one extra scalar degree of freedom on stellar structure. We use a perturbative technique capable of solving in a linearized regime. We have selected our main scenario concerning the application to lowdensity stellar objects. First, in a perturbative approach, we have obtained mass and radius solutions for low-mass solartype objects using a polytropic n=3 EoS. We have analyzed how internal temperature profiles and stellar luminosities are affected with respect to the reference case of GR. We provide a fit for masses and radii provided a central density value, i.e., M(ξ, ρ c ), R(ξ, ρ c ). We obtain that a change in the curvature of M(R) results when ξ changes from positive to negative values, connecting the latter to metastable stellar configurations. We also find that MG can affect stellar luminosity from its ξ dependence. Globally, the effective temperature from a µ L T eff 4 law results in a linear Log L-Log T eff behavior with a weak dependence on ξ. We anticipate that this could result in objects appearing brighter (as L is increased with any stable configuration ξ>0, but this seems hard to measure experimentally as internal stellar dynamics is complex. Furthermore, it seems challenging to disentangle this effect from other standard effects such as proper fluctuations of the star. Even the Sun is a weakly variable star, with its major source of fluctuation coming from the 11-year solar cycle and revealing a smaller periodic variation of about ∼±0.1%. It is nevertheless worth pointing out that both effects, involving standard and new physics, could indeed be present and need to be further studied. When a particular case of a high-density star is considered (WDs), we find that no constraint appears on the ξ parameter from calculated values of the Chandrasekhar mass as it is always larger than the GR value in the perturbative approach we use. Our treatment for both low-and high-density stellar case examples can help to understand the behavior of MG on small (not cosmological) scales.