Neutron Stars in Screened Modified Gravity: Chameleon vs Dilaton

TL;DR

This paper investigates how scalar fields in screened modified gravity models, specifically chameleon and dilaton types, influence neutron star structures and their observable mass-radius relationships.

Contribution

It provides a numerical analysis of scalar field profiles around neutron stars for chameleon and dilaton models, highlighting degeneracies with the equation of state.

Findings

Chameleon and dilaton behaviors are similar around neutron stars.

Large degeneracy exists between gravity parameters and the neutron star equation of state.

Modifications to the mass-radius relationship depend on model parameters.

Abstract

We consider the scalar field profile around relativistic compact objects such as neutron stars for a range of modified gravity models with screening mechanisms of the chameleon and Damour-Polyakov types. We focus primarily on inverse power law chameleons and the environmentally dependent dilaton as examples of both mechanisms. We discuss the modified Tolman-Oppenheimer-Volkoff equation and then implement a relaxation algorithm to solve for the scalar profiles numerically. We find that chameleons and dilatons behave in a similar manner and that there is a large degeneracy between the modified gravity parameters and the neutron star equation of state. This is exemplified by the modifications to the mass-radius relationship for a variety of model parameters.

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{M_{\mathrm{Pl}}^2}{2}R-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V(\varphi )\right]+S_m\left[{\Psi }_i,A^2(\varphi )g_{\mu \nu }\right].$$

$$A(\varphi )=1+{\int }^{\varphi }d\tilde{\varphi }\frac{\beta (\tilde{\varphi })}{M_{\mathrm{Pl}}}$$

$$G^{\mu \nu }=\frac{1}{M_{\text{pl}}^2}\left(T_m^{\mu \nu }+T_{\varphi }^{\mu \nu }\right)$$

$$T_{\varphi }^{\mu \nu }={\nabla }^{\mu }\varphi {\nabla }^{\nu }\varphi -g^{\mu \nu }\left(\frac{1}{2}{g}^{\alpha \beta }{\nabla }_{\alpha }\varphi {\nabla }_{\beta }\varphi +V(\varphi )\right)$$

$$\square \varphi =\frac{\partial V}{\partial \varphi }-\frac{\partial \ln A}{\partial \varphi }{T}_m.$$

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

$$T_{\nu }^{\mu }=A^4{\tilde{T}}_{\nu }^{\mu }$$

$$\rho =A^4\tilde{\rho },P=A^4\tilde{P}$$

$$V_{\mathrm{eff}}=V+\frac{1}{4}{A}^4(\tilde{\rho }-3\tilde{P})$$

$$m_{\mathrm{eff}}^2(\tilde{\rho },\tilde{P})\equiv \frac{d^2{V}_{\mathrm{eff}}}{d{\varphi }^2}=\frac{d^{2}V}{d{\varphi }^2}+(\frac{d\beta }{d\varphi }+4\frac{\beta }{M_{\mathrm{Pl}}})\frac{\beta }{M_{\mathrm{Pl}}}{A}^4(\varphi )(\tilde{\rho }-3\tilde{P})$$

$$d{s}^2=-e^{\nu (r)}d{t}^2+e^{\lambda (r)}d{r}^2+r^2\left(d{\theta }^2+{\sin }^2\theta d{\varphi }^2\right)$$

$$e^{-\lambda (r)}=1-\frac{2m(r)}{r}.$$

$$m^{\prime }$$

$${\nu }^{\prime }$$

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

$${\tilde{P}}^{\prime }=-\frac{1}{2}\left(\tilde{\rho }+\tilde{P}\right)\left({\nu }^{\prime }+2\frac{A^{\prime }}{A}{\varphi }^{\prime }\right).$$

$${\varphi }^{\prime \prime }+\left(\frac{2}{r}+\frac{1}{2}({\nu }^{\prime }-{\lambda }^{\prime })\right){\varphi }^{\prime }=e^{\lambda }\left[\frac{dV}{d\varphi }+A^3{A}^{\prime }(\tilde{\rho }-3\tilde{P})\right].$$

$$\frac{d{V}_{\mathrm{eff}}}{d\varphi }=\frac{dV}{d\varphi }+\frac{\beta }{M_{\mathrm{Pl}}}{A}^4(\varphi )(\tilde{\rho }-3\tilde{P})=0.$$

$$\delta \varphi =\sum \delta {\varphi }_{lm}(r)Y_{lm}(\Omega )e^{i\omega t}$$

$$-{\omega }^2\delta \varphi -e^{\nu -\lambda }\left[{\delta \varphi }^{\prime \prime }+\left(\frac{{\nu }^{\prime }-{\lambda }^{\prime }}{2}+\frac{2}{r}\right){\delta \varphi }^{\prime }\right]+e^{\nu }\left[\frac{l(l+1)}{r^2}+m_{\mathrm{eff}}^2\right]\delta \varphi =0.$$

