Equation of State of Neutron Stars with Junction Conditions in the Starobinsky Model

TL;DR

This paper investigates neutron star structure within the Starobinsky $f(R)$ gravity model, solving coupled differential equations with junction conditions, and finds minimal mass and radius constraints based on the model parameters.

Contribution

It provides a direct solution method for neutron star equations in $f(R)$ gravity using junction conditions and derives constraints on star mass and radius.

Findings

Minimal neutron star mass around 1.44 solar masses.

Mass and radius decrease with larger polytropic constant $ar{k}$.

Constraint on $ ext{alpha}$ parameter for viable solutions.

Abstract

We study the Starobinsky or $R^2$ model of $f(R)=R+\alpha R^2$ for neutron stars with the structure equations represented by the coupled differential equations and the \emph{polytropic} type of the matter equation of state. The junction conditions of $f(R)$ gravity are used as the boundary conditions to match the Schwarschild solution at the surface of the star. Based on these the conditions, we demonstrate that the coupled differential equations can be solved \emph{directly}. In particular, from the dimensionless equation of state $\bar{\rho} = \bar{k}\, \bar{p}^{\,\gamma}$ with $\bar{k}\sim5.0$ and $\gamma\sim0.75$ and the constraint of $\alpha\lesssim {1.47722}\times 10^{7}\, \text{m}^2$, we obtain the \emph{minimal} mass of the NS to be around 1.44 $M_{\odot}$. In addition, if $\bar{k}$ is larger than 5.0, the mass and radius of the NS would be smaller.

Tables

Table 1. Table 1: The results of the radius x s = r s / r ∗ subscript 𝑥 𝑠 subscript 𝑟 𝑠 subscript 𝑟 x_{s}=r_{s}/r_{*} and mass M ¯ ¯ 𝑀 \bar{M} with the polytropic exponent γ 𝛾 \gamma and central Ricci curvature R ¯ ​ ( 0 ) ¯ 𝑅 0 \bar{R}(0) for the R 2 superscript 𝑅 2 R^{2} model respect to the various α ¯ ¯ 𝛼 \bar{\alpha} with the fixed central pressure p ¯ ​ ( 0 ) = 1 ¯ 𝑝 0 1 \bar{p}(0)=1 and polytropic constant k ¯ = 5.0 ¯ 𝑘 5.0 \bar{k}=5.0 .

$\overline{\alpha }$	$x_s$	$\overline{M}$	$\gamma$	$\overline{R}\left(0\right)$
0.01	1.999	1.444	0.7525000000	8.95
0.0005	2.477	1.557	0.7503553926	35.00
GR ($\alpha =0$)	2.297	1.672	0.7503553926	16$\pi$

Table 2. Table 2: The results of the mass M ¯ ¯ 𝑀 \bar{M} , radius x s subscript 𝑥 𝑠 x_{s} and Ricci curvature R ¯ ​ ( 0 ) ¯ 𝑅 0 \bar{R}(0) at the center for different values of the polytropic constant k ¯ ¯ 𝑘 \bar{k} with α ¯ = 0.01 ¯ 𝛼 0.01 \bar{\alpha}=0.01 and γ ∼ 0.75 similar-to 𝛾 0.75 \gamma\sim 0.75 in the R 2 superscript 𝑅 2 R^{2} model.

$\overline{k}$	5.0	4.5	4.0	3.5	3.2
$\overline{M}$	1.444	1.673	1.959	2.394	2.716
$x_s$	1.999	2.228	2.504	3.032	3.513
$\overline{R}\left(0\right)$	8.95	8.40	7.54	6.15	4.93

Equations

$$S=\frac{1}{2\kappa }\int d^{4}x\sqrt{-g}f(R)+S_m,$$

$$f^{\prime }{R}^{\mu \nu }-\frac{1}{2}f{g}^{\mu \nu }-({\nabla }^{\mu }{\nabla }^{\nu }-g^{\mu \nu }\square )f^{\prime }=\kappa T^{\mu \nu },$$

$$f(R)=R+\alpha R^2.$$

