Stability of black holes with non-minimally coupled scalar hair to the Einstein tensor

Abstract

General relativity admits a plethora of exact compact object solutions. The augmentation of Einstein’s action with non-minimal coupling terms leads to modified theories with rich structure, which, in turn, provide non-trivial solutions with intriguing phenomenology. Thus, assessing their viability under generic fluctuations is of utmost importance for gravity theories. We consider static and spherically-symmetric solutions of a Horndeski subclass which includes a massless scalar field non-minimally coupled to the Einstein tensor. Such theory possesses second-order field equations and admits an exact black hole solution with scalar hair. Here, we study the stability of such solution under axial gravitational perturbations and find that it is linearly stable. The qualitative features of the ringdown waveform depend solely on the ratio of the two available parameters of spacetime, namely the black hole mass m and the non-minimal coupling strength $\ell _\eta$. Finally, we demonstrate the gravitational-wave ringdown transitions between three distinct patterns as the ratio $m/\ell _\eta$ increases; a state which is dominated by photon-sphere excitations and maintains a typical quasinormal ringdown, an intermediate long-lived state which exhibits gravitational-wave echoes and, finally, a state where the ringdown and echoes are depleted rapidly to turn to an exponential tail.

1 Introduction

Compact objects play a decisive role in contemporary astrophysics, as their relativistic collisions may provide crucial information concerning astrophysical processes in extreme-gravity conditions. The latest gravitational-wave (GW) detections by ground-based interferometers [1,2,3,4,5] have provided intelligence regarding the strong-field regime. The early stage of the gravitational ringdown of black hole (BH) mergers, described by quasinormal modes (QNMs) [6,7,8,9], further contributes to the understanding of their relaxation properties, as well as the governing theory of gravity. Nevertheless, a conclusive interpretation of the underlying gravitational theory has not yet been met. Therefore, it is expected that future ground and space-borne detectors will improve our perception of gravitational interactions, and in particular will shed light into the existence of exotic compact objects (ECOs) [10,11,12,13,14,15,16,17,18,19] which may possess unexpected multipolar and near-horizon structures that differ significantly from those of BHs.

ECOs are spacetime solutions of general relativity (GR) and modified gravity theories that describe compact objects with exotic properties and intriguing multipolar structure [20], such as BHs which evade the ‘no-hair’ theorem and give rise to additional spacetime parameters (besides the mass, spin and charge), wormholes that evade singularities and connect Universes [11] as well as horizonless compact objects which possess unexpected near-horizon structures that expel the event horizon and subsequent singularities through reflective centrifugal barriers [21]. The majority of ECOs, which possess a photon sphere, can naturally mimic BHs when perturbed and produce prompt ringdown waveforms in the time domain which are identical to those of BHs [22, 23]. This occurs due to the indifference of photon sphere excitations from external perturbations. The dominant effects of ECO ringdown only appear at late times in the form of successively damped repetitions of subsequent photon sphere excitation, know as echoes, which occur due to the entrapment of perturbations inside potential wells and the formation of quasibound states [24,25,26,27]. These modes represent the actual QNM content of the ECO, which in the frequency domain is dramatically different from the QNMs of BHs [28]. In what follows, we will loosely refer to the spectral content of the prompt ringdown as the QNM spectrum whenever the echo timescales are sufficiently large, even though these to do not necessarily correspond to the actual QNMs of the full eigenvalue problem.

Even though GR has withstood many experimental tests, remaining consistent with plenty of observations so far, modified theories of gravity attempt to describe phenomena where GR seems to fail, such as the construction of viable cosmological models for inflation and dark energy [29]. The most general scalar-tensor theory of gravity in four dimensions whose Lagrangian is constructed by the metric tensor and a non-minimally coupled scalar field is Horndeski’s theory [30]. This gravity theory contains subclasses that preserve a classical Galilean symmetry [31, 32], leads to second-order field equations and is free of ghost instabilities [33,34,35,36]. The most extensively studied subclass of Horndeski theory is represented by a Lagrangian with a scalar field non-minimally coupled to the Einstein tensor.

On large scales, the non-minimal coupling term of the aforementioned subclass possesses very intriguing effects on inflationary dynamics. From an inflationary model-building point of view, it allows for a very effective implementation of a slow-roll phase, due to the fact that it acts as a friction mechanism [37, 38], allowing potentials such as the Standard-Model Higgs [39], thus making it a very attractive term in Horndeski theory. A generalization of the non-minimal kinetic term and its application to inflation was recently analyzed in [40,41,42]. During the inflationary phase of the Universe in GR, scalar and tensor perturbations result in spectra which are red-shifted [43]. One then expects that the presence of a non-minimal kinetic term will magnify the red-shift behavior of the perturbation spectra because of the friction effect, which is subsequently related to the decrease of the Hubble parameter H during inflation (for a review of the effects of the non-minimal kinetic term in inflation see [44]). On the contrary, if the scalar fields are phantom fields with negative kinetic energy non-minimally coupled to the Einstein tensor, then the spectra of scalar and tensor perturbations produced during an inflationary phase are typically blue-shifted [45]. Their dynamics in cosmological setups, as well as the instabilities at which tachyons or ghosts appear in the infrared region around the present Hubble scale were discussed in [46]. Beyond inflation, the non-minimal kinetic coupling has also been utilized to construct cosmological models [47].

Apart from cosmological applications, the particular subclass of Horndeski theory allows the construction of BH solutions with scalar hair [48,49,50,51,52]. Consequently, an important aspect of these compact objects is their stability. Regarding the formation of stable hairy BHs, the ‘no-hair’ theorem should be evaded [53, 54], which translates to the existence of a balance mechanism to outweigh the gravitational force outside the BH event horizon. A typical example is offered by holography. A charged scalar field theory embedded into an anti-de Sitter (AdS) Lagrangian leads to the formation of horizon hair, as a result of the counterbalance between the attractive gravitational and the repulsive electromagnetic force [55, 56]. Then, according to the gauge/gravity duality, such mechanism allows a holographic phase transition which results to a conformal field theory describing a holographic superconductor on the AdS boundary [57,58,59], besides other interesting phenomena [60,61,62,63,64].

In the subclass of the Horndeski theory in which a scalar field is coupled kinetically to the Einstein tensor there is a direct coupling of matter to curvature and there exist local solutions in which this coupling appears as a primary charge in the metric functions of the resulting hairy BH, that may play the role of an effective negative cosmological constant, even though the action is absent of any cosmological constant term. An interesting aspect of such objects are their thermodynamical properties as well as their viability as novel compact objects. One of the most important requirements for the viability of these objects is their stability against perturbations. To that end, the stability of the hairy BH solution found in [49], against linear scalar perturbations, was recently assessed [65]. At the linearized level, the non-minimal coupling constant sources an effective asymptotic boundary where the effective potential of the wave equation that governs the propagation of scalar perturbations diverges. Such a boundary serves as a perfect reflector for incident scalar waves and generates a trapping region outside the photon sphere without the need of invoking a negative cosmological constant in the action of the theory. As a result, the ringdown signal of the BH exhibits successively damped echoes. Thus, scalarized compact objects in the particular Horndeski class can serve as alternatives to the standard echo sources which possess trapping regions beyond the photon sphere due to near-horizon structures [15, 21,22,23, 66,67,68,69,70,71].

Gravitational perturbations in modified theories of gravity provide information regarding the velocity with which GWs travel. The recent observations of GW170817 and GRB170817A, as well as its electromagnetic counterpart, imply that GWs travel at the speed of light, with deviations smaller than a few $10^{-15}$. The consequences of this experimental result for models of dark energy and modified gravity theories were discussed in [72, 73]. In particular these constraints on the speed of GWs were used to test some classes of Horndeski theory. A detailed discussion of the effects of the kinetic coupling on the speed of the GWs in the subclass of the Horndeski theory, in which the scalar field is coupled to the Einstein tensor, is provided in [74]. It was found that while the kinetic energy of a minimally coupled scalar field does not change under the cosmological evolution, the kinetic energy of the scalar field coupled to the Einstein tensor changes as the Universe expands. At the inflationary epoch it acts as a friction term and drives inflation with steep potentials, while as the Universe expands its contribution to the cosmological evolution is less important and at the late cosmological epoch is negligible, thus GWs propagate at the speed of light at late cosmological times.

