Neutron stars in $f(R)$ gravity and scalar-tensor theories

TL;DR

This paper investigates neutron star structures within $f(R)$ gravity and scalar-tensor theories, analyzing how modifications to gravity influence star properties and stability, with numerical solutions for different potentials and equations of state.

Contribution

It provides the first detailed numerical analysis of neutron stars in $f(R)$ and Brans-Dicke theories with various scalar potentials, highlighting stability issues and modifications to mass-radius relations.

Findings

Scalar coupling affects neutron star radius but not maximum mass.

Massive potentials pose stability challenges for neutron star solutions.

Certain scalar potentials allow stable neutron star solutions with modified mass-radius relations.

Abstract

In $f(R)$ gravity and Brans-Dicke theory with scalar potentials, we study the structure of neutron stars on a spherically symmetric and static background for two equations of state: SLy and FPS. In massless BD theory, the presence of a scalar coupling $Q$ with matter works to change the star radius in comparison to General Relativity, while the maximum allowed mass of neutron stars is hardly modified for both SLy and FPS equations of state. In Brans-Dicke theory with the massive potential $V(\phi)=m^2 \phi^2/2$, where $m^2$ is a positive constant, we show the difficulty of realizing neutron star solutions with a stable field profile due to the existence of an exponentially growing mode outside the star. As in $f(R)$ gravity with the $R^2$ term, this property is related to the requirement of extra boundary conditions of the field at the surface of star. For the self-coupling potential $V(\phi)=\lambda \phi^4/4$, this problem can be circumvented by the fact that the second derivative $V_{,\phi \phi}=3\lambda\phi^2$ approaches 0 at spatial infinity. In this case, we numerically show the existence of neutron star solutions for both SLy and FPS equations of state and discuss how the mass-radius relation is modified as compared to General Relativity.

Figures

Equations

$$\mathcal{S}=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[\frac{M_{\mathrm{pl}}^2}{2}F(\varphi )R+\left(1-6Q^2\right)F(\varphi )X-V(\varphi )\right]+\int {\mathrm{d}}^{4}x{\mathcal{L}}_m\left(g_{\mu \nu },{\Psi }_m\right),$$

$$F(\varphi )=e^{-2Q\varphi /M_{\mathrm{pl}}}.$$

$${\nabla }^{\mu }{T}_{\mu \nu }=0,$$

$${\mathcal{S}}_{f(R)}=\int {\mathrm{d}}^{4}x\sqrt{-g}\frac{M_{\mathrm{pl}}^2}{2}f(R)$$

$$Q=-\frac{1}{\sqrt{6}},V(\varphi )=\frac{M_{\mathrm{pl}}^2}{2}\left(FR-f\right),F=\frac{\partial f}{\partial R}=e^{-2Q\varphi /M_{\mathrm{pl}}}.$$

$${\mathcal{S}}_{\mathrm{dilaton}}=\int {\mathrm{d}}^{4}x\sqrt{-g}\frac{e^{-\Phi /M_{\mathrm{pl}}}}{2}\left(M_{\mathrm{pl}}^{2}R+g^{\mu \nu }{\partial }_{\mu }\Phi {\partial }_{\nu }\Phi \right).$$

$$\mathrm{d}{s}^2=-f(r)\mathrm{d}{t}^2+h^{-1}(r)\mathrm{d}{r}^2+r^2\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\phi }^2\right),$$

$$P^{\prime }+\frac{f^{\prime }}{2f}\left(\rho +P\right)=0,$$

$$\xi ={\log }_{10}(\rho /{\mathrm{g}\cdot \mathrm{cm}}^{-3}),\zeta ={\log }_{10}(P/{\mathrm{dyn}\cdot \mathrm{cm}}^{-2}),$$

$$\zeta (\xi )$$

$$f_0(x)={\left(e^x+1\right)}^{-1},$$

$$h(r)=1-\frac{2G\mathcal{M}(r)}{r},$$

$$M\equiv\lim_{r\to\infty}{\cal M}(r)=\frac{r}{2G}\left(1-h\right)\biggl{|}_{r\to\infty}\,.$$

$$P(r_s)=0.$$

$$\frac{f^{\prime }}{f}$$

$${\mathcal{M}}^{\prime }$$

$${\varphi }^{\prime \prime }$$

$$R=\frac{1}{M_{\mathrm{pl}}^2}\left\{(1-6Q^2)\left[h{\varphi }^{\prime 2}+\left(4V+\rho -3P\right)e^{2Q\varphi /M_{\mathrm{pl}}}\right]-6M_{\mathrm{pl}}Q{V}_{,\varphi }{e}^{2Q\varphi /M_{\mathrm{pl}}}\right\},$$

$$f^{\prime }(r=0)=0,h^{\prime }(r=0)=0,{\varphi }^{\prime }(r=0)=0,{\rho }^{\prime }(r=0)=0.$$

$$f(r)=f_0+\sum _{n=2}^{\infty }{f}_n{r}^n,h(r)=1+\sum _{n=2}^{\infty }{h}_n{r}^n,\varphi (r)={\varphi }_0+\sum _{n=2}^{\infty }{\varphi }_n{r}^n,\rho (r)={\rho }_0+\sum _{n=2}^{\infty }{\rho }_n{r}^n,$$

$$V(\phi)=V(\phi_{0})+\sum_{n=1}^{\infty}\frac{1}{n!}\frac{{\rm d}^{n}V}{{\rm d}\phi^{n}}\biggl{|}_{\phi=\phi_{0}}(\phi-\phi_{0})^{n}\,.$$

$$f(r)$$

$$h(r)$$

$$\varphi (r)$$

$$P(r)$$

$$f(r\to \infty )=1,h(r\to \infty )=1,{\varphi }^{\prime }(r\to \infty )=0,V({\varphi }_{\infty })=0,$$

$${\tilde{\rho }}_0$$

$$r_0$$

$$y\equiv \frac{\rho }{{\tilde{\rho }}_0},z\equiv \frac{P}{{\tilde{\rho }}_0},v\equiv \frac{V}{{\tilde{\rho }}_0},v_{,\phi }\equiv \frac{M_{\mathrm{pl}}{V}_{,\varphi }}{{\tilde{\rho }}_0},m\equiv \frac{3\mathcal{M}}{4\pi r_0^3{\tilde{\rho }}_0},\phi \equiv \frac{\varphi }{M_{\mathrm{pl}}},s\equiv \ln \frac{r}{r_0}.$$

$$\xi ={\alpha }_1+{\alpha }_2\ln y,\zeta ={\alpha }_3+{\alpha }_2\ln z,$$

$$z=\exp \left[\frac{\zeta (\xi )-{\alpha }_3}{{\alpha }_2}\right],$$

$$y_{,s}=-\frac{y(y+z)}{2z}{\left(\frac{\mathrm{d}\zeta }{\mathrm{d}\xi }\right)}^{-1}\frac{f_{,s}}{f}.$$

$$\frac{f_{,s}}{f}=-\frac{2(h-1)-16\pi e^{2s+2Q\phi }(z-v)+h[(6Q^2-1){\phi }_{,s}-8Q]{\phi }_{,s}}{2h(1-Q{\phi }_{,s})},$$

Full text

I Introduction

The dawn of gravitational-wave (GW) astronomy opened up a new possibility for probing the physics in the strong-field regime Abbott:2016blz . The accuracy of General Relativity (GR) is well confirmed on the weak gravitational background Will:2014kxa , but the theory can be subject to modifications in the region of high density with the large scalar curvature $R$. Now, we are entering the golden era in which the deviation from GR can be tested from the GW observations of strong gravitational sources such as black holes (BHs) TheLIGOScientific:2016src and neutron stars (NSs) TheLIGOScientific:2017qsa .