$$G_{\mu \nu }(1+2\alpha R)+\frac{\alpha }{2}{g}_{\mu \nu }{R}^2-2\alpha ({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square )R=\kappa T_{\mu \nu }$$

$$-R+6\alpha \square R=\kappa T.$$

$$d{s}^2=-e^{2\Phi (r)}d{t}^2+e^{2\Lambda (r)}d{r}^2+r^{2}d{\Omega }^2,$$

$$G_{tt}$$

$$G_{rr}$$

$$G_{\theta \theta }$$

$$G_{\phi \phi }$$

$$\Phi {}^{\prime }=\frac{-p{}^{\prime }}{\rho +p}.$$

$$\Lambda {}^{\prime }=\frac{r{m}^{\prime }-m}{r(r-2m)}$$

$$T_{tt}=\rho e^{2\Phi },T_{rr}=p{e}^{2\Lambda },T_{\theta \theta }=p{r}^2,T_{\phi \phi }=p{r}^2{\sin }^2\theta ,$$

$$\square R$$

$${\nabla }_t{\nabla }_{t}R$$

$${\nabla }_r{\nabla }_{r}R$$

$$m^{\prime }=$$

$$-\frac{\alpha (16\pi p{r}^3+4m(1+2\alpha R)-\alpha r^3{R}^2-8\alpha R^{\prime }r(r-2m))R^{\prime }}{4(1+2\alpha R)(1+2\alpha R+\alpha r{R}^{\prime })},$$

$$p^{\prime }=$$

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

$$\rho =k{p}^{\gamma }.$$

$$m_{\ast }$$

$$r_{\ast }$$

$${\rho }_{\ast }$$

$$p_{\ast }$$

$$R_{\ast }$$

$$\bar{m}{}^{\prime }=$$

$$-\frac{\bar{\alpha }((16\pi \bar{p}-\bar{\alpha }{\bar{R}}^2)x^3(R_{\ast }{r}_{\ast }^2)+4\bar{m}(1+2\bar{\alpha }\bar{R})(\frac{m_{\ast }}{r_{\ast }})-8\bar{\alpha }x(x-2\bar{m}(\frac{m_{\ast }}{r_{\ast }}))\bar{R}{}^{\prime })\bar{R}{}^{\prime }}{4(1+2\bar{\alpha }\bar{R})(1+2\bar{\alpha }\bar{R}+\bar{\alpha }x{\bar{R}}^{\prime })(\frac{m_{\ast }}{r_{\ast }})},$$

$$\bar{p}{}^{\prime }=$$

$$\bar{R}{}^{\prime \prime }=$$

$$+\frac{x^2((1+3\bar{\alpha }\bar{R})+16\pi \bar{\rho }){\bar{R}}^{\prime }(R_{\ast }{r}_{\ast }^2)}{(x-2\bar{m}(\frac{m_{\ast }}{r_{\ast }}))(1+2\bar{\alpha }\bar{R})}+\frac{2\bar{\alpha }\bar{R}{{}^{\prime }}^2}{1+2\bar{\alpha }\bar{R}},$$

$$ds^{2}=-\bigg{(}1-\frac{2\tilde{M}}{r}\bigg{)}\,dt^{2}+\bigg{(}1-\frac{2\tilde{M}}{r}\bigg{)}^{-1}\,dr^{2}+r^{2}\,d\Omega^{2},$$

$$[R]$$

$$[{\nabla }_{\mu }R]$$

$$\bar{R}(x_s)=0,\bar{R}{}^{\prime }(x_s)=0,\bar{p}(x_s)=0,\bar{m}(0)=0,\bar{R}{}^{\prime }(0)=0.$$

$$\frac{d\Phi }{dr}=\frac{16\pi p{r}^3+4m(1+2\alpha R)-\alpha r^3{R}^2-8\alpha r(r-2m)R^{\prime }}{4r(1+2\alpha R+\alpha r{R}^{\prime })(r-2m)}.$$

$$m^{\prime }=4\pi r^2\rho$$

$$-\frac{\alpha (16\pi p{r}^3+4m(1+2\alpha R)-\alpha r^3{R}^2-8\alpha r(r-2m)R^{\prime })R^{\prime }}{4r(1+2\alpha R)(1+2\alpha R+\alpha r{R}^{\prime })}\equiv 4\pi r^2{\rho }_{\text{eff}},$$

$${\rho }_{\text{eff}}=\rho$$

$$-\frac{\alpha (16\pi p{r}^3+4m(1+2\alpha R)-\alpha r^3{R}^2-8\alpha r(r-2m)R^{\prime })R^{\prime }}{16\pi r^3(1+2\alpha R)(1+2\alpha R+\alpha r{R}^{\prime })}.$$

Full text

I Introduction

The astrophysical observations from the Type Ia Supernovae Riess:1998cb ; Perlmutter:1998np , large scale structure Tegmark:2003ud ; Seljak:2004xh and baryon acoustic oscillations Eisenstein:2005su as well as cosmic microwave background Spergel:2003cb ; Spergel:2006hy ; Komatsu:2008hk indicate the necessity of new physics beyond the Einstein’s general relativity (GR). The modified theories of gravity Nojiri:2006ri ; Nojiri:2017ncd become more significant in order to explain the accelerated expansion phenomenon not only inflation Starobinsky:1980te ; Guth:1980zm ; Linde:1981mu ; Albrecht:1982wi ; Starobinsky:1987zz in the early epoch but also dark energy Weinberg:1988cp ; Sahni:1999gb ; Carroll:2000fy ; Peebles:2002gy ; Padmanabhan:2002ji ; Copeland:2006wr in the recent stage of the universe. A class of alternative theories of the modification from the geometric point of view is the so-called $f(R)$ gravity theories Barrow:1983rx ; Sotiriou:2008rp ; DeFelice:2010aj ; Nojiri:2010wj . In these theories, the Lagrangian density is modified by using an arbitrary function $f(R)$ instead of the scalar curvature $R$ of the Einstein-Hilbert term. The most well-known $f(R)$ model is the Starobinsky or $R^{2}$ model with $f(R)=R+\alpha R^{2}$, originally proposed to obtain the quasi-de Sitter solution for inflation Starobinsky:1987zz . Furthermore, several viable $f(R)$ gravity theories Hu:2007nk ; Starobinsky:2007hu ; Nojiri:2007as ; Tsujikawa:2007xu ; Cognola:2007zu ; ExpGra have been used to explain the cosmic acceleration problems.

In order to realize the structure of a compact star, one needs to know its equation of state (EoS), which characterizes the thermodynamic relation between the density $\rho$, pressure $p$ and temperature $T$ of the dense matter. Under the adiabatic assumption, the EoS is reduced to a polytropic relation $\rho=k\,p^{\,\gamma}$. This assumption has been discussed for the neutron stars (NSs) in the literature Lattimer:2012nd ; Chen:2015zpa ; Gandolfi:2015jma ; Reisenegger:2015crq ; Hebeler:2010jx ; Barausse:2007pn ; Santos:2009nu ; Cooney:2009rr ; Babichev:2009fi ; Reijonen:2009hi ; Orellana:2013gn ; Alavirad:2013paa ; Ganguly:2013taa ; Yazadjiev:2014cza . In particular, the allowed region of the polytopes has been shown in Hebeler:2013nza .

The compact relativistic star was first studied by Chandrasekhar Chandrasekhar , who assumed that a white dwarf is supported only by the completely degenerate electron gas, and then obtained the so-called Chandrasekhar limit of a white dwarf with the maximal mass of 1.44 $M_{\odot}$. Subsequently, Oppenheimer and Volkoff Oppenheimer:1939ne proposed a limit of $0.7\,M_{\odot}$ of a NS by considering a completely degenerate neutron gas. However, this approach is inappropriate due to the strong nuclear repulsive forces of neutrons and other strong interaction of the heavy hadrons in dense matter.

In the scenario of GR, the structure of the relativistic stars is determined by EoS of matter inside the stars without an explicit constraint, whereas it is expected that the $f(R)$ theories do provide some constraints with singularity problems Briscese:2006xu ; Nojiri:2008fk ; Frolov:2008uf . The relativistic stars in the modified gravities have been studied in the literature Kobayashi:2008tq ; Babichev:2009td ; Upadhye:2009kt ; Santos:2009nu ; Cooney:2009rr ; Babichev:2009fi ; Reijonen:2009hi ; Arapoglu:2010rz ; Jaime:2010kn ; Santos:2011ye ; Orellana:2013gn ; Cheoun:2013tsa ; Alavirad:2013paa ; Astashenok:2013vza ; Ganguly:2013taa ; Yazadjiev:2014cza ; Staykov:2014mwa ; Yazadjiev:2015zia ; Goswami:2015dma ; Staykov:2015mma ; Capozziello:2015yza ; Yazadjiev:2015xsj ; Hendi:2015pua ; Hendi:2015vta ; Bordbar:2015wva ; Hendi:2017ibm . It has been argued that the compact relativistic stars are difficult to exist due to the curvature scalar $R$ divergence inside the star in $f(R)$ Kobayashi:2008tq . However, the realistic EoS in the Starobinsky’s dark energy model Starobinsky:2007hu has been constructed in Ref. Babichev:2009td , in which $R$ does not diverge inside the star, so that the relativistic stars could occur in $f(R)$. The pure geometric study is formulated in Ref. Ganguly:2013taa , which imposes the junction conditions in $f(R)$ Senovilla:2013vra ; Deruelle:2007pt as the additional conditions to solve the coupled structure equations and obtain the final result indirectly.

In this study, we consider the $R^{2}$ model by performing the calculation only in the Jordan frame. In our discussions, we solve the coupled structure equations by the junction conditions approach directly rather than the perturbation methods Cooney:2009rr ; Arapoglu:2010rz ; Orellana:2013gn ; Cheoun:2013tsa 111Some other non-perturbative methods have been addressed in Refs. Astashenok:2014dja ; Capozziello:2011nr ; Astashenok:2017dpo . We show that the NSs can exist in the $R^{2}$ model under the polytrope assumption of EoS. The possible dimensionless EoS $\bar{\rho}\sim 5.0\,\bar{p}^{\,0.75}$ is concluded by the analysis of the various values of the dimensionless parameter $\bar{\alpha}$ in the $R^{2}$ model, where the bars represent the dimensionless quantities. The theoretical constraint on the coefficient $\alpha$ of the $R^{2}$ term in the model is given by $\alpha\lesssim 1.47722\times 10^{7}\,\text{m}^{2}$. By applying the resultant EoS and critical value of $\alpha$, the minimal mass of the NSs is obtained about 1.44 $M_{\odot}$ which is the same as the Chandrasekhar limit of the white dwarf Chandrasekhar . For a fixed parameter $\alpha$, we observe that the mass and the radius get larger when $k$ decreases, while the maximal value of $\bar{k}=5.0$ can be illustrated.

This paper is organized as follows. In Sec. II, we derive the coupled differential equations and show the boundary conditions for the spherically symmetric compact stars in the $R^{2}$ model. In Sec. III, we analyze the model parameter $\alpha$ and explore its reasonable value from the typical units in the neutron star system. We discuss our result of EoS under the specific choice of the initial conditions. Finally, we give conclusions in Sec. IV.

II Spherically Symmetric Solution of the $R^{2}$ model

The action of the $f(R)$ theories with matter is given by

[TABLE]

with $\kappa=8\,\pi$ and the conventional units of $G=c=1$. By the variation with respect to the metric $g_{\mu\nu}$, we have the modified Einstein equations

[TABLE]

with $T^{\mu\nu}$ the energy-momentum tensor and $\square=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ the D’Alembertian operator. In addition, “$\,{}^{\prime}\,$” in this paper denotes the differentiation with respect to its argument, e.g. $f^{\prime}(R)=df(R)/dR$. We will focus on the Starobinsky or $R^{2}$ model with the function of the Lagrangian density

[TABLE]

As a result, we obtain the following field equation

[TABLE]

In order to study the system of a compact star, we will study the solution with an ansatz given by the static spherical symmetric metric

[TABLE]

where $\quad d\Omega^{2}=d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}$ and $\exp(2\Lambda(r))=(1-2m(r)/r)^{-1}$ with $m(r)$ the mass function characterizing the mass enclosed within the radius $r$. In GR, $m(r)=\int^{r}_{0}4\pi\bar{r}^{2}\,\rho(\bar{r})\,d\bar{r}$ with $\rho(\bar{r})$ the density function. For the radius of the star $r_{s}$, $m(r_{s})=M$ can be identified as the total mass in the Newtonian limit. In the $R^{2}$ model, the mass function should be modified with some correction terms. However, it cannot be integrated by the density function $\rho$ directly. The function $\Phi(r)$ can be regarded as the effective relativistic gravitational potential. Subsequently, we can obtain the Einstein tensor from (5), given by