In any case, Horndeski theories do predict a modified speed of GW propagation. Even so, recent studies [75,76,77] demonstrate that with analogue versions of Horndenki gravity, which are based on teleparallel gravity constructed with a nonvanishing torsion tensor, one can device a more general Horndeski theory where GWs propagate with the speed of light without eliminating the coupling functions $G_4(\Phi ,X)$ and $G_5(\Phi ,X)$ that were highly constrained in standard Horndeski theory. Hence, in the teleparallel approach one is able to restore these terms, creating an interesting way to revive this theory of gravity. Even though our analysis still lies in a curvature-based formulation of gravity, it is still very interesting that there are ways of evading the tight constraints of Horndeski theory.

The purpose of this work is twofold. First, we investigate the effect of the kinetic coupling of the scalar field to the Einstein tensor on the stability of local solutions of the particular subclass of Horndeski theories. We will work with the hairy BH solution [49] for which scalar perturbations have been analyzed recently [65, 78] for a wide range of the kinetic coupling. Under these analyses the hairy BH was found to be linearly stable with echoes being present at late times on the ringdown waveform. Here, we perform a first step towards gravitational modal stability, thus extending our previous test scalar field analysis to axial gravitational perturbations. Since the kinetic coupling appears as a primary charge in the metric functions, we expect to get a better understanding on the stability of such objects, although a complete picture of gravitational modal stability can only be discerned when one considers not only the axial but also the polar sector of fluctuations which generally couple to the scalar hair present in scalar-tensor theories. Our second goal is to investigate how the kinetic coupling affects the ringdown waveform and attempt to assess the appearance of GW echoes in the parametric space of such geometries.

The work is organized as follows. In Sect. 2 we review the BH solution of the Horndeski theory with a scalar field kinetically coupled to the Einstein tensor. In Sect. 3 we discuss the general framework of axial gravitational perturbations and we derive the effective potential of the considered BH solution. In Sect. 4 we demonstrate the numerical scheme of time-domain integration. In Sect. 5 we study the evolution of the axial gravitational perturbations and finally in Sect. 6 we conclude this work.

2 Black hole solution with a scalar field kinetically coupled to Einstein tensor

In what follows, we will consider static solutions of a scalar-tensor theory in which the scalar field is kinetically coupled to the Einstein tensor. This is part of the most general scalar-tensor theory which yields second-order field equations, namely the Horndeski theory. The full Lagrangian is given by:

$$\begin{aligned} \mathcal {L}&=\sum _{i=2}^{i=5}\mathcal {L}_i~,\nonumber \\ \mathcal {L}_2&=K(\Phi ,X)~,\nonumber \\ \mathcal {L}_3&=-G_3(\Phi ,X)\square \Phi ~,\nonumber \\ \mathcal {L}_4&=G_4(\Phi ,X)R+G_{4,X}\left[ (\square \Phi )^2-(\nabla _\mu \nabla _\nu \Phi )^2\right] ~.\nonumber \\ \mathcal {L}_5&=G_5(\Phi ,X)G_{\mu \nu }\nabla ^\mu \nabla ^\nu \Phi -\frac{1}{6}G_{5,X}\left[ (\square \Phi )^3-3\square \Phi (\nabla _\mu \nabla _\nu \Phi )^2 +2(\nabla _\mu \nabla _\nu \Phi )^3\right] ~, \end{aligned}$$

(2.1)

where $X=-\frac{1}{2}\nabla _\mu \Phi \nabla ^\mu \Phi$. We consider a particular subset of Horndeski theory with non-trivial $\mathcal {L}_2=K(\Phi ,X)=2\varepsilon X$ and $G_4(\Phi ,X)=(8\pi )^{-1}-\eta X$ terms. The theory is described by the following action,

$$\begin{aligned} S=\int d^4x \sqrt{-g}\left[ \frac{R}{8\pi }-(\varepsilon g_{\mu \nu }+\eta G_{\mu \nu })\partial ^{\mu }\Phi \partial ^{\nu }\Phi \right] , \end{aligned}$$

(2.2)

where $g_{\mu \nu }$ is the metric tensor, $g=\mathrm{det}(g_{\mu \nu })$, R is the scalar curvature, $G_{\mu \nu }$ is the Einstein tensor, $\Phi$ is a real massless scalar field and $\eta$ is the non-minimal kinetic coupling parameter with dimensions of length-squared. The $\varepsilon$ parameter equals $\pm 1$, where in the case $\varepsilon =+1$ we have a canonical scalar field with positive kinetic term, while the case $\varepsilon =-1$ corresponds to a phantom scalar field with negative kinetic energy. Even though in the original Hordenski theory the kinetic energy of the scalar field is positive, in this work we will also considered the case where the scalar field’s kinetic energy is negative.

Varying the action (2.2) with respect to the metric tensor $g_{\mu \nu }$ and scalar field $\Phi$ provides the following field equations

$$\begin{aligned}&G_{\mu \nu }=8\pi \big [\varepsilon T_{\mu \nu } +\eta \Theta _{\mu \nu }\big ]~, \end{aligned}$$

(2.3a)

$$\begin{aligned}&[\varepsilon g^{\mu \nu }+\eta G^{\mu \nu }]\nabla _{\mu }\nabla _\nu \Phi =0~, \end{aligned}$$

(2.3b)

where

$$\begin{aligned} T_{\mu \nu }&=\nabla _\mu \phi \nabla _\nu \Phi - \frac{1}{2}g_{\mu \nu }(\nabla \Phi )^2~, \end{aligned}$$

(2.4)

$$\begin{aligned} \Theta _{\mu \nu }&=-\frac{1}{2}R\nabla _\mu \Phi \,\nabla _\nu \Phi \, +2\nabla _\alpha \Phi \,\nabla _{(\mu }\Phi R^\alpha _{\nu )} +\nabla ^\alpha \Phi \,\nabla ^\beta \Phi \,R_{\mu \alpha \nu \beta } +\nabla _\mu \nabla ^\alpha \Phi \,\nabla _\nu \nabla _\alpha \Phi \nonumber \\&\quad -\,\nabla _\mu \nabla _\nu \Phi \,\square \Phi -\frac{1}{2}(\nabla \Phi )^2 G_{\mu \nu }\nonumber \\&\quad +\,g_{\mu \nu }\left[ -\frac{1}{2}\nabla ^\alpha \nabla ^\beta \Phi \, \nabla _\alpha \nabla _\beta \Phi +\frac{1}{2}(\square \Phi )^2 -\nabla _\alpha \Phi \,\nabla _\beta \Phi \,R^{\alpha \beta } \right] ~. \end{aligned}$$

(2.5)

A static and spherically-symmetric BH solution to the aforementioned theory has been found in [49], where the scalar field of the theory depends only on the radial coordinate. The solution yields the constraint $\varepsilon \eta <0$, which leads to the definition of the following coupling parameter

$$\begin{aligned} \ell _\eta =\left| \varepsilon \eta \right| ^{1/2}~. \end{aligned}$$

(2.6)

In terms of the line element

$$\begin{aligned} ds^2=-g_{tt}(r)dt^2+g_{rr}(r)dr^2+g_{\theta \theta }(r)d\Omega ^2~, \end{aligned}$$

(2.7)

the BH solution corresponds to $g_{\theta \theta }(r)=r^2$ with $r\in (0,+\infty )$ and yields the following metric components

$$\begin{aligned} g_{tt}(r)&=-\frac{1}{4}F(r)~, \end{aligned}$$

(2.8a)

$$\begin{aligned} g_{rr}(r)&=\frac{(r^2+2\ell _\eta ^2)^2}{(r^2+\ell _\eta ^2)^2F(r)}~, \end{aligned}$$

(2.8b)

$$\begin{aligned} F(r)&=\left[ 3+\frac{r^2}{3 \ell ^2_\eta }-\frac{8 m}{r}+\frac{\ell _\eta }{r}\arctan \left( \frac{r}{\ell _\eta }\right) \right] , \end{aligned}$$

(2.8c)

while the scalar hair of the theory reads

