Stable spontaneously-scalarized black holes in generalized scalar-tensor theories

It has been shown that the synergy of a scalar field coupling with both the Ricci scalar and the Gauss-Bonnet invariant significantly affects the properties of scalarized black holes and neutron stars, including their domain of existence and the amount of scalar hair they carry. Here we study the radial stability of scalarized black-hole solutions. We demonstrate that they are stable against radial perturbations for Ricci couplings consistent with both a late-time cosmological attractor and the evasion of binary pulsar constraints. In addition, we investigate the effect of the Ricci coupling on the hyperbolicity of the equation governing linear, radial perturbations and show that it significantly reduces the region over which hyperbolicity is lost.


I. INTRODUCTION
The first direct detection of gravitational waves in 2015 [1] by the LIGO-Virgo collaboration signaled the start of a new era in astrophysics. Observations of gravitational waves generated by merging compact objects, namely black holes (BH) and neutron stars (NS), provide a powerful tool to probe the strong field, dynamical regime of gravity for the first time in history. Despite the success of General Relativity (GR) in the weak field limit, deviations from GR could become relevant when gravitational effects are extreme. Indeed, gravitational-wave observations have already been used to place constraints on theories that seek to modify GR and look beyond the Standard-Model framework [2][3][4][5][6].
Considering the overwhelming success of GR in the weak-field regime, noteworthy theories are those which exhibit novel gravitational effects that only manifest in the vicinity of black holes and compact neutron stars but lie dormant elsewhere. To that end, spontaneous scalarization is a strong gravity effect in which a phase transition endows black holes and neutron stars with an extra scalar configuration. The phase transition is induced by a coupling of the scalar field with curvature invariants, and the details of this coupling control the transition threshold. At a linear level, such couplings trigger a tachyonic instability, and consequently the development of scalar hair. The phenomenon was originally identified in Ref. [7] for a class of theories where it only takes place for compact stars. However, these theories are unable to induce scalarization on black hole spacetimes as they are subject to the no-hair theorems, [8][9][10][11], and so cannot support non-trivial scalar profiles. Recently, it been shown that another class of theories can exhibit black hole scalarization, [12,13], by including a coupling with the Gauss-Bonnet (GB) invariant, G = R µνλσ R µνλσ − 4R µν R µν + R 2 . Here R µ νλσ denotes the Riemann tensor, R µν the Ricci tensor, and R the Ricci scalar.
In particular, Ref. [13] considered the action where κ = 8πG/c 4 , X = −(∂φ) 2 /2, and showed that stationary and asymptotically flat black holes are described by the Kerr metric provided f (φ 0 ) = 0 for some constant φ 0 and f (φ 0 )G < 0. Further, it was found that −f (φ 0 )G acts as an effective mass squared for the scalar perturbation; hence, when f (φ 0 )G becomes positive and sufficiently large it triggers a tachyonic instability. A f (φ) ∝ φ 2 coupling was considered, being the simplest model that satisfied the conditions for developing the instability. Indeed, scalarization was found to occur and the scalarized black hole solutions exist in the same region of parameter space as the tachyonic instability.
Ref. [12] focused instead on the case where f (φ) ∝ e φ 2 , also showing that tachyonic instabilities occur and that scalarized black holes again exist in the same region of the parameter space.
Extensions of these models have examined the effects of additional terms in the action, including a bare mass and higher order corrections, different fields, or different types of instabilities as triggers of scalarization, e.g. [14][15][16][17][18][19][20]. In more recent work, it has been shown that a tachyonic instability leading to scalarization can be triggered by spin [21] and spacetimes describing such scalarized black holes have been generated [22,23]. The onset of scalarization is controlled by terms that contribute to linear perturbations around a GR background in all of these scenarios. The minimal action that contains all such terms for scalar-tensor theories, up to field redefinitions, and leads to second order equations upon variation was identified in Ref. [24]: arXiv:2204.01684v1 [gr-qc] 4 Apr 2022 where κ = 8πG/c 4 , X = −(∂φ) 2 /2, m φ is the bare mass of the scalar field, β is a dimensionless parameter, and α has dimensions of length squared. Since vacuum spacetimes in GR satisfy G µν = 0, the derivative coupling to G µν and the coupling between the scalar and R will not contribute to the onset of the tachyonic instability for black holes. These terms will however affect scalarization thresholds for neutron stars [25] and also affect the end state of scalarization for both neutron stars [26] and black holes [27]. Indeed, the tachyonic instability is quenched by non-linearities. Therefore, these terms affect the final properties of the scalarized configurations, as do any other terms that have been neglected in action (2) because they do not contribute to the linear perturbation equation. Understanding how the various (self)interactions beyond the scalar-Gauss-Bonnet coupling affect the scalar profile of a scalarized compact object is essential from an observational perspective. It has been shown that scalarized black holes are unstable under radial perturbations in the simplest, quadratic coupling scenario [15]. This issue can be overcome if one considers an additional quartic interaction in the Gauss-Bonnet coupling function, provided that the sign of the quartic coupling coefficient is opposite to the quadratic one [16]. However, addressing the instability with quartic, or exponential couplings is not entirely appealing from an effective field theory (EFT) perspective. This is because these terms have a higher mass dimension than other terms that could contribute non-linearly, e.g. a simple φ 4 self interaction. It was, indeed, shown in Ref. [17] that including self interactions for the scalar can lead to radially stable scalarized solutions.
More recently, Ref. [27] provided strong indications that a coupling of the scalar field with the Ricci scalar, already present in action (2), can also lead to radially stable solutions. This would be particularly interesting if proven to be true as this same coupling has already been shown to address observational viability issues for black hole scalarization models when β > 0: it dominates the scalar dynamics in the late-time cosmology and turns GR into a cosmological attractor [28]; it can also suppress scalarization of neutron stars [26] and hence, remove binary pulsar constraints. In this paper, we perform a radial perturbation analysis for scalarized solutions and fully explore the role of the coupling with the Ricci scalar in the stability of the solutions.
The structure of the paper is as follows: in Sec. II we, first, introduce our model and present the field equations. We reproduce the results already derived in previous work, concerning the background BH solutions, that describe static and spherically symmetric configurations. Then, we consider radial perturbations to the static background and discuss the properties of our model as a Schrödinger-like problem. In Sec. III we discuss our numerical results regarding the stability of the solutions and the hyperbolicity of the scalar perturbation equation. Finally, in Sec. IV we present our conclusions.

