The realistic models of relativistic stars in f(R) = R + alpha R^2 gravity

TL;DR

This paper investigates relativistic star models in f(R)=R+alpha R^2 gravity, revealing how positive alpha affects star structure, mass, and radius, with implications for gravitational mass and redshift, and discusses the problematic features for negative alpha.

Contribution

It provides a detailed analysis of neutron and quark stars in f(R)=R+alpha R^2 gravity without approximations, highlighting effects of alpha on star properties and stability.

Findings

Positive alpha causes non-zero scalar curvature at star surface.

Mass bounded by the star surface decreases with increasing alpha.

Negative alpha leads to oscillations and indefinite growth of gravitational sphere.

Abstract

In the context of f(R)=R + alpha R^2 gravity, we study the existence of neutron and quark stars with no intermediate approximations in the generalised system of Tolman-Oppenheimer-Volkov equations. Analysis shows that for positive alpha's the scalar curvature does not drop to zero at the star surface (as in General Relativity) but exponentially decreases with distance. Also the stellar mass bounded by star surface decreases when the value alpha increases. Nonetheless distant observers would observe a gravitational mass due to appearance of a so-called gravitational sphere around the star. The non-zero curvature contribution to the gravitational mass eventually is shown to compensate the stellar mass decrease for growing alpha's. We perform our analysis for several equations of state including purely hadronic configurations as well as hyperons and quark stars. In all cases, we assess that the relation between the parameter $\alpha$ and the gravitational mass weakly depend upon the chosen equation of state. Another interesting feature is the increase of the star radius in comparison to General Relativity for stars with masses close to maximal, whereas for intermediate masses around 1.4-1.6 solar masses, the radius of star depends upon alpha very weakly. Also the decrease in the mass bounded by star surface may cause the surface redshift to decrease in R^2-gravity when compared to Einsteinian predictions. This effect is shown to hardly depend upon the observed gravitational mass. Finally, for negative values of alpha our analysis shows that outside the star the scalar curvature has damped oscillations but the contribution of the gravitational sphere into the gravitational mass increases indefinitely with radial distance putting into question the very existence of such relativistic stars.

Tables

Table 1. Table 1: Parameters of compact stars for various equations of state in General Relativity (i.e., α = 0 𝛼 0 \alpha=0 ) and for several values of α 𝛼 \alpha in f ​ ( R ) = R + α ​ R 2 𝑓 𝑅 𝑅 𝛼 superscript 𝑅 2 f(R)=R+\alpha R^{2} gravity theories. The quantity Δ ​ M m ​ a ​ x Δ subscript 𝑀 𝑚 𝑎 𝑥 \Delta M_{max} holds for the maximal contribution of the gravitational sphere outside the star in the total value of the gravitational mass (having included the stellar mass). On the fourth column, the value for the corresponding gravitational mass is given in brackets. In the last column we show the maximum difference in the star radius as compared with the General Relativity counterpart. This maximum difference takes place for a star with a mass equal to the maximal mass in General Relativity for a given equation of state.

EoS	$\alpha$,	$M_{max}$,	$M_{s,max}$,	$\Delta M_{max},$	$\Delta r_{max}$,
	${10}^{10}$ cm2	$M_{\odot }$	$M_{\odot }$	$M_{\odot }$	km
	0	2.23	2.23	0	0
	0.5	2.22	2.18	0.045 (1.94)	0.05
AP4	5	2.24	2.12	0.182 (1.82)	0.45
	10	2.26	2.12	0.232 (1.88)	0.63
	20	2.29	2.11	0.296 (1.91)	0.84
	0	2.05	2.05	0	0
	0.5	2.05	2.02	0.042 (1.92)	0.08
SLY	5	2.10	1.97	0.176 (1.74)	0.58
	10	2.13	1.95	0.237(1.84)	0.80
	20	2.16	1.94	0.294 (1.86)	0.96
	0	2.49	2.49	0	0
	0.5	2.49	2.45	0.045 (2.29)	0.045
MPA1	5	2.51	2.38	0.188 (1.99)	0.55
	10	2.54	2.37	0.256 (2.01)	0.78
	20	2.57	2.34	0.320 (2.02)	1.0
	0	1.93	1.93	0	0
	0.5	1.93	1.90	0.035(1.88)	0.06
GM1nph	5	1.96	1.80	0.169 (1.85)	0.32
	10	1.98	1.76	0.237 (1.89)	0.47
	20	2.01	1.73	0.296 (1.82)	0.60
	0	1.88	1.88	0	0
	0.5	1.89	1.82	0.072 (1.72)	0.075
QEoS	5	1.92	1.73	0.209 (1.71)	0.31
	10	1.94	1.70	0.264 (1.71)	0.41
	20	1.96	1.67	0.313 (1.82)	0.51