$$\begin{aligned} (\partial _r\Phi )^2=\Psi ^2=-\frac{\varepsilon }{8\pi \ell _\eta ^2}\frac{r^2 g_{rr}}{r^2+\ell _\eta ^2}~, \end{aligned}$$

(2.9)

which implies that

$$\begin{aligned} 4\pi \eta \Psi ^2=\frac{1}{2}\frac{r^2 }{r^2+\ell _\eta ^2}g_{rr}> 0~, \quad \forall r> r_h \end{aligned}$$

(2.10)

where $r=r_h$ is the event horizon radius. It is important to note that the asymptotic behavior of the lapse function (2.8a) when $r\rightarrow +\infty$, becomes $F(r)\sim r^2/\ell _\eta ^2$, where the term $\ell _\eta ^2$ assumes a form of an effective cosmological scale, similar to that of an actual cosmological radius, with dimensionality length squared in geometrized units. Even so, we have to stress that the action does not contain any negative cosmological constant term and the emergence of this effective scale is solely due to the non-minimal coupling of the scalar field to the Einstein tensor.

Another important note is that the equations of motion of the scalar field, (2.3b), can be expressed as the conservation of the Noether current that corresponds to the shift symmetry of the Galileon, i.e. $\Phi \rightarrow \Phi +\delta \Phi$, where $\delta \Phi$ is constant. It can straightforwardly be found that the current is defined as

$$\begin{aligned} J^\mu =(\varepsilon g^{\mu \nu }+\eta G^{\mu \nu })\nabla _\nu \Phi ~. \end{aligned}$$

(2.11)

The BH solution satisfies the physical requirement that the norm of this current does not diverge at the horizon, by virtue of (2.6). The scalar hair, however, diverges at the horizon as one can readily see from (2.9). One can also deduce from (2.9) that the metric components can be expressed in terms of the scalar hair. As such, the scalar hair of the theory can be understood as an intrinsic part of the geometry. Finally, we note that due to the Bianchi identity, $\nabla ^\mu G_{\mu \nu }=0$, the equation (2.3a) leads to a differential consequence

$$\begin{aligned} \nabla ^\mu \big [\varepsilon T_{\mu \nu }+\eta \Theta _{\mu \nu }\big ]=0~. \end{aligned}$$

(2.12)

The substitution of expressions (2.4) and (2.5) into the Bianchi identity yields Eq. (2.3b). In other words, the equation of motion of scalar fields (2.3b) is the differential consequence of the system (2.3a) due to the general covariance and the absence of further degrees of freedom. Let us also note that the solution reproduces the Schwarzschild BH in the limit of $\ell _\eta \rightarrow +\infty$, therefore the metric can be understood as a hairy BH generalization of the Schwarzschild spacetime with effective AdS-asymptotics, when the spin-0 degree of freedom also acquires dynamics from the kinetic mixing with the graviton, i.e. the $G^{\mu \nu }\partial _\mu \Phi \partial _\nu \Phi$ term.

As one can observe, the metric functions (2.8) of the BH solution [49] depend only on the parameter $\ell _\eta$, besides the mass parameter m. Therefore, they do not contain the information of whether the produced compact object is made of normal or phantom matter as long as $\epsilon \eta <0$. Even though in [49] the scalar field is assumed to be canonical ($\epsilon >0$) with the non-minimal coupling being negative ($\eta <0$), we have checked that the same solution is obtained when the scalar field is phantom ($\epsilon <0$) and the non-minimal coupling is positive ($\eta >0$).

In finality, let us mention that in [38] it was argued that the subclass of the Horndeski theory under consideration yields wormhole solutions as well. We wish to note here that the wormhole solution derived there is just a coordinate artifact of the BH and does not correspond to a traversable wormhole. Let us consider the BH metric we previously described

$$\begin{aligned} ds^2=-\frac{1}{4}F(r)dt^2+\frac{(r^2+2\ell _\eta ^2)^2}{(r^2+\ell _\eta ^2)^2F(r)}dr^2+r^2d\Omega ^2. \end{aligned}$$

(2.13)

We perform the following coordinate transformation

$$\begin{aligned} u^2=r^2-a^2, \quad u\in (-\infty ,+\infty ) \end{aligned}$$

(2.14)

and the coordinate redefinition

$$\begin{aligned} dt^2=\frac{4}{\left( 3-\frac{8m}{a}+\frac{a^2}{3\ell _\eta ^2}+\frac{\ell _\eta }{a}\arctan \left( \frac{a}{\ell _\eta }\right) \right) }dT^2=CdT^2\,. \end{aligned}$$

(2.15)

Note that this coordinate transformation covers the BH region for $r>a$ twice. Fixing $a>r_h$, the corresponding geometry will cover only the region $r>r_h$ twice. Performing the coordinate transformation (2.14), one finds:

$$\begin{aligned} ds^2&=-CF(\sqrt{u^2+a^2})dT^2\nonumber \\&\quad +\,\frac{u^2(u^2+a^2+2\ell _\eta ^2)^2}{(u^2+a^2)(u^2+a^2+\ell _\eta ^2)^2F(\sqrt{u^2+a^2)}}du^2\nonumber \\&\quad +\,(u^2+a^2)d\Omega ^2, \end{aligned}$$

(2.16)

$$\begin{aligned} F(\sqrt{u^2+a^2)}&=\left( 3-\frac{8m}{\sqrt{u^2+a^2}}+\frac{u^2+a^2}{3\ell _\eta ^2}+\frac{\ell _\eta }{\sqrt{u^2+a^2}}\arctan \left( \frac{\sqrt{u^2+a^2}}{\ell _\eta }\right) \right) . \end{aligned}$$

(2.17)

The corresponding metric is the solution derived in [79], modulo the form of the $g_{tt}$ component, which was left as an indefinite integral. In particular, the integral factor in the $g_{tt}$ component found in [79] can be solved exactly to yield:

$$\begin{aligned} g_{tt}&= \frac{a}{\sqrt{u^2+a^2}}\exp \left[ \int _0^u\frac{u(u^2+a^2+2\ell _\eta ^2)^2}{\ell _\eta ^2(u^2+a^2)(u^2+a^2+\ell _\eta ^2)F(\sqrt{u^2+a^2})}dr\right] \nonumber \\&=\frac{a}{\sqrt{u^2+a^2}}\exp \left[ \int _0^u\left( \frac{d(\ln [\sqrt{u^2+a^2}F(\sqrt{u^2+a^2})])}{du}\right) du\right] \nonumber \\&=\frac{F(\sqrt{u^2+a^2})}{F(a)}\overset{(2.15)}{=} \frac{C}{4}F(\sqrt{u^2+a^2}) \end{aligned}$$

(2.18)

Obviously, a compact object cannot change nature due to a coordinate transformation, thus, the metric (2.16) is just the BH solution (2.8) written in a “bad” coordinate system and was falsely identified as a wormhole. This result is in accordance with the findings in [80] regarding the absence of static and spherically-symmetric wormhole solutions in the particular subclass of Horndeski theory.

3 Axial gravitational perturbations: general analysis

In this section we will undergo an analysis of axial perturbations of the BH solution using Chandrasekhar’s method [81]. The most general metric for an axisymmetric non-stationary spacetime is given by

$$\begin{aligned} ds^2=-e^{2f_0}dt^2+e^{2f_3}(d\phi -q_0dt-q_1 dr-q_2 d\theta )^2+e^{2f_1}dr^2+e^{2f_2}d\theta ^2~. \end{aligned}$$

(3.1)

This result is found by use of the Cotton–Darboux theorem, which states that any three-dimensional metric, $g=g^{ij}\partial _i\partial _j$, can always be brought to a diagonal form by a local coordinate transformation. It is clear that the background metric of our solutions can be described by $q_i=0$. In this gauge, axial perturbations are described by the non-vanishing values of $q_i$, while polar perturbations are described by $f_i\rightarrow f_i+\delta f_i$ and $q_i=0$.

For the purposes of writing down the explicit form of the equations (2.3a) for the most general form of the metric of (3.1), we shall obtain the components of the curvature tensors via Cartan’s structure equations. We choose the following tetrads to work with

