On the stability of scalarized black hole solutions in scalar-Gauss-Bonnet gravity

TL;DR

This paper investigates the stability of scalarized black hole solutions in scalar-Gauss-Bonnet gravity, revealing that the form of the scalar coupling determines whether solutions are stable or unstable.

Contribution

It clarifies how higher-order scalar-Gauss-Bonnet couplings influence the nonlinear stability of scalarized black holes.

Findings

Quadratic couplings lead to unstable scalarized black holes.

Exponential couplings can produce stable scalarized black hole solutions.

Higher-order couplings control the nonlinear quenching of tachyonic instabilities.

Abstract

Scalar-tensor theories of gravity where a new scalar degree of freedom couples to the Gauss-Bonnet invariant can exhibit the phenomenon of spontaneous black hole scalarization. These theories admit both the classic black holes predicted by general relativity as well as novel hairy black hole solutions. The stability of the hairy black holes is strongly dependent on the precise form of the scalar-gravity coupling. A radial stability investigation revealed that all scalarized black hole solutions are unstable when the coupling between the scalar field and the Gauss-Bonnet invariant is quadratic in the scalar, whereas stable solutions exist for exponential couplings. Here we elucidate this behavior. We demonstrate that, while the quadratic term controls the onset of the tachyonic instability that gives rise to the black hole hair, the higher-order coupling terms control the nonlinearities that quench that instability, and hence also control the stability of the hairy black hole solutions.

Equations

$$S=\frac{1}{2}\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R-\frac{1}{2}{g}^{ab}{\phi }_{;a}{\phi }_{;b}+f(\phi )\mathcal{G}\right],$$

$$\square \phi =-f_{,\phi }(\phi )\mathcal{G},$$

$$R_{ab}-\frac{1}{2}{g}_{ab}R=T_{ab}$$

$$f(\phi )\equiv \frac{1}{8}\bar{\eta }{\phi }^2+\frac{1}{16}\bar{\zeta }{\phi }^4,$$

$$f(\phi )=\frac{{\lambda }^2}{12}\left[1-\exp (-3{\phi }^2/2)\right]$$

$$\bar{\eta }={\lambda }^2,\bar{\zeta }=-\frac{3}{2}{\lambda }^2,$$

$$f_{,\phi }=\frac{\phi }{4}\left(\bar{\eta }+\bar{\zeta }{\phi }^2\right)=0.$$

$${\phi }_{\pm }=\pm \sqrt{|\bar{\eta }/\bar{\zeta }|}.$$

$$\square \phi =-f_{,\phi }(\phi )\mathcal{G}=-\frac{\phi }{4}\left(\bar{\eta }+\bar{\zeta }{\phi }^2\right)\mathcal{G},$$

$$\phi =\frac{1}{r}\sum _{\ell m}{\bar{\sigma }}_{\ell m}(r)Y_{\ell m}(\theta ,\varphi ).$$

$$\sigma$$

$$\eta$$

$${\sigma }^{\prime \prime }$$

$$\sigma ={\sigma }_{\mathrm{h}}+({\sigma }_{\mathrm{h}}-3\eta {\sigma }_{\mathrm{h}}-3\zeta {\sigma }_{\mathrm{h}}^3)(\rho -1).$$

$$\sigma ={\sigma }_0+\delta \sigma ,$$

$$\frac{{\mathrm{d}}^2\delta \sigma }{\mathrm{d}{\rho }_{\ast }^2}-V_{\mathrm{eff}}(\rho )\delta \sigma =0,$$

$$V_{\mathrm{eff}}\equiv \left(1-\frac{1}{\rho }\right)\left[\frac{\ell (\ell +1)}{{\rho }^2}+\frac{1}{{\rho }^3}-\left(\frac{3\eta }{{\rho }^6}+\frac{9\zeta }{{\rho }^8}{\sigma }_0^2\right)\right].$$

$${\int }_1^{\infty }\frac{V_{\mathrm{eff}}}{1-1/\rho }\mathrm{d}\rho$$

$$\mathrm{d}{s}^2=$$

$$+\exp [2\Lambda (r)+\epsilon F_r(t,r)]\mathrm{d}{r}^2+r^2\mathrm{d}{\Omega }^2,$$

$$\phi ={\phi }_0(r)+\epsilon \frac{{\phi }_1(t,r)}{r}.$$

$$\exp (2\Phi )$$

$$\exp (-2\Lambda )$$

$${\phi }_0$$

$${\phi }_{0,\mathrm{h}}^{\prime }=-\frac{r_{\mathrm{h}}}{{\phi }_{0,\mathrm{h}}\left(\eta +\zeta {\phi }_{0,\mathrm{h}}^2\right)}\left[1-\sqrt{1-\frac{6{\phi }_{0,\mathrm{h}}^2}{r_{\mathrm{h}}^4}{\left(\eta +\zeta {\phi }_{0,\mathrm{h}}^2\right)}^2}\right].$$