Table 2. Table 2: Calculation of B ​ ( r i ) ≡ B i 𝐵 subscript 𝑟 𝑖 subscript 𝐵 𝑖 B(r_{i})\equiv B_{i} at some r i subscript 𝑟 𝑖 r_{i} for B ​ ( 0 ) = 0.1 𝐵 0 0.1 B(0)=0.1 and MPA1 equation of state. The fourth column shows m ​ ( r i ) 𝑚 subscript 𝑟 𝑖 m(r_{i}) . Radial distances have been considered in the interval r i − 5 ⋅ 10 3 ​ km < r < r i + 5 ⋅ 10 3 subscript 𝑟 𝑖 ⋅ 5 superscript 10 3 km 𝑟 subscript 𝑟 𝑖 ⋅ 5 superscript 10 3 r_{i}-5\cdot 10^{3}\,{\rm km}<r<r_{i}+5\cdot 10^{3} km. As observed, masses m ​ ( r i ) 𝑚 subscript 𝑟 𝑖 m(r_{i}) grow with r i subscript 𝑟 𝑖 r_{i} . As expected, numerical calculations in General Relativity give the well-known result, i.e., the function ( B ∞ − B ​ ( r ) ) ​ r / 2 subscript 𝐵 𝐵 𝑟 𝑟 2 (B_{\infty}-B(r))r/2 remains constant being the gravitational mass is 2.48 ​ M ⊙ 2.48 subscript 𝑀 direct-product 2.48M_{\odot} for this case.

$r_i$, ${10}^4$ km	$B_i$	${\overline{B}}_i\left(0\right)$	$m\left(r_i\right)$	$B_{i,\mathrm{GR}}$
1	1.487088	0.067178	2.64	1.482622
2	1.487687	0.067151	2.83	1.483165
3	1.487900	0.067141	3.02	1.483346
4	1.488013	0.067136	3.19	1.483437
5	1.488085	0.067133	3.37	1.483491
6	1.488135	0.067131	3.55	1.483528
7	1.488173	0.067129	3.72	1.483553
8	1.488203	0.067128	3.90	1.483573
9	1.488227	0.067127	4.07	1.483588
10	1.488248	0.067126	4.25	1.483600
12	1.488280	0.067124	4.60	1.483618
14	1.488305	0.067123	4.95	1.483631

Table 3. Table 3: Results for B ∞ subscript 𝐵 B_{\infty} in the interval r i < r < r i + 10 4 subscript 𝑟 𝑖 𝑟 subscript 𝑟 𝑖 superscript 10 4 r_{i}<r<r_{i}+10^{4} km for B ​ ( 0 ) = 0.1 𝐵 0 0.1 B(0)=0.1 for the MPA1 equation of state and corresponding of B ¯ ​ ( 0 ) ¯ 𝐵 0 \bar{B}(0) in same interval. We give the maximum value of m ​ ( r ) 𝑚 𝑟 {m}(r) ( m m ​ a ​ x subscript 𝑚 𝑚 𝑎 𝑥 {m}_{max} ) for corresponding the B ¯ ¯ 𝐵 \bar{B} and B ∞ subscript 𝐵 B_{\infty} .

$r_i$, ${10}^4$ km	$B_{\infty }$	$\overline{B}\left(0\right)$	$m_{max}$
1	1.488284	0.067192	2.74
2	1.488325	0.067190	2.92
3	1.488352	0.067189	3.05
4	1.488372	0.067188	3.27
5	1.488387	0.067187	3.44
6	1.488400	0.067186	3.62
	1.488271	0.067919	2.71
	1.488362	0.067188	3.18

Equations

$$S=\frac{1}{16\pi }\int {\mathrm{d}}^{4}x\sqrt{-g}f(R),$$

$$f_R{G}_{\mu \nu }+({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }{\nabla }^{\alpha }{\nabla }_{\alpha })f_R-\frac{1}{2}(f(R)-R{f}_R)g_{\mu \nu }$$

$$=8\pi T_{\mu \nu },$$

$$T_{\mu \nu }=(\rho +p)u_{\nu }{u}_{\mu }-p{g}_{\mu \nu }.$$

$$\mathrm{d}{s}^2=-B(r)\mathrm{d}{t}^2+A(r)\mathrm{d}{r}^2+r^2(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2),$$

$$A^{\prime }$$

$${\psi }^{\prime }$$

$$p^{\prime }$$

$$R^{\prime \prime }$$

$$A(0)=1,\psi (0)=0,R^{\prime }(0)=0.$$

$$R\to 0,A\to 1,r\to \infty .$$

$$B(0)\to \frac{B(0)}{B_{\infty }}.$$

$$A(r)={\left(1-\frac{2m(r)}{r}\right)}^{-1}$$

$$m(r)=\frac{r}{2}\left(1-A(r{)}^{-1}\right).$$

$$B(r)=1-\frac{2m(r)}{r},$$

$$M\equiv \underset{r\to \infty }{\lim }r(1-B(r))/2.$$

$$B(r)=B_{\infty }\left(1-\frac{2m(r)}{r}\right),$$

$$M=\frac{1}{B_{\infty }}\underset{r\to \infty }{\lim }\frac{r(B_{\infty }-B(r))}{2},$$

$$\zeta$$

$$\zeta ={\zeta }_{low}K(a_1(\xi -c_{11}))+{\zeta }_{high}K(a_2(c_{12}-\xi )),$$

$${\zeta }_{low}$$

$${\zeta }_{high}$$

$$p=b(\rho -4B),$$

$$z_s=(1-2M_s/r_s{)}^{-1/2}-1,$$

$$B(r)=B_{\infty }\left(1-\frac{2m}{r}\right)$$

Full text

modified gravity; neutron stars; quark stars.

I Introduction