$${\omega }^2\simeq k^2+m_{\mathrm{eff}}^2$$

$$(1-3\tilde{P}/\tilde{\rho }){\Phi }_{\star }\gtrsim -\frac{1}{24{\beta }_0^2}$$

$$\tilde{P}=\tilde{P}(\tilde{\rho })$$

$$\tilde{\rho }={\left(\frac{\tilde{P}}{K}\right)}^{1/\Gamma }+\frac{\tilde{P}}{\Gamma -1}$$

$$d{s}^2=-e^{\nu (r)}d{t}^2+{\left(1-\frac{2m(r)}{r}\right)}^{-1}d{r}^2+r^{2}d{\Omega }^2$$

$${\tilde{P}}^{\prime }$$

$$m^{\prime }$$

$${\nu }^{\prime }$$

$$m(r)=\frac{{\tilde{\rho }}_0{r}^3}{6M_{\mathrm{Pl}}^2}$$

$$\tilde{P}(r)={\tilde{\rho }}_0\frac{{\left(1-\frac{2G_N{M}_{\ast }}{r_{\star }}\right)}^{1/2}-{\left(1-\frac{2G_N{M}_{\ast }{r}^2}{r_{\star }^3}\right)}^{1/2}}{{\left(1-\frac{2G_N{M}_{\ast }{r}^2}{r_{\star }^3}\right)}^{1/2}-3{\left(1-\frac{2G_N{M}_{\ast }}{r_{\star }}\right)}^{1/2}}.$$

$$\tilde{\rho }-3\tilde{P}={\tilde{\rho }}_0\frac{4{\left(1-\frac{2G_N{M}_{\ast }{r}^2}{r_{\star }^3}\right)}^{1/2}-6{\left(1-\frac{2G_N{M}_{\ast }}{r_{\star }}\right)}^{1/2}}{{\left(1-\frac{2G_N{M}_{\ast }{r}^2}{r_{\star }^3}\right)}^{1/2}-3{\left(1-\frac{2G_N{M}_{\ast }}{r_{\star }}\right)}^{1/2}}.$$

$${\tilde{\rho }}_c-3{\tilde{P}}_c={\tilde{\rho }}_0\frac{4-6{\left(1-2{\Phi }_{\star }\right)}^{1/2}}{1-3{\left(1-2{\Phi }_{\star }\right)}^{1/2}}$$

Full text

TOV equation, scalar fields, modified gravity, neutron star

1 Introduction

General Relativity is at the heart of our understanding of compact objects such as neutron stars [1]. This could be challenged by the observation of the current accelerated expansion of the universe [2]. Indeed, the cosmological acceleration has now been corroborated using multiple probes [3] since the original observations with supernovae [4, 5] and whilst the addition of a cosmological constant term to Einstein’s equations fits all the data extremely well there are appealing theoretical reasons to investigate alternative explanations [6]. One of these possibilities could be a modification of the laws of gravity [7]. As argued by Weinberg [8] any modification to the framework of general relativity necessarily introduces additional degrees of freedom. This leads to many infrared modified gravity theories which, in principle, could incorporate the late time acceleration of the universe. They all modify the gravitational interaction and this is bound by stringent experimental constraints in our local environment at least [9, 10]. Generically, these alternative theories are valid when they include a screening mechanism that “shields” the effects of such modifications to general relativity (see [7] for a detailed review). Many proposals exist and they have all been extensively studied using laboratory [11, 12, 13, 14, 15, 16, 17, 18, 20, 19, 21, 22] solar system [23], cosmological [24], and astrophysical [25] tests. These investigations all probe gravity in the weak field regime. With advances in experimental techniques, attention has increasingly focussed on tests of gravity in the strong regime for compact objects such as neutron stars and black holes, [26, 27]. The effect of modified gravity in the vicinity of astrophysical black holes has been studied in [28, 29]

Many studies of the effects of modifying gravity in the strong gravity regime have originally been limited to models subject to the chameleon screening mechanism. For instance, there is a substantial number of works on the nature of compact objects such as neutron stars in $f(R)$ gravity. For a long time it was claimed [30, 31] that relativistic stars do not exist in these models due to the presence of an easily accessible singularity and other divergences springing from the functional form of $f(R)$ [30, 32, 33]. Several exceptions are now known to exist, e.g. [34] gives the first numerical construction of static relativistic stars in $f(R)$ gravity. Other examples are known [35, 36] including numerical solutions corresponding to static stellar configurations with a strong gravitational field [37, 38]. More recently relativistic stars have been studied in Horndeski [39, 40] and beyond Horndeski models [41, 42, 43].