[TABLE]

II.1 Coupled Differential Equations

Considering a static perfect fluid with the energy-momentum tensor $T^{\mu\nu}=(\rho+p)u^{\mu}u^{\nu}+pg^{\mu\nu}$ with $u^{\mu}$, $\rho$, and $p$ denoting the 4-velocity, the density, and the pressure of the fluid respectively. The $\nu=r$ component of the conservation equation $\nabla_{\mu}T^{\mu\nu}=0$ gives

[TABLE]

In addition, we can obtain the identity

[TABLE]

via the definition of $\Lambda(r)$. In the local rest frame, $u_{t}=-e^{\Phi}$ and $u_{i}=0$ with $i=r,\theta$ and $\varphi$ the spatial coordinates, we have

[TABLE]

due to $u_{\mu}u^{\mu}=-1$. In addition, we can obtain the following identities for convenience

[TABLE]

Consequently, with the metric given by (5) and the energy momentum tensor given by the static perfect fluid, we can write the field equations as a set of differential equations

[TABLE]

Note that we have written $m^{\prime}=m^{\prime}(R,R^{\prime},p,m)$, $p^{\prime}=p^{\prime}(R,R^{\prime},p,m)$ and $R^{\prime\prime}=R^{\prime\prime}(R,R^{\prime},p,m)$ as algebraic functionals of $R,R^{\prime},p,m$. Here Eq. (13) is derived from the $tt$-component and the trace equation of the field equation (4a). In addition, Eq. (13b) is derived from the $rr$-component of the field equation (4a), and is also known as the modified Tolman-Oppenheimer-Volkoff (mTOV) equation Oppenheimer:1939ne . Finally, Eq. (13c) is derived from the trace equation (4b).