The discovery of the accelerated expansion of the universe Riess ; Perlmutter has recently led to an intensive study of cosmological applications of extended (also dubbed modified) theories of gravity which aim to overcome the limitations of the cosmological Concordance $\Lambda$CDM model, made of General Relativity (GR), a cosmological constant $\Lambda$ and Cold (non-baryonic) Dark Matter. Many of these extensions aim to avoid the usually invoked need of dark components which is required, for instance both in the $\Lambda$CDM model, where dark energy is nothing but a cosmological constant, and in phantom/quintessence-like models where a scalar field is needed, among others. Throughout the years, a plethora of modified gravity theories have been shown to be able to preserve the positive results of Einstein theory and obtain the required cosmological evolution with the most relevant cosmological epochs Capozziello2 ; Odintsov1 ; Turner . Among those attempts, one of the most successful ones is the so-called $f(R)$ theories where the gravitational Lagrangian includes powers of the Ricci scalar $R$ encompassed in an arbitrary function of $R$ (c.f. Odintsov-3 ; Capozziello_book ; Capozziello4 for extensive reviews on the subject).

Since in principle the number of viable modified gravities, and classes of models therein, is sufficiently large, the comparison between theoretical predictions and different catalogues of cosmological dataset, extracted from both the cosmic background and relativistic perturbations (mainly scalar and tensor) evolutions, remains mandatory. Usually, such comparisons suffer from the so-called degeneracy problem, meaning that several competitive gravitational theories are able to explain the same phenomena with roughly speaking the same statistical precision. Other approaches recently explored have considered model-independent (cosmographic) techniques in order to reconstruct the underlying theory capable of causing the observed cosmological evolution once the Copernican principle is assumed (c.f. Cruz-Cosmography ; Bamba:2012 and references therein).

Thus in order to further constrain the number of experimentally viable theories beyond Einsteinian gravity and obtain more rigid constraints on the parameters space of every theory, one possibility consists of studying the consequences of these proposals at the astrophysical level. Understandably, these high-energy configurations are the ideal arena to study the behaviour of extended theories of gravity in strong gravitational field regimes and may complement the picture as provided by the low-curvature late-time cosmological evolution. First of all, it is interesting to investigate the possibility of the existence and features of black holes delaCruzDombriz:2012xy and relativistic stars, such as white dwarfs and neutron stars, the latter being the objective of the research in this paper. Thus for viable classes of $f(R)$ Lagrangians, one needs to obtain the relations between the paradigmatic parameters of relativistic stars such as masses, radii, momenta of inertia, quadrupole momenta and Love number, among others. Unfortunately, at the present, fully precise astrophysical measurements cannot be provided for the aforementioned quantities, the exception being the mass measurements for neutron stars.

As extensively studied, once GR is assumed as the correct theory of gravity, there are certainly upper limits for realistic neutron star masses, being the theoretical mass limit itself increased along the years, from $\mathcal{O}(0.6M_{\odot})$ that Tolman-Oppenheimer-Volkoff (TOV) Oppenheimer:1939ne found for a free neutron gas equation of state (EoS), to $\mathcal{O}(2.2M_{\odot})$ limits when stiffer chiral interactions are allowed Dobado:2011gd . Notwithstanding, the discovery of a neutron star with mass $1.974M_{\odot}$ Demorest:2010bx by the Shapiro delay method, confirmed some years after by a precise sighting Antoniadis:2013pzd through binary system measurements, have indeed put into question the validity of GR in such high-energy gravitational environment and have been able to rule out exotic EoS proposals. Further observations have certainly confirmed the existence of neutron stars heavier than the GR naive prediction others_masses as well as double neutron stars and pulsar systems also violating the standard GR predictions (c.f. Dexheimer:2007mt for an exhaustive list of the latter) . In that realm, some attempts have also tried to reconcile these results with GR, assuming for instance strong magnetic fields Astashenok-1 or electric charges Jing:2015ota inside the star.

Neutron stars in the context of $f(R)$ theories were primarily studied in Kobayashi-Maeda ; Upadhye-Hu and then in abundant bibliography Feng:2017hje ; Pannia:2016qbj ; Wojnar:2016bzk ; Arapoglu:2016ozr ; Katsuragawa:2015lbl ; Fiziev:2015xpa ; Hendi:2015pua ; Momeni:2015vwa ; Zubair:2016kov ; Bakirova:2016ffk ; Resco:2016upv ; Moraes:2015uxq ; Sharif:2015jaa ; Sotani:2017pfj for both scalar-tensor $f(R)$ theories, as well as some closely-related theories. Also the possibility of emergence of new types of astrophysical objects in the context of $f(R)$ theories, such as stable stars with large magnetic fields, supermassive stars, among others (see for instance Laurentis ; Laurentis2 for some proposals).

Perturbative solutions of TOV equations for simple $f(R)$ models were investigated in Arapoglu ; Alavirad ; Astashenok-1 . Obviously, limitations of a perturbative approach lie mainly in the impossibility of comparing the (unknown) exact solution with the determined perturbative solution. Also one would expect that in strong gravity regimes, non-GR gravitational effects to be dominant and consequently, the GR limit cannot be simply considered as the leading contribution in the solution. Thus, for a careful description one needs to solve exact higher-order differential equations as we shall do herein.

Another approach in the literature has consisted of resorting to the equivalence between $f(R)$ theories and Brans-Dicke scalar-tensor theory in order to study realistic models of neutron and strange stars Yazadjiev2014 ; Doneva2014a ; Henttunen . Within this approach the description in terms of $f(R)=R+\alpha R^{2}$ gravity was proposed in Astashenok2015 and illustrated for quark stars. There, it was shown that for distant observers the gravitational mass of star increases with increasing $\alpha$ ($\alpha>0$), though the interpretation of this fact depends upon the frame in which calculations are performed.

