Exact theory for the Rezzolla-Zhidenko metric and self-consistent calculation of quasinormal modes

As demonstrated in Ref. [22], the configuration space of the theory defined by (1) is so large that practically any given metric (spherically symmetric or otherwise) can be admitted as an exact solution for some choice of the functions $f,F,\mathcal{V}$, and $\chi$. In particular, given a metric, the dynamics of the scalar field can be constrained in such a way that that particular $g$ and $\phi$ pair form an exact solution to equations (3) and (4) provided that $f$ belongs to some particular class of functions.

Specifically, given some metric $g$, if there exists a scalar-field solution $\phi$ to the kinematic constraint equation

$$X_{0}=F(\phi)R+\mathcal{V}(\phi)-\chi(\phi)\nabla_{\alpha}\phi\nabla^{\alpha}\phi,$$

(5)

for some constant $X_{0}$, then, for any $f$ such that $f(X_{0})=f^{\prime}(X_{0})=0$, that particular $g$ is an exact solution to (3) and (4). In general, equation (5) must be solved numerically.

As a concrete example, if $\phi$ is chosen such that $X_{0}=0$, then $g$ is a solution to the field equations for

$$f(X)=X^{1+\sigma},$$

(6)

for any $\sigma>0$. Theories of the form (6) represent a potentially infinitesimal deviation from standard scalar-tensor gravity at the level of the action, and generalise the $f(R)=R^{1+\sigma}$ theories considered by Buchdahl [35] and others. These latter theories are of potential astrophysical relevance since they have been successful in explaining the flatness of galactic rotation curves for $\sigma\sim 10^{-6}$ without the invocation of dark matter [36, 37]. Such values for $\sigma$ are also consistent with cosmological constraints coming from primordial nucleosynthesis [38] (see also Ref. [39]). For $\sigma\lesssim 10^{-4}$, stellar structure is largely unchanged relative to GR [40], and thus the theory is capable of accommodating massive neutron stars, such as J0740$+$6620 (with mass $M\sim 2.17M_{\odot}$ [41]), for stiff equations of state. For weakly-varying scalar fields, these features are likely to persist since the theory (6) reduces to the Buchdahl one exactly for $\nabla_{\mu}\phi=0$ and $\mathcal{V}=0$. A curious feature of theories with $\sigma\ll 1$ is that the function $f$ is not analytic, which means that one cannot employ Taylor expansions about small scalar curvatures to investigate the Newtonian limit of the theory. Strong-field tests (coming from QNMs, for instance) are therefore especially germane to theories such as (6). Note, however, that analytic models can also be constructed; for example, the theory with $f(X)=X+\alpha X^{2}+\alpha^{2}X^{3}/4$ for any $\alpha\neq 0$ admits exact solutions with $X_{0}=-2\alpha^{-1}$.

In any case, we emphasise here that we are not necessarily advocating for this particular theory as a realistic description for gravitational phenomena, especially since higher-order theories such as (1) are generally susceptible to the Ostrogradsky instability (though see Ref. [42]). The techniques presented in this work are only meant to illustrate that self-consistent studies of QNMs for parameterised spacetimes are possible, and that their properties are sensitive to the particulars of the gravitational action, be it (1) or otherwise. In particular, the scalar field solution to equation (5) depends on the functional forms of the potential $(F,\mathcal{V})$ and kinetic $(\chi)$ terms, and therefore behaves differently in different branches of the general theory (1). The QNMs are thus sensitive to the particulars of the dynamics in the scalar sector (see Sec. III. B).

To give a concrete example of how the QNMs of a black hole can be self-consistently studied and compared between different theories, we focus on axial perturbations. Axial perturbations have the benefit that scalar degrees of freedom do not couple to the tensor sector [49], which further simplifies the analysis. The background scalar field still plays a role in the behavior of the perturbations however [49, 50] and modifies the form of the Regge-Wheeler potential [43]. Additionally, we note that even if the perturbed scalar field $\delta\phi$ does not couple to the metric, equation (4) predicts the existence of scalar GWs, the amplitudes of which may be non-negligible in some cases [51].

Given some astrophysical measurements of axial QNMs, it has been shown that it is possible to reconstruct the spacetime metric using statistical techniques (e.g., the Bayesian methods detailed in Ref. [20]) if one assumes some fixed set of EOM. However, metric and action parameters may be intertwined in the sense that the functional structure of the EOM may itself be dependent on the non-Schwarzschild parameters, which complicates this procedure. The full track of the inverse problem is therefore a non-linear one: one attempts to astrophysically reconstruct the metric by matching QNM measurements to eigenvalues that arise from a metric-dependent operator. As a first step towards a general solution of this difficult problem, we show how one can build a theory around a given spacetime and study its axial perturbations self consistently. This could be then be combined with statistical approaches, such as those detailed in [20], in a future study.