For the perfect fluid, we assume the EoS is polytrope, i.e.,

[TABLE]

In order to simplify the calculations, we can choose the typical values $r_{*}$, $m_{*}$, $p_{*}$, $\rho_{*}$ and $R_{*}$ for the compact star system and express $r\equiv xr_{*}$, $m\equiv\bar{m}m_{*}$, $p\equiv\bar{p}p_{*}$, $\rho\equiv\bar{\rho}\rho_{*}$, $R\equiv\bar{R}R_{*}$ and $\alpha\equiv\bar{\alpha}\alpha_{*}\equiv\bar{\alpha}(1/R_{*})$ in terms of the dimensionless quantities $x$, $\bar{m}$, $\bar{p}$, $\bar{\rho}$, $\bar{R}$ and $\bar{\alpha}$, while the derivatives of $p$, $m$ and $R$ can be written as $p^{\prime}=\bar{p}^{\prime}(p_{*}/r_{*})$, $m^{\prime}=\bar{m}^{\prime}(m_{*}/r_{*})$, $R^{\prime}=\bar{R}^{\prime}(R_{*}/r_{*})$ and $R^{\prime\prime}=\bar{R}^{\prime\prime}(R_{*}/r_{*}^{2})$, respectively, where the prime of the dimensionless quantities denotes the derivative with respect to $x$. The polytropic type of EoS in terms of the dimensionless quantities can be given as $\bar{\rho}=\bar{k}\,\bar{p}^{\,\gamma}$ with $\bar{k}=k{\rho^{-1}_{*}}p_{*}^{\gamma}$. Since we are interested in the NSs in the $R^{2}$ model, it is convenient for us to define the following typical values in SI units,