Considering $f(R)$ theory directly, one can see that the matching conditions at the edge of the star, do not impose the scalar curvature to vanish Clifton:2012ry , so $R\neq 0$ is in principle expected outside the star. As a consequence, the spacetime region outside the star also contributes to the total gravitational mass as perceived by a distant observer (c.f. Resco:2016upv for a thorough discussion on this issue). By analogy with scalar-tensor theory of gravity in which that region outside the star is referred to as the dilaton sphere (or disphere) Fiziev , in the following we shall denote this domain as gravitational sphere (or gravisphere).

In this paper we shall consider in detail the possibility of existence and main features of neutron stars and quark (also dubbed strange) stars in the context of $f(R)=R+\alpha R^{2}$ gravity. The TOV set of equations for $f(R)$ theories is presented in the Section II. Then we shall briefly discuss the various equations of state (EoS) for nuclear matter in Sec. III. As mentioned above, for modified gravity theories, in principle one needs to investigate the solution of field equations outside the star since it may happen that the exterior spacetime is not Ricci flat and Schwarzschild is not the only permitted solution. As we shall show, the obtained exterior solution crucially depends upon the chosen EoS because the latter fixes the scalar curvature value at the star surface. The next two following section will be devoted to obtaining solutions for various EoS. In Sec. IV we shall show how for positive $\alpha$ one needs to choose the initial condition for $R(0)$ to satisfy the asymptotic (Schwarzschild) flatness requirement. The scalar curvature outside the star thus sharply decreases. It will be shown that the stellar mass, for a given central density (i.e., density at the centre of the star) is smaller than the GR counterpart. Nonetheless the emergence of the gravitational sphere around the star, provides an extra contribution to the gravitational mass as measured by a distant observer. As a result the gravitational mass overcomes the GR counterpart for central density values higher than a given critical central density. The value of this critical central density will be shown to be determined by both $\alpha$ and the chosen EoS. This effect takes place for various EoS, including the exotic case of quark stars. It is interesting to note that the values of radii for intermediate masses $1.4M_{\odot}<M<1.6M_{\odot}$ very weakly depend upon $\alpha$. Finally, in Sec. V we shall show how for $\alpha<0$ cases, although outside the star the scalar curvature possess damped oscillations, the gravitational mass contribution when taking into account the increases with the radial distance. We conclude the paper by showing our main conclusions in Sec. VI. Figures for all the studied cases are presented at the end of the paper and are intended to help the reader to understand the explained trends in the bulk of the paper.

II Equations for compact stars in $f(R)$ gravity

Let us consider the gravitational action for $f(R)$ gravity where $R$ denotes the scalar curvature

[TABLE]

where $G=c=1$ units will be used throughout the paper. The corresponding field equations can be obtained by varying action (1) with respect to the metric $g_{\mu\nu}$, yielding

[TABLE]

where $f_{R}\equiv{\rm d}f(R)/{\rm d}R$, $G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}\,R\,g_{\mu\nu}$ is Einstein tensor and $T_{\mu\nu}$ is energy-momentum tensor, which for a perfect fluid in a comoving frame can be written as

[TABLE]

As easily deduced taking the trace of Eq. (II), in GR ($f(R)\equiv R$) the scalar curvature is univocally defined by the energy-momentum content through an algebraic equation. On the contrary, in the context of $f(R)$ theories of gravity in the metric formalism, given the higher order of the field equations, $R$ may be thought of as an independent dynamical function.

For a spherically symmetric metric of the form

[TABLE]

the corresponding field equations as obtained from Eq. (II) for functions $A$ and $\psi=B^{\prime}/B$ (prime denotes derivative with respect to $r$) become

[TABLE]

and the trace of Eq. (II) gives the equation for the scalar curvature:

[TABLE]

where $f_{2R}={\rm d}^{2}f(R)/{\rm d}R^{2}$, $f_{3R}={\rm d}^{3}f(R)/{\rm d}R^{3}$. Finally one needs to specify the $f(R)$ model under consideration (in our case $f(R)=R+\,\alpha R^{2}$, with $\alpha$ is constant parameter) and add the EoS for completeness of the system, i.e., $p=f(\rho)$. Several EoS will be considered and explained in Sec. III.

Initial and boundary conditions: The following conditions at the centre of star ($r=0$) should be imposed in order to guarantee regularity and finiteness both of density and pressure everywhere. This provides

[TABLE]

Therefore there is only one free initial condition for the system of equations Eqs. (5)-(8), namely the value of the scalar curvature $R(0)$. At the surface of star, both pressure and density (for most EoS) cancel, i.e., $\rho_{s}=p_{s}=0$. Thus, from the star radius onwards, one needs to integrate our system of Eqs. (5)-(8) with $\rho=p=0$. The boundary conditions at spatial infinity can be determined from the reasonable requirement that the star exterior solution should be matched to the Schwarzschild one at infinity. Therefore

[TABLE]

From the equations above, one can see that solutions for $A(r)$, $\rho(r)$, $p(r)$ and $\psi(r)$ do not depend on the chosen value $B(0)$. However, in order to achieve, as desirable, Schwarzschild solution at the asymptotic spatial infinity, one must require $B(r)\rightarrow 1$ for $r\rightarrow\infty$. If $B(r)\rightarrow B_{\infty}\neq 1$ at $r\rightarrow\infty$, one can simply rescale the initial condition for $B(0)$ as

[TABLE]

On the other hand, on usually assumes $A(r)$ to be of the form

[TABLE]

and therefore the mass function $m(r)$ as a function of radial coordinate can be expressed as :