$$6{\phi }_{0,\mathrm{h}}^2{\left(\eta +\zeta {\phi }_{0,\mathrm{h}}^2\right)}^2<r_{\mathrm{h}}^4.$$

$$\exp (2\Phi )$$

$$\exp (-2\Lambda )$$

$${\phi }_0$$

$$h(r)\frac{{\partial }^2{\phi }_1}{\partial t^2}-\frac{{\partial }^2{\phi }_1}{\partial r^2}+k(r)\frac{\partial {\phi }_1}{\partial r}+p(r){\phi }_1=0.$$

$$x\equiv 1-r_{\mathrm{h}}/r,$$

$$\varphi_{1}(x)=\left\{\begin{array}[]{ll}e^{i\omega r_{*}}&r_{*}\to\infty~{}~{}(x\to 1)\\ e^{-i\omega r_{*}}&r_{*}\to-\infty~{}~{}(x\to 0)\end{array}\right.,$$

$$W={\left[{\phi }_1^{(-)}\frac{\mathrm{d}{\phi }_1^{(+)}}{\mathrm{d}x}-{\phi }_1^{(+)}\frac{\mathrm{d}{\phi }_1^{(-)}}{\mathrm{d}x}\right]}_{x=x_m}$$

Full text

I Introduction

A century after the inception of general relativity (GR), we have entered the era of gravitational wave astronomy. The LIGO/Virgo Collaboration has already observed ten binary black hole (BH) mergers and one binary neutron star merger Abbott et al. (2018a), and this number is only expected to grow as the sensitivity is increased and future detectors (such as LIGO-India and KAGRA) come online. This new observational window offers us the unprecedented opportunity to test gravity on new distance and energy scales using some of the most extreme objects in the Universe Berti et al. (2015); Yunes et al. (2016); Berti et al. (2018a, b); Barack et al. (2018). Indeed, the events already observed have provided new bounds on modified gravity theories and important consistency tests Abbott et al. (2016); Sakstein and Jain (2017); Creminelli and Vernizzi (2017); Baker et al. (2017); Ezquiaga and Zumalacárregui (2017); Crisostomi and Koyama (2018); Langlois et al. (2018); Dima and Vernizzi (2018); Abbott et al. (2018b).

Astrophysical BHs in GR are simple objects characterized by two numbers: their mass and spin. Because of their simplicity, they are attractive probes for testing GR. Any observational test of gravity must compare the predictions of GR with those coming from competing theories. Without alternative predictions, one cannot verify whether GR is correct, or even quantify the amount by which GR is the preferred theory. This is where the challenges typically arise: the vast majority of well-motivated modifications of GR are subject to no-hair theorems that preclude the existence of new BH charges. In particular, many of these models contain an additional scalar degree of freedom $\varphi$ Berti et al. (2015); Sotiriou (2015a), and the no-hair theorems preclude the existence of BHs with some new scalar charge $Q$ Hawking (1972); Mayo and Bekenstein (1996); Sotiriou and Faraoni (2012); Hui and Nicolis (2013); Sotiriou and Zhou (2014a); Sotiriou (2015b); Herdeiro and Radu (2015). A notable exception to this rule are theories that include a coupling of the scalar field to the Gauss-Bonnet invariant ${\cal G}=R^{2}-4R_{ab}R^{ab}+R_{abcd}R^{abcd}$, where $R$, $R_{ab}$ and $R_{abcd}$ are the Ricci scalar, the Ricci tensor and the Riemann tensor, respectively. Such scalar-Gauss-Bonnet (sGB) couplings arise in the low-energy effective field theory (EFT) derived from string theory, and it has long been known that they give rise to non-Schwarzschild and non-Kerr BHs (see e.g. Mignemi and Stewart (1993); Kanti et al. (1996)). A linear coupling between a scalar field and the Gauss-Bonnet invariant features prominently in the EFT of shift-symmetric scalars Sotiriou and Zhou (2014a), since the Gauss-Bonnet term is a topological invariant and a total divergence in four dimensions. Again, studies have found new hairy BH solutions in these theories Sotiriou and Zhou (2014a, b); Ayzenberg and Yunes (2014); Maselli et al. (2015a); Antoniou et al. (2018a, b). The fact that LIGO/Virgo has observed BH merger events consistent with GR implies that such couplings are necessarily small Witek et al. (2018).

Recently, a new possibility was pointed out and named “spontaneous BH scalarization”: certain sGB theories admit both the BH solutions of GR and hairy BH solutions Doneva and Yazadjiev (2018a); Silva et al. (2018). This phenomenon can occur in theories where a scalar field is coupled to the Gauss-Bonnet invariant and the coupling function respects $\mathbb{Z}_{2}$ symmetry and vanishes for some constant $\varphi_{0}$. The last condition guarantees that GR BHs are admissible solutions. The coupling with the Gauss-Bonnet invariant acts as an effective mass term for the scalar perturbations around these solutions. When the BH mass lies within a certain interval, this effective mass term is negative in parts of the BH exterior, triggering a tachyonic instability and producing a nonzero scalar charge. The scalarized solutions exist in BH mass bands whose onset coincides with the tachyonic instability, and whose termination is due to regularity conditions on the horizon that arise from nonlinear effects Doneva and Yazadjiev (2018a); Silva et al. (2018). Similar scalarization phenomena have been studied also for neutron stars Silva et al. (2018); Doneva and Yazadjiev (2018b), Reissner-Nordström and Kerr BHs Cardoso et al. (2013a, b); Doneva et al. (2018); Herdeiro et al. (2018); Brihaye et al. (2019), and scalar-tensor gravity coupled with Born-Infeld electrodynamics Stefanov et al. (2008); Doneva et al. (2010).

There is no a priori guidance for the functional dependence of the coupling function $f(\varphi)$. Reference Silva et al. (2018) focused on the quadratic coupling $f(\varphi)\sim\varphi^{2}$, as this is the simplest case where the tachyonic instability should be present and the leading-order term is expected to control the onset of the instability. Reference Doneva et al. (2018) focused on the exponential coupling $f(\varphi)\sim\exp(\beta\varphi^{2})$ instead. Recent studies suggest that scalarized BH solutions are unstable under radial perturbations for the quadratic coupling, while solutions within the exponential coupling model have better stability properties Blázquez-Salcedo et al. (2018).

This paper is concerned with understanding the nature of the instability of the quadratic coupling function. In particular, we will show that the radial instability of the quadratic model is directly linked to the fact that, in this model, the scalar field equation is linear in the scalar. This implies that gravitational backreaction is crucial for quenching the tachyonic instability that leads to scalarization, and that backreaction determines the properties of the scalarized solution in this model. Here we find that introducing nonlinearity in the scalar provides a different quenching mechanism for the tachyonic instability, changes the properties of the scalarized solutions, and removes the radial instability. The simplest setup to demonstrate these points is a theory where the quadratic coupling is augmented by a quartic term111The $\varphi^{4}$ correction to the coupling function introduces nonlinearity in $\varphi$ in the scalar’s equation and is the leading order term with this property in the exponential coupling case. Note that from an EFT perspective, augmenting the coupling function is not the most natural method for stabilizing the solutions. However, our goal here is to understand the role and properties of the coupling function. We will discuss EFT considerations in a forthcoming publication Macedo et al. (2019)..

The plan of the paper is as follows. In Sec. II we briefly review sGB gravity, the necessary conditions for the existence of scalarized BH solutions, and our chosen coupling functions. In Sec. III we study quartic sGB gravity in the decoupling limit, compute scalar field bound states and investigate their stability. In Sec. IV we obtain full nonlinear BH solutions in this theory and discuss their stability under radial perturbations. In Sec. V we summarize our findings.

II Scalar-Gauss-Bonnet gravity

In sGB gravity, a real, massless scalar field is coupled to gravity through the Gauss-Bonnet invariant ${\cal G}$. The action of sGB gravity is

[TABLE]

where $\varphi$ is the scalar field and $g_{ab}$ is the spacetime metric. We use geometrical units, such that $8\pi G=c=1$. The field equation for the scalar field in sGB gravity is

[TABLE]

while the equation for the spacetime metric is

[TABLE]

where $T_{ab}$ is the sum of the matter stress-energy tensor (which is vanishing for BH solutions) and an effective stress-energy tensor which depends on $f_{,\varphi}(\varphi)$ Antoniou et al. (2018a). Different choices of the function $f(\varphi)$ correspond to different sGB gravity theories. In particular, $f(\varphi)\sim\exp(\alpha\varphi)$ (where $\alpha$ can be different depending on the specific stringy scenario) corresponds to Einstein-dilaton Gauss-Bonnet (EdGB) gravity Mignemi and Stewart (1993); Kanti et al. (1996); Maeda et al. (2009); Torii and Maeda (1998); Pani and Cardoso (2009); Yunes and Stein (2011); Pani et al. (2011); Yagi et al. (2012); Yagi (2012); Ayzenberg and Yunes (2014); Maselli et al. (2015b); Kleihaus et al. (2016); Blázquez-Salcedo et al. (2016), which can arise in the low-energy effective actions of some string theories Gross and Sloan (1987); Metsaev and Tseytlin (1987); $f(\varphi)\sim\varphi$ corresponds to shift-symmetric sGB gravity, which is invariant under $\varphi\rightarrow\varphi+{\rm constant}$, and is so far the only known shift-symmetric scalar-tensor theory with second-order field equations to allow for asymptotically flat, hairy BH solutions Sotiriou and Zhou (2014a, b); Maselli et al. (2015a).

Remarkably, EdGB and shift-symmetric sGB gravity do not admit Schwarzschild BH solutions: all static, spherically symmetric BH solutions in this theory have nontrivial scalar field configurations. This feature, however, is not shared by all sGB gravity theories. As shown in Doneva and Yazadjiev (2018a); Silva et al. (2018), sGB gravity admits the BH solutions of GR (with a constant scalar field) if $f_{,\varphi}(\varphi_{0})=0$ for some constant $\varphi_{0}$, and it also admits BH solutions with nontrivial scalar field configurations if $f_{,\varphi\varphi}{\cal G}<0$. In these theories, the $f_{,\varphi\varphi}{\cal G}$ term in the field equations acts as a negative mass term, triggering a tachyonic instability222See also Myung and Zou (2018) for an alternative interpretation in terms of the Gregory-Laflamme instability., which can in principle lead to the development of “scalar hair”. This process of spontaneous scalarization is analogous to that studied in compact stars in scalar-tensor gravity Damour and Esposito-Farèse (1993, 1996), but, crucially, it does not rely on any coupling with matter, and therefore it could potentially be tested through the observation of gravitational waves from binary BH mergers.

Two examples of sGB gravity theories satisfying the conditions for spontaneous scalarization have been studied: quadratic sGB gravity Silva et al. (2018), in which $f(\varphi)\sim\varphi^{2}$, and exponential sGB gravity Doneva and Yazadjiev (2018a), where $f(\varphi)\sim 1-\exp(-3\varphi^{2}/2)$ 333Note that our notation is consistent with that of Silva et al. (2018) but different from the notation of Doneva and Yazadjiev (2018a), since the scalar field defined in Doneva and Yazadjiev (2018a) differs from our scalar field by a factor $2$..

In this work we will study quartic sGB gravity, with a coupling term of the form

[TABLE]

where $\bar{\eta}$, $\bar{\zeta}$ are coupling constants with dimensions of [length]2, and the numerical factors are chosen for convenience. The scalar field is dimensionless, while ${\cal G}$ has units of [length]-4. When $\bar{\zeta}=0$, we obtain the quadratic sGB gravity theory considered in Silva et al. (2018), which allows for spontaneous scalarization when $\bar{\eta}>0$. For small values of the scalar field, the exponential sGB gravity studied in Doneva and Yazadjiev (2018a)

[TABLE]

reduces to quartic sGB gravity with

[TABLE]

plus ${\cal O}(\varphi^{6})$ terms.

Since the field equations (3) reduce to Einstein’s equations when $f_{,\varphi}(\varphi)=0$, quartic sGB gravity admits the GR solutions provided that

[TABLE]

Therefore, it always admits the GR solutions with $\varphi\equiv 0$, and if $\bar{\eta}\bar{\zeta}<0$ it admits two additional, real-valued, constant scalar field solutions $\varphi\equiv\varphi_{+}$ and $\varphi\equiv\varphi_{-}$, where

[TABLE]

III Decoupling limit of quartic scalar-Gauss-Bonnet gravity

Our goal in this section is to understand whether a quartic term in the coupling function $f(\varphi)$ i.e., a nonzero value of $\bar{\zeta}$ can stabilize static, spherically symmetric BH solutions, which are known to be unstable in the quadratic case. Since quadratic sGB gravity only admits scalarized solutions if ${\bar{\eta}}\sim f_{,\varphi\varphi}>0$, in the following we shall assume ${\bar{\eta}}>0$, as in Silva et al. (2018); Blázquez-Salcedo et al. (2018).

To begin, we consider the decoupling limit in which we neglect the backreaction from the metric. Thus, we study the scalar field equation (2)

[TABLE]

on a fixed Schwarzschild background with mass $M$. As discussed in Doneva and Yazadjiev (2018a); Silva et al. (2018), this class of theories admits two kinds of solutions: the GR solutions, i.e. a constant scalar field with $\varphi\equiv 0,\varphi_{\pm}$; and the scalarized solutions, in which the scalar field has a nontrivial configuration.

III.1 Static bound-state solutions

As a first step, we look for static, bound-state solutions of quartic sGB gravity in the decoupling limit. As in Ref. Silva et al. (2018), we consider a time-independent scalar field and we expand it in spherical harmonics in standard Schwarzschild coordinates:

[TABLE]

Due to the nonlinearity introduced by the quartic term in Eq. (4), the wave equation is only separable for spherically symmetric configurations ($\ell=0$). Defining $\bar{\sigma}=\bar{\sigma}_{00}$ and introducing the dimensionless variables

[TABLE]

we obtain the following nonlinear differential equation:

[TABLE]

where primes denote derivatives with respect to $\rho$.

Introducing the tortoise coordinate $\rho_{\ast}=\rho+\log(\rho-1)$, this equation becomes a Schrödinger-like equation with a nonlinear potential. We cannot straightforwardly apply results from quantum mechanics to determine when Eq. (11) admits bound-state solutions, as done e.g. in Refs. Dotti and Gleiser (2005); Cardoso et al. (2009, 2013b) using criteria derived in Buell and Shadwick (1995), and therefore we must study Eq. (11) numerically. We start the integrations close to the event horizon $\rho=1$ (typically at $\rho=1+10^{-5}$) using a power series solution for $\sigma$ valid in the near-horizon region. To leading order we have

[TABLE]

We then integrate out to a large $\rho$ (typically $\approx 10^{5}$) for given values of $(\eta,\,\zeta)$ and an arbitrary value $\sigma_{\rm h}\equiv\sigma(1)$. For each choice of the coupling constants $\eta$, $\zeta$, and for each choice of $\sigma_{\rm h}$, Eq. (11) admits a unique solution with $\sigma^{\prime}_{\rm h}\equiv\sigma^{\prime}(1)=\sigma_{\rm h}(1-3\eta-3\zeta\sigma_{\rm h}^{2})$, $\varphi_{\rm h}=\sigma_{\rm h}$ and $\varphi^{\prime}_{\rm h}=-3\sigma_{\rm h}(\eta+\zeta\sigma_{\rm h}^{2})$. We find that the scalar field diverges as $\rho\rightarrow\infty$ if and only if $\varphi_{\rm h}^{\prime}/\varphi_{\rm h}>0$, i.e. if $\eta+\zeta\sigma_{\rm h}^{2}>0$. Since we assume that $\eta>0$, this condition is always satisfied if $\zeta>0$, while if $\zeta<0$ the condition is satisfied for $\varphi_{-}<\varphi_{\rm h}<\varphi_{+}$ [see Eq. (8)].

Since we are interested in scalarized solutions with the same asymptotic behavior as in GR, we require that the scalar field vanishes at infinity. We find that this condition can only be enforced for $\eta$ larger than a threshold value $\eta_{\rm thr}=0.726$; for each $\eta>\eta_{\rm thr}$, there is a discrete set of values of $\sigma_{\rm h}$ which correspond to the solutions satisfying the boundary condition at infinity. This value, i.e. $\bar{\eta}_{\rm thr}/M^{2}=4\eta_{\rm thr}=2.904$, coincides with the first (zero-node) eigenvalue in the quadratic theory of Ref. Silva et al. (2018). We focus on the nodeless solution because previous work Blázquez-Salcedo et al. (2018) showed that this is the only stable scalarized solution in the exponential theory. Therefore it is natural to ask whether the nodeless solution is stable in the (simpler) quartic theory.

The results of the integration of Eq. (11) are shown in Fig. 1, where we plot the scalar field at the horizon, $\varphi_{\rm h}$, compatible with the boundary conditions, in the case of $\zeta=-(3/2)\eta$ [corresponding to the choice in Eq. (6)], for the solution with no nodes. Note that quartic sGB gravity is symmetric under $\varphi\rightarrow-\varphi$, so we only show solutions with $\varphi>0$. For $\eta<\eta_{\rm thr}$ the only solution is $\varphi_{\rm h}=0$, and the scalar field is zero everywhere. A scalarized solution appears for $\eta>\eta_{\rm thr}$. As $\eta$ increases, the value of $\varphi_{\rm h}$ for the scalarized solution also increases, and it tends to the limit $\varphi_{+}$ as $\eta\rightarrow\infty$. The same qualitative behavior occurs for different (negative) values of $\zeta$.

In the same figure we show the corresponding curve for the quadratic theory Silva et al. (2018). In this case the scalar field equation in the decoupling limit is linear, therefore only a discrete set of values of $\eta$ fulfills the boundary conditions, and the zero-node solution corresponds to $\eta=\eta_{\rm thr}$. For this value all choices of $\varphi_{\rm h}$ are equivalent, since the solution of a linear equation is defined modulo an overall multiplicative constant.

III.2 Linear stability analysis in the decoupling limit

Let us now analyze the stability of the static, spherically symmetric solutions discussed in Sec. III.1. Let

[TABLE]

were $\sigma_{0}$ is a solution of Eq. (11) (e.g. found numerically as in Sec. III.1). Substituting this perturbation in Eq. (11) and neglecting ${\cal O}(\delta\sigma^{2})$ terms, we find the linear equation

[TABLE]

with effective potential

[TABLE]

Following Buell and Shadwick (1995), a sufficient (but not necessary) condition for instability is

[TABLE]

where $\varphi_{0}=\delta\sigma_{0}/\rho$ and the integral must be computed numerically.

In the quadratic sGB theory ($\zeta=0$) the solution is unstable for $\eta>5/6\simeq 0.83$. Remarkably, in this case the effective potential does not depend on $\varphi_{0}$. Therefore, when the Schwarzschild solution is unstable – the instability eventually leading to spontaneous scalarization – the bound state solution is also unstable. This simple qualitative reasoning suggests that, as shown in Blázquez-Salcedo et al. (2018), scalarized solutions in the quadratic sGB theory are always unstable.

Eq. (16) also shows that when $\zeta<0$ the contribution of the integral is positive, and therefore it tends to stabilize BH solutions. Moreover the integral term vanishes as $\varphi_{0}\to 0$. This suggests that in the quartic theory, for certain values of the coupling constants $(\eta,\,\zeta)$ the Schwarzschild solution is unstable, while the scalarized solution is not. This is consistent with the results found in Blázquez-Salcedo et al. (2018) for the exponential theory, which is equivalent to the quartic theory with $\zeta=-(3/2)\eta$ if we ignore terms of order ${\cal O}(\varphi^{6})$.

III.3 Numerical results

To assess the stability of the scalarized solutions under radial perturbations, we computed the effective potential and the integral in Eq. (16) numerically in the case $\zeta=-(3/2)\eta$ using the bound-state solutions $\varphi_{0}$ corresponding to the curve in Fig. 1. We find that the condition (16) (which is just a sufficient condition for instability) is never satisfied for the bound-state solutions. Looking at the minimum of the effective potential, we find that it is negative for $\eta_{\rm thr}=0.726<\eta<0.86$, while it is positive for $\eta>0.86$. This is an indication that the bound state solutions, at least for $\eta>0.86$, are linearly stable under radial perturbations. In Fig. 1 we mark solutions for which $V_{\rm eff}$ is not positive-definite by a dashed line, and those for which $V_{\rm eff}$ is positive everywhere by a solid line. In the GR limit ($\eta=0$) the theory reduces to a scalar field minimally coupled to gravity, and no-hair theorems Hawking (1972); Mayo and Bekenstein (1996); Sotiriou and Faraoni (2012) require $\varphi$ to be (any) constant, as indicated by the wiggly line.

These results are in qualitative agreement with those for the exponential theory, as can be seen from a comparison with Fig. 2 of Blázquez-Salcedo et al. (2018). A scalarized solution exists when $1/(2\sqrt{\eta})=M/\bar{\eta}^{1/2}$ (which is the same as $M/\lambda$ in their notation) is smaller than $0.59$, but for $0.54<M/\bar{\eta}^{1/2}<0.59$ the effective potential is not positive definite.

We now turn to a study of the perturbations of the fully coupled field equations [Eqs. (2) and (3)].

IV Nonlinear black hole solutions and their radial (in)stability

In GR, radial perturbations describe nonradiative fields. The perturbation equations can be solved analytically and correspond to a change in mass of the Schwarzschild BH solution, as expected from Birkhoff’s theorem Zerilli (1970). In modified theories of gravity radial perturbations can be radiative, with important consequences for the stability of the spacetime.

As mentioned in the introduction, Ref. Blázquez-Salcedo et al. (2018) studied radial perturbations of a static, spherically symmetric BH for a generic coupling function $f(\varphi)$ in sGB gravity. We mostly follow their treatment. For brevity, we will only outline the procedure to obtain the perturbation equations and our numerical calculation of the radial oscillation modes.

The spherically symmetric, radially perturbed spacetime up to first order in the perturbations has line element

[TABLE]

where $\varepsilon$ is a small bookkeeping parameter and ${\rm d}\Omega^{2}={\rm d}\theta^{2}+\sin^{2}\theta\,{\rm d}\phi^{2}$ is the line element of the unit-sphere. To the same order, the scalar field is given by

[TABLE]

By inserting Eqs. (17) and (18) into the field equations (2)-(3) and expanding in powers of $\varepsilon$ we get equations for the background metric functions $(\Phi,\,\Lambda)$ at zeroth order in $\varepsilon$, and for the radial perturbations $(F_{t},\,F_{r},\,\varphi_{1})$ at first order in $\varepsilon$. Let us first discuss the background equations and boundary conditions.

IV.1 Black hole scalarization in quartic sGB gravity

The zeroth-order equations for $(\Phi,\,\Lambda,\,\varphi_{0})$ can be cast as a coupled system of two first-order equations for $\Phi$ and $\Lambda$ and a second-order equation for $\varphi_{0}$ Kanti et al. (1996); Doneva and Yazadjiev (2018a); Blázquez-Salcedo et al. (2018). BH solutions are obtained by imposing that the metric functions $\exp(2\Phi)$ and $\exp(-2\Lambda)$ vanish and that the scalar field $\varphi_{0}$ be regular at the horizon:

[TABLE]

where $r_{\rm h}$ is the horizon radius, while $\varphi_{0,{\rm h}}$ and $\varphi_{0,{\rm h}}^{\prime}$ denote the scalar field and its first derivative at the horizon. Using a near-horizon expansion of the field equations, we find that BH solutions correspond to the condition (cf. Antoniou et al. (2018a); Silva et al. (2018); Doneva and Yazadjiev (2018a) or more details)

[TABLE]

When Eq. (LABEL:eq:der_field) is not satisfied, $\varphi^{\prime\prime}_{0}$ diverges at the horizon. By requiring that the first derivative of the scalar field at the horizon be real and using Eq. (LABEL:eq:der_field) we find a condition for the existence of the solutions:

[TABLE]

At large distances, an expansion of the background equations in powers of $r^{-1}$ leads to

[TABLE]

where $M$ is the total (ADM) mass and $Q$ is the scalar charge. Since we are not interested in cosmological effects and we require that the hairy solution has the same asymptotic properties as the GR solution, we will also impose $\varphi_{0,\infty}=0$.

To obtain scalarized BH solutions we proceed as follows. We integrate the differential equations for the background from the horizon outwards imposing the boundary conditions (19)–(21) with the constraint (23), using a guess value for the scalar field at the horizon $\varphi_{0,{\rm h}}$. Numerical solutions in the far-horizon region ($r\gg r_{\rm h}$) are then compared with the boundary conditions (24)–(25). Not all sets of $(\varphi_{0,\rm{h}},\eta,\zeta)$ allow for BH solutions satisfying both boundary conditions. This generates a boundary value problem that can be solved by a shooting method. In practice, we fix the values of $(\eta,\zeta)$, and we find the values of $\varphi_{0,{\rm h}}$ by shooting and requiring that the scalar field vanishes in the far region.

IV.2 Radial perturbations of scalarized black holes

Let us now consider the radial stability of the BH solutions found in the preceding section. By manipulating the first-order equations we can show that the functions $(F_{t},\,F_{r},\,\varphi_{1})$ are not independent: $F_{t}$ and $F_{r}$ can be written in terms of $\varphi_{1}$, where $\varphi_{1}$ obeys the differential equation Torii and Maeda (1998); Blázquez-Salcedo et al. (2018)

[TABLE]

Here $h$, $k$ and $p$ are functions of $r$ that depend on the background metric functions: cf. Eq. (14) of Blázquez-Salcedo et al. (2018).

By a suitable redefinition of the functions $(h,k,p)$ and using a harmonic-time decomposition $\varphi_{1}(t,r)=\varphi_{1}(r)e^{-i\omega t}$, we can write the above equation in a Schrödinger-like form. This is useful for analyzing the effective potential felt by the perturbations Blázquez-Salcedo et al. (2018). However here we will deal directly with the differential equation in the form (27), mainly because it is simpler to solve it numerically. We introduce a compactified dimensionless coordinate

[TABLE]

such that the horizon and spatial infinity are mapped to $x=0$ and $x=1$, respectively. To integrate Eq. (27) numerically we impose the standard boundary conditions at the horizon and at spatial infinity:

[TABLE]

where $r_{*}$ is the tortoise coordinate.

The differential equation (27) together with boundary conditions (29) yields a boundary-value problem for the complex eigenvalue $\omega=\omega_{R}+i\omega_{I}$. Stable modes have $\omega_{I}<0$, while unstable modes have $\omega_{I}>0$. Therefore, to study the radial stability of the solutions we can search for purely imaginary modes with $\omega_{I}>0$. To obtain these modes, we use again a shooting method. We perform two integrations starting at $x=0$ and at $x=1$. At each boundary we impose that the scalar field is zero and that its first derivative is constant. We can fix the scalar field amplitude to unity because Eq. (27) is linear. The integration of Eq. (27) from the horizon yields a first solution $\varphi_{1}^{(-)}$, and the integration from infinity yields a second solution $\varphi_{1}^{(+)}$. We match the two solutions at an intermediate point $x_{m}$. The eigenvalue $\omega=i\omega_{I}$ corresponds to the frequency at which the Wronskian

[TABLE]

vanishes. We checked that the modes are stable under variations of the numerically chosen values of the near-horizon radius, of the large radius representing spatial infinity, and of the matching point $x_{m}$. Additionally, we checked that our results reproduce those for the exponential and quadratic couplings presented in Ref. Blázquez-Salcedo et al. (2018).

IV.3 Numerical results

In Fig. 2 we plot the scalar charge-mass ($Q$–$M$) diagram corresponding to nodeless ($n=0$) scalarized BH solutions in the quartic theory for fixed values of $|\zeta/\eta|=0,\,0.5,\,1,\,1.5$ and $2.0$. For completeness we also show the corresponding diagram for the exponential coupling from Ref. Doneva and Yazadjiev (2018a). The vertical line represents the threshold $\eta=\eta_{\rm thr}$.

Numerically we found that solutions with $\zeta/\eta\lesssim-0.8$ bend to the left in the $Q-M$ diagram (and so do solutions corresponding to the exponential coupling), while solutions with $\zeta/\eta$ larger than this critical value bend to the right. Solutions illustrating this behavior near this critical value are shown in the inset of Fig. 2. This different behavior corresponds to different radial stability properties. When $\zeta/\eta=-0.5$ there is a gap in the parameter space in which BH solutions do not exist: the first derivative of the scalar field at the horizon is complex in this region. This happens because Eq. (23) cannot be satisfied for parameter choices that lie in the gap. We verified that the gap exists in the parameter range $-0.64<\zeta/\eta<0$, and it is indeed related to the polynomial form of the existence condition given by Eq. (LABEL:eq:der_field).

The two different branches of solutions that exist (e.g.) when $\zeta/\eta=-0.5$ present different stability properties, one being stable and the other unstable. We can understand this behavior qualitatively using intuition built from Sec. III. From (16), we know that not only the magnitude of the coupling parameter $\zeta$, but also the amplitude of the background solution $\varphi_{0}$ contributes to quenching the instability. In the unstable segment of the $\zeta/\eta=-0.5$ solutions, the scalar field has small amplitude (i.e., small scalar charge) and therefore it cannot quench the instability. This is not the case for more charged solutions, which are stable.

To analyze the stability, in Fig. 3 we plot the unstable mode frequencies for the same set of theories with $(\eta,\zeta)$ shown in Fig. 2, and also for Schwarzschild BH solutions in the same theory. Note that for Schwarzschild BH solutions $\varphi_{0}=0$, and therefore the frequency does not depend on $\zeta$. Hence, all unstable scalarized BHs branch out of the same threshold value of $\eta\approx 0.726$, in agreement with the value obtained in Sec. III. All BH solutions in quartic theories with $\zeta/\eta>0$ shown in Fig. 3 are unstable to radial perturbations, again in agreement with the analysis from the decoupling limit, and the instability time scale $\tau=|\omega_{I}^{-1}|$ of these modes decreases as $\zeta/\eta$ increases (as long as $\zeta/\eta>0$).

As noted when discussing Fig. 2, there is a gap in the parameter space of BH solutions with $\zeta/\eta=-0.5$. The two branches have different behavior in the $Q$–$M$ plane: the branch with small values of $Q/\bar{\eta}^{1/2}$ is more similar to solutions with $\zeta/\eta>-0.5$, and the branch with large values of $Q/\bar{\eta}^{1/2}$ is more similar to solutions with $\zeta/\eta<-0.5$. We performed a radial stability analysis searching for unstable modes in the two branches. We found no unstable modes in the large $Q/\bar{\eta}^{1/2}$ branch, but we found unstable modes in the small $Q/\bar{\eta}^{1/2}$ branch, and these are the $\zeta/\eta=-0.5$ modes shown in Fig. 3.

The case $\zeta/\eta=-0.5$, presenting stable and unstable solutions, is similar to the ones shown in the inset of Fig. 2. Note that the unstable modes for $\zeta/\eta=-0.5$ have small $M\omega_{I}$ (see inset of Fig. 3). For the solutions with values of $\zeta/\eta$ presented in the inset in Fig. 2 the mode frequency is even smaller, being challenging to find numerically. Our numerical findings suggest that in quartic sGB with $\zeta/\eta<-0.8$ scalarized BHs with $n=0$ are stable.

V Conclusions

In this work we have investigated the radial stability of scalarized BH solutions in sGB gravity. Motivated by the radial instability of quadratic sGB solutions found in Blázquez-Salcedo et al. (2018), we have shown that adding higher-order (quartic) corrections to the original quadratic sGB model of Silva et al. (2018) can stabilize the solutions.

Our analysis provided a clear picture for the physical interpretation of this results. At the linearized level, scalarization manifests as a tachyonic instability that triggers the growth of the scalar field. For the end-point of the instability to be a hairy solution, the tachyonic instability needs to be quenched by some nonlinear effects. In quadratic sGB gravity however, the field equation for the scalar field is linear in the scalar and hence the only quenching mechanism would be backreaction. This is nicely demonstrated by our decoupling limit analysis, where backreaction is entirely ignored and the tachyonic instability is always present. Nonetheless, within the same approximation, a higher-order coupling introduces nonlinearity in the scalar and provides strong quenching for the tachyonic instability. This highlights that the very existence of hairy solutions found in Silva et al. (2018) for the purely quadratic model relies on backreaction effects, and this is what renders them rather special. This seems to be reflected on their radial stability properties.

Indeed, a $\varphi^{4}$-term turns out to stabilize scalarized BH solutions. We computed unstable radial modes and we found none when the coupling parameters satisfy $\zeta/\eta<-0.8$. This suggests that scalarized BH solutions are stable in this region of the theory’s parameter space.

More generally, our results clearly demonstrate that the quadratic coupling between the scalar and the Gauss–Bonnet invariant controls the onset of the scalarization, whereas the higher-order corrections (in the scalar) control the end-point of the tachyonic instability that triggers scalarization, and hence they are crucial for the properties of the hairy black holes solutions.

Note Added – While this work was being completed, a preprint with similar conclusions appeared as an e-print Minamitsuji and Ikeda (2019). Where our works overlap, our conclusions agree with theirs.

Acknowledgments

This work was supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904 and by the COST action CA16104 “GWverse”. H.O.S was supported by NASA Grants No. NNX16AB98G and No. 80NSSC17M0041. T.P.S. acknowledges partial support from the STFC Consolidated Grant No. ST/P000703/1. J.S. was supported by funds provided to the Center for Particle Cosmology by the University of Pennsylvania. E.B. is supported by NSF Grant No. PHY-1841464, NSF Grant No. AST-1841358, NSF-XSEDE Grant No. PHY-090003, and NASA ATP Grant No. 17-ATP17-0225. C.F.B.M. would like to thank the Johns Hopkins University for kind hospitality during the preparation of this work and the American Physical Society which funded the visit through the International Research Travel Award Program.