A. Action and field equations
We will consider the following action motivated in part by simplicity and in part by the fact that we seek to understand the influence of the coupling between the scalar and R in stability considerations. One can think of this action as part of an EFT in which the scalar enjoys Z 2 symmetry and shift symmetry is broken by the couplings to curvature. Note that the model considered in Refs. [7,29] is equivalent to action (3) with α = 0 in linear perturbation theory [24] and our definition of β matches that of Refs. [7,29]. The field equations for the metric that one derives from action (3) are where The scalar field equation reads where the effective scalar mass is given by

B. Spherical Black Hole Solutions
The background considered is static and spherically symmetric, and described by the metric with a scalar field that depends only on the radial coordinate, φ = φ(r). The differential equations describing the unknown functions (A(r), B(r), φ(r)) can be obtained from (4) and (6) which are explicitly shown in Appendix A.
In [27] BH solutions that are asymptotically flat, and continuously connected to GR were found and their characteristics were studied. From these solutions one can infer the ADM mass M and scalar charge Q from the large r behaviour of B(r) and φ g rr For the expansions including terms up to O(r −6 ) refer to [27]. The solutions depend on the scaling of the BH mass and scalar charge with respect to the parameter α, controlling the GB coupling. When analyzing the stability of the solutions, a dimensionless re-parametrization proves particularly useful: The domain of existence for black holes with non-zero scalar charge is non-trivial. To analyze the parameter space in which these solutions exist, we scan (α/r 2 h , β) for BHs with a non-trivial scalar field configuration that vanishes asymptotically, in accordance with [27]. For α = 0 the only solution that is regular at the horizon and asymptotically flat is the Schwarzschild BH [10], while for β = 0, the allowed values for α are in agreement with [12,13] The left panel of Fig. 1 shows the existence domain in the (α/r 2 h , β) plane. Scalarized solutions with the desired properties, i.e. regular everywhere and asymptotically vanishing, exist in the shaded regions. We can classify the space of solutions existing in three domains bounded by a seemingly vertical line given by the critical value for scalarization α/r 2 h ≈ 0.18 and a parabolalike curve defined by the existence condition presented in Ref. [27].
To further analyze the solutions, in the center panel of Fig. 1, we show the domain of existence in the (M/r h , Q/r h ) plane, which is suggestive of several properties of these solutions. First, it appears that the map (α/r 2 h , β) → (M/r h , Q/r h ) is invertible offering a direct connection between the asymptotes of the solution and the underlying gravity model. The lines of constant β approaching the vertical line in the left figure, merge at (M/r h , Q/r h ) = (0.5, 0), which corresponds to the GR solution. These same lines appear to be bounded from above as they approach the parabola in (α/r 2 h , β). The BH solutions and their stability are better understood when characterized by the quantities (M ,Q) [17,30], as in the rightmost panel of Fig. 1. For clarity, in Fig. 2 we show the same domain but for select values of β.
Schwarzschild BHs are radially unstable forM <M c ≈ 1.174 [15,17]. It can be seen thatM c is the point where the curves converge in Fig. 2. WhenM <M c , scalar perturbations grow spontaneously and form a non-trivial scalar profile, and so are the energetically favourable solutions. We again stress that the final scalar field profile is determined by the full non-linear equations and not only the terms that trigger the instability.
In contrast, Schwarzschild BHs are stable in the region M >M c , and potential scalarized solutions would decay back to GR, as seen from purely energetic arguments. This reasoning can be confirmed by performing a radial linear stability analysis of the scalarized BH solutions.
As we increase the value of β, the curves start to occupy the regionM <M c , and their gradients become negative. In accordance with [17,27], we expect a negative gradient to correspond to radially stable solutions. From the right panel of Fig. 2, we see that the critical value of β for which part of the curves starts existing in the regionM <M c is approximately β ≈ 1.15, while the window of potentially stable solutions widens when β increases. Moreover, beyond β c ≈ 1.2 we do not find any solutions that we would expect to be radially unstable. Note that while the right panel in Fig. 1 shows that there are potentially stable solutions with β < 0, it is the β > 0 case that renders GR a cosmic attractor [28], so we will focus on β > 0 henceforth. While the arguments above provide a preliminary indication of the scalarized solutions' stability, it is necessary to investigate the problem by performing a full stability analysis before reaching a conclusion.