$$\begin{aligned}&\epsilon ^{0}_{\mu }=(e^{f_0},0,0,0)~, \end{aligned}$$

(3.2a)

$$\begin{aligned}&\epsilon ^{1}_{\mu }=(0,e^{f_1},0,0)~, \end{aligned}$$

(3.2b)

$$\begin{aligned}&\epsilon ^{2}_{\mu }=(0,0,e^{f_2},0)~, \end{aligned}$$

(3.2c)

$$\begin{aligned}&\epsilon ^{3}_{\mu }=(-e^{f_3}q_0,-e^{f_3}q_1,-e^{f_3}q_2,e^{f_3})~. \end{aligned}$$

(3.2d)

Under these tetrads, the basis is found to be

$$\begin{aligned} \epsilon ^{0}&=e^{f_0}dt\rightarrow dt=e^{-f_0}\epsilon ^0~, \end{aligned}$$

(3.3a)

$$\begin{aligned} \epsilon ^{1}&=e^{f_1}dt\rightarrow dr=e^{-f_1}\epsilon ^1~, \end{aligned}$$

(3.3b)

$$\begin{aligned} \epsilon ^{2}&=e^{f_2}d\theta \rightarrow d\theta =e^{-f_2}\epsilon ^2~, \end{aligned}$$

(3.3c)

$$\begin{aligned} \epsilon ^{3}&=-e^{f_3}q_0dt-e^{f_3}q_1dr-e^{f_3}q_2d\theta \nonumber \\&\quad +\,e^{f_3}d\phi \rightarrow d\phi =e^{-f_3}\epsilon ^3+q_0e^{-f_0}\epsilon ^0+q_1e^{-f_1}\epsilon ^1+q_2e^{-f_2}\epsilon ^2~. \end{aligned}$$

(3.3d)

The reasoning behind these tetrads is that they associate the perturbations to a single tetrad and thus allow for a decoupling of the equations of motion at first order. The spin connections can be derived from the tetrad postulate, where we associate the zero torsion condition with the Levi–Civita connection as

$$\begin{aligned} \omega ^a_{\mu b}=\epsilon ^a_{\nu }\epsilon ^{\lambda }_b\Gamma ^{\nu }_{\mu \lambda }-\epsilon ^{\lambda }_b\partial _{\mu }\epsilon ^a_{\lambda }~. \end{aligned}$$

(3.4)

Using Cartan’s second structure equation, we can derive the Riemann tensor,

$$\begin{aligned} R^{a}_{b\mu \nu }=2(\partial _{[\mu }\omega ^a_{\nu ] b}+\omega ^{a}_{c[\mu }\omega ^c_{\nu ] b})~, \end{aligned}$$

(3.5)

and from (3.5) all the necessary tensors for the equations of motion of the underlying gravity theory. From (2.3a), we know that

$$\begin{aligned} G_{ab}=\hat{T}_{ab}~, \end{aligned}$$

(3.6)

where $\hat{T}_{ab}=8\pi \big [\varepsilon T_{ab}+\eta \Theta _{ab}\big ]$ (note that indices a, b are Lorentz and not spacetime indices). However, the Einstein and stress-energy tensors acquire contributions from the perturbations. From the linearization of the equations of motion we find that only the $G_{03},G_{13},G_{23}$ terms are important at first order. In fact, equation $\delta G_{03}=\delta \hat{T}_{03}$ is degenerate, i.e. it is automatically satisfied by the other two equations. In particular, we find the following results

$$\begin{aligned} \delta G_{13}&=\frac{1}{2r^3\sin ^2\theta }\frac{1}{\sqrt{g_{tt}}}\frac{\partial }{\partial \theta }Q+\frac{r\sin \theta }{2g_{tt}\sqrt{grr}}\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial r}q_a -\frac{\partial }{\partial t}q_b\right] ~, \end{aligned}$$

(3.7a)

$$\begin{aligned} \delta \hat{T}_{13}&=\frac{4\pi \eta \Psi ^2}{g_{rr}}\delta G_{13}~, \end{aligned}$$

(3.7b)

$$\begin{aligned} \delta G_{23}&=-\frac{1}{2r^2\sin ^2\theta }\frac{1}{\sqrt{g_{tt}g_{rr}}}\frac{\partial }{\partial r}Q+\frac{\sin \theta }{2g_{tt}}\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial \theta }q_a -\frac{\partial }{\partial t}q_c\right] ~, \end{aligned}$$

(3.7c)

$$\begin{aligned} \delta \hat{T}_{23}&=\frac{4\pi \eta \Psi ^2}{g_{rr}}\left( \frac{1}{2r^2\sin ^2\theta }\frac{-1}{\sqrt{g_{tt}g_{rr}}}\frac{\partial }{\partial r}Q -\frac{\sin \theta }{2g_{tt}}\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial \theta }q_a-\frac{\partial }{\partial t}q_c\right] \right) \nonumber \\&\quad -\,4\pi \eta \frac{\partial }{\partial r}\left( \frac{\Psi ^2}{2g_{rr}}\right) \frac{1}{r^2\sin ^2\theta \sqrt{g_{tt}}\sqrt{g_{rr}}}Q~, \end{aligned}$$

(3.7d)

where ${Q=\left[ r^2\sin ^3\theta \frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\frac{\partial }{\partial \theta }q_b-r^2\sin ^3\theta \frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\frac{\partial }{\partial r}q_c\right] }$. Using the redefinition

$$\begin{aligned} \mathcal {F}=\left( 1-\frac{4\pi \eta \Psi ^2}{grr}\right) Q~, \end{aligned}$$

(3.8)

the differential equations we get from $\delta G_{ab}=\delta \hat{T}_{ab}$, using (3.7a)–(3.7d), read

$$\begin{aligned}&\frac{1}{r^4\sin ^3\theta }\sqrt{g_{tt}}\sqrt{g_{rr}}\left( \frac{g_{rr}}{grr-4\pi \eta \Psi ^2}\right) \frac{\partial }{\partial \theta }\mathcal {F}=\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial t}q_b -\frac{\partial }{\partial r}q_a\right] ~, \end{aligned}$$

(3.9a)

$$\begin{aligned}&\frac{1}{r^2\sin ^3\theta }\frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\left( \frac{g_{rr}}{grr+4\pi \eta \Psi ^2}\right) \frac{\partial }{\partial r}\mathcal {F}=\frac{\partial }{\partial t}\left[ \frac{\partial }{\partial \theta }q_a -\frac{\partial }{\partial t}q_c\right] ~. \end{aligned}$$

(3.9b)

Differentiating (3.9a) and (3.9b) by $\theta$ and r respectively, and adding them together yields the differential equation that governs axial perturbations

$$\begin{aligned}&\frac{r^4}{\sqrt{g_{tt}}\sqrt{g_{rr}}}\left( \frac{g_{rr} -4\pi \eta \Psi ^2}{g_{rr}}\right) \frac{\partial }{\partial r}\left[ \frac{1}{r^2}\frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\left( \frac{g_{rr}}{grr +4\pi \eta \Psi ^2}\right) \nonumber \right. \\&\quad \left. \frac{\partial }{\partial r}\mathcal {F}\right] -\frac{r^2}{g_{tt}}\frac{\partial ^2\mathcal {F}}{\partial t^2} +\sin ^3\theta \frac{\partial }{\partial \theta }\left[ \frac{1}{\sin ^3\theta }\frac{\partial }{\partial \theta }\mathcal {F}\right] =0~. \end{aligned}$$

(3.10)

Using a separation of variables, the angular component of equation (3.10) can be understood as the known ultraspherical differential equation with solutions the Gegenbauer polynomials, i.e. the angular component of the perturbation is the same as in the Schwarzschild case. This is to be expected, since both spacetimes are spherically symmetric. Thus, setting $\mathcal {F}=\mathcal {R}(r,t)\mathcal {S}(\theta )$, where in the Appendix A we explicitly show that $\mathcal {S}(\theta )=C^{-3/2}_{\ell +2}(\theta )$, (3.10) can be rewritten as