[TABLE]

According to the typical units, we can rewrite (13), (13b) and (13c) as the dimensionless equations:

[TABLE]

respectively. The dimensionless parameters $m_{*}/r_{*}=0.147688$ and $R_{*}r_{*}^{2}=0.074215$ characterize the compactness of the NS. In order to discuss the structure of the NS, we have to solve the three coupled equations (15), (15b) and (15) numerically with the EoS $\bar{\rho}=\bar{k}\,\bar{p}^{\,\gamma}$.

II.2 Boundary Conditions

In GR, the Birkhoff’s theorem states that the spherically symmetric vacuum solution must be given by the Schwarzschild metric. On the other hand, even though the absence of the Birkhoff’s theorem in $f(R)$ theories might lead to the non-uniqueness of this vacuum solution, the Schwarzschild metric can serve as a vacuum solution in $f(R)$ under some circumstances. It has been shown that the conditions of $R=0$ with $f(0)=0$ and $f^{\prime}(0)\neq 0$ for the existence of the Schwarzschild metric are satisfied in the Starobinsky model Nzioki:2013lca . As a result, we introduce the Schwarzschild vacuum solution for the exterior region. In this way, we can obtain the mass and radius of the star from the Schwarzschild metric once (15) is solved with proper boundary conditions.

In the following, we consider the star without thin shells. In order to match the solution at the surface of the star, we use the Schwarzschild solution for the exterior region ($r>2\tilde{M}$)

[TABLE]

where $\tilde{M}$ is the mass parameter in GR. The junction conditions for the $f(R)$ theories should be more restrictive as discussed in Refs. Senovilla:2013vra ; Ganguly:2013taa ; Deruelle:2007pt . The first and the second fundamental forms of the conditions are $[h_{\mu\nu}]=0$ and $[K_{\mu\nu}]=0$, respectively, where $[\,]$ denotes the jump at the boundary surface of the star. We can identify $\tilde{M}$ with $M=m(r_{s})$ only when the first fundamental form matches. However, there are two additional conditions for the scalar curvature across the surface Ganguly:2013taa , given by

[TABLE]