As mentioned earlier, we work with static and spherically³³3Theories of the form (1) generally allow for BHs with topological horizon structures also [52], though such objects are likely ruled out by electromagnetic observations of reflection spectra in accreting black holes [53]. symmetric spacetimes in this work for simplicity. The background line element, in Boyer-Lindquist coordinates $(t,r,\theta,\varphi)$, is taken to be

$$ds^{2}=-A(r)dt^{2}+B(r)dr^{2}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\varphi^{2}.$$

(7)

Making use of the Regge-Wheeler gauge [43], a generic axial perturbation can be written as

$$h^{\text{axial}}_{\mu\nu,\ell m}=\left(\begin{array}[]{cccc}0&0&h_{0}(t,r)S^{\ell m}_{\theta}(\theta,\varphi)&h_{0}(t,r)S^{\ell m}_{\varphi}(\theta,\varphi)\\ 0&0&h_{1}(t,r)S^{\ell m}_{\theta}(\theta,\varphi)&h_{1}(t,r)S^{\ell m}_{\varphi}(\theta,\varphi)\\ \ast&\ast&0&0\\ \ast&\ast&0&0\\ \end{array}\right),$$

(8)

where $S^{\ell m}_{\theta}(\theta,\varphi)=-\csc\theta\partial_{\varphi}Y_{\ell m}(\theta,\varphi)$ and $S^{\ell m}_{\varphi}(\theta,\varphi)=\sin\theta\partial_{\theta}Y_{\ell m}(\theta,\varphi)$ for spherical harmonic $Y_{\ell m}$, and the asterisk denotes a symmetric entry. The (in general complex) QNM frequency $\omega$ is introduced by taking a Fourier transform of the components $h_{0}$ and $h_{1}$ through

$$h_{0,1}(t,r)=\frac{1}{2\pi}\int^{\infty}_{-\infty}d\omega e^{-i\omega t}\tilde{h}_{0,1}(\omega,r).$$

(9)

The explicit dependence on $\omega$ and the tildes will henceforth be dropped for ease of presentation.

The methodology used to derive the EOM for $h_{0}$ and $h_{1}$ in the mixed scalar-$f(R)$ theory is remarkably similar to that of GR. In particular, the $\theta\varphi$-component of the $\mathcal{O}(\delta)$-field equations (3) allows for a direct expression of $h_{0}$ in terms of $h_{1}$, viz. [$f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})\neq 0$]

	$\displaystyle h_{0}(r)=$	$\displaystyle\frac{i}{2\omega B}\Bigg{[}h_{1}A^{\prime}+2Ah_{1}^{\prime}-\frac{Ah_{1}B^{\prime}}{B}$		(10)
		$\displaystyle+\frac{2Ah_{1}F^{\prime}(\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}})\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime}}{F(\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}})}+\frac{2Ah_{1}\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}f^{\prime\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})}{f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})}\Bigg{]},$

where we note that $\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}$ is a function of radius only for the spacetime described by (7). For the special case of $f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})=0$, a similar expression to that of (10) can be imposed, though without the final term on the second line. We note, however, that both expressions for $h_{0}(r)$ agree for $\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}(r)=0$ (relevant for the inverse problem discussed in Sec. II. A), i.e., when the constraint equation (5) is satisfied.

Moving forward, we introduce the tortoise coordinate $r^{\star}(r)\equiv\int dr\sqrt{\frac{B(r)}{A(r)}}$, and a new variable $Z$ through

$$h_{1}(r)=rZ(r)\sqrt{\frac{B(r)}{A(r)f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})F(\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}(r))}}$$

(11)

for $f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})\neq 0$ and $\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}(r)\neq 0$, and

$$h_{1}(r)=rZ(r)\sqrt{\frac{B(r)}{A(r)F(\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}(r))}}$$

(12)

otherwise, which allows us to rewrite the $r\varphi$-component of the field equations (3) as a Schrödinger-like equation. Algebraic manipulations eventually lead to the generalised Regge-Wheeler [43] equation in mixed scalar-$f(R)$ theory, which takes the familiar form