C. Radial perturbations
We can investigate the radial stability of the BH solutions by employing a perturbative approach. This not only allows us to further elucidate the timescale of the instability, but also to analyze possible modifications to the oscillatory spectrum of the BHs. It is worth noting that since a scalar degree of freedom is absent in GR, the radial modes contribute only to a shift in the mass in that case, and hence are not radiative in nature [31].
We start by considering time-dependent radial perturbations of the metric tensor and scalar field over a static and spherically symmetric background: where A 0 , B 0 and φ 0 are the time-independent background solutions, while A 1 , B 2 and φ 1 are the timedependent perturbations. We can then write down a system of equations for (A 1 , B 1 , φ 1 ), by substituting the metric and scalar perturbations in Eqs. (4) and (6). This system can be reduced to a second order partial differential equation system for φ 1 [15], namely where the coefficients depend only on the background solution.
When considering a quadratic coupling function with the Gauss-Bonnet term, Ref. [15] pointed out that for some values of the coupling, the equation describing the perturbations is not hyperbolic. While this can hinder the investigation of linear stability as a Cauchy problem, we can still examine the mode structure of the spacetime by looking into the frequencies ω.
To perform a mode analysis of the spacetime, we search for the natural frequencies of the system ω, such that φ 1 (t, r) = φ 1 (r)e −iωt . The perturbation equation (14) can be manipulated to the more familiar Schrödinger form: where φ 1 (r) = F (r)ψ(r) and the tortoise coordinate is defined through dr * = g(r) dr. We also define The real part of the frequency ω describes the resonant modes of the system, i.e., for which frequencies an initial perturbation would respond. The imaginary part of the frequency indicates the system's (modal) stability. For modes with a negative imaginary part, an initial perturbation decays exponentially, while when positive the perturbation grows and the system is rendered unstable. Generally, the effective potential identified from Eq. (16) is useful when one attempts to verify the presence of unstable modes as a negative value for the integral of the effective potential with respect to the tortoise coordinate indicates the existence of unstable modes [32], namely It is worth noting that while the condition (18) is sufficient in indicating the presence of unstable modes, it is not a necessary condition.
In what follows we will be making use of the compact- ified coordinate which maps all of spacetime, from the black hole horizon up to infinity, to a finite region, namely x ∈ [0, 1]. In Fig. 3 we plot the effective potential for two random normalized massesM corresponding to the left and right parts of the left panel of Fig. 2 To explore the modal structure of the spacetime, we have to impose proper boundary conditions in order to obtain the modes. These correspond to an ingoing wave at the horizon and an outgoing one at infinity We can see see that, for modes with ω I > 0 (unstable), they simplify to φ 1 (x → 0, 1) = 0. To finish this section, let us mention that the equations describing the radial perturbations in the Schwarzschild spacetime can be directly obtained from (14) by setting φ 0 = 0 and A = B = 1 − r h /r which specify g, C and U . The resulting scalar wave equation is where dr * = dr/(1 − r/r h ) is the tortoise coordinate of the Schwarzschild spacetime. The effective potential is identified as where the index d stands for decoupling. In Fig. 4 we plot the potential in the decoupling limit for some values of α, using the tortoise coordinate to improve visualization. Just as before, the equation describing radial perturbations can be used to access the stability properties of the Schwarzschild spacetime in sGB. From a direct integration of the potential, it is straightforward to see that (18) in this case yields α/r 2 h 0.208.

