Anisotropic compact stars in $f(R)$ gravity

We derive a new interior solution for stellar compact objects in $f\mathcal{(R)}$ gravity assuming a differential relation to constrain the Ricci curvature scalar. To this aim, we consider specific forms for the radial component of the metric and the first derivative of $f\mathcal{(R)}$. After, the time component of the metric potential and the form of $f(\mathcal R)$ function are derived. From these results, it is possible to obtain the radial and tangential components of pressure and the density. The resulting interior solution represents a physically motivated anisotropic neutron star model. It is possible to match it with a boundary exterior solution. From this matching, the components of metric potentials can be rewritten in terms of a compactness parameter $C$ which has to be $C=2GM/Rc^2<<0.5$ for physical consistency. Other physical conditions for real stellar objects are taken into account according to the solution. We show that the model accurately bypasses conditions like the finiteness of radial and tangential pressures, and energy density at the center of the star, the positivity of these components through the stellar structure, and the negativity of the gradients. These conditions are satisfied if the energy-conditions hold. Moreover, we study the stability of the model by showing that Tolman-Oppenheimer-Volkoff equation is at hydrostatic equilibrium. The solution is matched with observational data of millisecond pulsars with a withe dwarf companion and pulsars presenting thermonuclear bursts.


I. INTRODUCTION
Compact stars are the final stage of stellar evolution when the radial pressure, caused by the nuclear fusion in the inner core, is no longer able to oppose the gravitational forces. Gravitational compact objects involve white dwarfs, neutron stars, and black holes. The interior structure as well as the physical properties of compact stars are studied adopting the Einstein field equations of General Relativity (GR). Schwarzschild presented a first solution capable of describing the interior of the compact body in hydrostatic equilibrium [1].
The stellar core was first assumed to be a fluid with equivalent radial and tangential pressures [2]. However, the chances of generating anisotropies in compact stars are considerably higher due to the existence of strong interactions between particles. In this situation, the inner regions become non-uniform as discussed in [3]. Specifically, it is the fundamental nature of particles which can generate anisotropies. As a result, anisotropic compact stars can become more compact than isotropic stars. The consequences of local variations of relative inner fields were discussed using the equation of state and some interesting results were presented in [4]. Ruderman suggested that a nuclear system with very high densities indicates that anisotropic pressure can emerge in compact stars [5]. In such a case, the radial and tangential pressures are different within the stellar core. These problems can find room in the context of modified gravity.
Due to the issue that GR is not capable to address the whole gravitational phenomenology from ultraviolet to infrared scales, modified and extended theories of gravity have been proposed as effective approaches towards quantum gravity, from one side, and towards astrophysics and cosmology from the other side. These theories are considered very useful to investigate the expansion of the universe as well as the dark energy and dark matter issues which could be explained under the standard of geometrical effects considering curvature and torsion invariants into the gravitational action [6][7][8][9][10][11][12][13][14][15][16]. In particular, the so-called f (R) gravity has gained special interest due to the fact that it is a straightforward generalization of the Hilbert-Einstein action of GR. Various models of f (R) gravity have been studied to investigate problems like the late cosmic acceleration [17][18][19], the ghosts of gravitational field [20][21][22], the inflationary epoch [23] and the whole expansion of the universe [24][25][26]. In f (R) gravity, several basic results have been achieved like the categorizations of singularities, the causal structure and the classification of energy conditions [27][28][29][30][31].
Using such modified theories, it is possible to explain the collapse and stability of astrophysical objects which escape the standard picture of GR. For example the supermassive compact object resulting from the gravitational wave event GW190814 can be framed and interpreted in the context of f (R) gravity (see [32] and references therein).
In this framework, the study of compact stars has recently become an interesting research topic because observations are reporting several anomalous objects which cannot be explained by GR and standard equations of state. For example, anisotropic compact stars can be analyzed considering GR with a cosmological constant [33]. The compact structure can be discussed considering various forms of the energy density and radial pressure [34]. Singularity-free solution for compact stars can be derived in [35]. The homotopy perturbations procedure can be applied to obtain stable structures [36], and perturbative models of f (R) gravity have been applied in [37] to study neutron stars. Neutron stars with quarks cores can also be considered, and results can be physically consistent [38]. The quadratic and cubic forms of f (R) were discussed using numerical methods in [39], while anisotropic compact stars in the context of f (R) gravity were studied in [40]. Neutron and quark stars are studied in the framework of f (R) = R + αR 2 while the Brans-Dicke theory is taken into account in [41]. In [42], quark star models are derived with an equation of state in the non-perturbative form of f (R) gravity.
Besides these results, it is worth noticing that several modified gravities are receiving much attention to cure the shortcomings of GR like the cosmic accelerated expansion, the flat rotation curves of galaxies, possible wormhole solutions and other phenomena [43][44][45][46][47]. For example, a quadratic function of the Ricci scalar was studied by Starobinsky to address the inflationary behavior of early universe [48]. It was shown that higher-order f (R) gravity can solve issues related to supermassive neutron stars [37,42,[49][50][51]. The easiest way to prescribe f (R) gravity is assuming an analytic function whose first order is the Ricci scalar. In general, the equations of motion of f (R) are higher-order in derivatives and supply substantial classes of solutions that are different from GR. In the framework of f (R) gravity, the dynamical behavior of matter and curvature fields has been studied in several works. See, for example, [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69].
In particular, spherically symmetric vacuum black hole solutions in f (R) have been derived in [70][71][72][73]. Furthermore, from the existence of Noether symmetries, it is possible to find out other spherically symmetric solutions [74,75]. Using the same technique, it is possible to achieve also axially symmetric black hole solutions [76]. Non-trivial spherically symmetric black hole solutions for specific classes of f (R) are derived also in [77][78][79]. Due to the higherorder curvature terms in f (R), one can discuss strong gravitational fields in local objects. In this framework, many researches concentrate on the study of spherically symmetric black holes [80][81][82][83][84][85][86] and neutron stars .
Finally, it is worth noticing that f (R) gravity can be recast as scalar-tensor gravity [108] with a scalar potential of gravitational origin so the further degrees of freedom with respect to GR result as scalar field counterparts in the Einstein field equations [109][110][111][112].
The aim of the present paper is deriving a spherically symmetric inner solution assuming general forms of the Ricci scalar and the f (R) function. After, we want to test the physical reliability of this solution. According to the procedure, it is possible to construct self-consistent compact stellar models capable of describing anisotropic neutron stars.
The layout of the paper is the following: In Sec. II, a short summary of f (R) gravity is reported. In Sec. III, we apply the field equations of f (R) to a spherically symmetric metric with two different potentials for time and radial components. Assuming the form of radial component as well as of the first derivative of f (R), we are able to obtain the time component of the metric and the f (R) function. Starting from this information, we derive the density, the radial and the tangential pressures. In Sec. IV, we impose boundary conditions to match the interior with the exterior solution. We succeed in matching the constants characterizing the solution by a compactness parameter C which satisfies the constrain C << 0.5. In Sec. V, we discuss the physical conditions to show that the model is consistent with a realistic compact star. In Sec. VI, the stability, adopting the Tolman-Oppenheimer-Volkoff (TOV) approach and the adiabatic index, is discussed. We show that the model satisfies both these conditions. The results are matched with observational data in Sec. VII where we draw the conclusions.

II. A SUMMARY OF f (R) GRAVITY
In 4-dimension, the action of f (R) gravity is [113][114][115][116][117][118][119][120][121][122][123][124]: with the gravitational coupling κ = 8πG c 4 , where G is the Newtonian constant, c is the speed of light and g is the determinant of the metric. Here f (R) is an arbitrary function of the Ricci scalar R and L M (g µν , ξ) is the action of generic matter fields ξ which are minimally coupled to the metric g µν . Clearly, GR is recovered for f (R) = R.
The variations of action (1) w.r.t. the metric tensor g µν gives the non-vacuum field equations [125] Eµν where ✷ is the d'Alembert operator and f R = df dR . Here T µν is the energy-momentum tensor, defined as Field Eqs. (2) have the following trace: From Eq. (4), one can obtain f (R) in the form: Using Eq. (5) in Eq.
(2) we get [126] Eµν where The energy-momentum tensor T µν can be assumed with the following anisotropic form where u µ is a time-like vector defined as u µ = [c, 0, 0, 0] and ε µ is the unit space-like vector in the radial direction. In this study ǫ is the energy-density, p and p t are the radial and the tangential pressures respectively. The energymomentum tensor takes the diagonal form T µ ν = diag(−c 2 ǫ; p; p t ; p t ).

III. COMPACT STARS IN f (R) GRAVITY
In order to consider relativistic compact stars in f (R) gravity, let us assume the following spherically symmetric space-time: where a(r) and b(r) are unknown functions of the radius. The Ricci scalar of the metric (9) takes the form: with a ′ = da dr , a ′′ = d 2 a dr 2 and b ′ = db dr . For the line-element (9), the non-vanishing components of field Eqs. (4) and (6) have the form: Moreover, we introduce the characteristic density where R 0 is a characteristic radius. It can be used to re-scale the density, the radial and the tangential pressure respectively. We get the dimensionless variableŝ Eqs. (11) are four non-linear differential equations in six unknown variables a, b, F ǫ, p and p t . Using Er r and Eθ θ , one can define the anisotropy parameter of stellar structure as: which reduces to ∆(r) = 0 for isotropic systems. In order to find a physical solution of the above system of differential equations, let us assume the following forms for the radial component of metric potential g rr and the f (R) derivative function F as: This choice is physically motivated for the following reasons. Here y is a dimensionless variable and we assume the star can extend up to the maximal radius r = l. The parameters c 0 and c 1 are dimensionless.They will be determined later. Eqs. (15) show that when c 1 = 0, we recover GR. In other words, the constant c 1 is responsible for the deviation from GR. Using Eqs. (15) into (14), we get the metric potential g tt as where c 2 and c 3 are given by combinations of the parameters c 0 , c 1 and l. The anisotropy parameter ∆(y) has the form Using Eqs. (15) and (17) in Eq. (11), we get the exact form of energy, radial and tangential pressures that we report in Appendix A. Inserting Eq. (15) in the last equation of system (11), the trace of the field Eqs. (6), we get a very lengthy differential equation in the unknown f (R). We report the full equation in Appendix C. Developing such a differential equation up to the first order in O 1 y , we get all the information we need for our purposes. It is where m = 5c 0 2 c 2 2 + 7c 0 c 2 c 3 − c 3 2 , m 1 = c 3 + c 0 c 2 and m 2 = 173c 0 3 c 2 3 + 138c 0 2 c 2 2 c 3 + 48c 0 c 2 c 3 2 − 8c 3 3 . The solution of the above differential equation takes the form where c 4 , c 5 and c 6 are integration constants. For c 5 , the solution (20) is compatible with F (y), given by Eq. (15), up to the leading order. Equations of density, radial and tangential pressures show that, for c 1 = 0, we return to GR. This is due to the fact that, for F = 1, we have f (R) = R. This means that differences with GR come from c 1 = 0. See Appendix A. Using Eqs. (15) and (17), we get the Ricci scalar, up to the leading order, in the form +· · · .

IV. MATCHING CONDITIONS
The solution given in Appendix A has a non-trivial Ricci scalar as shown by Eq. (21). Let us match it with the exterior solution. To this aim, we are going to match solution given in Appendix A with the one presented in [79] that we report here: where M is the total stellar mass and 4M G < c 2 r. We have to match the interior spacetime metric (9), using Eqs. (15) and (17), with the exterior Schwarzschild spacetime metric given by Eq. (23) at the boundary of the star, r = l. The continuity of the metric functions across the boundary r = l gives Using the above conditions, we get the constraints on the constants c 0 , c 2 . These constants take the form where C = 2GM c 2 l . Eq. (25) shows that the compactness parameter must be C < 0.5. This is also a consistency condition.

V. THE PHYSICAL VIABILITY OF SOLUTION
Now we are going to test if the interior solution, reported in Appendix A, can describe a realistic physical star. To this aim, we have to derive some necessary conditions which have to hold for any realistic star. They are the following.

A. Energy-momentum tensor
It is well-known that, for a realistic interior solution, the energy-density, the radial and transverse pressures must be positively defined. Moreover, all the components of the energy-momentum tensor should be finite at the center of the star then become decreasing in the direction of the surface of the stellar structure and the radial pressure should exceeds the tangential one 1 . The plots of the energy-momentum components, density, radial and transverse pressures are shown in Fig.1. The values of parameters used in the figure are:ǫ(y = 0) c1=0,c3=0.5 = 4.481840660, ǫ(y = 0) c1=0.5,c3=0.5 = 3.731840660,p(y = 0) c1=0,c3=0.5 =p t (y = 0) c1=0,c3=0.5 = 1.493946887,p(y = 0) c1=0.5,c3=0.5 = p t (y = 0) c1=0,c3=0.5 = 1.243946887. Fig. 1 shows that the density, radial and tangential pressures are decreasing towards the stellar surface. Moreover Fig.1 shows that the values of the components of the energy momentum at the center, in case c 1 = 0, are greater than those for c 1 = 0. We can also deduce, from Fig.1, that density, radial and tangential pressures converge to the surface of the star more rapidly in case c 1 = 0 than in the case of c 1 = 0.5. Finally, it is easy to see that, when c 1 = 0, one getsp =p t at the center however, as we reach the surface of the star, we can see thatp t ≥p.
In Fig.2, we show the behavior of the anisotropy parameter which is defined as ∆(y) =p t −p. Also Fig. 2 shows that the anisotropic force is positive. This means that a repulsive force, due to p t ≥ p, is present.

B. Causality
In order to show the behavior of sound velocities, we must calculate the gradient of energy-density, radial and transverse pressures that take the forms given in Appendix B. Let us analyze the gradient of energy-momentum tensor components reported in Appendix B. The asymptotic forms arê which give a negative gradient of the energy-momentum components provided that c 0 > c 1 ≥ c 2 > c 3 . This condition is satisfied according to the choice in footnote 1. The behavior of density, radial and tangential pressure gradients are shown in Fig. 3. To ensure that causality conditions, either for the radial and transverse sound speeds v r 2 and v ⊥ 2 , have values less than the speed of light, we use equations reported in given in Appendix B and get Eqs. (27) show that both the conditions 1 > v r 2 > 0 and 1 > v t 2 > 0 hold using the choice in footnote 1. The behaviour of the radial and tangential sound speeds are shown in Fig. 4 (a) and (b).
The appearance of non-vanishing total radial force with different signs in different regions of the fluid is called gravitational cracking when this radial force is directed inward in the inner part of the sphere for all values of the radial coordinate r between the center and some value beyond which the force reverses its direction [127]. In Ref. [128], it is stated that a simple requirement to avoid gravitational cracking is 0 < v r 2 − v t 2 < c 2 . In Fig. 4 (c), we show that the solution given in Appendix A is stable against cracking for c 1 = 0 and c 1 = 0.5.

C. Energy conditions
Energy conditions are considered an important test for non-vacuum solutions. To satisfy the dominant energy condition (DEC), we have to prove thatǫ −p > 0 &ǫ −p t > 0. As shown in Fig. 5 (a) and (b), DEC is satisfied for suitable choices of parameter c 1 . Furthermore, the weak energy condition (WEC),ǫ +p > 0ǫ +p t > 0 and the strong energy condition (SEC),ǫ −p − 2p t > 0 are satisfied as shown in Figs. 6 (a), (b) and (c).

D. Mass-radius relation
The compactification factor u(r) is defined as the ratio between the mass and radius. It has a key role to understand physical properties of compact objects. Starting from the solution given in Appendix A, we can define the gravitational mass as where erf (x) is the error function that is defined as The compactification factor u(r) is then defined as Using Eq. (28) into (30), one can get the explicit form of the compactification factor. The behavior of the gravitational mass and the compactification factor are plotted in Fig. 7. Fig. 7 (a) shows that the gravitational mass increases when the radial coordinate increases while Fig. 7 (b) shows that the compactification factor decreases for increasing radial coordinate.

E. Equation of state
The study of equation of state (EoS) for neutral compact stellar objects is reported in [129] . In that case, EoS is linear. In the present case, we will show that EoS cannot non-linear. To investigate this issue, we define the radial and transverse EoS as: where ω r and ω t are the radial and transverse EoS. Using Eqs. in Appendix A, one can calculate the analytical form of ω r and ω t whose behaviors are plotted in Fig.8. As Figs. 8 (a) and (b) show, EoS are non-linear because of the contribution of the higher-order curvature terms.

VI. STABILITY OF THE MODEL
In order to study the stability problem we are going to use two different procedures; the first one is the Tolman-Oppenheimer-Volkoff (TOV) equation and the second is the study of the adiabatic index.

A. The Tolman-Oppenheimer-Volkoff equation
To study the stability of solution in Appendix A, let us assume the hydrostatic equilibrium by using the TOV equation [130,131] as presented in [132]. It is where M g (y) is the gravitational mass calculated at the radius y. It has the following equation Inserting Eq. (33) into (32), we get where , F a = 2(pt−p) y and F h = − dp(y) dy are the gravitational, the anisotropic and the hydrostatic forces respectively. The behavior of TOV equation, for the model given in Appendix A, is shown in Fig.9. This figure shows that hydrostatic and anisotropic forces are positive and dominated on the gravitational force which is negative to keep the system in static equilibrium either when c 1 = 0 or c 1 = 0.5.

B. The adiabatic index
The adiabatic index γ is defined as [133,134] γ = ǫ +p p dp dǫ .
It allows to link the structure of a spherically symmetric static object and the EoS of the interior solution. It can be used to investigate the stability [135]. Any interior solution is stable if its adiabatic index is greater than 4/3 [136]. If γ = 4 3 , then the isotropic sphere is in neutral equilibrium. Following Chan et al. [137], the stability condition of a relativistic anisotropic sphere, γ > Γ, must be satisfied. Here it is The behaviour of the adiabatic index is shown in Fig.10 which ensures the stability condition of our solution.

VII. DISCUSSION AND CONCLUSIONS
In this study, we derived an anisotropic spherically symmetric solution in the framework of f (R) gravity. To this aim, we used a line element containing two unknown functions and applied it to the field equations of f (R) deriving four non-linear differential equations: three of them are the field equations and the fourth one is the trace of the field equations. This system involves six unknown functions, two of them are the metric potentials, one is the form of f (R) and the other are the density, and radial and tangential pressures. We solved these differential equations assuming the forms of the radial component of the metric and of the derivative of f (R). Substituting these quantities in the anisotropy parameter, ∆ =p t −p r , we get the time component of the metric. From these results, we obtained the form of f (R), the density, radial, and tangential pressures. It is interesting to mention that our solution is characterized by a constant c 1 which parameterize the deviation with respect to GR. Specifically, for c 1 = 0, we return to GR. (10 15 gr/cm 3 ) dp dρ (c 2 ) | c 1 =0 dp dρ (c 2 ) | c 1 =0  Table I. The boundary density ρ l and slope of the radial pressure reported for various observed objects. The density and the derivative of radial pressure are derived for the values c1 = 0 and c1 = 0.5 at the stellar boundary l assuming a nuclear density of the order 10 15 gr/cm 3 .
Due to the fact that the obtained interior solution has a non-trivial form for the Ricci scalar, we match it with a spherically-symmetric exterior solution having also a non-trivial form of Ricci scalar. From this matching, it is possible to relate the mass and the maximal radius of the compact object with the compactness parameter, that is C = 2GM c 2 l which must be C << 0.5 for the physical consistency of the model. The further step is to see if this solution represents a realistic star. To this aim, we discussed the necessary conditions required by any physically consistent compact stellar object and show that our model satisfies all these conditions. In particular, we considered two cases: one for c 1 = 0 and another for c 1 = 0.5.
Finally, we checked if the model satisfies the stability conditions. To this aim, we derived the TOV equation and showed that it is in static equilibrium either when c 1 = 0 or c 1 = 0.5. We also studied the adiabatic index to complete the picture of the stability and showed that the model satisfies the adiabatic index γ > 4/3.
These results can be applied to the observational data of different millisecond pulsars with white dwarf companions and pulsars presenting thermonuclear bursts as reported in Table I. For these objects, it is possible to calculate the density at the center and at the surface and then derive the EoS at the center and at the surface. It is interesting to see that, in the case of millisecond pulsars with a white dwarf companion, and, in the case of pulsars presenting thermonuclear bursts, our solution is compatible with observations. From the observational data in literature, as reported in the Table I, it is possible to derive the total mass and the radius of the given pulsar.
To conclude, we derived an interior solution in the framework of f (R) assuming physically motivated conditions to obtain the form of f (R). This solution has a non-trivial form of the Ricci scalar. The internal solution can be matched with the external one considering compacteness and anisotropy parameters. In a forthcoming paper, we will consider a similar situation assuming a charged interior solution.
The same as above for gradients.
whereǫ ′ = dǫ dy ,p ′ = dp dy andp ′ t = dpt dy . Appendix C: Derivation of Eq. (19) Eqs. (15) can be inserted in the last equation of system (11), the trace (6). We get a lengthy differential equation in y which we report here for the sake of completeness. It is where N (y), N 1 (y), N 2 (y), N 3 (y) are complicated functions of y obtained by summing up the coefficients of the same derivatives. The asymptotic form of this differential equation, up to O 1 y , gives Eq. (19).