In this paper, we focus on screened modified gravity [44, 45, 46] with two different screening methods: the chameleon [23, 47, 48] and the Damour-Polyakov mechanisms [49, 50, 51]. We concentrate on two typical examples. The first one is the inverse chameleon model which has already been analysed in the strong gravity regime and serves as a benchmark for our numerical results [38]. The second one is the environmentally dependent dilaton [50] which is motivated by string theory in the strong coupling regime and obeys the conjectured least coupling principle. It provides a simple and concrete example where the Damour-Polyakov mechanism is at play. We then apply the modified Tolman-Oppenheimer-Volkov equations due to the presence of the screened scalar field to these two models and study constant density and polytropic models of neutron stars within the chameleon and dilaton theories. We also numerically compute the scalar profiles and mass-radius relationships for chameleons and dilatons [52, 53, 54], and comment on parameter degeneracy for both models.

The outline of this paper is as follows. In section 2 we briefly review the screened modified gravity models that will concern us and restrict ourselves to static, spherically symmetric spacetimes. In section 3 we discuss the Tolman-Oppenheimer-Volkoff equation and our numerical procedure. We present our numerical results in section 4. We conclude with possible implications in section 5.

2 Screened Modified Gravity

2.1 The models

Our starting point is the Einstein-Hilbert action with the conformally coupled scalar field in the Einstein frame,

[TABLE]

Here $M_{\rm Pl}^{2}=1/8\pi G_{N}$ is the reduced Planck mass and $S_{m}$ is the action for matter fields (denoted generically as $\Psi_{i}$), which couple minimally to the Jordan frame metric, $\tilde{g}_{\mu\nu}=A^{2}(\phi)g_{\mu\nu}$. We assume $A(\phi)$ is close to 1 for values of $\phi$ in the allowed range below the Planck scale, i.e. we have

[TABLE]

where $A(\phi)$ and $\beta(\phi)$ will be specified later explicitly for both chameleons and dilatons. As we are always in the low energy regime, we take $\phi\ll M_{\rm Pl}$ throughout.

The Einstein equation can be derived as

[TABLE]

where the energy momentum tensor of the scalar field is given by

[TABLE]

and $T^{\mu\nu}_{m}$ is the Einstein frame energy-momentum of matter. Defining the trace of the energy momentum of matter $T_{m}\equiv g^{\mu\nu}T^{m}_{\mu\nu}$, the scalar field equation can then be written as

[TABLE]

In particular we need to pay attention to the conformal transformation between the Einstein frame metric $g_{\mu\nu}$ and the Jordan frame metric $\tilde{g}_{\mu\nu}$. The perfect fluid energy momentum tensor in the Einstein frame can be written as

[TABLE]

where $u^{\mu}$ is the 4-velocity of the fluid elements, $\rho$ the energy density and $P$ the pressure. The Jordan frame energy-momentum tensor ${\tilde{T}}_{\mu\nu}$ and the Einstein frame energy momentum tensor, $T_{\mu\nu}$, are related by

[TABLE]

so that, for a perfect fluid, we have

[TABLE]

The trace of the energy momentum tensor is $T_{m}\equiv g^{\mu\nu}T^{m}_{\mu\nu}=-\rho+3P$ in the Einstein frame. The numerical values for the pressure and density are essentially the same in both frames given our constraint $\phi\ll M_{\rm Pl}$ and $A\sim 1$.

In the following we shall concentrate on two types of models subject to the chameleon and Damour-Polyakov mechanisms. Both types of scenarios require the existence of a minimum of the effective potential $\phi({\tilde{\rho}},{\tilde{P}})$

[TABLE]

in terms of the Jordan matter density and pressure, taken as conserved in the Jordan frame and scalar field independent. For chameleons the effective mass of the scalar field at the minimum of the effective potential

[TABLE]

must be large enough in dense environments, such as the interior of a star, in order to Yukawa suppress any effect of the scalar field on outside bodies. On the other hand, the Damour-Polyakov mechanism operates when $\beta(\phi({\tilde{\rho}},{\tilde{P}}))$ becomes extremely small in dense environments and effectively suppresses all fifth force effects.

2.2 Static and spherically symmetric configurations

We are interested in the effect of screened modified gravity in the context of compact and relativistic objects. In such a case, we can restrict ourselves to a static, spherically symmetric background geometry with

[TABLE]

where the scalar field is such that $\phi=\phi(r)$ and we parameterise

[TABLE]

The $t$ and $r$ components of Einstein’s equations then become,

[TABLE]

where the fluid parameters are defined in the Jordan frame, i.e. we analyse the Einstein equations in the Einstein frame, where the gravitational constant is a true constant, as a function of the conserved fluid parameters as defined by the Jordan frame. Therefore we use the conservation equation in the Jordan frame

[TABLE]

which gives us the relation

[TABLE]

The Klein-Gordon equation also simplifies to

[TABLE]