III. NUMERICAL RESULTS
In order to find the modes of scalarized BHs, we follow the direct integration method presented in previous works on the same subject [15,16]. We briefly summarize the method here.
After picking a value for ω, we integrate Eq. (14) from the horizon and infinity, using in-going and-outgoing waves as boundary conditions respectively 1 . In practice, the integration starts from finite values very close and very far away from the horizon's position, such that the potential is small. This changes the boundary conditions which are no longer purely in-going and out-going waves, but are rather given by the Taylor expansion of the field at the horizon and infinity. Using the two separate solutions, we can demand that they are linearly dependent on a given frequency ω. This is done by examining the Wronskian, given by where φ (−) 1 represents the solution obtained by integrating from the horizon and φ (+) 1 the solution obtained from infinity. The Wronskian vanishes when the value of ω is a QNM frequency. An alternative approach is to integrate from the horizon using the in-going boundary condition, up to a large distance r ∞ /r h . Then one can decompose the solution at infinity onto in-going and out-going waves, and should the value of ω be a QNM frequency, the ingoing amplitude is zero. For both methods, an initial choice for ω is made, and root-finding algorithms are employed to solve for the QNM frequency. This is usually called the shooting method, as one starts at one end, "shooting" for the value of ω for which the boundary condition is satisfied at the other end.
Both methods are suitable for finding stable quasinormal modes with large quality factors, i.e., with large |ω r /ω i |. This usually means that the fundamental mode is easily found. The reason is that the function describing radial perturbations grows exponentially in r, as ∼ e ωir * . The method also works remarkably well with unstable modes, as the perturbations decrease exponentially with r. As such, this method is reliable in identifying regions where BHs are linearly unstable.
BHs solutions of the theory (2) were already studied to some extent in Ref. [27] and we revisited some in Figs. 1 and 2. Now we present new insights considering the stability and modes of the Schwarzschild and scalarized BHs.

A. Radial oscillations and the existence of purely imaginary modes
While the Schwarzschild BH is a solution in both GR and sGB theories, its dynamical response to perturbations can be completely different. In fact, it is precisely this difference that allows for spontaneous scalarized BHs, where the radial perturbations instabilities lead to the scalar hair. Here we investigate the radial mode structure of BHs in sGB, elucidating some major points considering Schwarzschild and scalarized solutions.

Scalarized BHs
Let us look into the solutions presented in Fig. 2. As discussed, we expect to find unstable modes for the re-gionM >M c , where the Schwarzschild BH is stable and energetically favourable. We performed a search for the modes using the shooting method. In Fig. 5 we plot the unstable mode frequencies for the scalarized solutions considering different values for β, some of them presented in the right panel of Fig. 2. These frequencies are purely imaginary. The imaginary part of the Schwarzschild fundamental mode is also plotted, and it is independent of our choice of β. We notice that all scalarized solutions withM >M c ≈ 1.175 are unstable.
In agreement with the predictions made by observing Fig. 2, for β > β c ≈ 1.15, the behavior of the curves changes. We notice that the curve for which β = 1.2 begins not fromM c , but from someM ≈ 1.168 <M c . This means that the minimum BH mass for this parameter is no longerM c , but smaller. Once again, the reasoning for this can be understood qualitatively from purely energetic arguments. For β = 1.2, we can have two scalarized solutions for a given massM , presenting different charges. Solutions with higher charges are unstable, decaying to the scalarized BH with a smaller charge. A similar feature was observed for the case of scalarized BHs with self-interaction [17]. Overall, the timescale of the instability τ = |ω −1 I | increases as β increases, indicating the shift from unstable to stable solutions.

