Unscreening scalarons with a black hole

TL;DR

This paper challenges the common belief that scalar degrees of freedom are always suppressed near black holes, showing that low accretion rate black holes can unscreen scalar fields, with new techniques developed for their analysis.

Contribution

The authors develop a novel, accurate method to study scalar field dynamics around black holes, including static profiles and stability, applicable to various modified gravity models.

Findings

Scalar fields can remain unscreened around low accretion rate black holes.

New computational tools for analyzing scalar field stability near black holes.

Application to multiple modified gravity models, including f(R) and symmetron.

Abstract

It is typically believed that the additional degrees of freedom in any modification of gravity are completely suppressed by the large energy densities coexisting with an astrophysical black hole. In this paper, we find that this might not always be the case. This belief holds for black holes formed via gravitational collapse in very dense environments, whereas the black holes with sufficiently low accretion rates that have low matter densities inside innermost stable circular orbit will generally unscreen chameleons. We develop a novel technique to study the dynamics of accretion of a scalar field onto a Schwarzschild-like black hole which is accurate on both short and long time scales. In particular, we study the behavior of the extra scalar degree of freedom in the Starobinsky and Hu-Sawicki $f(R)$ theories, for the symmetron model, and for the Ratra-Peebles model. Aside from calculating non-trivial static field profiles outside the black hole, we provide the tools to study the (in)stability and evolution towards the equilibrium solution for any generic well behaved set of parameters and initial conditions. Our code is made publicly available for further research and modifications to study other models.

Equations

$$L=\left(1-\sqrt{\frac{8}{9}}\right)\dot{M}\sim 0.06\dot{M}.$$