One of the simplest modifications from GR is known as $f(R)$ theories, in which the Lagrangian contains nonlinear functions of $R$ Bergmann ; Ruz . In the Starobinsky $f(R)$ model characterized by the Lagrangian $f(R)=R+R^{2}/(6m^{2})$ Staro , where $m^{2}$ is a positive mass squared, the existence of the $R^{2}$ term can drive cosmic inflation in the early Universe. In the same manner, $f(R)$ theories have been extensively applied to the physics of late-time cosmic acceleration fRearly1 ; fRearly1b ; fRearly2 ; fRearly3 ; fR1 ; fR2 ; fR3 ; fR4 ; fR5 . The existence of higher-curvature term in $f(R)$ theories can also be potentially important in local objects on the strong gravitational background. In this vein, many papers devoted to the study of spherically symmetric and static BH delaCruzDombriz:2009et ; Nelson:2010ig ; Kehagias:2015ata ; Canate:2015dda ; Yu:2017uyd ; Canate:2017bao ; Sultana:2018fkw and NS Cooney:2009rr ; Arapoglu:2010rz ; Orellana:2013gn ; Astashenok:2013vza ; Ganguly:2013taa ; Yazadjiev ; Capozziello:2015yza ; Resco:2016upv ; Feng:2017hje solutions in the Starobinsky model and pure $R^{2}$ gravity. We note that $f(R)$ theories are equivalent to Brans-Dicke (BD) theories Brans with a scalar potential of the gravitational origin Ohanlon ; Chiba03 (see also Ref. Chakraborty ). In BD theories with the potential of a positive mass squared, there is the “no-hair” theorem of BHs forbidding the existence of a nontrivial scalar hair Hawking:1972qk ; Bekenstein:1995un ; Sotiriou:2011dz . This does not allow the presence of hairy BH solutions in the Starobinsky model delaCruzDombriz:2009et ; Nelson:2010ig ; Yu:2017uyd ; Canate:2017bao .

In BD theories, the nonminimal coupling to gravity mediates the fifth force between the scalar field $\phi$ and matter Fujii . This property is particularly transparent in the Einstein frame where the field $\phi$ directly interacts with matter with a universal coupling constant $Q$ chame1 ; chame2 ; Amendola:2006kh ; Yoko . For example, the $f(R)$ gravity in the metric formalism corresponds to $Q=-1/\sqrt{6}$ fRreview . This scalar-matter coupling plays an important role for studying the existence of NS solutions in BD theories and $f(R)$ gravity. In the Starobinsky $f(R)$ model mentioned above, the potential in the Einstein frame is given by $V_{{\rm E}}(\phi)=(3/4)m^{2}M_{\rm pl}^{2}[1-e^{-\sqrt{6}\phi/(3M_{\rm pl})}]^{2}$, where $\phi=(\sqrt{6}M_{\rm pl}/2)\ln[1+R/(3m^{2})]$ and $M_{\rm pl}$ is the reduced Planck mass fRreview . In the regime $|\phi|\ll M_{\rm pl}$, this potential reduces to $V_{{\rm E}}(\phi)\simeq m^{2}\phi^{2}/2$, so the scalar field has a constant mass $m$ with the matter coupling $Q=-1/\sqrt{6}$.

If we apply $f(R)$ theories to NSs on the spherically symmetric and static background, both the potential and matter coupling contribute to the scalar-field equation inside the star. For the vacuum exterior, the $\phi$ derivative of Einstein-frame potential $V_{{\rm E}}(\phi)$ mostly determines the field profile outside the star. For the massive potential $V_{{\rm E}}(\phi)=m^{2}\phi^{2}/2$, the field equation contains a growing-mode solution of the form $\phi\propto e^{mr}/r$ outside the body, where $r$ is the distance from the center of symmetry. In this case, we do not realize the asymptotic flat boundary conditions $\phi\to 0$ and ${\rm d}\phi/{\rm d}r\to 0$ at spatial infinity. This exponential growth of $\phi$ can be avoided by imposing the boundary conditions $\phi=0$ and ${\rm d}\phi/{\rm d}r=0$ at the surface of star, which amounts to considering the Schwarzschild exterior without any scalar-field contribution to the metric. However, as claimed in Ref. Ganguly:2013taa , the existence of such additional conditions does not allow the natural realization of NS solutions for arbitrary equations of state (EOSs). This property holds not only for $f(R)$ gravity with the massive potential $V_{{\rm E}}(\phi)=m^{2}\phi^{2}/2$ but also for the Starobinsky $f(R)$ model.

On the other hand, the $f(R)$ models of late-time cosmic acceleration fR1 ; fR2 ; fR3 ; fR4 are constructed to have a density-dependent effective scalar mass $m_{\phi}$ to accommodate the chameleon mechanism chame1 ; chame2 ; Faulkner:2006ub ; Capozziello:2007eu in over-density regions with the nonrelativistic background (see also Ref. Brax ). In such cases the scalar field is heavy inside the NS, but it can be practically massless outside the star. In spite of accessible curvature singularities in those $f(R)$ dark energy models fR2 ; Frolov:2008uf ; Kobayashi:2008tq ; Kobayashi:2008wc , the existence of relativistic stars was shown for a constant density profile Upadhye:2009kt and for a polytropic EOS Babichev:2009td ; Babichev:2009fi . In BD theories with the inverse power-law potential $V_{{\rm E}}=M^{4+n}\phi^{-n}$ ($n>0$), there are also chameleon-like solutions for relativistic stars with the constant density Tsujikawa:2009yf . As in the standard chameleon solution on the Minkowski background chame1 ; chame2 , the field $\phi$ and its $r$ derivative do not vanish outside the star. Hence the boundary conditions $\phi=0$ and ${\rm d}\phi/{\rm d}r=0$ at the surface of star ($r=r_{s}$) are not mandatory for the nearly massless scalar field at the distance $r>r_{s}$.

In this paper, we study the NS solutions in BD theories with the scalar potential and the general coupling $Q$ for two realistic EOSs: SLy SLY and FPS FPS . For this purpose, we use the analytic representations of these two EOSs presented in Ref. Haensel:2004nu . Our analysis is sufficiently general in that it covers massless BD theories and $f(R)$ gravity as special cases. In particular, we would like to clarify the difference of NS solutions between the potential $V(\phi)$ with a constant mass and the potential allowing the asymptotic behavior $V_{,\phi\phi}\equiv{\rm d^{2}}V/{\rm d}\phi^{2}\to 0$ at spatial infinity. For this purpose, we consider BD theories with the self-coupling potential $V(\phi)=\lambda\phi^{4}/4$ in the Jordan frame. The same potential was also introduced in the Einstein frame for accommodating the chameleon mechanism on the nonrelativistic background Gubser:2004uf . It is not yet clear whether the similar chameleon-like solutions arise on the relativistic background with realistic EOSs. We show the existence of new NS solutions for both SLy and FPS EOSs, despite a different property of the field profile from that on the weak gravitational background.

The NS solution with the potential $V(\phi)=\lambda\phi^{4}/4$ is in contrast to that in massive BD theories with $V(\phi)=m^{2}\phi^{2}/2$. While the scalar field in the latter case is subject to exponential growth by the constant mass $m$ outside the body, the former potential evades this problem due to the property that the mass squared $V_{,\phi\phi}$ approaches 0 as $r\to\infty$. We compute the mass $M$ and radius $r_{s}$ of NSs in BD theories with/without the potential $V(\phi)=\lambda\phi^{4}/4$ in order to see the signature for the modification of gravity from GR. For the purely massless case ($\lambda=0$), the modification to the radius $r_{s}$ tends to be significant for increasing $|Q|$ of order 0.1, while the maximum NS mass is hardly changed. As $\lambda$ increases, the theoretical curve of the mass-radius relation approaches that in GR by reflecting the fact that the field tends to be heavy inside the star.

This paper is organized as follows. In Sec. II, we derive the full equations of motion on the spherically symmetric and static background in BD theories with the scalar potential $V(\phi)$ in both Jordan and Einstein frames. We also discuss the boundary conditions at the center of star and at spatial infinity. In Sec. III, we study the mass-radius relation of NSs in BD theories with $V(\phi)=0$ and investigate how much modification from GR arises for different couplings $Q$. In Sec. IV, we consider the massive potential $V(\phi)=m^{2}\phi^{2}/2$ and show the difficulty of obtaining NS solutions consistent with the boundary conditions at spatial inifnity. In Sec. V, we investigate how the NS solutions can be realized by the self-coupling potential $V(\phi)=\lambda\phi^{4}/4$ and compare the mass-radius relation with that in GR. Sec. VI is devoted to conclusions.