In our assumption with the static and spherically symmetric metric, the curvature $R$ is only a function of $r$. By matching of the second fundamental form to make the pressure vanishing at the boundary surface Senovilla:2013vra , the boundary conditions are reduced to $R(r_{s})=0$, $R^{\prime}(r_{s})=0$ and $p(r_{s})=0$. Inside the star, we have to determine the boundary conditions at the center of the star. There are two first-order and one second-order differential equations in Eq. (13). Hence, only four boundary conditions are required to solve these coupled ordinary differential equations. To satisfy the regularity conditions at the center of the star, we must have $m(0)=0$, $p^{\prime}(0)=0$, $\rho{{}^{\prime}}(0)=0$ and $R^{\prime}(0)=0$ Ganguly:2013taa , in which two of them are redundant. According to Eq. (13b), $p^{\prime}(0)=0$ is automatically satisfied as long as $m(0)=0$ and $R^{\prime}(0)=0$ as $r\rightarrow 0$. In addition, $\rho$ and $p$ are related by EoS in (14), leading to $p^{\prime}(0)=0$ and $\rho{{}^{\prime}}(0)=0$, so that only conditions $m(0)=0$ and $R^{\prime}(0)=0$ are left.

Consequently, we have three boundary conditions at the surface and two boundary ones at the center written in the dimensionless forms, given by

[TABLE]

These boundary conditions are referred to as the Schwarzschild boundary conditions. Mathematically, since there are four undetermined integration constants $c_{1},c_{2},c_{3}$ and $c_{4}$ in (15), only four in (18) are enough to solve it. However, these integration constants should be associated with the model parameter $\alpha$ and $(\gamma,\bar{k})$ in the EoS. The fifth one in (18) can be used to constrain the parameter space of $(\alpha,\gamma,\bar{k})$. For example, if we choose $\bar{m}(0)=\bar{R}{{}^{\prime}}(0)=\bar{p}(x_{s})=\bar{R}(x_{s})=0$, then we have to determine whether $\bar{R}{{}^{\prime}}(\alpha,\gamma,\bar{k};x)|_{x=x_{s}}$ satisfies $\bar{R}{{}^{\prime}}(x_{s})=0$ for fixed values of $\alpha$, $\gamma$ and $\bar{k}$.

According to the mTOV equation in (13) and conservation equation in (7), we have

[TABLE]

In the region outside of the star ($r\geq r_{s}$), the pressure and scalar curvature as well as the derivative of the scalar curvature should be continuous, resulting in $p(r)=0$, $R(r)=0$ and $R^{\prime}(r)=0$ by (18). It can be checked that the exterior solution of (19) coincides with the Schwarzschild solution $e^{2\Phi(r)}=1-2\tilde{M}/r$.

III Analysis and Results

III.1 Determination of $\alpha$

In principle, Eq. (13) can be regarded as the GR results with $\alpha R^{2}$ as the modification term. For example, Eq. (13b) corresponds to the TOV equation in GR Oppenheimer:1939ne when $\alpha\rightarrow 0$. Similarly, we can recover $m^{\prime}$ and $R^{\prime\prime}$ equations in GR for (13) and (13c) with $\alpha\rightarrow 0$. By separating the GR contribution, Eq. (13) can be rewritten as

[TABLE]

where

[TABLE]

In the limits of $\alpha\rightarrow 0$ and $R\rightarrow 8\pi(\rho-3p)$, we have $m^{\prime}\rightarrow 4\pi r^{2}\rho$, which is the same as result in GR.

However, in the numerical analysis, there are problems of choosing $\alpha$ for the system. On one hand, the main numerical difficulty arises from (13c), in which

[TABLE]

by taking $\alpha\rightarrow 0$. Furthermore, we have the boundary conditions $R^{\prime}(0)=0$ and $m(0)=0$ as $r\rightarrow 0$, and obtain

[TABLE]

which implies the singularity of $\mathcal{R}^{\prime\prime}$ as $\alpha\rightarrow 0$ under the numerical calculation. As a result, we encounter the fine-tuning problem of $p(0)$ and $R(0)$. On the other hand, we would like to discuss the upper bound for $\bar{\alpha}$. In the dimensionless form $x=r/r_{*}$, (15) with (III.1) and (23) in $x\rightarrow 0$ can be read as

[TABLE]

and

[TABLE]

respectively, where $(\rho_{*}\,r_{*}^{3})/m_{*}={0.502513}$ and $R_{*}r_{*}^{2}=\rho_{*}r_{*}^{2}=7.42\times 10^{-2}$, which characterizes the compactness of a star. In order to determine the proper value of $\bar{\alpha}$, we use (25) to rewrite (24) as

[TABLE]

Since $\bar{R}$ is convex upward around $x=0$, we expect that $\bar{\mathcal{R}}{{}^{\prime\prime}}\leq 0$. Then, we have $8\pi(\bar{\rho}-3\bar{p})-\bar{R}\geq 0$ for $\bar{\alpha}\geq 0$, which can be seen from (25). We can choose $\bar{\alpha}\lesssim(m_{*}/r_{*})(R_{*}r_{*}^{2})={0.010961}$ and obtain the inequality