$$\rho =\{\begin{array}{ll}\sigma {\rho }_0;& r_g<r<r_{\text{ISCO}}\\ {\rho }_0;& r\ge r_{\text{ISCO}},\end{array}$$

$$\rho =(1-\sigma ){\rho }_0\frac{1+\tanh a_0(r-r_{\text{ISCO}})}{2}+\sigma {\rho }_0,$$

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

$$d{s}^2=\left(1-\frac{r_g}{r(x)}\right)\left[-d{t}^2+d{x}^2\right]+r^2(x)d{\Omega }^2.$$

$$\begin{split}\biggl{[}-\frac{\partial^{2}}{\partial t^{2}}+&\frac{\partial^{2}}{\partial x^{2}}+\frac{2}{r(x)}\left(1-\frac{r_{g}}{r(x)}\right)\frac{\partial}{\partial x}\biggr{]}\phi(t,x)\\ &-\left(1-\frac{r_{g}}{r(x)}\right)V^{\prime}_{\text{eff}}(\phi)=0.\end{split}$$

$$S=\int \frac{f(R)}{16\pi G}\sqrt{-g}{d}^{4}x+S_m[g_{\mu \nu },\psi ],$$

$$\begin{array}{rl}{f}_R{R}_{\mu \nu }-\frac{1}{2}f& g_{\mu \nu }=8\pi G{T}_{\mu \nu }+\\ & {\nabla }_{\mu }{\nabla }_{\nu }{f}_R-g_{\mu \nu }\square f_R,\end{array}$$

$$\square f_R=\frac{1}{3}(2f-f_{R}R)+\frac{8\pi G}{3}T,$$

$$\square \varphi =V^{\prime }(\varphi )-\mathcal{F},$$

$$\frac{dV}{dR}=\frac{dV}{d\varphi }\frac{d\varphi }{dR}=\frac{1}{3}\left(2f-f_{R}R\right)f_{RR},$$

$$\varphi (R)$$

$$V(R)$$

$$f_{\text{S}}=R+\lambda \left[\frac{1}{{\left(1+(R/R_0{)}^2\right)}^n}-1\right]R_0,$$

$$f_{\text{HS}}=R-\frac{\alpha (R/R_0{)}^n{R}_0}{1+\beta (R/R_0{)}^n},$$

$$a-(n-1)c=0,$$

$$b-nc=0,$$

$$V_{\text{eff}}(\varphi )=V(\varphi )-\mathcal{F}\varphi .$$

$$V^{\prime }(\varphi )=\mathcal{F}.$$

$${\tilde{f}}_{\text{S}}=R+\lambda \left[\frac{1}{{\left(1+(R/R_0{)}^2\right)}^n}-1\right]R_0+\frac{\mu }{R_0}{R}^2,$$

$${\tilde{f}}_{\text{HS}}=R-\frac{\alpha (R/R_0{)}^n{R}_0}{1+\beta (R/R_0{)}^n}+\frac{\mu }{R_0}{R}^2,$$

$$\begin{split}S=\int\biggl{[}\frac{R}{16\pi G}&-\frac{1}{2}(\nabla\phi)^{2}-V(\phi)\biggr{]}\sqrt{-g}\,d^{4}x\\ &+S_{m}[A^{2}(\phi)g_{\mu\nu},\psi],\end{split}$$

$$V(\varphi )=V_0-\frac{\mu }{2}{\varphi }^2+\frac{\lambda }{4}{\varphi }^4,$$

$$A(\varphi )=1+\frac{\hslash }{2}\frac{{\varphi }^2}{m_{\text{S}}^2},$$

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi G\left(T_{\mu \nu }^{[\varphi ]}+T_{\mu \nu }^{[m]}\right),$$

$$\square \varphi =V_{\text{eff}}^{\prime }(\varphi ).$$

$$V_{\text{eff}}(\varphi )={\tilde{V}}_0+\frac{1}{2}\left(-\frac{\hslash T}{m_{\text{S}}^2}-\mu \right){\varphi }^2+\frac{\lambda }{4}{\varphi }^4,$$

$$\varphi (r,t=t_0)=\{\begin{array}{ll}{\varphi }_{\ast },& r<r_{\text{ISCO}}\\ 0,& r\ge r_{\text{ISCO}},\end{array}$$

$$V(\varphi )=V_0+\frac{{\gamma }^{n+4}}{{\varphi }^n}$$

$$A(\varphi )\approx 1+\frac{\epsilon \varphi }{m_{\text{C}}},$$

$$V_{\text{eff}}(\varphi )=V_0+\frac{{\gamma }^{n+4}}{{\varphi }^n}+\frac{\epsilon \varphi }{m_{\text{C}}}(-T).$$

$${\bar{\varphi }}_{\text{C}}={\left(n{\gamma }^{n+4}\frac{m_{\text{C}}}{-\epsilon T}\right)}^{\frac{1}{n+1}},$$

Full text

I Introduction

Einstein’s theory of general relativity (GR) is one of the greatest achievements of modern science, due to its ability to describe many of the gravitational phenomena with impressive level of detail. Nevertheless, there are many motivations to look for alternative theories: the observation of cosmic acceleration, the unknown nature of dark matter and dark energy and the lack of an ultraviolet completion of general relativity related with the unavoidable presence of spacetime singularities in this theory. Any successful modification of Einstein’s theory should solve some of these conundrums, while remaining consistent with the local tests of general relativity and do no harm to the standard behaviour of matter as discussed in Refs. Will:2014kxa ; Burgess:2013ara .

For example, $f(R)$ gravity was postulated in Ref. Stelle:1976gc as an attempt to produce a renormalisable theory of quantum gravity. It introduces a ghost-free functional of curvature, analogous to the notion of free energy in thermodynamics. Since its first incarnation, however, different forms of $f(R)$ have been used to solve issues at cosmological scales (see Refs. Starobinsky:1980te ; Hu:2007nk ; Starobinsky:2007hu ). Notwithstanding, either considering a functional source of gravitational interaction different from the Ricci scalar, or simply adding extra gauge degrees of freedom in the gravitational sector has consequences at all energy scales. Gravity is not only defined by the metric, but also by other fields. The force exerted by any extra degree of freedom typically contradicts the stringent local constraints compatible with general relativity. Screening mechanisms such as theories in Refs. Khoury:2003aq ; Pietroni:2005pv ; Vainshtein:1972sx are designed to circumvent this issue by diluting the sources of this force.

In the case of an astrophysical black hole, one naively expects that every extra force will be screened by large environmental energy densities. In the particular case of a force sourced by a scalar, screening is expected to be a consequence of the no-hair theorem, as stated in Ref. Bekenstein:1995un , applicable in the case of a source with positive energy density. Nonetheless, as noticed in Ref. Davis:2014tea , the presence of a non-trivial distribution of accreting matter and the way this is coupled to the scalar field allows scalar hair outside the event horizon.

In this paper, we simulate the accretion of the additional scalar degree of freedom in the $f(R)$ models presented in Refs. Hu:2007nk ; Starobinsky:2007hu , the symmetron model in Ref. Hinterbichler:2010es and the Ratra-Peebles chameleon in Ref. Ratra:1987rm in a dense environment. We assume that in such a region we can form a stable Schwarzschild black hole via gravitational collapse. This spherically symmetric black hole interacts with the screened extra real scalar field in a non-trivial matter distribution and accretes around the event horizon. The accreting matter density profile proposed here is a slight generalization to the suggested in Ref. Davis:2014tea . Even when the accretion in astrophysical black holes takes place in rotating spacetime solutions, it is interesting to explore the test case of a spherically symmetric system.

The purpose of this study is to visualize the process of accretion of the scalar field for various choices of model parameters, and discuss the time scales involved in the convergence to a non-trivial static profile. We designed a compact and efficient spectral code, flexible to modification and capable of producing time-dependent solutions of a scalar field in a Schwarzschild-like background. It is available for further application to other models at https://github.com/andrei-v-frolov/accretion.

The layout of this paper is as follows. In section II, we define the matter environment and the spacetime geometry in which the extra scalars propagate. In section III, we review the equations of motion and particular features of the scalar degrees of freedom in the chosen $f(R)$ models and study how these fields accrete. In section IV, we modify these models by adding terms proportional to $R^{2}$ and show the dynamical resolution of the curvature singularities in the Hu-Sawicki and the Starobinsky model. We extend our treatment for the symmetron model and the Ratra-Peebles chameleon in sections V and VI, respectively. Finally, we present summary of our results and conclusions in section VII. Numerical implementation is discussed in detail in Appendix A.

II Matter distribution for spherically symmetric black holes

Although astrophysical black holes are rotating systems, spherically symmetric solutions remain interesting when exploring some of the dynamical subtleties of the extra scalars in modified gravity. From Refs. Scheel:1994yn ; Frolov:2004rz , we learn that the fate of a collapsing spherically symmetric system in standard scalar-tensor theories is to become a Schwarzschild black hole.

It is then sensible to ask whether the Schwarzschild solution remains a valid description of spacetime even in the presence of accreting matter. To answer this question, we briefly review the dynamics of a steady flow of matter to estimate the accretion rate of the black hole. We calculate the luminosity assuming the power lost by a generic inflow of particles traveling from infinity to the innermost stable circular orbit (ISCO) at $r_{\text{ISCO}}=6GM/c^{2}$

[TABLE]

Using the table of luminosities of active galactic nuclei included in Ref. Khorunzhev:2012ga , we note that the accretion rate of a supermassive black hole with $M_{\text{BH}}=10^{9}M_{\odot}$ is $\dot{M}\lesssim 1M_{\odot}/\text{yr}$. The accretion rate in a rotating black hole of similar luminosity could be even smaller, since $L\sim 0.42\,\dot{M}c^{2}$ for near-extremal rotation. This implies that it would take $\sim 10^{7}$ years to change the mass of the black hole by 1%. Therefore, it is reasonable to assume that spacetime is static within time scales we are interested in.

At this point of the discussion, it is necessary to provide an approximate expression of the matter density distribution outside the black hole’s horizon. This is not by any means a full discussion of the radial structure equations for accretion disks. However, we provide sufficient arguments to justify our choice of a matter distribution. It must remain nearly static within the time scale estimated previously. For that purpose, we consider that radial matter density at a given radius $r$ is proportional to some positive power of the time that particles spend in orbits passing through that $r$. Stable orbits are possible only when $r\geq r_{\text{ISCO}}=6GM/c^{2}$. Therefore, the dynamics of a test particle moving in Schwarzschild spacetime only provides two scenarios in which matter can be found at $r<r_{\text{ISCO}}$: (I) these are on “no-return” trajectories towards the horizon or (II) particles are travelling on eccentric trajectories with a minimal radius smaller than $r_{\text{ISCO}}$. The latter case is highly unlikely as viscous forces spread anisotropies all along the accretion disk via diffusion. In Ref. Frank:2002bk , one can find a simplified linear model explaining that the cause of viscosity is the radial propagation of angular momentum from one orbit to another in the disk of matter. In this case, a final state of radial homogeneity is reached very rapidly. In a more realistic scenario, diffusion is driven by non-linear viscous forces studied in magnetohydrodynamical (MHD) simulations of rotating (and slightly magnetized) systems, as explained in Ref. Abramowicz:2011xu . In any case, the probability of finding matter on orbits within innermost stable circular orbit (ISCO) is greatly reduced, and in consequence, so is the matter density in those regions. For black holes with low accretion rates, this density drop could be sufficient to unscreen the scalarons within immediate vicinity of the black hole, as we will show in the following sections.

Using the arguments aforementioned and considering $\sigma\in[0,1)$ as the density contrast parameter, our rough prescription for the pressureless matter distribution is

[TABLE]

where $r_{g}=2GM/c^{2}$ is gravitational radius where black hole horizon is located. This coincides with the distribution suggested in Ref. Davis:2014tea when $\sigma=0$. For numerical reasons, the use of a smooth matter profile approximated by a hyperbolic tangent is more convenient

[TABLE]

where $a_{0}\gg 1$ and $\rho_{0}>0$. The shape of this matter distribution is a crude approximation of the results presented in Ref. Abramowicz:2011xu when $\sigma<10^{-4}$. As a consequence, we will setup all the wave equations for the scalar field using the Schwarzschild geometry

[TABLE]

and the matter density profile suggested in (3). However, it is convenient to change from Schwarzschild to tortoise coordinates $dr=\left(1-r_{g}/r\right)dx$

[TABLE]

The coordinate change from $r$ to $x$ and its inversion is discussed in detail in Appendix A.6. We neglect back-reaction of the scalaron dynamics on the background geometry for all the scalar-tensor theories explored in this project, and treat the metric as static in the scalaron equations of motion. Using (5), we can write the equation of motion of the scalar field $\Box\phi=V_{\text{eff}}^{\prime}(\phi)$ which appears in (10) and (25) as a spatially damped wave equation in 1+1 dimensions

[TABLE]

This choice of coordinates is sufficient for our purposes since the scalar solutions we seek do not need to cover the black hole’s interior $(r<r_{g})$. For simplicity of the implementation, we use units of $r_{g}=1$ in the code, but other choices can also be considered without difficulty.

III Scalar accretion in $f(R)$ theories

In this section we describe the two examples of $f(R)$ theories considered in this work: the Starobinsky and Hu-Sawicki models, where the Ricci scalar $R$ in the Einstein-Hilbert action is replaced by a function of $R$. We will consider equations of motion derived in the metric formalism. Alternative formalisms like the Palatini or the metric-affine mentioned in Refs. Sotiriou:2008rp ; Olmo:2005hd ; Sotiriou:2006qn ; BeltranJimenez:2017doy can change the number and/or the nature of the degrees of freedom that emerge. In particular, in Ref. Sotiriou:2008rp , we see that the Palatini formulation of $f(R)$ is equivalent to a $\omega_{0}=-3/2$ Brans-Dicke theory, and no new dynamical degrees of freedom appear. In all cases, one needs to ensure that the model remains ghost-free in the gravity sector and that there are no tachyonic modes for it to be viable. A recent discussion on ghosts in various formulations appeared in Ref. Koivisto:2013kwa .

III.1 Chameleons in Starobinsky and Hu-Sawicki model

We first consider $f(R)$ theories described, in the Jordan frame, by the action

[TABLE]

where $g_{\mu\nu}$ is the Jordan frame metric and $\psi$ are the matter fields. Varying the action with respect to the metric $g_{\mu\nu}$ we obtain equations of motion which replace Einstein’s equation in the $f(R)$ models; they are

[TABLE]

where $f_{R}\equiv\partial f/\partial R$. The two terms in the second line of the equation (8) contain fourth-order derivatives of the metric, a signal that a new degree of freedom emerges in the theory. This can be seen explicitly by taking the trace of the equation above

[TABLE]

which yields a second order equation of motion for the real field $f_{R}$ with a canonical kinetic term under the influence of an effective potential $V^{\prime}\equiv(2f-f_{R}R)/3$ and an external force term ${\cal F}\equiv-8\pi G\,T/3$ with $T\equiv T^{\mu}_{\mu}$. By defining $\phi\equiv f_{R}-1$ we can rewrite equation (9) simply as

[TABLE]

where prime denotes derivative with respect to $\phi$. Alternatively, one can explicitly see the emergence of the extra degree of freedom, usually dubbed “scalaron”, by mapping the action (7) into the Einstein frame, as described in Ref. Maeda:1988ab ; DeFelice:2010aj . Solving the equation (10) requires the knowledge of the potential $V(\phi)$ which is defined in a parametric form via

[TABLE]

or, integrating with a choice $V\left|{}_{R=0}\right.=0$ for a constant,

[TABLE]

When plotting the scalar potentials in Fig. 1, we observe that they are generally multi-valued, with turning points at field values where $f^{\prime\prime}(R)=0$. One must be aware of the branch choice when determining the curvature value $R$ and the effective force $V^{\prime}(\phi)$ for the field $\phi$ in the wave equation (10). The branch we are interested in is the one connected to the large curvature $R\rightarrow+\infty$ where Einstein gravity is recovered by screening. The particular models we study are defined by specific forms of $f(R)$

[TABLE]

which correspond to the Starobinsky – in Ref. Starobinsky:2007hu – and Hu-Sawicki – in Ref. Hu:2007nk – models respectively. Hence, by replacing the solution of (13) in $f_{\text{S}}$ or $f_{\text{HS}}$ (and their derivatives) the potentials are completely determined as functions of $\phi$ for every physical choice of parameters. A particular feature of (15) is that it has an apparent extra parameter compared to (14). However, it is entirely free, and we can reduce the number of parameters by considering the transformation $\alpha\rightarrow\lambda^{a}\alpha$, $\beta\rightarrow\lambda^{b}\beta$ and $R_{0}\rightarrow\lambda^{c}R_{0}$. $f_{\mathrm{HS}}$ is invariant under this transformation if

[TABLE]

for a given value of $n$. It is therefore possible to map any solution for one set of parameters to an equivalent one for a different set following (16), hence we fixed the parameter $\beta$ to be 1 throughout the rest of the paper and one can convert to other choices via these transformations. Our choices of $\beta$ and the crossover curvature scale $R_{0}$ are useful to compare these results with the solutions from the Starobinsky model. As summarized in Ref. Sotiriou:2008rp , $f(R)$ models of gravity in metric formalism must have $f^{\prime}(R)>0$ and $f^{\prime\prime}(R)>0$ to avoid ghost-like gravitons and tachyonic scalarons, respectively.

From equation (10), we can define an environmentally dependent effective potential $V_{\text{eff}}(\phi)$ that provides the same equations of motion in the regions of constant force ${\cal F}$ by

[TABLE]

One peculiar feature of the effective potential is that the extra term coming from the interaction with matter provides an external source term. As a result of this, the no-hair theorem in its usual form is in general not applicable. In this paper, the presence of matter with $T^{\mu}_{\mu}\neq 0$ (which excludes electromagnetic radiation) is not neglected.

One must consider any particular choice of model parameters for $f_{\text{HS}}$ and $f_{\text{S}}$ that could emerge from their corresponding renormalization group flows. Therefore, it is prudent to study the flow lines in parameter space by exploring the stability of the scalar wave equation for different choices of model parameters, even when we consider cases where the model does not match with current observations. The field solutions are screened in the same way as we described in section VI, henceforth $f_{\text{S}}$ and $f_{\text{HS}}$ scalarons can also be dubbed as “chameleons”.

III.2 Accreting chameleons in Hu-Sawicki and Starobinsky models

Now we present our results after evaluating the scalar equations of motion in (10). In order to proceed, we must first find the effective potentials corresponding to (14) with arbitrary densities of matter.

In the left panel of Fig. 1 we observe that the features of this potential correspond to what is usually called a chameleon field. Just as described in Ref. Brax:2008hh , the steepness and depth of the potential grows with the surrounding matter density in which the field propagates. The effective potential of the Hu-Sawicki model is plotted in the right panel of Fig. 1.

In both cases, we observe that the formation of curvature singularities do not require infinite energies to be achieved. In addition to that, we notice the existence of an equilibrium configuration for the field, the corresponding minimum in the effective potential $V_{\text{eff}}^{\prime}(\phi)=0$ is defined by

[TABLE]

In the presence of accreting matter outside $r_{\text{ISCO}}=3\,r_{g}$ the field is screened, which implies a very small value of it in this region. Inside $r_{\text{ISCO}}$, we can set the matter density in (3) to zero for now. Naively, one should not expect significant differences in our results calculated in the Jordan frame when compared with what is expected in the Einstein frame after conformal transformations: inside $r_{\text{ISCO}}$, the matter density is zero regardless of any value of the field; outside, the screening sets the conformal coupling to one. It is typically argued that there are changes in the metric and equations of motion of a test particle when these entities are rephrased in this frame: the derivatives of the conformal coupling might modify these entities in a non-negligible way. This is true in general, however these changes do not represent a significant contribution to the solutions we present since these modifications are always proportional to the first and second radial derivatives of the field, which is a smooth function outside the event horizon.

In order to find the screening value of $\rho_{0}$ in (3), we evaluate the equilibrium condition in (18) at the screened value of the field $\phi_{0}$ far away from $r_{\text{ISCO}}$ (which is very close to zero). Numerically, this is more convenient than (but still equivalent to) finding the equilibrium value of $\phi$ for a given value of $\rho_{0}$ from (18) since the dependence of $V^{\prime}(\phi)$ on the field is implicit.

The evolution of the system requires initial conditions. In what follows, we set the units of the field $\phi$ in terms of a pivot value $\phi_{*}$, which for $f(R)$ models is dimensionless and is chosen to be simply $\phi_{*}=1$. We consider a completely screened initial configuration of the field $\phi(r,t_{0})=-10^{-4}\phi_{*}$ which is originally spatially homogeneous, and assume $\dot{\phi}(r,t_{0})=0$ for the initial field velocity. Our choice of initial conditions is the same in our treatment of both $f(R)$ models we present here. In the left panel of Fig. 2, we see the accretion of the Starobinsky chameleon until it approaches to its static solution, choosing $n=1$, $\lambda=2.0$ and $R_{0}=10^{-2}/r_{g}^{2}$ as the model parameters to run the simulation.

The chameleon accretes around the horizon, then it oscillates slightly around the static solution. Gradients do not cancel outside the horizon for $r\leq r_{\text{ISCO}}$, which make the “hair” profile non-trivial. In addition to this, the screening outside $r_{\text{ISCO}}$ is not lost during accretion. Which shows that the equilibrium condition in (18) is being held.

Field evolution of the Hu-Sawicki chameleon is plotted in the right panel of Fig. 2, where we used the the same initial conditions. $n=1$, $\alpha=2.0$ and $R_{0}=10^{-2}/r_{g}^{2}$ were the model parameters chosen for the numerical evolution. Hu-Sawicki chameleon also accretes around the event horizon resulting in a non-trivial hair solution. We study the negative field branches of the potentials in Fig. 1 that roughly scale as $\phi^{1/\kappa}$ with $\kappa>2$, therefore the field becomes less massive in the regions where its amplitude deviates from the screened value.

The convergence into a hair solution as well as its shape are sensitive to choice of model parameters. Static solutions are found using the relaxation method described in Appendix A.5, where we also discuss all the details related to the numerical evolution. In Fig. 3, we evaluate the solutions of (6) in the static limit for different parameters of the Starobinsky model. We find the solutions for different values of $n$, while keeping $\lambda=1.0$ as a constant. Likewise, in Fig. 4 we evaluate the change of the static field profile when $\lambda$ varies and $n=1$ is kept as a constant. In the same way, we represent static solutions for different values of the density contrast parameter $\sigma$.

Different values of $\lambda$ and $R_{0}$ define how effective is the modification of gravity with respect to GR. In particular, $\lambda$ controls the depth and vertical extension of the effective potential and in consequence, it affects the existence and stability of the field solutions, while $R_{0}$ sets the crossover curvature scale. In particular, we chose $R_{0}=10^{-2}/r_{g}^{2}$ to have the same value throughout this paper.

As expected, we notice the reduction of the field amplitude in $r\in(r_{g},r_{\text{ISCO}}]$ for larger values of $\sigma$. However, it is not required to impose $\sigma=0$ to obtain a hairy solution. Furthermore, we calculate the static solutions corresponding to different choices parameters of the Hu-Sawicki model. The shape of the static solution for different values of $n$ can be found in Fig. 6. Additionally, steady-state solutions for different values of $\alpha$ and $\sigma$ are represented in Fig. 7 and Fig. 8, respectively. In all of these cases, the rest of the parameters were kept as constants. The dynamical results exposed show the accretion of the chameleon solutions and proof the existence of non-trivial stable solutions of the static version of (6) outside $r=r_{g}$ for the parameters we chose. The radial scalar flux ${\cal J}\equiv 4\pi r^{2}\,T[\phi]^{r}_{0}=4\pi r^{2}\,\phi_{,x}\phi_{,t}$ is represented in Fig. 9, showing no propagation outside $r_{\text{ISCO}}$.

We show an additional way to test that the static profiles in Figs. 3-7 are suitable representations of field configurations around the minimum of the effective potentials in Fig. 1. In the large curvature limit, one can notice that $V^{\prime}(\phi)\simeq R/3$, and so the equilibrium condition in (18) implies $R\simeq-8\pi G\,T$ just as it is in GR. As a consistency check, in Fig. 11 we tested the validity of the general relativistic limit in the case of a non-trivial value of the density contrast in the Hu-Sawicki model for a configuration that remains fully screened. This result remains valid for different values of the density contrast and further extends to the Starobinsky model.

So far, we discussed a few cases where we notice a smooth evolution into a non-trivial static solution. Nonetheless, the existence of static hair solutions is not enough to ensure smooth convergence to them. In Fig. 1, we find the field values where infinite curvature is reached are close to the minima of the potential for different choices of model parameters, which is consistent with the results in Refs. Frolov:2008uf ; Appleby:2009uf . Small field excursions from the minimum are sufficient to form curvature singularities outside the event horizon, and it can happen dynamically.

Equilibrium solutions for the field are defined by locating the “valleys” of $V_{\text{eff}}(\phi)$, which in these cases are not distant from reaching the point of infinite scalar curvature. Further modifications of the functional form of $f(R)$, for example the addition of a term proportional to $R^{2}$, can create a potential wall which will shield the field from reaching singularity. We will explore this alternative in the next section.

In Fig. 11, we show the dynamical formation of a curvature singularity outside the black hole’s event horizon. We chose $n=3$ and $\lambda=2.1$ as parameters for the Starobinsky model. Similarly, in Fig. 13, we pick $n=3$ and $\alpha=4.31$ to form a naked singularity in the Hu-Sawicki model. The profiles of field $\phi$ themselves do not diverge or show irregular behaviour while the singularities are formed, but the curvature $R$ goes singular under algebraic inversion $R\equiv R(\phi)$ at $\phi=0$. The singularity is rather weak, and most likely of integrable type, but it is nonetheless a naked singularity formed in evolution of the regular data, which shows cosmic censorship conjecture is violated in these models.

IV Resolving curvature singularities in $f(R)$ theories

In Figs. 11 and 13, we observed the formation of integrable curvature peaks in the Starobinsky and Hu-Sawicki models. In both cases, these were located in $r\in(r_{g};r_{\mathrm{ISCO}}]$ and appeared as a consequence of small field excursions from the potential minima reaching the field value corresponding to infinite curvature. In this section, we briefly discuss that the addition of an extra “mass” term in the functional form of $f(R)$

[TABLE]

is enough to remediate the divergencies appearing in both models for the same choices of model parameters chosen in the previous section, as suggested in Refs. Kobayashi:2008wc ; Appleby:2009uf . Here the singular point is avoided by adding an infinite barrier that regularizes the potential and its derivatives. In the same reference, it is possible to find constraints of the value of $\mu$ mostly related with the expected decay time of the scalarons in cosmological scales. Here we chose $\mu=10^{-6}$ to leave the low-curvature features of the model unaffected, while putting a cap at moderate curvature values to avoid numerical dynamical range issues which might cloud the discussion. Note that the extra mass term does not appear in the definition of $V^{\prime}(\phi)$ described in (9), regardless of any choice of $\mu$. In Fig. 14, we evaluated the dynamics of the scalar curvature in tortoise coordinates considering $n=3$ and $\lambda=2.1$ as parameters for the Starobinsky model and $n=3$, $\alpha=4.31$ for the Hu-Sawicki model. These choices produced unstable evolution truncated in the large curvature regime when using $f_{\mathrm{S}}$ and $f_{\mathrm{HS}}$. After adding the corrections in (19) and (20), the curvature peaks are limited and absorbed without reaching infinite values, with evolution towards an equilibrium configuration continuing without further inconveniences. However, even when the evolution of the curvature peaks is more benign these are still formed outside the event horizon. Thus, if such cusps are not observed, their existence in the model will pose constraints not only for specific choices of initial conditions or parameters, but for the entire subspace of the model parameters connected to the troublesome region by the renormalization flow. For the Hu-Sawicki model, the parameter scaling (16) that leaves action invariant will involve $\mu$, nonetheless it is still possible to find a family of parameters with similar curvature features starting from only one solution.

V Scalar accretion for the symmetron model

V.1 Action and equations of motion

We now turn to the symmetron model described by a scalar-tensor action of the form

[TABLE]

where $g_{\mu\nu}$ is now in the Einstein frame metric and $\psi$ are the matter fields minimally coupled to the Jordan metric $A^{2}(\phi)g_{\mu\nu}$. For simplicity, we will describe the dynamics of the scalar $\phi$ in the Einstein frame as in Ref. Davis:2014tea . The symmetron is then modelled with a potential

[TABLE]

and a coupling function

[TABLE]

where $\lambda$ and $\mu$ are positive coupling constants and $m_{\text{S}}^{2}$ is a high mass scale that suppresses any contributions higher or equal than $\mathcal{O}(\phi^{4}/m_{\text{S}}^{4})$. The equation of motion for the Einstein frame metric are

[TABLE]

while for the real scalar $\phi$ the equation of motion reads

[TABLE]

The effective potential $V_{\text{eff}}\equiv V(\phi)-TA(\phi)$ is defined as

[TABLE]

and, as in the case of $f(R)$ theories, it is dependent on the environment through the trace of stress-energy tensor $T\equiv T^{\mu}_{\mu}$. The coefficient of $\phi^{2}$ in equation (26) changes sign depending on the magnitude of $T$, and determines the shape of the symmetron effective potential. In case of an environment made solely of dust ($T=-\rho$), we can define $\rho_{\text{crit}}=\mu m_{\text{S}}^{2}/\hbar$ such that for energy densities $\rho<\rho_{\text{crit}}$ the effective potential becomes a shape of a mexican hat. For vanishing density, the two minima are at $\phi=\pm\phi_{*}$ with $\phi_{*}\equiv\sqrt{\mu/\lambda}$. In regions of high density, $\rho>\rho_{\text{crit}}$, there is one single minimum at $\phi=0$. Hence $A(\phi)=1$ for high densities and the field decouples from matter. This mechanism allows the symmetron model, with the proper parameters, to pass the solar system tests of GR. In the case of this model, the difference between the field dynamics in high density regions described in the Jordan and in the Einstein frame is not substantial because of the quadratic dependence of the conformal factor on $\phi/m_{\text{S}}$.

V.2 Accreting symmetrons

In this section we evaluate the dynamics of the accreting symmetrons. Previously, we discussed a simplified model equipped with spontaneously broken $\mathbb{Z}_{2}$ symmetry and an environmentally dependent mass. The shape of the potential and the effects of the coupling with matter are shown in Fig. 13. In regions of high matter overdensities, the only stable field configuration is $\phi=0$; however, as we discussed in the previous section, the field takes a non-zero vacuum expectation value in zones of lower energy density. In consequence, considering the mass distribution in (3), the initial radial profile used is

[TABLE]

which can be smoothed by using an hyperbolic tangent in the same manner as in the matter profile. An equivalent but more physically motivated way to choose the initial field configuration is to place it in equilibrium $V_{\text{eff}}^{\prime}(\phi)=0$ for a given matter distribution. Once again we will assume $\dot{\phi}(r,t=t_{0})=0$ for initial field velocity. With these initial conditions and choosing $\sigma=0$, we compute the evolution of the field profile for two different sets of model parameters: (a) $\lambda=10^{-2}/\hbar$, $\mu=8\times 10^{-4}/r_{g}^{2}$, $m_{\text{S}}^{2}=10^{-3}\hbar^{2}/r_{g}^{2}$ and (b) $\lambda=10^{2}/\hbar$, $\mu=8\times 10^{-1}/r_{g}^{2}$, $m_{\text{S}}^{2}=\hbar^{2}/r_{g}^{2}$. In Fig. 15, we calculated the evolution of the field towards equilibrium for both choices. Depending on parameters chosen, symmetrons do not always form non-trivial static hair solutions. The flux of the symmetron field ${\cal J}\equiv 4\pi r^{2}\,\phi_{,x}\phi_{,t}$ is shown in Fig. 16. All the ingoing scalar fluxes calculated here settle to zero smoothly after initial transient.

The difference in evolution lies in the contribution of the matter source to the effective potential: In the left panel of Fig. 16, we supressed the effects of the “external force” driven by the static matter density. Therefore this limit case is consistent with the standard no-hair theorem. However, that is not the case for the model depicted in the right panel, where the mass of the black hole is the same as the coupling parameter $m_{\text{S}}$. In accordance with our description in the previous subsection, we will explore the cases in which we can find non-trivial static solutions. From Fig. 13, we can recognize at least one of equilibrium field configurations which correspond to the different vacua for a given shape of the deformed “Mexican hat” potential. In Figs. 17-19, we compute the static solutions for different values of the model parameters.

VI Accretion of the Ratra-Peebles chameleons

VI.1 Equations of motion and setup

We now study the prototypical example of the chameleon screening, defined in the Einstein frame by the same scalar-tensor action as in (21), with the Ratra-Peebles potential

[TABLE]

and an approximate coupling function

[TABLE]

where $\gamma$ and $\varepsilon$ are positive model parameters and $m_{\text{C}}$ plays the role of a high mass scale where the screening is effective. Here we also choose to work in the Einstein frame being consistent with the procedures followed in Ref. Davis:2014tea . The equations of motion for both the metric and the scalar field are the same as in (24) and (25), but now the effective potential is given by

[TABLE]

This is another case of an effective potential dependent on the environment. Due to the term proportional to $T$, it is possible to find an equilibrium configuration for the field from $V^{\prime}_{\text{eff}}(\phi)=0$

[TABLE]

in the case of pressureless matter $(T=-\rho)$, we see that the field is suppressed in dense regions (with energies larger than $m_{\text{C}}$) and unscreened in regions with lower densities. The effective mass of the field is given by

[TABLE]

which becomes larger in dense environments. It is widely known that chameleon screening occurs in a different way as in the case of the symmetron as a consequence of the reduction of the field correlation length due to a larger effective mass of the field. In a low density environment, the minimum can be found at $\phi=\phi_{*}$ with $\phi_{*}\equiv\left(n\gamma^{n+4}\,m_{\text{C}}/(\varepsilon\sigma\rho_{0})\right)^{\frac{1}{n+1}}$, where $\sigma\rho_{0}$ is the density at $r<r_{\text{ISCO}}$.

VI.2 Accreting chameleons

Here we evaluate the accretion of the Ratra-Peebles chameleons in the background described in section II. In Fig. 21 we observe the shape of the potential found as a function of the field $\phi$, and the effects of the coupling with matter. The effective mass – i.e. the concavity of the effective potential – increases with the environmental matter density of dust.

We derive initial conditions compatible with (18) by replacing the matter distribution proposed in (3) into (31), avoiding the value of $\sigma=0$ to not have any divergencies in the initial field profile at $r<r_{\text{ISCO}}$. The constitution of the Ratra-Peebles model offers the possibility of finding the corresponding initial field configuration as a function of the surrounding matter density, the reverse process can also be coded without inconveniences. Assuming static initial conditions and $\sigma=0.01$, $\varepsilon=10^{2}$, $\gamma=0.3m_{\text{C}}$ and a surrounding dust density of $-T=\rho_{0}=3\times 10^{-2}m_{\text{C}}/r_{g}^{3}$, we find the evolution of the field from initial conditions towards a static solution in the left panel of Fig. 20. From this figure, we can notice that the evolution of the field is not significantly different from our results for the Starobinsky and the Hu-Sawicki models in Fig. 2, which is consistent with the effects of the chameleon screening in these models. Additionally, the ingoing field flux is depicted in the right panel of the same figure. Here the flux is regular and converges to the limit where there is no other source apart from the static matter distribution. Ingoing flux lines are represented around $x\left(r_{\text{ISCO}}\right)\approx 3.89$. We also present static solutions for different choices of the model parameters. In Figs. 22 - 24, we show the changes in shape of the hair profiles for different parameter choices. In this case, the tuning of the parameter in the runaway potential is sensitive to changes in the orders of magnitude: it is simple to suppress the whole contribution this part of the potential by accident, due to the $(n+4)$ power of $\gamma$.

VII Discussion

In this paper, we describe the dynamics of scalar accretion onto a Schwarzschild black hole in the presence of a static matter distribution modelled by (3). In particular, we studied the accretion of the extra scalar degrees of freedom appearing in two models of $f(R)$ gravity, the Starobinsky and Hu-Sawicki model, as well as in the symmetron model. Stable convergence to static scalar hair profiles results from varying parameters for each specific model. In the cases of the Starobinsky and Hu-Sawicki model, we can obtain dynamical chameleon solutions with singular curvature outside horizon without requiring an infinite energy budget. In the case of accreting symmetrons, it is not always possible to form a non-trivial static solution since it depends on the strength of the coupling with matter. More concretely, it depends on how large is the energy scale $m_{\text{S}}$ compared to the mass of the black hole. Our results for the field fluxes are included for all the cases we studied, along with the static solutions for different model parameters.

Even when the simulations of astrophysical rotating black holes suggest a large density contrast, we noticed that the formation of non-trivial static solutions does not require an absolute vacuum environment close to the black hole’s event horizon, where $r\in(r_{g},r_{\text{ISCO}}]$. We did not consider the effects of backreaction of the field in the spacetime solutions since these are small even during the formation of integrable naked singularities.

We acknowledge the progress made in Refs. Davis:2014tea ; Davis:2016avf , wherein approximate analytic expressions for scalar hair solutions in the case of rotating black holes are found. Additionally, this paper discusses the possibility of a non-negligible ratio between the radiated power from extra scalar sources and the quadrupole gravitational radiation in GR, which might be testable by the future generation of gravitational-wave detectors. In our approach, apart from calculating static solutions in different circumstances, we evaluated the scalar accretion dynamically in such a way that it is possible to converge to a hairy or a “bald” solution, depending on the model and its parameters. These results also motivate further explorations on the effect of fifth forces confined by screening, surrounded by a non-trivial matter profile for merging binary systems.

The existence of non-trivial field profiles has also been studied for scalars with non-canonical kinetic terms, such as in the Galileons studied in Refs. Babichev:2013cya ; Babichev:2016fbg . Those solutions are described in vacuum environments and do not accrete into the black hole’s horizon. The form of the equations of motion of Galileon is not semi-linear, and is harder to study numerically, as Galileons propagate with speeds that vary at different locations. Numerical implementation of a dynamical code designed to compute evolution of such scalar fields will be covered in a future project. We are also considering to extend the techniques developed here to full 3D scalar field scattering by a black hole.

Appendix A Numerical Implementation

A.1 Scalar field equations of motion

Equations of motion describing evolution of a scalar field $\phi$ with a (non-linear) self-interaction potential $V(\phi)$ and an external force term ${\cal F}$ propagating on a fixed background spacetime are described by a (semi-linear) PDE

[TABLE]

where $\Box$ denotes a covariant d’Alembert operator. For a spherically symmetric black hole described by the Schwarzschild metric

[TABLE]

where $d\Omega^{2}$ is the metric on a unit sphere, and the metric function $g(r)$ is

[TABLE]

the left hand side of the equation of motion is simply

[TABLE]

This can be reduced to a one-dimensional wave equation with constant propagation speed by introducing the tortoise coordinate $x$ by $\partial_{x}=g(r)\,\partial_{r}$. With this re-definition, the scalar field equation of motion reads

[TABLE]

Explicit form of the tortoise coordinate $x$ for the Schwarzschild spacetime can be obtained by integrating

[TABLE]

Tortoise coordinate $x$ is vastly preferable for numerical integration of the wave equation over areal coordinate $r$ since the characteristic speed is constant on the sampled time slice, but the added difficulty with this choice is that accurate $r(x)$ inversion is quite not-trivial numerically, as detailed in Appendix A.6.

The standard numerical evolution scheme would involve first-order Hamiltonian dynamical system

[TABLE]

However, as we will see in the following Section, it is easier to handle absorbing boundary conditions if we rewrite equations of motion in a flux-conservative form by introducing auxiliary variables $u\equiv\partial_{t}\phi$ and $v\equiv r^{2}\partial_{x}\phi$, so that equations of motion become

[TABLE]

The first equation is the identical rewrite of the wave equation (37), while the second is the integrability condition requiring that the partial derivatives of $\phi$ commute.

A.2 Absorbing boundary conditions

The scalar field degree of freedom $\phi$ asymptotes to a free field evolution near horizon (where $g\rightarrow 0$), and a massive field evolution far away from the black hole (where $V^{\prime}(\phi)\rightarrow{\cal F}$). Physically, excitations in $\phi$ take infinite amount of time $t$ to reach both boundaries, yet truncating or compactifying the evolution domain for numerical purposes will inevitably lead to spurious reflections unless special care is taken. The best technique to avoid spurious reflections is to introduce absorbing boundary conditions via Perfectly Matched Layers (PMLs) as described in Ref. Johnson:wt , which damp the solution at the boundaries while guaranteeing identically vanishing reflection coefficient at the absorption layer. This is achieved by analytic continuation of the equations of motion into the complex domain

[TABLE]

which turns the oscillatory travelling waves $e^{ikx-i\omega t}$ into exponentially decaying functions of $x$ instead. To make attenuation length independent of $\omega$, frequency dependent contour deformation $f^{\prime}=\gamma(x)/\omega$ is chosen and $i/\omega$ is transformed back into explicit integration operator in the time domain. Applying this idea to the scalar field equations of motion in flux-conservative form (A.1) for an arbitrary damping function $\gamma(x)$, we obtain

[TABLE]

To turn the inverse time evolution operator $\partial_{t}^{-1}$ into a differential equation form, introduction of a third auxiliary variable $w$ is in order. With re-definition $u\rightarrow u+w$, the non-reflecting PML equations of motion then become

[TABLE]

The damping function $\gamma(x)$ can be quite arbitrary, but it should have compact support near the boundaries to not affect the evolution in the interior, and have sufficient support and magnitude to absorb the impinging waves which hit the boundary during the expected evolution.

A.3 Spectral basis

As the scalar field is usually quite stiff and does not form shocks in the course of evolution, the method of choice to evaluate derivative operators is spectral, as described in Ref. Boyd . Compactifying the tortoise coordinate $x$ on a scale $\ell$

[TABLE]

and introducing a Chebyshev basis on interval $y\in[-1,1]$

[TABLE]

we arrive at the spectral representation of the solution

[TABLE]

truncated to a finite number of modes. While Galerkin method to discretize equations of motion can be employed, the simplest method to evaluate derivative operators is pseudo-spectral, where equations of motion are solved on a Gauss-Lobatto grid

[TABLE]

One does not have to explicitly find coefficients $c_{n}$ to evaluate the derivative operators of a function $\phi(x)$ sampled on a collocation grid $x_{i}$. Instead, derivative operators like ${\cal D}_{ij}$ and ${\cal L}_{ij}$ can be found in advance by solving linear matrix equations

[TABLE]

and so on for every basis function $T_{n}$ evaluated at all nodes $x_{i}$.

A.4 Gauss-Legendre integrator

Packing the scalar field variables $\phi,u,v,w$ evaluated at the collocation grid points $x_{i}$ into a state vector $\mathbf{y}\equiv\{\phi(x_{i}),u(x_{i}),v(x_{i}),w(x_{i})\}$, the wave equation (43) reduces to an autonomous dynamical system

[TABLE]

which can be integrated by an implicit Runge-Kutta method, as presented in Ref. Butcher:1964fx

[TABLE]

where the trial directions $\mathbf{g}^{(i)}$ are defined by

[TABLE]

Particularly accurate choice of coefficients for a time integrator corresponds to a Gauss-Legendre quadrature, where the trial directions are evaluated at the zeroes of the (shifted) Legendre polynomial

[TABLE]

with coefficients $a^{i}_{j}$ and $b_{j}$ set by

[TABLE]

The resulting time integration method is A-stable and symplectic for Hamiltonian problems, and is extremely easy to implement using a simple iterative scheme.

A.5 Static solver

Static configurations of the field $\phi$ have $\partial_{t}\phi=0$ and can be found by solving a (semi-linear) elliptical problem

[TABLE]

One can improve a trial solution $\bar{\phi}$ using Newton’s method by linearizing $\phi=\bar{\phi}+\delta\phi$ and solving

[TABLE]

which translates the residual {\cal R}=-{\cal L}\bar{\phi}+g\Big{(}V^{\prime}(\bar{\phi})-{\cal F}\Big{)} into a correction $\delta\phi$ by solving a set of linear equations

[TABLE]

With the basis as chosen in the last section, this scheme converges to machine precision in about 16 iterations or so for most of the potentials.

A.6 Inverting tortoise coordinate

Accurately inverting Schwarzschild tortoise coordinate

[TABLE]

to yield areal coordinate $r$ as a function of $x$ turns out to be a rather non-trivial task, despite appearances. The problem is that asymptotic for large positive $x$, where $r\simeq x-2M\ln\left(x/2M-1\right)$, and for large negative $x$, where $r\simeq 2M$ with exponentially suppressed metric function $\ln g\simeq x/2M-1$, have vastly different derivatives with respect to $x$ (which hampers numerical schemes like Newton’s method), and no closed form algebraic inverse.

A trick that works for the entire usable range of $x$ is to solve for an approximation variable $q\simeq x-2M$ instead

[TABLE]

which (unlike $x$) is easily invertible to yield $r$

[TABLE]

and can be readily found by Newton’s method iterating $q\rightarrow q+\delta q$ with

[TABLE]

as the derivative

[TABLE]

is of order one on the entire domain of definition of $x$. One still has to be careful to avoid numerical overflows in the exponents or catastrophic loss of precision when taking logarithms of one plus a small number, which can be achieved by evaluating

[TABLE]

in different limits.

Acknowledgements.

We would like to thank Andrew DeBenedictis, Michael Desrochers, Claire Maulit, Justine Munich, Levon Pogosian, Misao Sasaki and Alexander Vikman for their valuable comments and discussions on the previous drafts of this paper. This project was partly funded by the Discovery Grants program of the Natural Sciences and Engineering Research Council of Canada and it was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293. AZ is supported in part by the Bert Henry Memorial Entrance Scholarship at SFU.