[TABLE]

Mass extraction for negative $\alpha$: For the performed analysis with $\alpha<0$, we have seen that $A(r)$ oscillates around unity at sufficiently large distances and therefore one cannot extract the gravitational mass from this parametrization correctly. For $r\rightarrow\infty$ we have $B\sim{\rm e}^{2\phi}$ where $\phi$ can be interpreted as the Newtonian potential for the gravitational field at distance $r$, in other words, $\phi(r)=-m(r)/r$. Therefore $B(r)\approx 1-\frac{2m(r)}{r}$ at sufficiently large distance $r$. Provided the gravitational mass has a well-defined limit one can assume that $m(r)\approx M=const$ from some $r=r_{c}$ onwards ($r_{c}$ being a sufficiently large distance) and thus, the gravitational mass $M$ can be extracted from $B(r)$. Consequently for an initial value of $B(0)$ so that $B(r)\rightarrow 1$ at $r\rightarrow\infty$ we have

[TABLE]

and therefore the gravitational mass $M$ can be calculated as the limit

[TABLE]

Mass extraction for positive $\alpha$: As we shall see in Sec. IV, in the event $\alpha>0$, the limit $m(r)\rightarrow M$ at spatial infinity is well defined. In the event that the value of $B(0)$ is arbitrarily chosen, one first needs to obtain the asymptotic value $B_{\infty}=\lim_{r\rightarrow\infty}B(r)$, so that function $B(r)$ becomes

[TABLE]

so a simple relation for gravitational mass

[TABLE]

is obtained.

We need carefully investigate the possibility of appearance of contributions with $r^{\gamma}$ ($\gamma>0$) in $B(r)$. In this case we have no finite limit for mass. The behaviour of the resulting solution for the metric functions $A$ and $B$ in the outer region depends values $A(r)$, $B(r)$ and its derivatives on star surface ($r=r_{s}$). For the latter, specific values depend upon the chosen EoS.

III Realistic equations of state

For our investigation we have considered a representative set of four realistic EoS.

Neutron stars: Below we provide the main features and assumptions in the EoS under study herein.

• $1)$

the well-known SLy4 EoS SLy ; SLy-4 is obtained from many body calculations with simple two-nucleon potential. For completeness, we have also considered the AP4 EoS APR which hosts a three-nucleon potential and Argonne 18 potential with UIX potential;

• $2)$

as an example of EoS based upon the relativistic mean-field (RMF) calculations we take the generalisation of a model firstly considered by Glendenning and Moszkowski (GM) GM . This generalisation, dubbed GM1nph, has been proposed for cold neutron star matter in $\beta$-equilibrium containing the baryon octet and electrons GM-1 ;

• $3)$

relatively stiff MPA1 EoS MPA obtained from relativistic Dirac-Brueckner-Hartree-Fock formalism. This EoS accounts for the energetic contributions originated by the exchange between $\pi$ and $\rho$ mesons.

On Fig. 1 we plot dependence of pressure from energy density for these four EoS. One can see that for intermediate densities these EoS seem indistinguishable. For GM1 EoS there is a softening effect at high densities due to appearance of hyperons (mainly $\Sigma^{-}$, $\Lambda^{0}$ and $\Sigma^{0}$). In our calculations we have used analytical representations of these EoS. Thus, the SLy EoS can be represented as

[TABLE]

where $\zeta\equiv\log(P/{\rm dyn}\,{\rm cm}^{-2})$, $\xi\equiv\log(\rho/\mbox{g}\,\mbox{cm}^{-3})$, and $K(x)\equiv 1/({\rm e}^{x}+1)$. The coefficients $a_{i}$ can be found in Camenzind . For MPA1, AP4 and GM1nph EoS take the more complicated dependence shown in Eksi , namely

[TABLE]

with

[TABLE]

where $a_{i}$ and $c_{i}$ are constant coefficients provided in Eksi .

Quark stars: Apart from the EoS described above, we have also considered the case of quark (or so called strange) stars Itoh ; Witten . It is assumed that quark stars consist of $u$, $d$ and $s$ quarks and electrons. The deconfined quarks can form a colour superconductor. This leads to a softer EoS with possible observable effects on the minimum allowed mass, radii, cooling behaviour and other observables. One of the simplest EoS for quark matter is provided by the so-called MIT bag model Jaffe ; Simonov as follows:

[TABLE]

where $B$ is the bag constant. The value of parameter $b$ depends on both the chosen mass for the strange quark $m_{s}$ and QCD coupling constant and usually varies from $b=1/3$ ($m_{s}=0$) to $b=0.28$ ($m_{s}=250$ MeV). The value of $B$ lies in interval $0.98<B<1.52$ in units of $B_{0}=60$ MeV/fm3 Sterg . In the following, we shall consider the case $B=1$ and $b=0.31$.

IV Results for $\alpha>0$

Following the procedure described in Sec. II we have calculated the mass-radius and mass-central density relations for some values of parameter $\alpha$ for all EoS explained in Sec. III above. Results are depicted in Figs. 2, 4, 6, 8 and 10. More specifically, for a given initial pressure, we have integrated the system of Eqns. (5)-(8) using the EoS described in Sec. III. As soon as $p(r=r_{s})=0$, we label that radial coordinate value as the star radius $r_{s}$. From this value onwards, obviously we are in vacuum and the aforementioned system of equations is integrated with zero density and zero pressure. The determination of the asymptotic mass $M$ follows the procedure described in III, so a pair $\{M_{,}\,r_{s}\}$ is found, which gives rise to results in Figs. 2, 4, 6, 8 and 10. One needs to be aware that as in the GR case, not every initial pressure value is able to host a stable star, so in this sense every pair $\{M_{,}\,r_{s}\}$ corresponds to a different initial pressure.