In this paper we adopt the natural units $c=\hbar=1$, where $c$ is the speed of light and $\hbar$ is reduced Planck constant. When these fundamental constants are needed in numerical computations, we recover them and use their concrete values $c=2.9979\times 10^{10}~{}{\rm cm\cdot s^{-1}}$ and $\hbar=1.0546\times 10^{-27}~{}{\rm erg\cdot s}$, together with the Newton gravitational constant $G=6.6743\times 10^{-8}~{}{\rm g^{-1}\cdot cm^{3}\cdot s^{-2}}$.

II Equations of motion

We begin with the action of scalar-tensor theories accommodating BD theories with a scalar potential $V(\phi)$,

[TABLE]

where $g$ is the determinant of metric tensor $g_{\mu\nu}$, $F(\phi)$ is a function of a scalar field $\phi$, $R$ is the Ricci scalar, $Q$ is a constant, $X=-g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi/2$, and ${\cal L}_{m}$ is the action of matter fields $\Psi_{m}$. In Ref. Yoko , it was shown that BD theories Brans with the potential corresponds to the nonminimal coupling:

[TABLE]

In the limit that $Q\to 0$, the action (1) reduces to that of a canonical scalar field. The constant $Q$ characterizes the coupling between the field $\phi$ and the gravity sector. This coupling constant is related to the BD parameter $\omega_{\rm BD}$ as $2Q^{2}=1/(3+2\omega_{\rm BD})$ Yoko . The matter energy-momentum tensor is defined by $T_{\mu\nu}=-(2/\sqrt{-g})\delta{\cal L}_{m}/\delta g^{\mu\nu}$. Assuming that the matter fields are minimally coupled to gravity, there is the continuity equation

[TABLE]

where $\nabla^{\mu}$ is the covariant derivative operator.

The metric $f(R)$ gravity given by the action

[TABLE]

belongs to a subclass of the graviton-scalar action in Eq. (1), with the correspondence fRreview

[TABLE]

For $f(R)$ containing nonlinear functions in $R$, the scalar degree of freedom $\phi$ arises from the gravity sector. In this case the field potential $V(\phi)$ does not vanish, so the field $\phi$ generally has a nonvanishing effective mass.

In string theory, the low-energy effective action contains the so-called dilaton field $\Phi$ coupled to gravity Gas1 ; Gas2 . The lowest-order graviton-dilaton action in 4-dimensional spacetime takes the form

[TABLE]

After the field redefinition $\Phi\to 2Q\phi$, the action (6) reduces to the graviton-scalar action in Eq. (1), with the correspondences $Q^{2}=1/2$ and $V(\phi)=0$. Thus, the massless dilaton can be also accommodated in BD theory with the specific coupling $Q^{2}=1/2$.

In this paper, we will consider BD theories with general couplings $Q$ including metric $f(R)$ theories and dilaton gravity in the absence/presence of $V(\phi)$. We deal with the Jordan frame given by the action (1) as a physical frame and derive the equations of motion on the spherically symmetric and static background. We also discuss the field configuration in the Einstein frame in which the matter sector is directly coupled to $\phi$.

II.1 Jordan frame

We study the NS solutions on the spherically symmetric and static background given by the line element

[TABLE]

where $f$ and $h$ are functions of the distance $r$ from the center of symmetry. For the matter sector, we consider a perfect fluid whose energy-momentum tensor is given by $T^{\mu}_{\nu}={\rm diag}\,(-\rho(r),P(r),P(r),P(r))$, where $\rho(r)$ is the energy density and $P(r)$ is the pressure. From the continuity equation (3), we obtain

[TABLE]

where a prime represents the derivative with respect to $r$. To relate $P$ with $\rho$ for realistic NSs, we resort to the analytic representations of SLy and FPS EOSs Haensel:2004nu . Introducing the notations

[TABLE]

the two EOSs can be parameterized as

[TABLE]

where

[TABLE]

and the coefficients $a_{1,\cdots,18}$ for the SLy and FPS are given in Table 1 of Ref. Haensel:2004nu . Taking account of additional functions to Eq. (10), it is also possible to accommodate other EOSs like BSk19, BSk20, and BSk21 Potekhin:2013qqa .

Instead of the metric $h$, it is convenient to introduce the mass function ${\cal M}(r)$ defined by

[TABLE]

where the gravitational constant $G$ is related to $M_{\rm pl}$ as $G=(8\pi M_{\rm pl}^{2})^{-1}$. We define the ADM mass $M$ of the star, as

[TABLE]

The star radius $r_{s}$ is determined by the condition

[TABLE]

Varying the action (1) with respect to $g_{\mu\nu}$ and $\phi$, the equations of motion on the background (7) read

[TABLE]

where $V_{,\phi}\equiv{\rm d}V/{\rm d}\phi$. On using these equations, the Ricci scalar is expressed as

[TABLE]

which shows that not only the matter density and pressure but also the field kinetic energy and potential generally contribute to $R$.

The regularities of solutions at the center of NSs demands the following boundary conditions

[TABLE]

As long as the mass function has the dependence ${\cal M}(r)\propto r^{3}$ at leading order, we also have $h(r=0)=1$ from Eq. (12). For the consistency with Eq. (19), we expand $f,h,\phi,\rho$ around $r=0$, as

[TABLE]

where $f_{0},f_{n},h_{n},\phi_{0},\phi_{n},\rho_{0},\rho_{n}$ are constants. On using Eq. (10), the pressure can be written in the form $P(r)=P_{0}+\sum_{n=2}^{\infty}P_{n}r^{n}$, where $P_{0},P_{n}$ are constants. We also expand the potential in terms of the Taylor series, as

[TABLE]

Then, the solutions consistent with Eqs. (8) and (15)-(17) around $r=0$ are given by

[TABLE]

Both the coupling $Q$ and the $\phi$ derivative of $V(\phi)$ lead to the variation of $\phi$ around the center of body. They also give rise to modifications to $f$, $h$, $P$ in comparison to the theories with $Q=0$ and $V(\phi)=0$.

The asymptotic flatness at spatial infinity requires that

[TABLE]

where $\phi_{\infty}\equiv\phi(r\to\infty)$. For the power-law potential $V(\phi)=\lambda_{n}\phi^{n}$, the field value $\phi_{\infty}$ is equivalent to 0. In this case, the nonminimal coupling (2) approaches the value 1 of GR in the limit $r\to\infty$. For the massless scalar field without the potential, we impose the boundary condition $\phi_{\infty}=0$ besides the first three of (26). Since only the ratio between $f^{\prime}$ and $f$ appears in Eqs. (8) and (15)-(17), the constant $f_{0}$ in the expansion of Eq. (22) can be chosen as any arbitrary constant. The asymptotic value of $f$ at spatial infinity is generally different from 1, but it can be shifted to 1 by the time reparametrization. For $M={\rm constant}$, Eq. (12) shows that the function $h$ approaches 1 as $r\to\infty$. The field value $\phi_{0}$ at $r=0$ can be determined by a shooting method to satisfy the boundary conditions (26) at spatial infinity.

For the numerical purpose, we introduce the density $\tilde{\rho}_{0}$ and the distance $r_{0}$, as

[TABLE]

where $m_{n}=1.6749\times 10^{-24}$ g is the neutron mass and $n_{0}=0.1~{}{\rm(fm)}^{-3}$ is the typical number density of NSs. It is convenient to define the following dimensionless variables

[TABLE]

Then, the quantities $\xi$ and $\zeta$ in Eq. (9) are expressed, respectively, as

[TABLE]

where $\alpha_{1}=\ln(\tilde{\rho}_{0}/{\rm g\cdot cm^{-3}})/\ln 10$, $\alpha_{2}=(\ln 10)^{-1}$, and $\alpha_{3}=\ln(\tilde{\rho}_{0}\,c^{2}/{\rm dyn\cdot cm^{-2}})/\ln 10$. Then, the EOS translates to the form

[TABLE]

where $\zeta(\xi)$ is the function on the right hand side of Eq. (10). From the continuity Eq. (8), the derivative $y_{,s}\equiv{\rm d}y/{\rm d}s$ is expressed as

[TABLE]

From Eqs. (15)-(17), we have

[TABLE]

with

[TABLE]

Solving Eqs. (31)-(36) with the boundary conditions (22)-(25) outwards, we know the values of $y$, $z$, $f$, $m$, and $\varphi$ inside the star. Outside the star, we can simply set $y=0=z$ and solve Eqs. (33)-(36) for $f$, $m$, and $\varphi$. Defining $m_{\infty}\equiv m(r\to\infty)$, the ADM mass $M$ of star can be computed as