As already stated, the source term on the right hand side can be seen as an effective potential for the scalar field up to the metric component $e^{\lambda}$. A local extremum of the potential is characterised by

[TABLE]

The effective mass of the scalar field $\frac{d^{2}V_{\rm eff}}{d\phi^{2}}$ at the minima is always positive if the potential $V(\phi)$ is convex and $({\tilde{\rho}}-3{\tilde{P}})>0$. As it turns out, the latter is an important condition for the stability of the solution.

2.3 Stability

A stability analysis of the scalar equation of motion 17 provides further insight. Following  [55] we note that since we are in a static, spherically symmetric geometry, we can perturb the scalar into spherical harmonics ($Y_{lm}(\Omega)$) and ignore all other perturbations to leading order. This is enough for the qualitative argument presented here and this implies that

[TABLE]

and the mode equation becomes

[TABLE]

Focussing on the radial part, we get an approximate dispersion relation of the form

[TABLE]

where $k$ is proportional to the inverse of the typical length scale in the problem, the stellar radius $r_{\star}$. For stability we require $\omega^{2}>0$. Using ${\tilde{\rho}}_{\star}\sim M_{\star}/r_{\star}^{3}$ and $\Phi_{\star}\sim M_{\star}/(M_{\rm Pl}^{2}r_{\star})$, we end up with the condition

[TABLE]

where $\beta_{0}$ is the coupling in vacuum outside the object, $M_{\star}$ the stellar mass and $\Phi_{\star}$ the Newtonian potential at its surface. For neutron stars with $\Phi_{\star}\sim 0.1$ and $\beta_{0}\sim 1$, we find that as soon as $(1-3{\tilde{P}}/{\tilde{\rho}})\lesssim-1$ we would generically expect instabilities to appear. Eventually, this would back react and prevent the existence of the object itself.

The nature of this instability depends essentially on one more ingredient which characterises the compact object: the fluid equation of state, or more precisely the relation between pressure and matter density

[TABLE]

which we specify in the Jordan frame and which closes the system of equations. The equation of state deep inside a neutron star is not very well known. Several approximations exist and the uncertainties are reasonably large. For our purposes it will be sufficient to focus on a polytropic equation of state for the interior region of the neutron star. See [56]–[57] for a detailed discussion. We thus use

[TABLE]

where $\Gamma$ is the polytropic index and $K$, of mass dimension $4(1-\Gamma)$, is a scale to be determined below. Realistic neutron stars can be parameterised by piecewise polytropic equations of state with $\Gamma\sim 1-3$. As it turns out, the lower the polytropic index, the stiffer the mass radius relationship. We use values close to $2$ for our numerical work.

3 The scalar profile for chameleons and dilatons

3.1 Tolman-Oppenheimer-Volkoff equation

We have to simultaneously solve equations 13–23 and given their complexity we will have to resort to numerics. We will come to the numerical implementation of our relaxation algorithm in the next section, but before we do that it is helpful to recall the solution in GR, i.e. when neglecting the scalar back reaction.

The GR problem was solved by Tolman, Oppenheimer and Volkoff via what is generally known as the Tolman-Oppenheimer-Volkoff (TOV) equation. It constrains the structure of a spherically symmetric body of isotropic material which is in static gravitational equilibrium, as modelled by GR. This equation is derived by solving Einstein’s equations for a general time-invariant, spherically symmetric metric. Let us consider an interior Schwarzschild metric which gives the following line element

[TABLE]

where $m(r)/G_{N}$ has the simple interpretation of being the mass inside the radius $r$ whilst $\nu(r)$ or more precisely $\nu(r)/2$ is known as the metric potential. In vacuum, outside the star, $m=G_{N}M_{\star}$, $T_{m}^{\mu\nu}=0$ and $\nu(r)=\ln(1-2G_{N}M_{\star}/r)$ where $M_{\star}$ is the total stellar mass. Now assuming a barotropic equation of state (i.e. the pressure is a function of density alone), the system is governed by the set of equations which can be obtained by discarding the scalar field from the equations derived previously

[TABLE]

where conventionally eq 26 is known as the TOV equation.

Using for the density a constant $\tilde{\rho}=\tilde{\rho}_{0}$, this leads to

[TABLE]

and one gets for the pressure

[TABLE]

The trace of the energy momentum tensor can be computed as,

[TABLE]

This quantity is of great interest as it directly couples to the scalar field which is here considered as a test field. As $\tilde{\rho}-3\tilde{P}$ is an increasing function of $r$ , we find its minimal value at the centre of the star

[TABLE]