$$\begin{aligned}&\frac{r^4}{\sqrt{g_{tt}}\sqrt{g_{rr}}}\left( \frac{g_{rr} -4\pi \eta \Psi ^2}{g_{rr}}\right) \frac{\partial }{\partial r}\left[ \frac{1}{r^2}\frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\left( \frac{g_{rr}}{g_{rr} +4\pi \eta \Psi ^2}\right) \nonumber \right. \\&\quad \left. \frac{\partial }{ \partial r}\mathcal {R}\right] -\frac{r^2}{g_{tt}}\frac{\partial ^2\mathcal {R}}{\partial t^2}-(\ell +2)(\ell -1)\mathcal {R}=0~. \end{aligned}$$

(3.11)

where $\ell$ is the angular quantum number of the perturbation. In order to continue, it will prove useful to set

$$\begin{aligned} h=\frac{\sqrt{g_{tt}}}{\sqrt{g_{rr}}}\, , \qquad A=\left( \frac{g_{rr}}{g_{rr}+4\pi \eta \Psi ^2}\right) \, ,\qquad B=\left( \frac{g_{rr}-4\pi \eta \Psi ^2}{g_{rr}}\right) ~,\qquad \end{aligned}$$

(3.12)

thus, simplifying equation (3.11) to

$$\begin{aligned} B r^2 h\frac{\partial }{\partial r}\left[ \frac{1}{r^2}Ah\frac{\partial }{\partial r}\mathcal {R}\right] -\frac{\partial ^2\mathcal {R}}{\partial t^2}-\frac{(\ell +2)(\ell -1)}{r^2}g_{tt}\mathcal {R}=0~. \end{aligned}$$

(3.13)

By further defining a scalar field

$$\begin{aligned} u=\frac{\mathcal {R}}{r}\sqrt{A}~, \end{aligned}$$

(3.14)

and using the tortoise coordinate transformation $h=dr/dr_*$, we find that the equation that governs the axial perturbations reads

$$\begin{aligned} h\frac{\partial }{\partial r}\left[ h\frac{\partial }{\partial r}u\right] -\frac{1}{AB}\frac{\partial ^2u}{\partial t^2}-\left[ \frac{(\ell +2)(\ell -1)}{r^2}\frac{g_{tt}}{AB} +\frac{r}{\sqrt{A}}h\frac{\partial }{\partial r}\left( h\frac{\partial }{\partial r}\left( \frac{\sqrt{A}}{r}\right) \right) \right] u=0~. \end{aligned}$$

(3.15)

If we consider a time dependence of the form $u(r,t)\sim u(r)\exp (i \omega t)$, then (3.15) yields a non-trivial wave equation of the following form

$$\begin{aligned} h\frac{d}{dr}\left[ h\frac{d}{dr}u\right] +\frac{\omega ^2}{AB}u -\left[ \frac{(\ell +2)(\ell -1)}{r^2}\frac{g_{tt}}{AB}+\frac{r}{\sqrt{A}}h\frac{d}{dr}\left( h\frac{d}{dr}\left( \frac{\sqrt{A}}{r}\right) \right) \right] u=0~. \end{aligned}$$

(3.16)

As such, the corresponding Regge–Wheeler-like potential reads

$$\begin{aligned} V(r)=\frac{(\ell +2)(\ell -1)}{r^2}\frac{g_{tt}}{AB}+\frac{r}{\sqrt{A}}h\frac{d}{dr}\left( h\frac{d}{dr}\left( \frac{\sqrt{A}}{r}\right) \right) ~. \end{aligned}$$

(3.17)

It is important to note that by fixing the metric components to those of the BH solution, Eq. (3.17) asymptotes to the standard Schwarzschild Regge–Wheeler effective potential in the limit where $\ell _\eta \rightarrow +\infty$. Another crucial aspect of Eq. (3.16) is the presence of a multiplication factor 1/AB to the gravitational perturbation frequency $\omega$, which signifies the existence of a modified speed of GW propagation [74].

Equation (3.16) demonstrates that one can reduce the problem of axial gravitational perturbations around the compact object under consideration into a single one-dimensional scattering problem with an effective potential. By applying the procedure outlined above, on the BH solution (2.8)–(2.9) we find the corresponding effective potentials a gravitational perturbation induces. An illustration of these potentials for various choices of the angular index of the perturbation $\ell$ is given in Fig. 1. The effective potential possesses a peak for sufficiently large $\ell _\eta$, at the Schwarzschild limit, which is located arbitrary close to the photon sphere at $r=3m$. For a myriad of BH solutions, Ref. [82] demonstrates that the angular frequency and instability timescale of null geodesics that are trapped in unstable circular orbits at the photon sphere are associated with the oscillation frequency and decay rate of eikonal QNMs. In turn, the existence of such centrifugal potential barrier is responsible for the prompt ringdown and photon sphere QNMs found in the response of a plethora of perturbed BHs [83,84,85,86,87,88,89]. In what follows, we will show that the aforementioned analogy only holds at the GR limit and away from it such duality is broken.

Asymptotically, the potential approaches a constant positive value, a behavior very different from the case of scalar perturbations [65], but still, one which encodes the effective non-asymptotically-flat nature of spacetime. A similar behavior was also observed in [90]. More specifically, at spatial infinity, the BH potential at zeroth order yields $V(r\rightarrow \infty ) \sim \mathcal {C}+ \mathcal {O}\left( 1/r\right)$ where the constant $\mathcal {C}$ depends on the parameters $\ell$ and $\ell _\eta$ and can be calculated only numerically. A similar asymptotic behavior was also observed in [91] for the case of vector perturbations in the scalarized BH considered here.

We must note here the dimensionality of various quantities appearing in the following figures, in order to avoid repetition and cluttering on our discussion. According to the geometrized units utilized here, the BH mass m and coupling $\ell _\eta$ have dimensions of length, while the perturbation u is dimensionless, as well as the ratio $m/\ell _\eta$, which makes it an appropriate scale for our analysis. In turn, the effective potential V(r) has dimensionality length to the power of $-2$, as expected, while the frequency $\omega$ has length dimensions.

Fig. 1

Left: BH effective potential with $\ell _\eta = 100\,,\,m = 0.1$ and different values of angular numbers $\ell$. Right: BH effective potential with $m = 0.5\,,\,\ell = 2$ and various coupling constants $\ell _\eta$

4 Time-domain integration scheme

Here, we briefly demonstrate the numerical scheme of time-domain integration, first proposed in [92], which yields the temporal response of a metric perturbation as it propagates on a fixed background spacetime. By defining $u(r_*,t)=u(i\Delta r_*,j\Delta t)=u_{i,j}$, $V(r(r_*))=V(r_*)=V(i\Delta r_*)=V_{i}$, $A(r_*)=A(i\Delta r_*)=A_i$ and $B(r_*)=B(i\Delta r_*)=B_i$ Eq. (3.16) takes the form

$$\begin{aligned} \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j} }{\Delta r^2_*} - \frac{1}{A_i B_i} \frac{ u_{i,j+1} - 2u_{i,j} + u_{i,j-1} }{\Delta t^2} - V_i \,u_{i,j} = 0\;. \end{aligned}$$

(4.1)

Then, by using as initial condition a Gaussian wave-packet of the form $\psi (r_*,t) = \exp \left[ -(r_*-c)^2/(2\sigma ^2) \right]$ and $\psi (\rho _*,t<0) = 0$, where c and $\sigma$ correspond to the median and width of the wave-packet, we can derive the time evolution of u

$$\begin{aligned} u_{i,j+1} = A_i B_i \,\Delta t^2 \left( \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j} }{\Delta r^2_*} - V_i u_{i,j} \right) + 2 u_{i,j} - u_{i,j-1}~, \end{aligned}$$

(4.2)

where the Courant–Friedrichs–Lewy (CFL) condition requires $\Delta t/\Delta r_*< 1/(A_i B_i)$. To calculate the precise values of the potential $V_i$, we integrate numerically the equation for the tortoise coordinate and then solve with respect to the corresponding radial coordinate. Moreover we require the vanishing of perturbations at radial infinity by imposing reflective boundary conditions $u_{i_{max}, j} = 0$, since our solution is effectively asymptotically AdS [8, 93]. For further details regarding the numerical scheme and its convergence see Appendix B.

5 Evolution of axial gravitational perturbations