[TABLE]

where $M_{\odot}=1.9884\times 10^{33}$ g is the solar mass.

II.2 Einstein frame

Under the so-called conformal transformation

[TABLE]

the action (1) can be transformed to that in the Einstein frame without the nonminimal coupling Fujii . Here and in the following, we use the roman subscript “${\rm E}$” to represent quantities in the Einstein frame. The Ricci scalars in two frames are related to each other, as $R=(R_{{\rm E}}-6g^{\mu\nu}_{{\rm E}}\partial_{\mu}\omega\partial_{\nu}\omega+6\square_{{\rm E}}\,\omega)F$, where $\omega=(1/2)\ln F=-Q\phi/M_{\rm pl}$ and $\square_{{\rm E}}\omega=(1/\sqrt{-g_{{\rm E}}})\partial_{\mu}(\sqrt{-g_{{\rm E}}}\,g^{\mu\nu}_{{\rm E}}\partial_{\nu}\omega)$. On using the property $\sqrt{-g}=F^{-2}\sqrt{-g_{{\rm E}}}$ and dropping a boundary term associated with $\square_{{\rm E}}\omega$, the action (1) reduces to

[TABLE]

where

[TABLE]

From Eq. (39), it is clear that the canonical scalar field $\phi$ is directly coupled to matter fields in the Einstein frame. The matter energy-momentum tensor in the Einstein frame, which is defined by $(T_{\mu\nu})_{{\rm E}}=-(2/\sqrt{-g_{{\rm E}}})\delta{\cal L}_{m}/\delta g^{\mu\nu}_{{\rm E}}$, is related to that in the Jordan frame, as $(T_{\mu\nu})_{{\rm E}}=T_{\mu\nu}/F$. The energy density $\rho_{{\rm E}}$ and the pressure $P_{{\rm E}}$ of perfect fluids in the Einstein frame is given by $(T^{\mu}_{\nu})_{{\rm E}}={\rm diag}(-\rho_{{\rm E}},P_{{\rm E}},P_{{\rm E}},P_{{\rm E}})$, so there are the following relations

[TABLE]

Varying the action (39) with respect to $\phi$ and using the relation $\partial{\cal L}_{m}/\partial\phi=-\sqrt{-g_{{\rm E}}}\,T_{{\rm E}}\,F_{,\phi}/(2F)$, where $T_{{\rm E}}=-\rho_{{\rm E}}+3P_{{\rm E}}$ is the trace of energy-momentum tensor, it follows that

[TABLE]

This shows that the matter coupling $Q$ modifies the dynamics of $\phi$.

In the Einstein frame, we consider the spherically symmetric and static background given by the line element

[TABLE]

Since ${\rm d}s_{{\rm E}}^{2}=F{\rm d}s^{2}$, the distance $r$ and the metrics $f,h$ in the Jordan frame are related to those in the Einstein frame, as

[TABLE]

We introduce the mass function ${\cal M}_{{\rm E}}(r_{{\rm E}})$ in the Einstein frame, as

[TABLE]

together with the asymptotic mass

[TABLE]

In the Einstein frame, the star radius $r_{s}$ corresponds to

[TABLE]

where $\phi_{s}$ is the field value at the surface of star.

On using the correspondence (46) with Eqs. (12) and (47), it follows that

[TABLE]

The existence of terms $e^{Q\phi/M_{\rm pl}}$ and ${\rm d}\phi/{\rm d}r_{{\rm E}}$ lead to the difference between ${\cal M}(r)$ and ${\cal M}_{{\rm E}}(r_{{\rm E}})$. If ${\rm d}\phi/{\rm d}r_{{\rm E}}$ decreases faster than $1/r_{{\rm E}}^{2}$ and $\phi$ approaches 0 as $r_{{\rm E}}\to\infty$, then we have $M=M_{{\rm E}}$. On the other hand, if the field at large distances has the radial dependence

[TABLE]

where $\alpha$ is a constant, it follows that

[TABLE]

Even when $\phi_{\infty}=0$, the nonvanishing radial derivative (51) leads to the difference between $M$ and $M_{{\rm E}}$. As we will discuss in Sec. III, this difference appears for BD theories with $V(\phi)=0$.

In the Einstein frame, the matter continuity Eq. (8) reads

[TABLE]

while the scalar-field Eq. (42) reduces to

[TABLE]

Varying the action (39) with respect to $(g_{\mu\nu})_{{\rm E}}$, the metric $f_{{\rm E}}$ and the mass function ${\cal M}_{{\rm E}}$ obey

[TABLE]

which show that the field potential $V_{{\rm E}}$ and the kinetic energy $({\rm d}\phi/{\rm d}r_{{\rm E}})^{2}$ modify the values of $f_{{\rm E}}$ and ${\cal M}_{{\rm E}}$ in GR. From Eqs. (53) and (54), we find that the field $\phi$ and matter interact with each other through the coupling $Q$. While the equations of motion in the Einstein frame are simpler than those in the Jordan frame, the EOS (10) needs to be transformed to the relation between $P_{{\rm E}}$ and $\rho_{{\rm E}}$. The boundary conditions at $r_{{\rm E}}=0$ and $r_{{\rm E}}\to\infty$ are similar to those in the Jordan frame, i.e.,

[TABLE]

and

[TABLE]

The analytic solutions to $f_{{\rm E}},h_{{\rm E}},\phi,P_{{\rm E}}$ expanded around $r_{{\rm E}}=0$ can be obtained in a similar way to those derived in Eqs. (22)-(25) in the Jordan frame. The resulting solutions consistent with Eqs. (53)-(56) are given by

[TABLE]

where $f_{{\rm E}0},\rho_{{\rm E}0},P_{{\rm E}0}$ are constants corresponding to $f_{0},\rho_{0},P_{0}$ in Eq. (20), respectively. The field value $\phi_{0}$ is iteratively known to satisfy the boundary conditions (58) at spatial infinity. Solving Eqs. (53)-(56) numerically and transforming the solutions back to the Jordan frame, the physical observables like $r_{s}$ and $M$ should coincide with those computed directly in the Jordan frame for given model parameters. We will address this issue in Sec. III.

III Massless Brans-Dicke theories

We first consider massless BD theories without the scalar potential, i.e.,

[TABLE]

In this case, we can set $v=0$ and $v_{,\varphi}=0$ in Eqs. (33)-(35). We recall that the EOS is written as the form (32), with $\zeta(\xi)$ and $z$ given by Eqs. (10) and (31), respectively. Numerically, we solve Eqs. (33)-(35) with Eqs. (31), (32), and (36) by using the boundary conditions (22)-(25) around $r=0$.

From Eqs. (24) and (25), the scalar field and pressure around $r=0$ are given, respectively, by

[TABLE]

If $Q=0$, then the scalar field stays constant ($\phi=\phi_{0}$) around the center of star. Indeed, this can be also confirmed by the field Eq. (54) in the Einstein frame. For $Q=0$, the general solution to Eq. (54) is expressed in the integrated form

[TABLE]

where $\alpha$ is a constant. To avoid the divergence of the integral in Eq. (66) around $r_{{\rm E}}=0$, we require that $\alpha=0$ and hence $\phi$ is constant at any radial distance $r_{{\rm E}}$.

For $Q\neq 0$, the scalar field varies with the increase of $r$. In the following, we consider the negative coupling

[TABLE]

without loss of generality. We also focus on the coupling in the range $Q^{2}\leq 1/2$, under which the pressure (65) decreases with the growth of $r$. Our analysis covers dilaton gravity ($Q=-1/\sqrt{2}$) as a special case. The massless BD theory with $Q=-1/\sqrt{6}$, which we discuss in this section, is different from $f(R)$ gravity, in that the latter contains the nonvanishing potential $V(\phi)$. We will study the BD theory with $V(\phi)\neq 0$ in Secs. IV and V.