$$0=\frac{d^{2}}{d{r^{\star}}^{2}}Z+\left[\omega^{2}-V_{\ell}(r)\right]Z.$$

(13)

In the general case where $f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})\neq 0$, the potential $V_{\ell}$ has the lengthy but simple form

$$V_{\ell}(r)=V_{0,\ell}(r)-\frac{fA}{Ff^{\prime}}+\frac{3f^{\prime\prime}\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}}{Bf^{\prime}}\left[\frac{4A+rA^{\prime}}{4r}-\frac{AB^{\prime}}{4B}+\frac{3AF^{\prime}\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime}}{2F}+\frac{A\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime\prime}}{2\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}}+\frac{A\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}f^{\prime\prime}}{4f^{\prime}}\right]+\frac{3f^{\prime\prime\prime}A(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime})^{2}}{2Bf^{\prime}}.$$

(14)

The term $V_{0,\ell}(r)$, which reads

$$V_{0,\ell}(r)=\frac{\ell\left(\ell+1\right)A}{r^{2}}+\frac{3\left(AB^{\prime}-BA^{\prime}\right)}{2rB^{2}}+\frac{3F^{\prime}\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime}}{BF}\left[\frac{4A+rA^{\prime}}{4r}-\frac{AB^{\prime}}{4B}+\frac{AF^{\prime}\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime}}{4F}\right]+\frac{3A(\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime})^{2}F^{\prime\prime}}{2BF}+\frac{3AF^{\prime}\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime\prime}}{2BF},$$

(15)

is the piece that can be considered to ‘survive’ for $f(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})=f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})=0$ and $\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}}^{\prime}(r)=0$ [cf. equations (11) and (12)]. The latter term (15) in particular is the one most relevant for solutions of the inverse problem. In this case, given some background metric, the scalar field is chosen such that the constraint equation (5) is satisfied, and (15) differs from the GR form explicitly through the presence of $\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}$ and $F$ and implicitly through $A$ and $B$. It is worth pointing out that the only $\ell$-dependent part within $V_{\ell}(r)$ comes from the $\ell(\ell+1)A/r^{2}$ term, as expected of theories which predict that (helicity-2) gravitational waves travel at the speed of light; in the geometric-optics limit $\ell\rightarrow\infty$, the perturbation equation (13) reduces to the null geodesic equation (see Ref. [20] for a discussion).

Note that we have considered background spacetimes that are vacuum in this work. This is important to mention since, strictly speaking, any given spacetime can also arise as a solution to the Einstein equations exactly if one allows for arbitrary stress-energy tensors. Perturbations can be studied self-consistently in this approach also. However, in this latter instance, it is unlikely that the matter sector will be associated with a physical Lagrangian (e.g., arising from a fluid) or abide by well-motivated energy restrictions (e.g., dominant energy condition) for a generic background [54]. Such an approach is therefore not entirely physical in the study of disturbed BHs, since it would be difficult to explain the presence of exotic matter following a vacuum merger event, for example. Nevertheless, a self-consistent approach to the study of perturbations of arbitrary spacetimes within GR is formally possible and the EOM are simply given by $\delta R_{\mu\nu}=0$ if one assumes that the matter sector is unperturbed by the ringing.

This latter scheme, which we refer to as a ‘GR-like’ approximation throughout, could also be used to study the perturbations of a given spacetime (7). In the case of constant scalar field, the ‘GR-like’ Regge-Wheeler equation (see, e.g., equations (14)–(16) in Ref. [20]) is recovered exactly from (14) for either linear $f$ or for those $f$ with the property that $f^{\prime}(\overset{{\scriptscriptstyle 0}}{X}{\vphantom{\scriptstyle X}})=0$ [expression (15)]. In many cases of interest therefore, one might expect that this and similar schemes provide a fair representation for the QNMs, even if they are inexact (see the discussion in Ref. [20] and below). One of the aims of this present work is to provide a quantification for its accuracy, at least in the simple case of a truncated RZ metric (see Sec. IV).

Expressions (14) and (15) imply that QNMs are, in general, sensitive to the particular choices of $f$, $F$, $\mathcal{V}$, and $\chi$. As a demonstration, consider (6) for the Brans-Dicke choices with vanishing potential, $F(\phi)=\phi$ and $\chi(\phi)=\chi_{0}/\phi$, where the quantity $\chi_{0}$ is akin to the Brans-Dicke coupling constant. The $\chi_{0}\rightarrow\infty$ limit therefore corresponds to a more GR-like theory; when $f$ is linear, the theory reduces to Brans-Dicke and therefore to GR in this limit. In any case, the constraint equation (5) depends on the value of $\chi_{0}$, and therefore the scalar field does too when considering a solution to the inverse problem. Since the axial potentials $V_{\ell}$ depend on $\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}$ they also implicitly depend on $\chi_{0}$, and the QNMs shift depending on the branch of the general theory (1) under consideration.