where $\tilde{\rho}_{c}$ and $\tilde{P}_{c}$ are respectively the density and pressure at the stellar centre and $\Phi_{\star}=G_{N}M_{*}/r_{*}$ is the gravitational potential at the stellar surface. Thus if $\Phi_{\star}$ is greater that $5/18$, we can see that $\tilde{\rho}-3\tilde{P}$ is negative in the stellar centre and this potentially causes an instability in the scalar equation of motion [37, 38]. When this is not the case, one needs to resort to numerical results in order to understand the full dynamics. In the following, we shall concentrate on inverse power law chameleons and environmentally dependent dilatons for which the dynamics will be analysed.

3.2 Chameleons

The inverse power law chameleons combine an inverse power law potential for a scalar field with a constant coupling to matter [23]. It is well known that the chameleon mechanism works for these models and that the local tests of gravity impose a constraint on the scale $M$ appearing in the scalar potential, which can be taken of the order of the dark energy scale [23]. The coupling function is simply given by

[TABLE]

where $\beta$ is constant. The interaction potential is

[TABLE]

where $n\geq 1$ is an integer of order one, and we define $V_{0}\equiv M^{4+n}$ to simplify notation. Keeping only the leading order term from the coupling function, we see that the effective potential is

[TABLE]

Numerically, we first focus on the simple case of a constant density star with an equation of state $\tilde{\rho}=\tilde{\rho}_{0}$. For numerical ease we use a smoothed out ($1-\tanh$) transition to the asymptotic matter distribution surrounding the star. This can either be vacuum or a constant density Schwarzschild de Sitter geometry. Realistic stars should require matching to the time dependent cosmological background geometry, however our solution should capture all the salient features of this more general case.

3.3 Dilatons

We follow a similar procedure for the dilaton fields. The environmentally dependent dilaton [50] is such that the coupling function has a minimum at $\phi=\phi_{d}$ as prescribed by the least coupling principle. The potential is taken to be exponentially decreasing as conjectured for the string dilaton in the strong coupling regime. The screening behaviour of dilatons is a consequence of the fact that their coupling function has a minimum at $\phi_{d}$,

[TABLE]

or equivalently we have for the coupling

[TABLE]

The dilaton self-interaction potential is

[TABLE]

where again we define $V_{0}\equiv M^{4}e^{-\phi_{d}/M_{\rm Pl}}$ to simplify notation. The effective potential is thus

[TABLE]

which can be expanded in $\phi/M_{\rm Pl}$. Solar system tests of the dilaton models imply that typically one must have $a_{2}\gtrsim 10^{6}$ [50] guaranteeing that the mass of the dilaton, on large cosmological scales, is larger than $10^{3}H_{0}$ [46], where $H_{0}$ is the Hubble rate now.

3.4 Numerical procedure for neutron stars

The solution of the equations of motion given a specified equation of state requires a numerical procedure. We solve the above set of equations 13–17 together with the equation of state numerically. We follow the general framework of  [58] where the authors study relativistic stars in $f(R)$ gravity. In our case we expect the scalar field to exhibit a smoothed out kink-like behaviour extrapolating between the minimum inside the neutron star to the minimum of the ambient matter density outside the stellar surface. The boundary conditions are set at $r=0$ and not at the surface of the star $r=r_{\star}$. For numerical ease we begin by rescaling all the variables

[TABLE]

where $r_{0}$, $M_{0}$ are parameters specific to the star/modified gravity theory being considered. We will make them explicit later.

Substituting the above variables into the Einstein and scalar equations of motion we obtain,

[TABLE]

The parameter $M_{0}$ essentially sets a mass scale for the scalar field and will depend on the modified gravity theory considered. To close this system of equations, we also need the equation of state for the fluid making up the star and written out in terms of the rescaled variables this reads,

[TABLE]

We use a polytropic equation of state and in our rescaled variables we get

[TABLE]

$\hat{K}$ here is the overall normalisation of the polytropic equation of state and for us it is just a numerical constant. We pick the value of $0.39$ which is consistent with what is used in the literature [60]. We now need to pick a value for the scale $r_{0}$ used to normalise all length scales in our problem. A convenient choice is of the order of

[TABLE]

so that for $\Gamma=2$ and $\hat{K}=0.39$ we get ${\hat{\rho}}_{\rm c}=1$ and ${\hat{P}}_{c}\simeq 0.23$. With these choices $r_{*}\simeq r_{0}\times\mathcal{O}(1)$ for highly relativistic stars, where $r_{*}$ is the stellar radius.

The choice of the various numerical parameters is motivated by fixing the normalisation scale to a specific, stable, stellar solution of the Tolman-Oppenheimer-Volkoff equations in GR. It can be shown [58] that such a solution corresponds to a star with central energy density of $\tilde{\rho}_{c}=6.47\times 10^{14}\;\rm{g/cm}^{3}$, a global mass of $M_{\star}=3.10M_{\odot}$ and a gravitational potential of $|\Phi_{\star}|\simeq 0.21$. The stellar configuration in chameleon and dilaton gravity theories is almost identical to this solution. We use a relaxation type algorithm and vary the relaxation parameter to ensure fast and smooth convergence. We vary the grid size by hand to zoom in on regions where the scalar field varies rapidly. A non-homogeneous grid would provide better resolution but is beyond the scope of this work. We do however note that for the theories and parameter ranges considered, there is no significant loss of resolution.