Under the condition $\rho_{0}>3P_{0}$, Eq. (64) shows that $\phi^{\prime}(r)<0$ around $r=0$. In the left panel of Fig. 1, we plot $\phi(r)$, $-\phi^{\prime}(r)$, and ${\cal M}(r)$ versus $r/r_{0}$ for the SLy EOS with the central density $\rho_{0}$ satisfying the condition $Q(\rho_{0}-3P_{0})<0$. In this case, $-\phi^{\prime}(r)$ is positive around $r=0$ and it linearly grows as $-\phi^{\prime}(r)\propto r$ for $r\lesssim 0.01r_{0}$. In this regime, the field $\phi$ slowly decreases from the central value $\phi_{0}\simeq 0.203M_{\rm pl}$ according to Eq. (64). The field derivative $-\phi^{\prime}(r)$ reaches the maximum around the surface of star ($r_{s}\simeq 0.13r_{0}$). In the left panel of Fig. 1, we observe that the large variation of $\phi$ starts to occur around $r=r_{s}$. In the right panel, we also find that both $\rho$ and $P$ rapidly drop down for $r>0.1r_{0}$.

Outside the star, both $\phi(r)$ and $-\phi^{\prime}(r)$ decrease with the increase of $r$. Since $\rho_{{\rm E}}=0=P_{{\rm E}}$ for $r>r_{s}$, the solution to Eq. (54) is the same as Eq. (66). For the large distance $r\gg r_{s}$, both $f_{{\rm E}}$ and $h_{{\rm E}}$ approach 1 with $\phi_{\infty}\to 0$ and $r_{{\rm E}}\to r$, so the field derivative is given by

[TABLE]

In Fig. 1, we can confirm that $-\phi^{\prime}(r)$ decreases in proportion to $1/r^{2}$ for $r\gg r_{s}$, with $\phi(r)$ approaching 0 at spatial infinity. For the massless scalar, the field value $\phi_{0}$ at $r=0$ satisfying the boundary condition $\phi_{\infty}=0$ can be identified in the following way. First, we perform the numerical integration by choosing $\phi_{0}=0$ and then find the asymptotic value $\phi_{\rm asy}$ at $r\gg r_{s}$ (say, at $r=10^{20}r_{s}$). Then, we run the code again with the value $\phi_{0}=-\phi_{\rm asy}$ at $r=0$. This second run leads to the asymptotic value $\phi_{\infty}$ converging to 0. It is important to identify the appropriate value of $\phi_{0}$ in this way because the mass $M$ and radius $r_{s}$ are affected by the nonminimal coupling term $e^{-2Q\phi_{0}/M_{\rm pl}}$. In the numerical simulation of Fig. 1, the mass function ${\cal M}$ quickly approaches the asymptotic value $M=1.83M_{\odot}$ for $r>r_{s}=0.13r_{0}=11.7$ km.

In Fig. 2, we show the mass $M$ in the Jordan frame versus the radius $r_{s}$ for the SLy (left) and FPS (right) EOSs with the coupling $Q=-1/\sqrt{6}$. As the central density $\rho_{0}$ grows from the value of order $10^{14}$ g $\cdot$ cm*-3*, the mass $M$ increases by reaching a maximum $M_{\rm max}$. For SLy, we have $M_{\rm max}\simeq 2.05M_{\odot}$ with the radius $r_{s}\simeq 10.5$ km around the density $\rho_{0}=2.7\times 10^{15}$ g $\cdot$ cm*-3*. For this maximum mass the condition $\rho_{0}>3P_{0}$ holds, but as $\rho_{0}$ increases further, the system eventually enters the region with the fully relativistic EOS satisfying $\rho_{0}<3P_{0}$. With the growth of $\rho_{0}$ in the fully relativistic region, the radius $r_{s}$ gets smaller, by reflecting the fact that the second term on the right hand side of Eq. (65) increases. This is also accompanied by the decrease of $M$. For $Q=-1/\sqrt{6}$, the maximum mass exceeds $2M_{\odot}$ with the SLy EOS, while, the FPS EOS gives rise to the maximum mass $M_{\rm max}\simeq 1.80M_{\odot}$ with $r_{s}\simeq$ 9.67 km and $\rho_{0}=3.5\times 10^{15}$ g $\cdot$ cm*-3*.

In Fig. 2, we also plot the mass $M_{{\rm E}}$ computed in the Einstein frame. As we estimated in Eq. (52), there is the relation $M=M_{{\rm E}}-8\pi M_{\rm pl}Q\alpha$ for $\phi_{\infty}=0$, where $\alpha$ corresponds to the coefficient in Eq. (68). The constant $\alpha$ is related to the field derivative $\phi^{\prime}(r)$ at $r=r_{s}$. If we extrapolate the solution (64) up to the radius of star, it follows that $\phi^{\prime}(r_{s})\approx Q(\rho_{0}-3P_{0})r_{s}/(3M_{\rm pl}e^{-2Q\phi_{0}/M_{\rm pl}})$. If the exterior solution (68) at spatial infinity is also extrapolated down to $r=r_{s}$, then the coefficient $\alpha$ can be estimated as $\alpha\approx Q(\rho_{0}-3P_{0})r_{s}^{3}/(3M_{\rm pl}e^{-2Q\phi_{0}/M_{\rm pl}})$. Although this is a crude estimation under which the coefficient $\alpha$ is inaccurate, we may generally express $\alpha$ in the form

[TABLE]

where $\beta$ is a constant at most of order 1. In this case, the two masses $M$ and $M_{{\rm E}}$ are related to each other, as

[TABLE]

where $M_{0}\equiv 4\pi r_{s}^{3}\rho_{0}/3$ corresponds to the mass of star with the constant density $\rho_{0}$. The difference between $M$ and $M_{{\rm E}}$ arises from the nonvanishing coupling $Q$. For $\beta>0$, $M$ is smaller than $M_{{\rm E}}$ with the difference $6\beta Q^{2}M_{0}$. In the numerical simulation of Fig. 2 ($Q=-1/\sqrt{6}$), even when $\rho_{0}-3P_{0}$ is negative in the large $\rho_{0}$ region, the constant $\beta$ in Eq. (69) is positive outside the star and hence $M<M_{{\rm E}}$. We also compute the mass $M$ from $M_{{\rm E}}$ by using the transformation relation (50). As we observe in Fig. 2 (bold dashed line), the mass $M$ obtained from the Einstein-frame mass $M_{{\rm E}}$ exactly coincides with the one directly computed in the Jordan frame. This shows the consistency of our calculations in both Jordan and Einstein frames.

In Fig. 3, we plot the mass $M$ versus $r_{s}$ for the SLy and FPS EOSs with four different values of $Q$. The solid line corresponds to $Q=0$, i.e., the massless scalar field in GR. As we see in case (b), the modification to $M$ and $r_{s}$ induced by the coupling $Q$ is small for $|Q|\leq 0.1$, but the difference from the $Q=0$ case arises for $|Q|>0.1$. The change of $r_{s}$ is particularly significant for large $|Q|$, like cases (c) and (d) in Fig. 3. The mass $M$ is also subject to modifications by the nonvanishing $Q$, but the maximum mass $M_{\rm max}$ does not exceed the corresponding value in GR for both SLy and FPS EOSs. The change of $r_{s}$ induced by the coupling in the range $|Q|~{}\gtrsim 0.1$ is the main signature of distinguishing between massless BD theories and GR from observations.

IV Brans-Dicke theories with constant scalar mass

In this section, we study NS solutions in BD theories with a constant scalar mass $m$. This is characterized by the potential

[TABLE]

in the Jordan frame.

The Starobinsky $f(R)$ model given by

[TABLE]

falls in this category in the regime $|R|\ll m^{2}$. To see this, we use the fact that the scalar degree of freedom $\phi$ is related to the Ricci scalar $R$, as Eq. (5), with $Q=-1/\sqrt{6}$. In the Starobinsky model, there is the correspondence

[TABLE]

In this case, the scalar potentials in Jordan and Einstein frames are given, respectively, by

[TABLE]