As has been explained in the Introduction, in the context of $f(R)$ extended theories of gravity, one needs to distinguish between stellar mass bounded by star surface ($M_{s}=m(r_{s})$) and gravitational mass as measured by a distant observer ($M=m(r\rightarrow\infty)$). Obviously, in GR these two values coincide, but in $R^{2}$-gravity the gravitational sphere around the star emerges and contributes to the gravitational mass defined as $M(r_{\infty})$. We have observed that for the same central density, $R^{2}$ stellar masses $M_{s}$ are smaller than their GR counterparts. This decrease is then partially compensated by energetic contribution of the gravitational sphere when the gravitational mass $M$ is computed. Consequently, the resulting gravitational mass-radius relation ends up to differ from the GR prediction. The effective radius of the gravitational sphere can be determined as the radius at which the metric solution is sufficiently close to the Schwarzschild solution (for example the radius at which the scalar curvature is of the order of magnitude $R\sim 10^{-5}$ in units of $r_{g}^{-2}\equiv c^{4}/G^{2}M^{2}_{\odot}$). We have observed that the effective radius turns out to be proportional to $\alpha^{1/2}$.

On the other hand, Figs. 3, 5, 7, 9 and 11 represent the Ricci scalar curvature evolution outside the star. For $R^{2}$ gravity, the curvature turns out to be positive on star surface and decreases outside the star so that $R\rightarrow 0$ at $r\rightarrow\infty$. In the aforementioned figures, one can see the following features: Firstly the maximal gravitational mass $M_{max}$ (defined as ${\rm d}M/{\rm d}r_{s}=0$) of star for each EoS slowly increases with $\alpha$. The value $\Delta M=M_{max}(\alpha)-M^{\rm GR}_{max}\sim\sqrt{\alpha}$. At first glance one cannot discriminate between $R^{2}$-gravity (for the considered values of $\alpha$) and GR because such mass differences lie within typical errors of neutron star mass measurements from pulsar timing. Indeed, from Table I in Lattimer2 one can see that 1-$\sigma$-uncertainties for neutron star masses for X/ray-optical binaries and neutron star-white dwarf binaries are $\sim 0.1-0.2M_{\odot}$ with few exceptions. Only for neutron stars with small masses precise measurements of masses with errors $<0.01M_{\odot}$ (neutron star-neutron star binaries) are available. In conclusion, we have no well-established data on radii of these low-mass neutron stars.

Another in principle observable astrophysical quantity for these gravitational configurations is the so-called the surface gravitational redshift $z_{s}$ Cottam , which is defined as

[TABLE]

and therefore $z_{S}$ is determined by stellar mass $M_{s}$. According to our results, in the context of $R^{2}$-gravity for $\alpha>0$, we conclude that $z_{s}<z_{s,\,{\rm GR}}$. Moreover, the expected absolute deviation of the surface redshift from the GR value for a given mass, i.e., $\Delta z=z_{s,\,{\rm GR}}-z_{s}$ turns out to be weakly dependent on the gravitational mass for typical masses of neutron stars. In principle this discrepancy could help to discriminate between GR and $R^{2}$-gravity, and by extension to establish the validity of different classes of $f(R)$ models and other extended theories of gravity. Indeed, in the case of $R^{2}$-gravity we have shown that $M_{s}<M$. Thus, the measure of the gravitational redshift would allow us to label the star by its stellar mass. On the other hand, as stated above, by means of binary systems one can measure the gravitational mass $M$. Accordingly, the discrepancy between these two quantities could in principle witnesses in favor of $R^{2}$-gravity, or eventually other extended theories of gravity for which $M\neq M_{s}$. Unfortunately, at the present moment, one have no well-defined data about neutron-star gravitational redshifts nor radii.

Another interesting feature is the increase of star radius for a given gravitational mass close to maximal mass $M_{max}$ when compared to GR counterparts (see Table 1). In principle this fact may open the possibility to obtain an upper limit on parameter $\alpha$ by using observational data for neutron stars radii. Unfortunately, at the present time we have no well-established such data.

From the gravitational mass-radius relation in Figs. 2, 4, 6, 8 and 10 we see that for purely hadron EoS, the radii of neutron stars with typical masses $\sim 1.4-1.5M_{\odot}$ hardly depends on the value of $\alpha$. This effect is more noticeable for SLy and MPA1 EoS. As an illustrative example, let us mention that for GM1 EoS with inclusion of hyperons we see in Fig 8 that for several values of $\alpha$ the mass-radius curves intersect in the vicinity of $1.6M_{\odot}$. Therefore the data corresponding to the mass-radius relation in the vicinity of either large or small masses can help to discriminate between GR and the gravitational modification under study here. Similar statements can be made for the other two neutron stars EoS.

Finally, for quark stars we have the following scenario: the density on star surface of strange stars is nonzero, namely $\rho_{s}=4B$. Therefore, the scalar curvature evaluated at the star’s surface, as predicted by GR, would be $R_{s}=32\pi B$ whereas outside the star, the scalar curvature drops to zero suddenly. On the contrary, in $R^{2}$-gravity one can construct solution with smooth decreasing of scalar curvature as can be seen in Fig. 11.

V Results for $\alpha<0$