We now turn to the boundary conditions. We have three first-order differential equations and one second-order differential equation and thus we require five boundary conditions.

Demanding regularity at the centre of the star ($a=0$) implies

[TABLE]

where the second condition ensures that the metric functions are well behaved at the stellar centre. This corresponds to a vanishing mass function $m(r)$ at the origin. 2. 2.

The condition on $\nu$ corresponds to defining the time coordinate as the proper time of a static observer at a point $a_{1}$ far from the stellar surface.

[TABLE] 3. 3.

As $a_{1}\rightarrow\infty$ we expect

[TABLE]

This “asymptotic” value ${\hat{\phi}}_{\infty}$ corresponds to the minimum of the effective potential $V_{\rm eff}$ far from the star where this behaviour can be solved for analytically. In the limit of $a\rightarrow\infty$ the scalar equation becomes

[TABLE]

and for a massive scalar field the solution looks like

[TABLE]

where $\hat{m}$ is the scalar mass and $B$ is a constant. Numerically it is easier to use this fall off rather than imposing a hard boundary condition. 4. 4.

The boundary conditions for the pressure are either

[TABLE]

For a constant density star we choose ${\hat{P}}(a_{1})=-{\hat{\rho}}_{\infty}$ so that the solution is asymptotically de Sitter, while for the polytropic equation of state it is easier to fix the pressure at the stellar centre and this once again gives an asymptotically de Sitter solution.

4 Numerical Results

In this section we present the results of our numerical integration. We carry out the integration for chameleon and dilaton modified gravity models.

4.1 Chameleons

For fixed values of the coupling $\beta\sim 1$ and the parameter $n=1$, the numerical solution changes with the central density ${\hat{\rho}}_{c}$, the mass scale $M_{0}$ and the asymptotic density $\tilde{\rho}_{\infty}$. More concretely, we take the density to vary according to

[TABLE]

Choosing $\epsilon$ to be sufficiently small we recover the step function profile where ${\hat{\rho}}\sim{\hat{\rho}}_{c}$ for $a\lesssim 1$ and ${\hat{\rho}}\sim{\hat{\rho}}_{\infty}$ for $a\gtrsim 1$ and thus the stellar radius is around $a\sim 1$. The ${\hat{\rho}}_{\infty}$ value should be the typical galactic density $\tilde{\rho}_{\infty}\sim 10^{-24}\rm{g/cm^{3}}$. This implies that $\hat{\rho}_{\infty}\sim 10^{-38}$. This is a more appropriate asymptotic density than typical cosmological densities as we expect neutron stars to be in galaxies and as is apparent from Eq. (51) for small enough values of $\epsilon$, ${\hat{\rho}}\sim{\hat{\rho}}_{\infty}$ for $a\gtrsim 1$. We do however note, that given the large density contrast between the neutron star interior and the asymptotic density used, we do not find our numerics to be very sensitive to this parameter. For numerical convenience we will use a larger value of $\hat{\rho}_{\infty}\sim 10^{-4}$. Realistic stars will require matching to the time dependent cosmological background geometry, however our solution should capture all the salient features of this more general case.

We always set $M_{0}$ such that the asymptotic value of the scalar field is normalised to a value $\mathcal{O}(1)$. This is done purely for numerical ease and in particular we use the convenient value of

[TABLE]

Thus the asymptotic value of the scalar field is given by

[TABLE]

so that we have ${\hat{\phi}}_{\infty}\sim 1.3$. We define the rescaled potential as $\hat{V}=V\frac{r^{2}_{0}}{M_{0}^{2}}$ expressed in terms of the rescaled variables as

[TABLE]

and the effective mass as

[TABLE]

For our numerics, we use the central density of the star to define the effective mass. Note that all of this is defined up to an overall numerical parameter $\hat{V}_{0}$. Realistic values for the chameleon with $n=1$ are such that $V_{0}\sim 10^{-15}$ eV. This is far too small a number for numerics and we instead use $\hat{V}_{0}\sim\mathcal{O}(0.01)$.

The effective rescaled mass is large in the stellar interior and small outside and thus we expect the scalar to track the local minima of the effective potential in the stellar interior and smoothly evolve towards ${\hat{\phi}}_{\infty}$ outside the star. We can further see that the sign of $({\hat{\rho}}-3{\hat{P}})$ is crucial in determining if this effective mass is well defined and consequently if the solution is stable.