both of which vanish at $\phi=0$. Expanding $V(\phi)$ and $V_{{\rm E}}(\phi)$ around $\phi=0$ in the regime $|\phi|\ll M_{\rm pl}$, it follows that $V(\phi)\simeq V_{{\rm E}}(\phi)\simeq m^{2}\phi^{2}/2$. Thus, the scalar field has a constant mass $m$ around the potential minimum. In this regime, we have $R\simeq\sqrt{6}m^{2}\phi/M_{\rm pl}$ from Eq. (73) and hence $|R|\ll m^{2}$. For $|\phi|\lesssim M_{\rm pl}$, the Einstein-Hilbert term $R$ dominates over $R^{2}/(6m^{2})$ in the Lagrangian (72). In Refs. Cooney:2009rr ; Arapoglu:2010rz ; Orellana:2013gn ; Astashenok:2013vza ; Ganguly:2013taa ; Yazadjiev ; Capozziello:2015yza ; Resco:2016upv , the NS solutions in the model (72) were studied in both Jordan and Einstein frames.

The parameter $a=1/(6m^{2})$ in the Starobinsky $f(R)$ model is constrained to be $a<5\times 10^{11}$ m2 from the binary pulsar data Naf:2010zy . This translates to the bound

[TABLE]

For general couplings $Q$, the bound on the mass $m$ is subject to modifications. As long as $|Q|={\cal O}(0.1)$, the order of the upper bound on $m^{-1}$ should be similar to that of Eq. (75).

The potential in the Einstein frame corresponding to Eq. (71) is given by

[TABLE]

Outside the star, the field Eq. (54) obeys

[TABLE]

For the potential (76), the boundary conditions of $\phi$ at spatial infinity correspond to

[TABLE]

In the asymptotic regime characterized by $r_{{\rm E}}\gg r_{s}$, we can employ the approximations that both $f_{\rm E}$ and $h_{\rm E}$ are close to 1 with $|Q\phi/M_{\rm pl}|\ll 1$. Then, the field Eq. (77) outside the star approximately reduces to

[TABLE]

This has the following solution

[TABLE]

where $c_{1}$ and $c_{2}$ are integration constants. The constant $c_{1}$ should vanish to satisfy the boundary conditions (78). In the exterior region of star close to its surface, the solution (80) is subject to modifications. Nevertheless, provided that $1+2Q\phi/M_{\rm pl}>0$, the last term on the left hand side of Eq. (77) leads to the growth of $|\phi|$ for the boundary conditions satisfying $\phi\neq 0$ or ${\rm d}\phi/{\rm d}r_{{\rm E}}\neq 0$ at the surface of star. In other words, we require that

[TABLE]

to avoid the increase of $|\phi|$ outside the body. Unless the boundary conditions (81) are satisfied, the exponential growth of $|\phi|$ starts to occur for the distance $r_{{\rm E}}\gtrsim 1/m\equiv r_{c}$. Under the bound (75), the critical distance $r_{c}$ corresponds to $r_{c}=19.3r_{0}=1.73\times 10^{3}$ km, which is about $10^{2}$ times as large as the typical radius of NSs ($\sim 10~{}$km). If $1+2Q\phi/M_{\rm pl}<0$ at some distance $r$, then the field $\phi$ exhibits damped oscillations with the decreasing amplitude ($|\phi|\propto 1/r_{{\rm E}}$). Then, the system enters the regime in which the condition $1+2Q\phi/M_{\rm pl}>0$ is satisfied, so the scalar field is eventually subject to exponential growth for the distance $r_{{\rm E}}\gtrsim 1/m$.

The boundary conditions (81) imply that the field does not contribute to the solution outside the body. Setting $\phi(r)=0$ and $\phi^{\prime}(r)=0$ in Eq. (18) for the vacuum exterior ($r\geq r_{s}$), we have $R=0$ and hence the external region of star corresponds to the Schwarzschild geometry. This fact was recognized in Ref. Ganguly:2013taa for the Starobinsky $f(R)$ model. From Eq. (18), there are the following particular relations in $f(R)$ gravity ($Q=-1/\sqrt{6}$):

[TABLE]

For the quadratic potential (71) and the Starobinsky potential given in Eq. (74), the boundary conditions (81) translate to

[TABLE]

which coincide with those derived in Ref. Ganguly:2013taa by using the junction conditions of $f(R)$ gravity Deruelle:2007pt .

Around the center of star, the scalar field $\phi$ has the radial dependence (61) in the Einstein frame. In comparison to massless BD theories studied in Sec. III, the potential-dependent term $V_{{{\rm E}},\phi}(\phi_{0})$ leads to the additional variation of $\phi$ around $r_{{\rm E}}=0$. We note that the mass squared $m^{2}$ explicitly appears as one of the coefficients of $r_{{\rm E}}^{4}$ in the expansion of $\phi(r_{{\rm E}})$ in Eq. (61). For the boundary conditions where the combination $M_{\rm pl}V_{{{\rm E}},\phi}(\phi_{0})+Q(\rho_{{\rm E}0}-3P_{{\rm E}0})$ is close to 0, the nonvanishing effective mass around the potential minimum leads to the variation of $\phi$ with respect to $r$. Indeed, this is the case for the chameleon scalar field where the variation of $\phi$ occurs mostly around the surface of body (“thin shell”) chame1 ; chame2 ; Tsujikawa:2009yf .

For the massive potential (71), the field value $\phi_{0}$ at the center of NSs needs to be chosen to satisfy the boundary conditions (81) at $r=r_{s}$. The EOS of NSs affects the field profile $\phi(r)$ not only around $r=0$ but also in the whole interior region of star ($0\leq r\leq r_{s}$). We numerically solve Eqs. (31)-(36) from the vicinity of $r=0$ for both SLy and FPS EOSs with general couplings $Q$. For the central matter density in the range $\rho_{0}>10^{14}$ g $\cdot$ cm*-3*, we could not find appropriate values of $\phi_{0}$ at $r=0$ avoiding the exponential growth outside the body. This means that the boundary conditions (78) are not consistently satisfied at spatial infinity for SLy and FPS EOSs.

Even if we try to fine-tune the field value $\phi_{0}$ to match with the Schwarzschild exterior for some other EOSs, it is difficult to realize the boundary conditions of $\phi$ and ${\rm d}\phi/{\rm d}r$ which exactly vanish at $r=r_{s}$. They are highly sensitive to a slight change of the EOS. Even when we find a value of $\phi_{0}$ compatible with the conditions (81) for a particular EOS, the same property no longer holds under a tiny change of the EOS. Since the EOSs are determined by the nuclear reaction inside NSs, a particular EOS chosen to match with the Schwarzschild exterior does not generally correspond to the realistic physical EOS Ganguly:2013taa . For the potential (71) with SLy and FPS EOSs, we did not numerically find regular solutions of the field and metrics consistent with all the boundary conditions at $r=0,r_{s},\infty$.

The above discussion is based on the scalar potential (71), but the same conclusion also persists for the Starobinsky $f(R)$ model given by the Lagrangian (72). In the Starobinsky model we require that $R/(6m^{2})\to 0$ as $r\to\infty$ for the asymptotic flatness, so the field needs to enter the region in which $|\phi/M_{\rm pl}|$ is smaller than 1 at some distance. In this regime the scalar potential reduces to $V(\phi)\simeq m^{2}\phi^{2}/2$, so the field is eventually subject to exponential growth around the distance $r\gtrsim 1/m$ for nonvanishing $\phi(r)$ or $\phi^{\prime}(r)$ outside the star. We performed numerical simulations in the Starobinsky model by varying $\phi_{0}$ at $r=0$ and did not find consistent solutions satisfying all the boundary conditions discussed above for both SLy and FPS EOSs. The analysis of Ref. Ganguly:2013taa based on the polytropic EOS $\rho=\kappa P^{9/5}$ also reached the same conclusion.

In Refs Arapoglu:2010rz ; Orellana:2013gn ; Ganguly:2013taa , the authors studied NS solutions in the $f(R)$ model $f(R)=R-aR^{2}$, where $a$ is a positive constant. In this case, the mass squared $m_{\phi}^{2}$ of the gravitational scalar field is negative ($m_{\phi}^{2}\simeq 1/(3f_{,RR})=-1/(6a)$). Then the field $\phi$ exhibits damped oscillations at large distances, so the exponential increase of $|\phi|$ can be avoided at the background level. However, the linearly perturbed version of Eq. (42) shows that the dynamical equation of motion for the field perturbation $\delta\phi$ (i.e., the equation associated with the second time derivative of $\delta\phi$) contains the negative mass squared, which induces the instability of perturbations. This property holds not only for the $f(R)$ model $f(R)=R-aR^{2}$ with $a>0$ but also for BD theories with the tachyonic mass squared.