Radial modes for the Schwarzschild spacetime
Since Schwarzschild BHs are stable forM >M c and the dynamical response is different from GR, it is natural to investigate the impact of sGB terms on the Schwarzschild BH spectrum. The fundamental mode was already analyzed in Ref. [30], elucidating how the transition from stable to unstable modes occurs. Here we take an additional step, looking into the first overtone. The results of this subsection are independent of β [cf. Eq. (21)].
To find the modes for the Schwarzschild BH we use the continued fraction (CF) method, as presented in Appendix B (see also Ref. [30]). In Fig. 6 we show the results. Starting from the rightmost part of the plots, we see that forM 1 (or equivalently α 2 M ) we recover the modes for a minimally coupled scalar field in GR, as expected. As we decreaseM , both the real and the imaginary part of the fundamental mode (n = 0) approach zero monotonically atM =M c . Beyond that point the frequency becomes purely imaginary and the corresponding mode unstable. For the first overtone (n = 1), however, we notice that the real part seems to approach zero before the instability occurs. This would imply that Schwarzschild BHs withM ≈ 1.87 in this theory have a purely imaginary mode, with the first one being ωr h ≈ −i0.55. The real part of the first overtone goes to zero again when this mode becomes unstable (purely imaginary), in the regionM <M c . Higher overtones present a similar behavior, suggesting that we might have many of these purely imaginary modes even for stable Schwarzschild BHs in sGB theories.
We highlight here that the usual CF method is not suitable to study purely imaginary modes [33,34]. As the real part of the modes decreases, more terms in the CF expansion are needed in order to compute reliable values for the QNM frequencies. Typically, near the purely imaginary mode, we consider N = 5 × 10 4 or more terms in the CF expansion. Curiously, only recently this class of modes was computed with a reliable precision for the Kerr spacetime [33,34]. We shall not attempt to generalize this method for sGB theories, but it would be interesting to further investigate the full spectrum of the Schwarzschild spacetime in the theory in light of such tools.
Purely imaginary modes for BHs are not an exclusive feature of sGB theories. As mentioned above, even within GR, the existence of such a class of modes has been known for quite some time. For instance, Schwarzschild BHs in GR have algebraically special modes that are purely imaginary [35], and many works have studied these modes for Kerr BHs. While the physical implications of these purely imaginary modes in BH spacetimes are still poorly understood, it is interesting to see them arising in the first overtone for Schwarzschild BHs in sGB. This feature, which seems to be an exception for GR BHs, seems to be common in BHs in sGB.

B. The hyperbolic nature of the equations
An important property of partial differential equations pertains to hyperbolicity. In physical theories, hyperbolic equations are necessary to describe the time evolution of initial data. In GR, where equations are quasilinear, proof of strong hyperbolicity (e.g. through the harmonic gauge choice) establishes the well-posedness of the theory (see Ref. [36] for a pedagogical discussion). In linear equations such as (14), the situation is fairly straightforward. Eq. (14) is hyperbolic provided that g 2 (x) > 0.
It has been shown that for the exponential coupling, the equation describing the radial perturbations is not hyperbolic for a variety of solutions [15]. Interestingly, even though the equations are not hyperbolic, we can still proceed to search for unstable modes. It was shown that for scalarized BHs in the exponential model with low mass, an unstable mode arises in the region where hyperbolicity is broken. Here we investigate the impact of the Ricci term on the hyperbolicity of the equation describing radial perturbations.
We start by analyzing the behavior of the coefficient g(r) 2 in the near-horizon regime. In this section, we shall replace the quadratic coupling with the Gauss-Bonnet term by a generic function of the form αφ 2 → αf (φ), in order to compare our results with other works in the literature. Note that the definition of the normalised charge and mass are unchained under this substitution. Using the expansion of the background near the event horizon we find that the coefficient g(r) 2 appearing in Eq. (14) behaves as where In order to investigate whether the Ricci term helps maintaining the hyperbolic nature of the equation, we look into large positive values of β that, from our previous analysis, we know correspond to stable scalarized solutions. We find which indicates that the Ricci term in the action acts in favor of the hyperbolic character of the equation. Note, however, that in order to verify whether the equation maintains its hyperbolicity we need to properly solve the BH solution and obtain the coefficient g(r) 2 numerically. Now, let us consider a scenario that is known to break hyperbolicity and check the influence of the Ricci term. We consider the exponential model presented in Ref. [15], where looking into solutions with a fixed φ h and varying β. In Fig. 7 we compare two scalarized solutions with φ h = 2 and 2.5. We observe that the radial domain in which the perturbation equation is non-hyperbolic decreases as β increases, as predicted by Eq. (27). Further, we observe that for some limiting value of β the region with g 2 < 0 seems to vanish and the solution is hyperbolic for all r > r h . We note, however, that as we approach this threshold the solutions for a given φ h are increasingly hard to find, and beyond the threshold solutions cease to exist. It seems that while the β-term improves the hyperbolicity of the perturbation equation, it still is not enough to ensure it for all values of (r, φ h ) for a given β. To further illustrate the hyperbolic properties of the equations as a function of the background solution, Fig. 8 shows curves of constant β and varying φ h in the normalised charge-mass plane for the exponential model. The dotted part of the curves corresponds to regions in which the radial perturbation equations are not hyperbolic. That the β = 2 curve ends is indicative of our inability to find any solutions past this point. It seems that while the additional term helps with hyperbolicity, the parameter space of the solutions is truncated.