[TABLE]

The last two terms in the square bracket represent the first-order and second-order corrections in the $R_{*}r_{*}^{2}$ unit, respectively. Therefore, we derive $\alpha=\bar{\alpha}/R_{*}\lesssim 1.47722\times 10^{7}\,\text{m}^{2}$. In addition, several constraints on $\alpha$ from the observational data have been derived in  Arapoglu:2010rz ; Cheoun:2013tsa ; Naf:2010zy . Moreover, Gravity Probe B Everitt:2009 gives $\alpha\lesssim 5\times 10^{11}\,\text{m}^{2}$; the precession measurement of the pulsar B in the PSR J0737-3039 system Breton:2008xy yields $\alpha\lesssim 2.3\times 10^{15}\,\text{m}^{2}$; and the strong magnetic NS Arapoglu:2010rz ; Cheoun:2013tsa results in $\alpha\lesssim 10^{5}\,\text{m}^{2}$. Furthermore, it has been shown that the ghost-free condition $f^{\prime\prime}(R)\geq 0$ Sotiriou:2008rp leads to $\alpha>0$. Here, it should be noted that only within the condition $\bar{R}\rightarrow 8\pi(\bar{\rho}-3\bar{p})$ can we have a finite $\bar{\mathcal{R}}{{}^{\prime\prime}}$ in the limit $\bar{\alpha}\rightarrow 0$. This condition assures that the $R^{2}$ model is consistent with GR in $\alpha\rightarrow 0$.

III.2 Numerical Results

By using the Runge-Kutta 4th-order (RK4) procedure, Eq. (15) can be solved by choosing $\bar{p}(0)$ and $\bar{R}(0)$ as the central values with boundary conditions $\bar{m}(0)=\bar{R}{{}^{\prime}}(0)=\bar{p}(x_{s})=0$. We can obtain $\bar{R}(x_{s})$ and $\bar{R}{{}^{\prime}}(x_{s})$ by applying random values of $\bar{p}(0)$ and $\bar{R}(0)$ numerically. In terms of the problem of (23), we have to find out the appropriate values of $\bar{p}(0)$ and $\bar{R}(0)$ to satisfy $\bar{R}(x_{s})=0$ and $\bar{R}{{}^{\prime}}(x_{s})=0$, which maintain the Schwarzschild boundary conditions(18).

The parameters $\bar{k}$ and $\gamma$ affect the behaviors of the coupled equations (15) as well as the boundary values at the surface. Clearly, they can be determined once our boundary conditions are fixed in the numerical calculations. All the results are given in the typical units $m_{*}$,222$m_{*}=M_{\odot}$ $r_{*}$,333$r_{*}=10$ km $\rho_{*}$, $p_{*}$, and $R_{*}$ as defined in Sec. IIA. We look for the reasonable EoS for $\bar{\alpha}=0.01$ and $0.0005$ and compare the results with GR ($\alpha=0$). For simplicity, we keep the high pressure at the center of the star to be $\bar{p}(0)=1$ initially. Then, we end up the calculation with $\bar{p}(x_{s})=10^{-6}$ at the surface of the star, corresponding to the density at the bottom of the NS’s outer crust around $10^{13}\sim 10^{14}\,\text{kg}\,\text{m}^{-3}$. We keep $\bar{k}=5.0$ and fine-tune the parameters $\gamma$ and $\bar{R}(0)$ in order to satisfy $\bar{R}(x_{s})=0$ and $\bar{R}^{\prime}(x_{s})=0$. The results are given in the TABLE 1. From this table, we find that for a smaller $\bar{\alpha}$, $\bar{R}(0)$ is larger, and the same goes for $\bar{M}$, which are the generic feature of the model. The behaviors of the growing $\bar{\alpha}$ and decreasing $\bar{M}$ have been also discussed in Ref. Astashenok:2017dpo with the realistic EoS instead of the polytropic one in this study. We note that the different choices of $\bar{k}$ will be shown in TABLE 2.

In FIG. 1, we illustrate the deviation of the interior region of the star in the $R^{2}$ model from GR. The profiles of the scalar curvature $\bar{R}$ and its derivative $\bar{R}^{\prime}$ are shown in FIG. 1a. Clearly, these two quantities satisfy the boundary conditions $\bar{R}(x_{s})=0$ and $\bar{R}^{\prime}(x_{s})=0$. The results of the negative trace of the energy momentum tensor in GR and the $R^{2}$ model in the interior of the NSs are displayed in FIG. 1b, illustrating similar behaviors. However, the conduct of the scalar curvature in the $R^{2}$ model is different from that of GR with $R=-\,8\pi T=\kappa\,(\rho-3p)$ due to the $R^{2}$ term.