Regardless of the fact that the spacetime utilized here does not describe a BH immersed in a Universe with a negative cosmological constant, it is meaningful to compare it with Schwarzschild-AdS BHs, since $\ell _\eta$ introduces an effective cosmological scale to the geometry considered. It what follows, we will adopt the categorization from Horowitz and Hubeny [94] regarding AdS BHs determined by two dimensionful parameters, namely the AdS radius $r=R_{AdS}$ and the event horizon radius $r=r_h$. The BH solution in the present study depends also on two parameters: m and $\ell _{\eta }$. The value of mass controls the position of the event horizon $r_{h}$ (and consequently of the photon sphere) and $\ell _{\eta }$ dictates the value of the effective cosmological radius $r=R_\mathrm{eff}$. However, one key difference between the two solutions is that the first is a ‘bald’ BH embedded in an AdS Universe, whereas the scalarized solution is ‘dressed’ with scalar hair whose existence creates the effective AdS-like asymptotics. In this sense, the parameter $\ell _\eta$ controls the strength of the scalar hair and as a consequence the value of $R_\mathrm{eff}$. In our case, the effective cosmological radius is given by $R_\mathrm{eff} = \sqrt{3}\,\ell _{\eta }$. One may categorize BHs in an AdS Universe [94] as (i) small size BHs with $r_h<< R_{AdS}$, (ii) intermediate size BHs with $r_h\sim R_{AdS}$ and (iii) large BHs with $r_h>> R_{AdS}$. We utilize a similar classification to distinguish between small ($r_h<< R_\mathrm{eff}$), intermediate ($r_h\sim R_\mathrm{eff}$) and large hairy BHs ($r_h>> R_\mathrm{eff}$) (see Fig. 2).

In what follows, we apply the numerical procedure outlined above, to calculate the temporal response of axial gravitational perturbations on the BHs of the above categories. In the following figures we obtain the perturbation response at a position arbitrarily close to the event horizon $r-r_h=10^{-5}$, though we have checked that the same results are obtained if we calculate the response at any position outside the event horizon. Furthermore, we have performed some typical tests to ensure the validity of the integration method. Specifically, we have calculated the response of gravitational perturbations on the BH considered here, in the limit where $\ell _\eta \rightarrow \infty$, where the effect of the scalar hair is suppressed, (we have chosen $\ell _\eta =1000$ though even for $\ell _\eta =10$ the potential V(r) converges to the Regge–Wheeler one). By using the Prony method [95] we can extract the spectral content from the temporal response at the large coupling limit and investigate if the modes extracted solely from the prompt ringdown agree with the standard axial gravitational QNMs of Schwarzschild BHs [81]. In Fig. 3 we show the prompt ringing of small BHs for $\ell =2$ gravitational perturbations. We only consider the case where the BH mass is $m=0.1$ in order to obtain a clear ringing phase, since for the range of couplings we considered from 4 to 100 the echo timescales are very large and the extraction of QNMs from the prompt ringdown is possible. Figure 3 indicates that decreasing $\ell _\eta$ has a menial effect on the spectrum while for $\ell _\eta =100$ the extracted mode asymptotes to the fundamental Schwarzschild QNM with accuracy $\sim 0.1\%$.

For completeness, we have further calculated the instability timescale of null geodesics (Lyapunov exponents) at the photon sphere [82] and found that at the GR limit the correspondence between null geodesics and eikonal (large $\ell$) QNMs still holds. This is expected since at this limit the BH approaches Schwarzschild and GWs propagate with the speed of light. Therefore this analysis serves as another validity test of our numerical results and justifies the existence of a modified propagation speed of GWs. In fact, for the case with $m=0.1$, $\ell _\eta =100$ ($m/\ell _\eta =10^{-3}$) the instability timescale of null geodesics and the extracted decay rate of the fundamental QNM for $\ell =10$ (approximately eikonal) axial perturbations have a percentage difference of less than 0.5%. On the other hand, when $\ell _\eta$ is not large enough then the propagation speed of GWs is modified, in accordance with Eq. (3.16), and this leads to a significant inconsistency between null geodesics and eikonal fundamental QNMs, as expected. For example, by choosing $m=5$, $\ell _\eta =1$ ($m/\ell _\eta =5$) the instability timescale of null geodesics and the extracted decay rate of the fundamental QNM for $\ell =10$ (approximately eikonal) axial perturbations have a percentage difference $\sim 50$%. Therefore we conclude that the fundamental eikonal QNMs are not always associated with null geodesics at the spacetime, as it was also shown in [96].

Fig. 2

The difference between the effective cosmological and event horizon radius as a function of m and $\ell _\eta$, where $R_\mathrm{eff} = \sqrt{3}\,\ell _{\eta }$. For $m<<\ell _\eta$ we have $r_h<< R_\mathrm{eff}$, for $m\simeq \ell _\eta$ we have $r_h\simeq R_\mathrm{eff}$, while for $m>>\ell _\eta$ we have $r_h>> R_\mathrm{eff}$. The minima of the curves occur at approximately $m/\ell _\eta \simeq 3$ and indicate the points in the parametric space $(m,\ell _\eta )$ for which $r_h>> R_\mathrm{eff}$

Fig. 3

Left: oscillatory prompt response of axial gravitational perturbations with $\ell =2$ of the BH considered with $m=0.1$ and varying $\ell _\eta$. Right: fundamental $\ell =2$ modes extracted from the prompt ringing phase of the BH considered with $m=0.1$ and $\ell _\eta =4,\,10,\,50,\,100$ from left to right

Figure 4 presents the temporal response of the BH against axial perturbations in the small ratio regime $m<<\ell _\eta$ ($m/\ell _\eta \sim \mathcal {O}(10^{-3})$). The initial quasinormal ringdown is quite similar to that of a Schwarzschild BH. Such behavior is expected and can be attributed to the high value of $\ell _\eta$ relative to the mass. Perturbations with higher angular number $\ell$ decay faster and with higher frequency since more energy is carried away from the photon sphere. This phenomenon is expected since similar behavior appears for gravitational perturbations and QNMs in Schwarzschild BHs [8]. The late time behavior however, shows that, instead of a power-law cutoff, the field settles to a constant value which is related to the asymptotic value that the effective potential acquires (see Fig. 1) and the expectancy of late-time echoes. The eventual late-time tail should be more evident for large BHs since echoes will be washed out rapidly at the event horizon.

In Fig. 5 the evolution of perturbations for $m<\ell _\eta$ ($m/\ell _\eta \sim \mathcal {O}(10^{-1})$) is displayed. The most obvious effect one observes is the emergence of echoes following the initial quasinormal ringdown. In this parametric region the relation between the mass of the BH and the coupling $\ell _\eta$ becomes more transparent. By keeping $\ell _\eta$ fixed and increasing the mass, perturbations will have to travel a shorter distance between the photon sphere and the effective AdS boundary induced by the scalar field leading to repetitions in the signal which appear in shorter timescales. Analogously, similar behavior is obtained when one keeps the mass fixed and decreases the coupling. This pattern was also observed in [65] for the case of scalar perturbations though test scalar fields travel with the speed of light, in contrast to axial gravitational waves in our analysis which have a variable propagation speed (see Eq. (3.16)). This means that a null geodesic analysis, similar to that in [21,22,23] where the echo timescales are approximated by the time that light takes to travel from a boundary to the photon sphere and back, is rendered pointless. Our case is much more intricate since one cannot consider null geodesics anymore but rather has to analyze waves traveling in a dispersive medium with varying propagation speed in different regimes. We have performed a trivial null geodesic analysis and the results we obtained are expected, that is for large $\ell _\eta$ the propagation speed of GW approaches the one of light and the echo timescales can be properly approximated, while as the coupling decreases the echo timescales predicted by null geodesics are completely inconsistent with the actual timescales of echoes obtained by our numerical integration. Nevertheless, such investigation reinforces our discussion regarding the existence of a modified GW speed of propagation.

Fig. 4

Left: effective potential of axial gravitational perturbations, with varying $\ell$, of the scalarized BH with $\ell _\eta = 100$, $m = 0.1$. The parametric region considered here corresponds to small BHs in the sense that $r_h<<R_\mathrm{eff}$ ($m/\ell _\eta \sim \mathcal {O}(10^{-3})$). Right: response of axial gravitational perturbations with respect to the effective potential of BHs considered at the left subfigure