This result is familiar from the study of Kerr or Schwarzschild BHs, which are known to ringdown differently in different theories of gravity [55, 47, 48], sometimes so dramatically that the same BH may be unstable in one theory and not another [56, 57] (see below). For the present case in the study of the inverse problem, this implies that ringdown analyses place constraints on the metric behaviour and on the theory simultaneously; the QNMs for any given quantum numbers depend on $\chi_{0}$ and the free parameters within the metric in a non-trivial way. This adds a layer of complication from a data analysis perspective, since there may be a degeneracy between the physical BH parameters and those parameters that appear within the action functional (see, e.g., Ref. [58]).

The dependence of the QNMs on $\chi_{0}$ implies that stability can be theory-dependent. Specifically, the appearance of a deep enough turning point for $V_{\ell}<0$ (in a ‘negative gap’) implies the existence of bound states (i.e., $\omega^{2}<0$ modes) that lead to (linear) instability [59, 60] (see also Ref. [61]). On the other hand, if the potential $V_{\ell}$ is positive-definite, then it is well known that the Schrödinger equation (13) admits no bound states. This condition is not necessary however: using the so-called S-deformation technique [62], Kimura [63] has presented evidence that if there exists a continuous and bounded function $S$ such that

$$0=V_{\ell}(r)+S^{\prime}(r)-S(r)^{2},$$

(16)

then stability is also assured. These techniques are not necessary for our analysis, since we consider cases with sufficiently small deformation parameters and large enough $\chi_{0}$ (i.e., more GR-like values) so that the negative gap is never too deep to allow for bound states. In general though, a given BH may be stable for $\chi_{0}$ values above some threshold value $\chi_{c}$ that depends on the deformation parameters, and unstable for theories with $\chi_{0}<\chi_{c}$. A thorough investigation of the (in)stability of BHs in different branches of the general theory (1) will be conducted elsewhere, since the WKB method used here is not well-suited for identifying bound states.

In general, the Schrödinger-like equation (13) is subject to some set of boundary conditions. These are chosen such that we have pure outgoing radiation at infinity and purely ingoing radiation at the horizon, i.e., that as $r^{\star}\rightarrow\pm\infty$, we have [68]

$$Z(\omega,r)\sim e^{\pm i\omega r^{\star}}.$$

(21)

Equation (13), subject to the boundary conditions (21), constitutes an eigenvalue problem for the QNMs. However, since we are mostly interested in the frequencies of the QNMs and not so much the functional form of $Z$ itself, it is not necessary to solve equation (13) formally (see Refs. [10, 11, 12] for some classical reviews on the topic).

Although there are a number of semi-analytic and numerical methods available, we operate with the WKB method [69, 70, 71] because it is especially suited for problems in which other standard approaches may require significant adjustments or extensions [72]. Other techniques involving continued fractions (e.g., Leaver’s method [73]), phase-integrals (e.g., [74]), or time-dependent integration (e.g., [75]) typically require an analytic understanding of the asymptotic properties of the potential, which may be absent when one only has numerical potentials defined up to some finite radius at their disposal, as in our case. In the order-$j$ WKB method, the QNM frequency is determined through

$\displaystyle\frac{iQ_{0}}{\sqrt{2Q^{\prime\prime}_{0}}}-\sum_{j=2}\Lambda_{j}(n)=n+\frac{1}{2},$

(22)

for overtone number⁴⁴4Throughout this work we only work with the least-damped (and therefore most relevant for astrophysical observation) modes for brevity, which have $n=0$ by definition. However, overtones are generally expected to be active for a newborn BH and may be important when reconstructing BH parameters from astrophysical data; see Ref. [76] for a discussion. $n$, where $Q(r^{\star})\equiv\omega_{n}^{2}-V_{\ell}(r^{\star})$ is evaluated at the maximum of the potential and primes denote derivatives with respect to the tortoise coordinate $r^{\star}$. The $\Lambda_{j}$ in equation (22) depend on $2j^{\text{th}}$-order derivatives of $V_{\ell}$ and have lengthy forms that we will not write out here; see, e.g., Ref. [71].