IV. CONCLUSIONS
We have explored the influence of a coupling between a scalar field and the Ricci scalar on linear perturbations around scalarized black holes. This coupling, which is expected to be present in scalarization models based on EFT considerations, has already been shown to be crucial for observational viability: it can make GR a cosmological attractor, thereby providing the right conditions for compact objects without cosmological finetuning [28], and it can also suppress neutron star scalarization, thus removing binary pulsar constraints [26]. It also affects the amount of scalar charge scalarized black holes can carry [27]. We performed a radial stability analysis and have numerically shown that this same term renders scalarized BH solutions stable, confirming the expectations of Ref. [27].
This happens for values of the coupling constants that are within the same range considered in the aforementioned papers. Indeed, we have not found unstable modes for β larger than some critical value, β ≈ 1.2. Choosing a β above this threshold is consistent with having a cosmological attractor [28] and could quench neutron star scalarization for a range of α that still leads to black hole scalarization [26,27]. It is worth noting that negative values of β are also capable of stabilizing the solutions in a similar manner, but since they do not give rise to the cosmological attractor feature, they seem to be less interesting.
We have also performed a radial mode analysis in the Schwarzschild spacetime and looked for QNM modes. We were able to illustrate an interesting property: Beyond the fundamental mode one can find purely imaginary modes in a region of the parameter space where Schwarzschild BHs are stable.
Finally, we analyzed the effects of the scalar-Ricci coupling on the hyperbolicity of the scalar perturbation equation, using the exponential GB coupling as an example. We demonstrated that it actually improves the hyperbolic nature of the problem, by reducing the region of the parameter space where hyperbolicity breaks down. This happens as the additional term changes the scalar field profile, as determined by the full non-linear field equations, which in turn changes the coefficients of the linear perturbation equations. It would be interesting to generalize the hyperbolicity analysis to more general perturbations and beyond the linear level, in order to check if the coupling with the Ricci scalar could have a positive effect when considering hair formation by collapse [37] or binary mergers [38]. plement the continued fraction (CF) method to understand how the Schwarzschild spacetime becomes unstable. The CF method relies on providing a semi-analytical approximation of the wave function through the Frobenius method [39]. In summary, the solution can be written in terms of coefficients that are computed by a recursive relation that takes the form of a continued fraction, justifying the name. Since the method requires the analytical form of the coefficients (at least for the expansion to be implemented), we focus mostly on the stability analysis of Schwarzschild BHs.
There are different methods used to solve the above equations. Usually, in BH mode analysis, one reduces the nterm recurrence relation to a three-term one by successive Gaussian elimination steps, finding the appropriate CF expression to obtain the mode [39,40]. For example, in the six-term recurrence relation, we have to use three Gaussian elimination steps. We shall use the above to compute the modes for the Schwarzschild spacetime in sGB theories. We usually truncate the series from a higher value for n = N moving backwards up to n = 1, solving the series [40]. The number N of necessary terms depends on the quality factor of the modes, beginning higher for low values of lower quality factors. We note that for α = 0 we recover the standard recurrence relation for a free scalar field in the Schwarzschild spacetime.