In summary, we showed that BD theories with the positive constant mass squared $m^{2}$ generally face the problem of realizing stable NS field configurations consistent with all the boundary conditions. Apart from Ref. Ganguly:2013taa , this fact was overlooked in most of the past works about NS solutions in $f(R)$ gravity with the positive constant $m^{2}$ outside the star. For the massless field studied in Sec. III, there is no need of satisfying the conditions (81) at $r=r_{s}$ due to the absence of the exponential growing term $e^{mr_{{\rm E}}}$ and hence the stable NS solution can be easily obtained for a given EOS.

V Brans-Dicke theories with self-coupling potential

For the NS solution discussed in Sec. IV, the field mass squared $m_{\phi}^{2}=V_{,\phi\phi}$ is a positive constant both inside and outside the star. Instead, we can consider other scalar potentials with the effective mass depending on the matter density. This is the case for a chameleon scalar field where the mass is large in the region of high density, whereas the field is light outside the compact object chame1 ; chame2 . For example, the $f(R)$ models of late-time cosmic acceleration fR1 ; fR2 ; fR3 ; fR4 are designed to have a heavy mass in large-curvature regimes to suppress the propagation of fifth forces around a compact body on the weak gravitational background, while the mass of gravitational scalar field is as light as today’s Hubble expansion rate $H_{0}$ fR1 ; Faulkner:2006ub ; Capozziello:2007eu . In such $f(R)$ models, the existence of relativistic stars was shown in Refs. Babichev:2009td ; Upadhye:2009kt by considering the constant-density profile or polytropic EOS.

The relativistic star can be also present in BD theories with the potential where $V_{,\phi\phi}$ depends on $\phi$. Indeed, the numerical simulation of Ref. Tsujikawa:2009yf confirmed the existence of relativistic stars with a constant-density profile for the inverse power-law potential $V(\phi)=M^{4+n}\phi^{-n}$ ($n>0$). This is analogous to the chameleon solution on the weak gravitational background, in that the scalar field is heavy inside the star and its variation occurs mostly around its surface (thin shell). In the exterior region, the scalar field becomes nearly massless due to the significant dropdown of matter density. In this case, it is possible to avoid the exponential increase of $\phi$ outside the star by sending $V_{,\phi\phi}$ to 0 as $r\to\infty$. Thus, we do not need to impose the boundary conditions (81) at $r=r_{s}$ for the models in which the scalar field is nearly massless outside the star.

In what follows, we study NS solutions in BD theories with the self-coupling potential

[TABLE]

where $\lambda$ is a positive constant. In this case, the second derivative of potential corresponding to the mass squared of scalar field is given by

[TABLE]

which depends on $\phi$. The $\phi$ dependence of $V_{,\phi\phi}$ allows a possibility for the chameleon mechanism to be at work. Indeed, this is the case for a compact body with the constant density in a weak-field limit Gubser:2004uf . We extend the analysis to the strong gravitational background by using the SLy and FPS EOSs.

If the chameleon mechanism works for a nonrelativistic compact body with the constant density $\rho_{{\rm E}}$, there is the balance $V_{{\rm E},\phi}+Q\rho_{{\rm E}}/M_{\rm pl}\simeq 0$ in Eq. (54) except for a thin-shell region near the surface. On the relativistic background, the pressure $P_{{\rm E}}$ also contributes to the field dynamics. Moreover, the realistic NSs have the $r$-dependent density and pressure. At spatial infinity, the scalar field needs to obey the boundary conditions (78). In general, the chameleon-like boundary conditions satisfying $V_{{\rm E},\phi}+Q(\rho_{{\rm E}}-3P_{{\rm E}})/M_{\rm pl}\simeq 0$ around the center of NSs do not give rise to the asymptotic solutions (78) on the relativistic background. In other words, we need to numerically identify the field value $\phi_{0}$ at $r=0$ consistent with the boundary conditions at spatial infinity.

For the potential (84), as long as $|\phi|$ gradually approaches 0 as $r\to\infty$, the second derivative of potential (85) also goes to 0. This is the way of avoiding the exponential growth of $\phi$ outside the body induced by the constant mass term $m$. Accordingly, it is not necessary to impose the boundary conditions (81) at the surface of NSs. Unlike the potential $V(\phi)=m^{2}\phi^{2}/2$, the absence of the conditions (81) does not restrict the forms of EOSs inside the star.

Using the dimensionless variables introduced in Eq. (29), the last term ${\cal F}_{\phi}\equiv[4QV+V_{,\phi}M_{\rm pl}+Q(\rho-3P)]/(M_{\rm pl}hF)$ in Eq. (17) can be expressed as

[TABLE]

where

[TABLE]

with $r_{0}M_{\rm pl}=1.107\times 10^{39}$. If $|Q\varphi|$ is smaller than the order 1, the scalar potential modifies the field dynamics for $|\tilde{\lambda}\varphi^{3}|\gtrsim|8\pi Q\left(y-3z\right)|$. Due to the largeness of the dimensionless quantity $r_{0}M_{\rm pl}$, even the values of $\lambda$ and $\varphi$ much smaller than 1 can give rise to the large contribution to Eq. (17). For example, we have $|\tilde{\lambda}\varphi^{3}|>1$ for $\lambda^{1/3}\varphi>10^{-26}$.

Numerically, we solve the background equations in the Jordan frame by randomly choosing the value of $\phi_{0}$ at $r=0$ and then identify its appropriate value consistent with the boundary conditions (78) by the shooting method. For the practical computation, we perform the integration up to the distance $r=10^{8}r_{0}$ by checking that both $\phi(r)$ and $\phi^{\prime}(r)$ sufficiently approach 0.

In Fig. 4, we plot $\phi_{0}$ versus the central density $\rho_{0}$ for $Q=-1/\sqrt{6}$ with four different values of $\tilde{\lambda}$. The EOS is chosen to be SLy (left) and FPS (right). The model with $\lambda\neq 0$ can be regarded as $f(R)$ gravity with the self-coupling potential (84). The vanishing self-coupling ($\tilde{\lambda}=0$) corresponds to the massless scalar field studied in Sec. III. For $\tilde{\lambda}\neq 0$ with a given central density $\rho_{0}$, there exists a unique value of $\phi_{0}$ consistent with the boundary conditions (78) for both SLy and FPS EOSs. When $\rho_{0}$ is of order $10^{14}$ g $\cdot$ cm*-3*, the condition $\rho_{{\rm E}}>3P_{{\rm E}}$ holds and hence the term $Q(\rho_{{\rm E}}-3P_{{\rm E}})/M_{\rm pl}$ in Eq. (54) is negative for $Q<0$. Provided that $\phi_{0}>0$, this matter-coupling term counteracts the increase of $\phi$ induced by the potential $V_{{\rm E},\phi}=3\lambda\phi^{2}$. In the full relativistic region where the opposite inequality $\rho_{{\rm E}}<3P_{{\rm E}}$ holds, the field value $\phi_{0}$ corresponding to the asymptotic solution (78) is typically negative. In Fig. 4, we observe that, for increasing $\tilde{\lambda}$, there is a tendency that $\phi_{0}$ decreases (apart from the density $\rho_{0}>4\times 10^{15}$ g $\cdot$ cm*-3* for the SLy EOS).

In Fig. 5, we show the second derivative of potential $\tilde{V}_{,\phi\phi}\equiv V_{,\phi\phi}r_{0}^{2}=3\lambda\phi^{2}r_{0}^{2}$ versus the distance $r$ for several different values of $\tilde{\lambda}$ with $Q=-1\sqrt{6}$. The field value at $r=0$ is chosen to be $\phi_{0}\simeq 5.6\times 10^{-2}M_{\rm pl}$ for SLy and $\phi_{0}\simeq 6.8\times 10^{-2}M_{\rm pl}$ for FPS. As $\tilde{\lambda}$ increases, $\tilde{V}_{,\phi\phi}$ tends to be larger. The variation of $\phi$ inside the NS is not significant except for the region around the surface of star. Still, the field profiles in the present model are different from that of the chameleon scalar satisfying the relation $V_{{\rm E},\phi}+Q(\rho_{{\rm E}}-3P_{{\rm E}})/M_{\rm pl}\simeq 0$ in most internal regions of a compact body. Outside the star, $\tilde{V}_{,\phi\phi}$ decreases toward the asymptotic value 0. As we observe in Fig. 5, the field $\phi$ does not vanish at the surface of star. Hence the boundary conditions (81) at $r=r_{s}$ do not hold for BD theories with the self-coupling potential. In other words, the field $\phi$ contributes to the geometry of the external region of star.