We plot the scalar field profile for a range of central densities ${\hat{\rho}}_{c}$ in the left panel of figure  1. For low central densities the field smoothly transitions from its minima inside the star to the asymptotic minimum far from the stellar surface. Increasing this scaled central density while keeping the radius fixed has the effect of increasing the physical density of the star and this in turn means that the field is very stiffly held on to its minimum value inside the star. There is now a sharp transition close to the surface of the star, much like the chameleon thin shell effect. The field then smoothly approaches its asymptotic global minimum outside the star. It is interesting to note that this behaviour is almost exactly what we get if we keep the central density of the star fixed and tune the effective mass of the chameleon scalar inside the star. This is shown in the right panel of figure 1. There are two complementary explanations for this behaviour. In the first case it is easier to think in terms of stellar parameters. Our relaxation algorithm tries to minimise the net energy (potential plus gradient) with our boundary conditions of ${\hat{\phi}}^{\prime}=0$ at $a=0$ and ${\hat{\phi}}$ fixed asymptotically. As the energy density inside the star increases it is increasingly unfavourable for the scalar to deviate significantly from its minimum value. Close to the stellar surface the smoothed out kink in $\rho$ forces the scalar to quickly evolve to its asymptotic value. Changing the mass of the scalar has a similar effect, a large mass, i.e. a small Compton wavelength, scalar switches between the two minima sharply whereas a small mass, large Compton wavelength, scalar has a sluggish transition.

For a polytropic equation of state for the neutron star we get a qualitatively similar behaviour. A slight modification to eq. 43 means that we can smoothly transition to the ambient background/cosmological density and this also helps us avoid certain numerical instabilities close to the stellar surface. This point is discussed in greater detail in [58]. We thus use

[TABLE]

where ${\hat{\rho}}_{\infty}$ is the ambient density outside the star and is equal to zero for vacuum. The left panel in figure  2 shows the evolution of the scalar field for a range of scalar masses respectively. The right panel shows that similar profiles for different values of the central masses can be obtained by changing $\Gamma$ respectively. Hence there is a strong degeneracy between the neutron star equation of state and the modified gravity model parameters.

Finally for the sake of completeness we plot the scalar field profile for the case of ${\hat{\rho}}-3{\hat{P}}<0$ in the centre of the star in the left panel of fig. 3. As mentioned previously, in this case the scalar field does not have a well defined minimum inside the star and diverges. This divergence obviously means that the system is no longer stable as the energy momentum tensor associated with the scalar also diverges implying that the background geometry is no longer stable.

4.2 Dilaton

In the case of dilatons, we get the following expansion in $\phi/M_{\rm Pl}$ for the effective potential

[TABLE]

The minimum of the effective potential is obtained at a value $\phi_{\rm min}$ defined by

[TABLE]

As for chameleons we rescale the field so that at infinity when ${\tilde{P}}_{\rm\infty}/{\tilde{\rho}}_{\rm\infty}=-1$ we get a rescaled field of order $\mathcal{O}(1)$. Moreover, it is convenient to shift and redefine $(\phi-\phi_{d})\to\phi$. Using the scale

[TABLE]

we obtain that

[TABLE]

and at infinity, ${\hat{\phi}}_{\infty}=0.5$. We can now express the rescaled effective potential $\hat{V}_{\rm eff}$ in terms of the rescaled variables

[TABLE]

and derive the corresponding mass. It should be noted that these expressions are defined up to an overall numerical constant $\hat{V}_{0}$ which we take to be $\mathcal{O}(0.01)$ for numerical ease.

We begin with a constant density star with equation of state $\rho=\rho_{0}$ and a smoothed out transition to the asymptotic matter distribution surrounding the star. We plot the scalar field profile for a range of central densities ${\hat{\rho}}_{c}$ in the left panel of figure 4 and for a range of effective dilaton masses in the right panel of figure 4. We can see that for low central densities the field smoothly transitions from its minimum inside the star to the asymptotic minimum far from the stellar surface. Increasing the central density or the mass of the dilaton have a similar effect on scalar profile. The transition becomes sharper and more kink like and occurs closer to the surface. The dilaton scalar profiles are however noticeably smoother than the chameleon profiles from the previous section.

For a polytropic neutron star we use the modified equation of state 56. In the panels of figure  5, we show the evolution of the scalar field for a range of scalar masses. Again and as for chameleons we show in the left panel that changing $\Gamma$ can be mimicked by a change of the central mass and vice versa, i.e. the equation of state is degenerate with a change of the dilaton parameters.

Finally in fig. 3 we exemplify the scalar instability that develops in the dilatonic case for ${\hat{\rho}}-3{\hat{P}}<0$. This criterion therefore appears to be a good model independent way to test the onset of such an instability for stellar structures in both the chameleon and dilaton cases.

5 Implications and conclusions