While the WKB method has been extended up to 13th order using Padé approximants [thereby requiring $26^{\text{th}}$-order derivatives of $V_{\ell}(r)$] [77], numerical stability and implementation considerations limit us to the $4^{\text{th}}$-order scheme. The relative errors of our results have been checked using both $3^{\text{rd}}$ and $4^{\text{th}}$ order approximations and are, at worst, at the percent level for a few cases (imaginary parts). For most cases however the errors are significantly smaller. An additional numerical check has been performed using the Pöschl-Teller approximation [78, 79] to verify the robustness of the WKB method. Note that for a few parameter combinations the $\ell=2$ potentials develop a small negative gap left of the maximum (see Figure 2 below), which is likely responsible for the percent-level errors mentioned previously. In any case, the various numerical checks performed in this work make us confident that the numerical routines for finding the potential [checking that the right-hand side of (19) vanishes to machine precision] and computing the QNMs (checking agreement between $3^{\text{rd}}$ and $4^{\text{th}}$ orders and the Pöschl-Teller approximation) give reliable results.

We graph the $\ell=2$ potential profiles $V_{0,2}(r)$ (solid curves) for two illustrative cases with $a_{1}=b_{1}=0$ and $M=\chi_{0}=1$ but $\epsilon=-0.1$ (red) and $\epsilon=-0.3$ (blue) in Fig. 2. For contrast, the corresponding potentials computed using the GR-like scheme $[\overset{\scriptscriptstyle{0}}{\phi}{\vphantom{\phi}}^{\prime}(r)=0]$ are represented by dashed curves. The respective domains of the potential functions vary with $\epsilon$, because the horizon location $r_{\mathcal{H}}$ (shown by vertical dashed lines) is sensitive to this quantity. In the case of low $|\epsilon|$ (i.e., more Schwarzschild-like), we see that the GR-RZ and exact potential functions overlap almost exactly: the greatest difference between the two curves is $\sim 1.5\%$ near the respective peaks (shown by solid diamonds). However, since the properties of the QNMs are directly tied to the location and value of the peak [cf. equation (22)], even a small difference can lead to a non-trivial disparity in the real and imaginary components (see Sec. IV. B). For the $\epsilon=-0.3$ case, the difference between the exact and approximate potentials is more noticeable, and the relative difference reaches $\sim 7\%$ near the respective peaks. Note that a negative gap develops near the horizon in the exact potential for the large $|\epsilon|$ case, though is not deep enough to produce bound states for the chosen kinetic value $\chi_{0}=1$. For $\chi_{0}\ll 1$ or $\epsilon\ll-0.3$, bound states may exist in principle. More generally, the disagreement between the GR and exact schemes scales inversely with $\chi_{0}$ and directly with $|\epsilon|$ (see Figures 3 and 4 below, respectively).

In this section we present computations of QNMs using the numerical methods detailed in the previous section for a variety of representative cases.

Fig. 3 shows the $\ell=2$ (blue circles) and $\ell=3$ (orange squares) QNM frequencies as functions of $\chi_{0}$ for the RZ parameters relevant to the scalar field solutions shown in Fig. 1 (i.e., $\epsilon=0.1,a_{1}=0.15$, and $b_{1}=-0.4$). For small values of $\chi_{0}$ in the range $0.1\leq\chi_{0}\lesssim 1.5$, we see a substantial and monotonic variation in both the real ($\approx 33\%)$ and imaginary ($\approx 15\%$) components of the $\ell=2$ QNMs. Variations on a similar scale are likewise observed in the $\ell=3$ case. Turning points occur in the graphs once a critical value of $\chi_{0}\approx 1.5$ is reached, and the real (imaginary) frequencies begin to increase (decrease) rather than decrease (increase). In the large $\chi_{0}$ limit, both sets of frequencies asymptote towards particular GR-like values (see below). In any case, we find that a continuous (with respect to $\chi_{0}$) spectrum of frequencies within some bounded range can be achieved for a fixed set of RZ parameters. This result has the interesting implication that, given only a single measurement of the fundamental frequency of a newborn object (as in the case of GW150914 [2]), many RZ parameter sets can be accommodated within some⁵⁵5Note, however, it may be the case that that particular branch is Ostrogradsky-unstable [42] or does not respect independent constraints coming from cosmology [36, 37], neutron-star astrophysics [40], or Solar-system dynamics [67]. branch of the general theory (1). This stresses the necessity of having multiple QNM measurements in placing constraints on BH behaviour and strong gravity simultaneously (see also Refs. [64, 14, 15]).