In Fig. 6, we plot the mass-radius relation of NSs for $Q=-1/\sqrt{6}$ with four different values of $\tilde{\lambda}$, together with the prediction of GR. While this is derived by the calculation in the Jordan frame, the same mass $M$ and radius $r_{s}$ can be obtained by integrating Eqs. (53)-(56) in the Einstein frame and transforming back to the Jordan frame. The case $\tilde{\lambda}=0$ corresponds to the massless scalar field plotted in Fig. 2 as the solid lines. If $\tilde{\lambda}$ is smaller than the order $10^{3}$, the values of $M$ and $r_{s}$ are similar to those in the massless case. For $\tilde{\lambda}\geq{\cal O}(10^{3})$, the difference from the $\tilde{\lambda}=0$ case starts to appear for both SLy and FPS EOSs.

As $\tilde{\lambda}$ increases further, the mass-radius relations tend to approach that in GR, by reflecting the fact that $V_{,\phi\phi}$ gets larger inside the star (see Fig. 5). For intermediate values of the self-coupling like $\tilde{\lambda}=10^{4}$, the maximum NS mass and the corresponding radius are slightly larger than those in GR. If we choose a larger coupling $|Q|$ than that in Fig. 6, the modification to $r_{s}$ tends to be more significant in comparison to the change of $M$. This is analogous to the mass-radius relation plotted in Fig. 3 for BD theories without the potential.

Finally, we should mention whether the potential (84) can follow from a specific model of $f(R)$ gravity. Let us consider the $f(R)$ Lagrangian

[TABLE]

where $a$ and $p$ are positive constants. In this case, the Jordan-frame potential (5) yields $V=M_{\rm pl}^{2}a(p-1)R^{p}/2$ with $R=[(e^{\sqrt{6}\phi/(3M_{\rm pl})}-1)/(ap)]^{1/(p-1)}$, i.e.,

[TABLE]

In the regime $|\phi|\ll M_{\rm pl}$, this potential approximately reduces to

[TABLE]

where $V_{0}$ is a positive constant. Hence the self-coupling potential (84) corresponds to the power

[TABLE]

For the potential (90), we have

[TABLE]

As long as the power $p$ is in the range

[TABLE]

$V_{,\phi\phi}$ is positive for $\phi>0$. Moreover, it has the asymptotic behavior $V_{,\phi\phi}\to 0$ as $\phi$ approaches 0 at spatial infinity. Then, as in the self-coupling potential (84), the $f(R)$ model (88) with $1<p<2$ can avoid the exponential growth of $\phi$ outside the star.

In summary, we presented models in BD theories with the potential (84) and in $f(R)$ gravity given by the Lagrangian (88), which can avoid the problem of exponential growth of $\phi$ outside the star. The existence of regular NS solutions was confirmed in the numerical simulation of Fig. 5. As $\lambda$ gets closer to 0, the mass-radius relation tends to deviate from that in GR as shown in Fig. 6.

VI Conclusions

In this paper, we studied NS solutions on the spherically symmetric and static background in $f(R)$ gravity and BD theories with/without the scalar potential $V(\phi)$. For this purpose, we used the SLy and FPS EOSs given by the parametrization (10) in the Jordan frame. In Sec. II, we obtained the full background equations in both Jordan and Einstein frames together with the solutions of metrics, field, and pressure expanded around the center of star. The explicit relation between the ADM masses in Jordan and Einstein frames is also derived in Eq. (50), which can be used for checking the consistency of calculations in two frames.

In Sec. III, we discussed NS solutions with the scalar field profile $\phi(r)$ in BD theories with $V(\phi)=0$. As we see in Eq. (64), the coupling $Q$ leads to the variation of $\phi(r)$ around $r=0$, with the growth $|\phi^{\prime}(r)|\propto r$. The field derivative $|\phi^{\prime}(r)|$ reaches a maximum value around the surface of star ($r=r_{s}$) and then it starts to decrease for $r>r_{s}$. For the large distance far away from the surface, $|\phi^{\prime}(r)|$ decreases in proportion to $1/r^{2}$. We found the way of identifying the field value $\phi_{0}$ at $r=0$ consistent with the boundary conditions $\phi(r)\to 0$ and $\phi^{\prime}(r)\to 0$ at spatial infinity. We numerically computed the mass-radius relation of NSs in the Jordan frame and showed that the calculation in the Einstein frame gives the same result after transforming back to the Jordan frame. The mass-radius relation exhibits the difference from that in GR for $|Q|>0.1$. As $|Q|$ increases, the radius $r_{s}$ is subject to large modifications, while the maximum mass reached in SLy and FPS EOSs is hardly changed in comparison to that in GR.

In Sec. IV, we considered BD theories with the quadratic potential $V(\phi)=m^{2}\phi^{2}/2$ and studied the effect of mass $m$ on the NS configuration. Far outside the surface of star, the field equation in the Einstein frame is of the form (79), which contains the growing-mode solution $e^{mr_{{\rm E}}}/r_{{\rm E}}$. To avoid the exponential increase of $|\phi|$ induced by the constant mass $m$, the scalar field is restricted to obey the boundary conditions (81) at $r=r_{s}$. This amounts to imposing the Schwarzschild geometry outside the star, without the scalar-field contribution to the metric. Such boundary conditions are not generally satisfied for arbitrary NS EOSs, in which case the field $\phi$ is subject to exponential growth for the distance $r\gtrsim 1/m$. We performed the numerical simulation for both SLy and FPS EOSs and did not find the NS configuration consistent with all the boundary conditions at $r=0,r_{s},\infty$. This is also the case for the Starobinsky $f(R)$ model given by Eq. (72).

In Sec. V, we extended the analysis to BD theories with the self-coupling potential $V(\phi)=\lambda\phi^{4}/4$. Since the second derivative of potential $V_{,\phi\phi}=3\lambda\phi^{2}$ goes to 0 for $\phi(r)$ approaching 0 at spatial infinity, it is possible to avoid the exponential growth of $\phi$ induced by the mass term without imposing the boundary conditions (81) at $r=r_{s}$. For given $\lambda$, $Q$, $\rho_{0}$, we identified the field value $\phi_{0}$ at $r=0$ leading to the appropriate boundary conditions $\phi(r)\to 0$ and $\phi^{\prime}(r)\to 0$ as $r\to\infty$. In general, we found that the chameleon-like boundary conditions at $r=0$ do not give rise to the appropriate NS solutions for both SLy and FPS EOSs and that, even for $\rho<3P$ around $r=0$, there are consistent NS configurations. We computed the mass-radius relation for $Q=-1/\sqrt{6}$ with several different values of $\lambda$ and showed that, with increasing $\lambda$, the theoretical curves tend to approach that of GR, see Fig. 6.

We have thus found the new type of NS solutions in BD theories in the presence of the self-coupling potential corresponding to $f(R)$ theories with the Lagrangian (88). Since the mass-radius relation is different from that in GR, the star’s compactness parameter $M/r_{s}$ is subject to modifications. This leads to the difference for the tidal Love number of NSs Flanagan:2007ix ; Damour:2009vw ; Binnington:2009bb ; Hinderer:2009ca , so that the deviation from GR can be potentially tested in the GW observations of NS coalescence like GW170817 TheLIGOScientific:2017qsa . It is of interest to clarify whether the same property holds in BD theories with a wide variety of potentials, e.g., the potential accommodating the chameleon mechanism proposed in the context of late-time cosmic acceleration fR1 ; fR2 ; fR3 ; fR4 , by using realistic EOSs, and put constraints on theories via tidal Love numbers of NSs. We leave the detailed analysis for the computation of tidal Love numbers in BD theories with the potential and $f(R)$ theories for a future work.

Acknowledgements.

RK is supported by the Grant-in-Aid for Young Scientists B of the JSPS No. 17K14297. ST is supported by the Grant-in-Aid for Scientific Research Fund of the JSPS No. 19K03854 and MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas “Cosmic Acceleration” (No. 15H05890).