Fig. 5

Left: effective potential of an $\ell =2$ gravitational perturbation of the scalarized BH with $\ell _\eta =5$ for varying mass m. The cases we consider here correspond to intermediate size BHs with $r_h\sim R_\mathrm{eff}$ ($m/\ell _\eta \sim \mathcal {O}(10^{-1})$). Right: response of axial gravitational perturbations with respect to the effective potential of BHs considered at the left subfigure

To obtain a complete picture regarding the effect of the ratio $m/\ell _\eta$ on the BH’s response to fluctuations we have plotted in Fig. 6 the time evolution of perturbations for a wide range of masses keeping $\ell _\eta$ fixed. As m grows the echoes are replaced by quasinormal oscillations, while further increment of the mass leads to a single quasinormal ringdown followed by a late-time tail. We conclude that this behavior stems from the shape of the effective potential which decreases in amplitude as m increases. This leads to an increasingly smaller region where trapped modes, which lead to echoes, can occur, and thus the quasinormal ringing of the BH dominates over the echoes which are quickly suppressed.

When the mass becomes proportional ($m/\ell _\eta \sim \mathcal {O}(10^0)$) or significantly larger than $\ell _\eta$ ($m/\ell _\eta \sim \mathcal {O}(10^2)$), negative wells develop in the effective potential in the vicinity of the event horizon (see Figs. 7, 8). Despite the negative well formation, the time-domain profiles show an exponential decay of the signal without any indication of a linear instability. On the contrary, more massive objects lead to signals with shorter quasinormal ringing stages, due to the absence of a photon sphere peak, and with faster decay rates even though the corresponding effective potentials develop even deeper negative wells. The exponential nature of the eventual late-time behavior of perturbations is related to the effective AdS asymptotics of our spacetime which requires the imposition of reflective boundary conditions at infinity and is in agreement with what occurs in perturbations of AdS BHs [97]. The tail in these cases appears because echoes are subdominant and vanish very rapidly at the event horizon, thus the asymptotic behavior is probed faster. We expect that even perturbations of the small BHs in study will eventually possess an exponential tail but at much later times which our numerical scheme cannot probe.

Fig. 6

Left: effective potential of an $\ell =2$ gravitational perturbation of the scalarized BH with $\ell _\eta =0.1$ for varying mass m. Right: response of axial gravitational perturbations with respect to the effective potential of BHs considered at the left subfigure

Fig. 7

Left: effective potential of an $\ell =2$ gravitational perturbation of the scalarized BH with $\ell _\eta =1$ for varying mass m. The cases we consider here correspond to intermediate size BHs with $r_h\sim R_\mathrm{eff}$ ($m/\ell _\eta \sim \mathcal {O}(10^{0})$). Right: response of axial gravitational perturbations with respect to the effective potential of BHs considered at the left subfigure. The responses were shifted in time for illustration purposes

Fig. 8

Left: effective potential of an $\ell =1$ gravitational perturbation of the scalarized BH with $\ell _\eta =5$ for varying mass m. The cases we consider here correspond to large BHs with $r_h>> R_\mathrm{eff}$ ($m/\ell _\eta \sim \mathcal {O}(10^{2})$). Right: response of axial gravitational perturbations with respect to the effective potential of BHs considered at the left subfigure. A shift in time was applied on the responses for illustration purposes

Fig. 9

Time evolution of an $\ell =2$ gravitational perturbation of the scalarized BH with fixed ratio $m/\ell _\eta \simeq 0.1$ (left), $m/\ell _\eta \simeq 0.3$ (middle) and $m/\ell _\eta =1$ (right)

From the above, we conclude that the scalarized BH spacetime is modally stable under axial gravitational perturbations, where the qualitative features of the response depend solely on the ratio $m/\ell _\eta$. In Fig. 9 we demonstrate the above statement for the case of intermediate size BHs. Our numerics show that a similar analogy occurs irregardless of the BH’s size. Even though the source of echoes in our BH is related to the asymptotics of spacetime, and not to the nature of the near-horizon structure, our results are in accordance with perturbations in wormholes with decreasing throat radii [71] and black-bounce models, which transition from regular BHs to wormholes [98].

6 Conclusions

In this work we studied static and spherically symmetric solutions of a Horndeski subclass which includes a massless scalar field non-minimally coupled to the Einstein tensor. Such theory admits an exact BH solution ‘dressed’ with scalar hair whose existence induces an effective negative cosmological constant even though the BH does not reside in an AdS Universe. We have studied the modal stability of such solutions under axial gravitational perturbations, with time evolution techniques, and complementary QNM extraction, that solve the linearized gravitational wave equation. Our results designate that the BH under study is linearly stable against axial perturbations, with decaying temporal responses akin to ringdown waveforms. The qualitative features of the ringdown waveform depend solely on ratio of the two available parameters of spacetime, namely the BH mass m and non-minimal coupling strength $\ell _\eta$. We have further demonstrated that as $m/\ell _\eta$ increases, we have gravitational-wave ringdown transitions between three distinct response patterns, namely a state with a typical quasinormal ringdown ($m/\ell _\eta \lesssim 10^{-2}$), an intermediate long-lived state which exhibits gravitational-wave echoes ($10^{-2}\lesssim m/\ell _\eta \lesssim 10^{-1}$) and a state where the ringdown and echoes are depleted rapidly to give turn to an exponential tail ($m/\ell _\eta \gtrsim 10^{-1}$).

Regardless that our findings point towards linear stability, we only considered the axial sector of gravitational fluctuations. In generality, one must investigate the polar sector of gravitational perturbations as well in order for a complete stability analysis to be established. This extension can be extremely challenging with what regards the achievement of writing the perturbation equation into a one-dimensional Zerilli-like equation and the stability of spacetime itself, since the polar degrees of freedom generically couple to the scalar hair in scalar-tensor theories. A first step towards the aforementioned direction is the consideration of radial perturbations which are a good proxy to polar ones [99,100,101,102]. Radial perturbations can also couple the scalar field with the metric components, thus can serve as more sensible probe to the overall linear stability of the hairy BHs under consideration.

Besides dealing with temporal evolution techniques, another interesting direction would be a complete frequency domain analysis of axial and polar gravitational QNMs which is still lacking in the particular family of BH solutions, in a similar manner as in Refs. [103, 104] where scalar QNMs have been discussed. Furthermore, since the BH geometry in study possesses a propagation speed for GWs that differs from that of light, it will be paramount to investigate potential observational imprints in order to disentangle possible degeneracies between GW phase modifications and environmental effects and avoid misinterpreting GWs in modified gravity with strongly-lensed GR GWs [105, 106].

Finally, in a recent analysis [107], a class of mechanical models were studied, where a canonical degree of freedom interacts with another one with a negative kinetic term, i.e. with a ghost. Surprisingly, it was shown that the classical motion of the system is completely stable for all initial conditions, even though one would expected that such system to be unstable due to the presence of a ghost field. In our case, we have dealt with a conceptually analogue system, consisting of a scalarized BH for which the kinetic energy of the scalar hair can be positive or negative (first degree of freedom) provided that the strength of the non-minimal coupling to the Einstein tensor has the opposite sign (second degree of freedom), being attractive of repulsive respectively. Regardless of the case, we find that the BH is stable under axial perturbations, thus providing an illustration that the classical mechanics analysis in [107] can potentially apply to BH physics.

Appendix

1.1 A. Solution of the angular differential equation

In this appendix, we present the solution of the angular part of the differential Eq. (3.10). The corresponding differential equation is:

$$\begin{aligned} \sin ^3\theta \frac{d}{d\theta }\left[ \frac{1}{\sin ^3\theta }\frac{d\mathcal {S}(\theta )}{d\theta }\right] +\mathcal {A}\, \mathcal {S}(\theta )=0~, \end{aligned}$$

(A1)

where $\mathcal {A}$ is the separation constant. By performing a change of variables of the form $x=\cos \theta$ we obtain the following differential equation:

$$\begin{aligned} (1-x^2)\mathcal {S}''+2x\mathcal {S}'+\mathcal {A}\mathcal {S}=0~. \end{aligned}$$

(A2)

Note that this is very similar to the Legendre differential equation albeit with one sign change. This differential equation is called the ultraspherical or Gegenbauer differential equation. There exist three alternate forms of the equation that yield the same result. We are going to show the two we are interested in here.

First form:

$$\begin{aligned} (1-x^2)\mathcal {S}''-2(m+1)x\mathcal {S}'+(\ell -m)(\ell -m+1)\mathcal {S}=0~. \end{aligned}$$

(A3)

The first form has the following solutions

$$\begin{aligned} \mathcal {S}=(x^2-1)^{-m/2}\left[ C_1P^m_{\ell }(x)+C_2\mathcal {Q}^m_{\ell }(x)\right] ~, \end{aligned}$$

(A4)

where $P^m_\ell (x)$ and $\mathcal {Q}^m_\ell (x)$ are the Legendre functions of the first and second kind, respectively. Note that $m=-2$.

Second form:

$$\begin{aligned} (1-x^2)\mathcal {S}''-(2n+1)x\mathcal {S}'+k(k+2n)\mathcal {S}=0~. \end{aligned}$$

(A5)

The second form has the following solutions

$$\begin{aligned} \mathcal {S}=(x^2-1)^{(1-2n)/4}\left[ C_1P^{n-1/2}_{-1/2+k+n}(x)+C_2\mathcal {Q}^{n-1/2}_{-1/2+k+n}(x)\right] ~. \end{aligned}$$

(A6)

Note that $n=-3/2$. If $-1/2+k+n$ is an integer, then this solution yields the Gegenbauer polynomials, $C^n_k$. Equating (A3) with (A4), yields $-1/2+k+n=l\rightarrow k=\ell +2$, and thus we are left with the solution of

$$\begin{aligned} \mathcal {S}=C^{-3/2}_{\ell +2}(\theta )~. \end{aligned}$$

(A7)

Therefore, the differential equation (A2) yields the separation constant $\mathcal {A}=(\ell +2)(\ell -1)$ .

1.2 B. Convergence tests

Here, we discuss in depth our numerical scheme which is briefly analyzed in Sect. 4. The essential equations in play are Eqs. (4.1) and (4.2), together with the CFL condition and the vanishing of perturbations at radial infinity. In terms of the tortoise coordinate $r_*$, we observe that when r tends to infinity, $r_*$ tends to a finite constant which we denote as $r_*^{max}$. The implications of the behavior of $r_*$ are twofold: firstly, the reflective boundary condition in terms of $r_*$ takes the form $u(r_*^{max},t) = u_{i_{max}, j} = 0$ and secondly, our region of interest in the $(r_* - t)$ diagram lies on the left of the vertical line $r = r_*^{max}$ as seen in Fig. 10.

Fig. 10

Diagram of the numerical grid in the $(r_* - t)$ plane. The points on the red lines are determined through the boundary and initial conditions. The points on the green line $u(r_*=0,t)$ are the results shown in this paper and they correspond to $r=1.00001\,r_{h}$. The rest of the grid points are calculated through the recursive relation (4.2) starting from the point $u_{i_{max}-1,1}$. The blue arrows demonstrate a graphical depiction of the time evolution Eq. (4.2) for the point $u_{i_{max}-1,1}$

It is important to note that the values of the finite constant $r_*^{max}$ are proportional to the value of the coupling $\ell _\eta$ i.e. $r_*^{max} \sim \ell _\eta$ (see Table 1). This means that the value of $\ell _\eta$ dictates the range of $r_*$ since $r_* \in (-\infty \,,\,r_*^{max}]$. A second important consequence of the above proportionality is that, as $\ell _\eta$ increases we also need to increase the number of grid points N in order to keep the value of $\Delta r_*$ sufficiently small. To better understand why this is occurring we need to delve into the technical details concerning the procedure executed by our code.

The first step is to find the function $r(r_*)$ by numerically solving the differential equation of the tortoise coordinate

$$\begin{aligned} \frac{dr(r_*)}{dr_*} = \sqrt{\frac{f(r(r_*))}{g(r(r_*))}}, \end{aligned}$$

(B1)

together with the condition $r(r_*=0)=1.00001\,r_{h}$ which fixes the integration constant. Hence, after the integration we have $r(r_*\rightarrow -\infty )\rightarrow r_{h}$, $r(r_*=0) = 1.00001\,r_{h}$ and $r(r_*\rightarrow r_*^{max})\rightarrow \infty$, meaning that $r_* \in (-\infty \,,\,r_*^{max}]$. However, in order to define a numerical grid with which we will perform the time evolution of u, we need to work on a finite interval of $r_*$. We do so by choosing a sufficiently large negative value¹ (which we denote by $r_*^{min}$) as the second end of the interval of $r_*$. Thus, in the context of the numerical integration we will work on the interval $r_*^{numerical} \in [r_*^{min}\,,\,r_*^{max}]$ even though in principle $r_* \in (-\infty \,,\,r_*^{max}]$.

The final step of our code which calculates the time domain profiles expects as inputs the values of $r_*^{min}\,,\,r_*^{max}$ and N in the $r_*$ direction. It then calculates the spatial step $\Delta r_*$ from the relation

$$\begin{aligned} \Delta r_* = \frac{r_*^{max} + |r_*^{min}|}{N} \end{aligned}$$

(B2)

and the time step from $\Delta t = c\, \Delta r_*$ where c is positive constant value satisfying the CFL condition that should not be confused with the speed of light. The fact that $r_*^{min}$ is constant throughout all of our evolutions and that $r_*^{max} \sim \ell _\eta$ implies, through Eq. (B2), that as $\ell _ \eta$ increases we also need to increase N in order to keep the value of $\Delta r_*$ sufficiently small (see Table 2).

Table 1 Reference values indicating the analogy $r_*^{max} \sim \ell _\eta$

Table 2 The corresponding grid spacing $\Delta r_*$ for various choices of grid points N

Fig. 11

Left: convergence curve for $m=0.1\,,\,\ell _\eta =100$ corresponding to a case where we obtain a clear prompt ringdown. As $u(r_*,t)|_{best}$ we choose the value of u for $N=12,000$ grid points, i.e. $u(r_*,t)|_{12,000}$, indicating the best approximation. All the points are extracted at $r_*=0$ and $t/m=225.228$. Right: Convergence curve for $m=0.5\,,\,\ell _\eta =5$ corresponding to a case where we obtain echoes after the initial ringdown. As $u(r_*,t)|_{best}$ we choose the value of u for $N=2800$ grid points, i.e. $u(r_*,t)|_{2800}$, indicating the best approximation. All the points are extracted at $r_*=0$ and $t/m=234.637$

Finally, we can produce convergence curves to provide some quantitative information regarding the accuracy of our numerical integration scheme. To produce these curves, we first calculate the values of $u(r_*,t)$ at a given point as the grid spacing $\Delta r_*$ is reduced by increasing N. We will denote these values by $u(r_*,t)|_N$. We then use the value $u(r_*,t)$ for the maximum number of grid points (i.e. the smaller grid spacing $\Delta r_*$) as a reference value indicating the best approximation to the true value of u at that point. We will denote that value by $u(r_*,t)|_{best}$. To calculate the error we subtract each value $u(r_*,t)|_N$ for every different N from the value of the best approximation and then take its absolute value i.e.

$$\begin{aligned} |\text {Error}|_{N} = \Big | u(r_*,t)|_{best}-u(r_*,t)|_{N} \Big |~. \end{aligned}$$

(B3)

The diagrams in Fig. 11 demonstrate that our code achieves numerical convergence irrespective of whether our compact object responds with a clear ringdown or with a signal with echoes i.e. in rather different regions of our parametric space $(m,\ell _\eta )$. Even though for the first case of Fig. 11 on the left the code converges rapidly, we expect that the same will occur for the second case depicted in Fig. 11 on the right if we further increase the number of grid points. As a final note, we stress the fact that even though the chosen values of the grid points N are very different for the convergence curves in Fig. 11, the corresponding grid spacing $\Delta r_*$ is of the same order of magnitude for both cases, as can be seen in Table 2.