Overplotted in Fig. 3 are the GR-like QNMs (black dots), which of course do not vary as a function of $\chi_{0}$. For small of values of $\chi_{0}$, the two schemes predict distinct behaviour, as noted previously. In the limit $\chi_{0}\rightarrow\infty$ however, the scalar dynamics are heavily suppressed (cf. Fig. 1) and the curves match as the theory approaches the Buchdahl one [35], much in the same way that the Brans-Dicke theory (i.e., linear $f$) approaches GR in this limit. In practice, values $\chi_{0}\gtrsim 10^{4}$ lead to numerical indistinguishability between the real and imaginary components for the GR-like and exact schemes for this particular set of RZ parameters. Note though that this does not mean that the frequencies approach the Schwarzschild values: the $\ell=2$ fundamental mode in the Schwarzschild case (not shown) has real value $\text{Re}(\omega_{\text{GR}})=0.374$ [75], for instance, and is still $\sim 8\%$ different from the RZ value in this particular case even when $\chi_{0}\rightarrow\infty$. Such a disparity does not, however, exceed the limits imposed by observations of GW150914, where constraints on the fundamental frequency were placed at roughly the $\sim 10\%$ level relative to the GR (Kerr) values at $90\%$ confidence [2].

As another example, we consider the case of $a_{1}=b_{1}=0$ with fixed $\chi_{0}=1$ but varying $\epsilon$. This case is illustrated in Fig. 4 in a similar style to Fig. 3. In this instance, the GR-like scheme is approached as $\epsilon\rightarrow 0$, where in fact the Schwarzschild QNMs are recovered exactly [75] in both schemes, as expected. Overall, we see that the GR-like approximation is a robust one: disagreements, relative to the exact case, in the real and imaginary parts of the QNM frequencies are at most $\sim 3\%$ for $|\epsilon|\lesssim 0.3$ and $\ell=2$. For larger values of $\chi_{0}$ and $\ell$, the disagreements fall even further. This adds strength to the claims made in Ref. [20], who made use of the GR-like approximation within a Bayesian scheme to show how QNM data can be used to reconstruct a spacetime metric to a high degree of accuracy. In either scheme, however, we see that between the $\epsilon=0$ (Schwarzschild value) and $\epsilon=-0.3$ cases, the real parts of the QNM frequencies change by $\approx 20\%$, which exceeds the limits imposed by GW150914 [2]. Demanding that the frequencies match to the Schwarzschild values within $\lesssim 10\%$ leads to the constraint $|\epsilon|\lesssim 0.16$ for $\chi_{0}=1$, though we note that such a direct comparison with GW data is imprecise at this stage because we do not model rotation.

Although the parameter space of RZ-QNMs is very large (different quantum numbers, RZ parameters, and theory variables) and an exhaustive study is impractical, it is instructive to consider a few additional cases. The $\ell\leq 5$ QNM frequencies for a few cases with scattered values of $\epsilon$, $a_{1}$, and $b_{1}$ are shown in Table 1 (2) for $\chi_{0}=1$ ($\chi_{0}=0.5$). For example, the second and third columns list the frequencies for cases where all RZ parameters vanish except for $b_{1}$. The real parts of the frequencies are largely similar to the Schwarzschild values in these cases for $b_{1}=-0.38$, where we find $\text{Re}(\omega)=0.38$ for $\chi_{0}=1$ while $\text{Re}(\omega)=0.368$ for $\chi_{0}=0.5$, which are marginally larger and smaller than the Schwarzschild value $\text{Re}(\omega_{\text{GR}})=0.374$, respectively. The imaginary components are substantially greater in either case ($\gtrsim 40\%$ difference) though, implying that the modes would be damped out faster [since the damping time $\propto-\text{Im}(\omega)^{-1}$] than the corresponding GR case. For more negative values of $b_{1}$, however, even the real values diverge significantly from the Schwarzschild values, and are thus likely ruled-out from GW observations [3]. The damping times for the cases shown in the final two columns, which have $\epsilon=0.1$ and $b_{1}=-0.19$ but $a_{1}=0.15$ and $a_{1}=0.3$, respectively, are very similar to each other, and only weakly depend on $\ell$. This implies that even high-$\ell$ modes may be important in the early characterisation of BH ringdown for some non-Schwarzschild metrics, since they may not necessarily be damped out faster than their low-$\ell$ counterparts.