We have studied scalar field profiles around neutron stars for different modified gravity models corresponding to two screening mechanisms, i.e. the chameleon and Damour-Polyakov mechanisms. This has been done for two classes of models: the inverse power law chameleons and the environmentally dependent dilaton. The latter is particularly interesting as it is expected to arise in the strong coupling regime of string theory. We have solved the dynamics of these models using a relaxation algorithm allowing us to study the modified Tolman-Oppenheimer-Volkoff equations for constant density stars and the slightly more realistic polytropic neutron stars for a range of model parameters.

We find non-trivial modifications to the stellar structure due to the addition of a scalar degree of freedom to gravity. In particular we find that deep inside the star, the scalar field settles down to the minimum of the effective potential which depends on the modified gravity model used and the details of the equation of state used to model the neutron star. Close to the stellar surface, the scalar field starts responding to the changing ambient density and this response is either sluggish or extremely rapid depending on the model parameters used, i.e. eventually on the central mass of the scalar field. Finally far away from the neutron star the scalar field settles down to the minimum of the effective potential with the ambient cosmological density. Thus the scalar profile on average looks like a smoothed out kink-like solution extrapolating between the two minima. We find that the scalar field develops an instability deep in the stellar centre for ${\tilde{\rho}}-3{\tilde{P}}<0$. This criterion appears to be robust for a range of parameter values and for both chameleon and dilaton theories. Indeed a crude theoretical estimate shows that we would expect the radial mode of spherically symmetric scalar perturbations to develop an instability for ${\tilde{\rho}}-3{\tilde{P}}\lesssim 0$.

For realistic models which pass the solar system tests, we expect the mass of the scalar field to be very large deep inside the star, due to the very large central density. Hence the kink-like solution should vary very abruptly when interpolating between the two minima, implying that the effects of the scalar field on the dynamics of star should be reduced. This can be inferred from the figure 5 where the profile is depicted and the figure 7 where the change of the mass-radius relationship is shown to be smaller when the central mass increases. Of course, taking into account the actual solar system constraints on the models and the resulting effects on the field profile and the mass radius relationship is beyond the scope of the present paper, as it very likely that different numerical methods would have to be used.

Another interesting aspect of the numerical solutions is the very high degree of degeneracy between parameters of the modified gravity theory and parameters of the neutron star equation of state. We find that by changing the effective mass parameter of chameleons and dilatons we can mimic any change in the scalar profile that results from varying the stellar central density or the polytropic index of the equation of state. This can be strengthened by analysing the mass-radius relationship for the chameleon, as in figures 6, and for the dilaton, as in figures 7. We have determined the radius $r_{\star}$ of the star as the point where the pressure vanishes $\tilde{P}(r_{\star})=0$ and the mass as $G_{N}M_{\star}=m(r_{\star})$. We have used the rescaled mass and radius defined by

[TABLE]

as $K$ has mass dimension $4(1-\Gamma)$, which has the advantage of making the overall normalisation of the polytropic equation of state $K$ not very relevant to our results [59].

It can be seen that varying the central mass of the scalar for instance has a similar effect on the mass-radius relationship as varying the polytropic index for the stellar equation of state. As the experimental constraints on the stellar equation of state for neutron stars are weak [60], it is likely that for dilatons and chameleons the neutron star observables that depend only on the mass and the radius of the star are not guaranteed to distinguish between small differences in the equation of state versus small modifications to gravity. It was argued in [62] that observables that depend on the effective surface gravity of the neutron star can in principle break this degeneracy. However it has been shown that this might not be the case for $f(R)$ theories [63]. It may be possible to break such a degeneracy by using the near independence on the equation of state of the relationship between the reduced moment of inertia $\bar{I}=I/G_{N}^{2}M_{\star}^{3}$ and the compactness $C=G_{N}M_{\star}/r_{\star}$ of the star [61]. Other observables might be useful, such as cooling rates of neutron stars. A neutron star cools via photon and neutrino emission. Whilst the photon luminosity depends only on the density of the photosphere surrounding the neutron star, neutrino cooling is sensitive to the central temperature. As a result, it would be interesting to analyse if observations of the surface temperatures of neutron stars could lead to useful constraints on the deviations from general relativity [64]. Such an analysis for chameleons and dilatons, however, is beyond the scope of this work. Finally it would be extremely interesting to study how our results change when the equation of state for the matter in the star is modified, i.e. not of the polytropic type anymore. We hope to come back to this point in another publication.

Acknowledgments.

We would like to thanks E. Babichev and D. Langlois for suggestions. RJ is supported by the Cambridge Commonwealth Trust and Trinity College, Cambridge. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 690575. This article is based upon work related to the COST Action CA15117 (CANTATA) supported by COST (European Cooperation in Science and Technology). ACD is supported in part by STFC UK under ST/L000385/1 and ST/L000636/1.