In this Section we have considered the case $\alpha=-0.05$ km2 as an illustrative example in order to get a better understanding of the behaviour of the system of Eqs. (5)-(8) when applied to $R^{2}$-gravity models with negative parameter $\alpha$. The considered central densities for the studied examples below correspond to stars possessing a GR maximal mass for each of the chosen EoS. The system (5)-(8) is thus integrated up to a large enough distance ($r\approx 3\cdot 10^{5}$ km).

We have analysed the $B(r)$ dependence with the radial coordinate and studied the behaviour of this metric function as $r\rightarrow\infty$ so by the use the Eq. (17) we have tried to extract the corresponding gravitational mass. Another (but completely equivalent) way we have followed consisted of trying to find such an initial value $B(0)$ for which the function ${B}(r)\rightarrow 1$ as $r\rightarrow\infty$. We shall designate this value as $\bar{B}(0)$. For the latter method, one can consider a sequence of $\bar{B}_{i}(0)$ values such that ${B}(r_{i})=1$. Then extrapolating this sequence for $r_{i}\rightarrow\infty$ gives the required $\bar{B}(0)$. The gravitational mass can then be calculated according to Eq. (15).

We have also performed the following calculations: after having assumed that $m(r)$ varies slowly at large distances within some interval $\Delta r$ (we take $\Delta r\approx 10^{4}$ km), we have approximated $B(r)$ by a functional dependence

[TABLE]

in this $\Delta r$ interval. Thus, as a rough estimation one can identify the corresponding mass for given $r$ with $m$. For GR this scheme gives correct value for gravitational mass, whereas for $R^{2}$ gravity our method shows that that the found mass grows with the radial distance.

In order to illustrate the reasoning, we have chosen $B(0)=0.1$. Table 2 below summarises the found results for the MPA1 EoS case. We have also calculated the values $\bar{{B}}_{i}(0)$ at which for $r_{i}$ the condition $1-{B}(r_{i})=\beta$, is accomplished with $\beta=10^{-4}$. In GR $B(r)$ exactly scales like $1/r$ and therefore one can extract the gravitational mass $M$ using dependence $B(r)$ for an arbitrary $B(0)$ value, by simply approximating this dependence to the well-known function $B_{\infty}(1-2M/r)$. Thus, for viable theories beyond GR the function $m(r)$ is expected to tend to constant value asymptotically, i.e., there is a plateau for $m(r)$.

In the context of $f(R)=R+\alpha R^{2}$ ($\alpha<0$) models, one can assume that the aforementioned plateau for the function $m(r)$ will be found at sufficiently large distances. Nonetheless, our method shows that the approximation of $\bar{B}(0)$ (or equivalently $B_{\infty}$ for an arbitrary $B(0)$) strongly depends on the considered interval of distances.

In particular we have considered the approximation of $B(r)$ by a function of the form $B_{\infty}(1-2m/r)$ in intervals $r_{i}<r<r_{i}+10^{4}$ km for some $r_{i}$ and results have been shown in Table 3). In this case the dependence $(1-B(r)/B_{\infty})r/2$ has a maximum near $r_{i}$. We have also considered the interval $10^{3}\,{\rm km}<r<5\cdot 10^{4}$ km for $B(r)$. In this case we have for $B_{\infty}=1.488271$ and the maximal value for ${m}(r)$ is $2.71M_{\odot}$. If for instance the consider distance interval is now $2\cdot 10^{4}\,{\rm km}<r<7\cdot 10^{4}$ km one obtains $B_{\infty}=1.488362$ and maximal of $m(r)$ is $3.15M_{\odot}$. Thus our method leads us to conclude that $B_{\infty}-B(r)\nsim r^{-1}$ even for large distances. Consequently, the main point of our analysis is that the function $B_{\infty}-B(r)$ (or $1-B(r)$ for calibrated initial condition for $B(0)$) tends to [math] more slowly than $r^{-1}$ at $r\rightarrow\infty$. Therefore we are led to conclude that we have no well-defined limit of mass since the function $m(r)$ seems to grow with radial distance. In fact the value of $\bar{B}(0)$ (or $B_{\infty}$) decreases (increases) with growing of maximal limit of radial interval under consideration.

Of course one can assume that the results above only take place for some EoS only. In order to see if this is true, we have considered all the other three-nucleon EoS from Sec. III, namely AP4, SLy and GM1. For those, we have obtained that the $m(r)$ growth rate with the radial coordinate slightly depends on the specific EoS, but there is no qualitative difference with respect to the trend fully explained above.

VI Conclusions

In this paper we have investigated the main features and existence of realistic models of compact neutron and quarks stars in one paradigmatic extension of General Relativity. Namely we have studied these static and spherically symmetric configurations for a class of fourth-order $f(R)$ extended theories of gravity, whose gravitational Lagrangian adopts the form $f(R)=R+\alpha R^{2}$. As explained in the bulk of the article, the resolution of the generalised system of Tolman-Oppenheimer-Volkoff equations for positive $\alpha$ has indeed provided the existence of solutions, whereas for negative $\alpha$ it has been impossible to recover a weak-field (Newtonian) limit far away from the star, concluding the impossibility of stable configurations in such a region of the parameter space.