For the density $\bar{\rho}$ and mass function $m/M_{\odot}$ profiles of the star, we exhibit $\bar{k}=5.0$ with $\bar{\alpha}=0.01$ and $0.0005$ in FIG. 2. We see that the deviation of the density in the $R^{2}$ model from GR is small in FIG. 2a, whereas that of the resultant mass is large in FIG. 2b. The end-points of the curves in FIG. 2b correspond to the mass $\bar{M}$ and radius $x_{s}$ shown in TABLE 1. In the $R^{2}$ model, the mass function in (III.1) is deviated from GR due to the geometric effect of the effective density $\rho_{\text{eff}}$.

In TABLE 1, the mass of the NS exceeds the Chandrasekhar limit ($1.44\,M_{\odot}$) of the white dwarf Chandrasekhar . Note that if the collapsing process is supplied only by gravity, the Chandrasekhar limit could be considered as a lower bound of the mass for a star whose ultimate destiny is a NS or black hole.

According to our analysis in the $R^{2}$ model, which allows a lighter NS than that in GR as shown in FIG. 2b. Furthermore, from the upper limit $\alpha\lesssim 1.47722\times 10^{7}\,\text{m}^{2}$, we find that the minimal mass of the NS is around 1.44 $M_{\odot}$ for $\gamma\sim 0.75$ and $\bar{k}=5.0$ (see solid line in FIG. 2b).

By fixing $\alpha$ equal to the critical value $\bar{\alpha}=0.01$, we can analyze the properties of the the NS in the minimal mass condition. The profiles of the mass function $m$ of the radial coordinate $r$ with $\gamma\sim 0.75$, $\bar{k}=5.0$, 4.5, 4.0, 3.5 and $3.2$ are shown in FIG. 3, respectively. In TABLE 2, we list the mass $\bar{M}$ and radius $x_{s}$ of the NSs and their corresponding Ricci curvature $\bar{R}(0)$ at the center. From the table, we observe that the mass becomes larger as $\bar{k}$ gets smaller, whereas $\bar{R}(0)$ becomes smaller. We note that the case of $\bar{k}=3.0$ due to $\bar{\rho}>3\bar{p}$ for ordinary matter inside the NS has been excluded in our discussion. On the other hand, we expect that the mass of the NS is not smaller than the Chandrasekhar limit and the value $\bar{k}$ can not be larger than 5.0. The reasonable maximal value of $\bar{k}$ can be determined as 5.0.

IV Conclusions

We have addressed the $R^{2}$ model on a compact star, especially on the NS through the junction conditions. We have solved the mTOV equation rather than the perturbation method in the literature. In order to satisfy the junction conditions (Schwarzschild conditions), the central pressure $p(0)$ and Ricci scalar $R(0)$ should be well-selected. In $f(R)$ gravity, more specifically, the $R^{2}$ model, the parameters $k$ and $\gamma$ in the polytropic EoS can be constrained by $p(0)$ and $R(0)$ due to the coupled structure equations. With the junction condition, in particular, we have shown that there exists the solution of EoS $\bar{\rho}=\bar{k}\,\bar{p}^{\,\gamma}$ with $\bar{k}\sim 5.0$ and $\gamma\sim 0.75$.

For the upper limit $\alpha=1.47722\times 10^{7}\,\text{m}^{2}$, we have obtained the minimal mass of the NS. Under $\bar{\rho}=5.0\,\bar{p}^{\,0.75}$, the typical value of the NS mass is around 1.44 $M_{\odot}$. We have shown that $\bar{k}$ has the maximal value of $\bar{k}=5.0$. In our discussion, we have only considered the ghost-free $f(R)$ theories ($\alpha>0$). One could have heavier NSs when taking a negative $\alpha$ into account under the polytrope assumption in Ref. Orellana:2013gn . For $\alpha>0$, our result of the polytropic EoS is consistent with that in Ref. Orellana:2013gn .

Finally, we remark that in our derivation, we have obtained the same coupled structure equations (15) as those in Ref. Ganguly:2013taa after some proper arrangements. However, by fine-tuning the EoS, we have solved (15) with the junction conditions (18) directly rather than the indirect method used in Ref. Ganguly:2013taa .

Acknowledgements.

We thank Professor Radouane Gannouji for communications and discussions. This work was partially supported by National Center for Theoretical Sciences and MoST (MoST-104-2112-M-007-003-MY3 and MoST-104-2112-M-009-020-MY3).