Note that our analysis has not considered any perturbative approximation whatsoever when dealing with the system of equations resolution nor a shooting method has been required owing to the chosen variables. Thus, we have found that for the same central density the stellar mass as bounded by the star surface decreases when the parameter $\alpha$ increases. This result can be understood as a proof that in extended theories of gravity, a more intense gravitational interaction leads to stable stars which for the same mass in Einsteinian gravity would either not exist or be condemned to be black holes. More specifically, for the parameters space of quadratic $f(R)$ models which are in agreement with theoretical constraints and cosmological data, the $f(R)$ predictions for the neutron stars maximal mass ranged from 2-2.6 $M_{\odot}$ (see Table 1) and therefore agreed with observed stars as those found in others_masses . Similar results were found for the so-called Hu-Sawicki $f(R)$ model Resco:2016upv .

Also, we have seen that, unlike previous works Capo which naively defined the mass as the matter integral over the star volume, we have shown how a more careful matching of interior and exterior regions causes the total gravitational mass to be unbounded for sufficiently large values of $\alpha$ (with specific values depending on the equation of state under consideration). Consequently this more general consideration of conditions to be imposed at the edge of the star is twofold.

First, these reasonable matching conditions imply that the scalar curvature is not identically null on the edge of the star, even so (for $\alpha>0$) it is exponentially suppressed and goes to zero and accordingly asymptotic flatness is eventually recovered. This can be thought of as being equivalent to the existence of an additional energetic content around the star which might be assigned to the additional extra scalar mode present in $f(R)$ theories. It is precisely this mode which can help to prevent the gravitational collapse and increase the gravitational mass as observed by a distant observer. Thus in the event of merging objects of this kind, the available energy would be higher than General Relativity counterparts and the emission of gravitational waves from such a configuration might be feasible and detectable by present and future interferometers Resco:2016upv ; Damour:1996ke ; Berti:2015itd ; Koyama:2011xz ; Sibandze:2016agp ; BeltranJimenez:2017doy .

Secondly, we have indeed corroborated the existence of exterior static and spherically symmetric physical solutions differing from the usual Schwarzschild solution (non Ricci flat) provided $f(R)$ theories are assumed to represent the underlying gravitational theory in high-curvature environments. This fact constitutes a straightforward violation of the Jebsen-Birkhoff theorem in $f(R)$ theories Clifton:2012ry ; Nzioki:2013lca . Whence it seems that the most natural spacetime solutions as generated by a spherical stable neutron or quark star is not the Schwarzschild solution once Einsteinian gravity is replaced by simple extensions of the standard Einstein-Hilbert Lagrangian.

Concerning the observational imprints to be expected from our results, let us claim that the net effect on the total mass is a faint, although $\alpha$-parameter, dependent increase in the mass as perceived by a distant observer with respect to the General Relativity prediction for the same central density and equation of state. Also, we have concluded the increase of star radii in vicinity of maximal masses as well as the radii reduction for stars with masses of the order of the Sun mass.

Then, we can think of at least two ways for discriminating between Einsteinian gravity and possible modifications from scalar-tensor theories, and by extension other extended theories of gravity with analogous effects. Firstly, high-accuracy measurements of surface redshift can shed some light on the underlying theory. Indeed, for the same asymptotic gravitational mass, the pure stellar mass bounded by surface as generated by $R^{2}$-gravity would be smaller than the Einsteinian value and therefore the measured surface redshift would be smaller than predicted. Secondly, another observable effect would be the determination of gravitational maximal masses, which we have calculated for competitive both neutron plasma and quark stars, using X/ray-optical binaries and neutron star-white dwarf binaries. Notwithstanding, the present precision in mass measurements is not enough to discriminate between competing theories yet. Both methods have been shown here as competitive proposals to break the degeneracy problem afflicting alternative theories of gravity.

Finally, for $\alpha<0$ our analysis has shown that outside the star the additional energetic effect explained above, i.e., the appearance of a gravitational sphere, is also present. Nevertheless, the scalar curvature suffers damped oscillations instead of an exponential attenuation. Moreover, in these cases the contribution of the gravitational sphere into the gravitational mass has been shown to increase indefinitely with radial distance. Therefore we conclude that $R^{2}$-gravity with $\alpha<0$ predictions are not consistent with the existence of realistic models of compact relativistic stars, at least for the studied equations of state. The inadmissibility of $R^{2}$ gravity with $\alpha<0$ was made earlier in cosmological contexts so our result shows that these classes of models would not be acceptable in the astrophysical context either in order to obtained finite asymptotic masses.

As a final comment let us mention that, the value for $\alpha$ in the Starobinsky inflationary model in good agreement with Planck latest results provides a much smaller value Riotto than those considered in this article. Thus, the $R^{2}$ models presented here should be understood as effective corrections to General Relativity appearing at curvatures of the order of those characterising the dynamics of neutron/quark stars.

Acknowledgements.

We would like to thank Miguel Aparicio Resco for useful and meaningful discussions on the topic. The work of AVA is supported by project 1.4539.2017/8.9 (MES, Russia). SDO is supported by MINECO (Spain), project FIS2016-76363-P and by CSIC I-LINK1019 Project. AdlCD acknowledges financial support from Consolider-Ingenio MULTIDARK CSD2009-00064, i-Link1019, Spanish Ministry of Economy and Science and FIS2016-78859-P, European Regional Development Fund and Spanish Research Agency (AEI), Management Committee Member for the Cosmology group - University of Cape Town, CANTATA COST action CA15117, EU Framework Programme Horizon 2020, University of Cape Town Launching Grants Programme and National Research Foundation grant 99077 2016-2018, Ref. No. CSUR150628121624 and NRF Incentive Funding for Rated Researchers (IPRR), Ref. No. IFR170131220846.