Constraining modifications of black hole perturbation potentials near the light ring with quasinormal modes

In modified theories of gravity, the potentials appearing in the Schr\"odinger-like equations that describe perturbations of non-rotating black holes are also modified. In this paper we ask: can these modifications be constrained with high-precision gravitational-wave measurements of the black hole's quasinormal mode frequencies? We expand the modifications in a small perturbative parameter regulating the deviation from the general-relativistic potential, and in powers of $M/r$. We compute the quasinormal modes of the modified potential up to quadratic order in the perturbative parameter. Then we use Markov-chain-Monte-Carlo (MCMC) methods to recover the coefficients in the $M/r$ expansion in an ``optimistic'' scenario where we vary them one at a time, and in a ``pessimistic'' scenario where we vary them all simultaneously. In both cases, we find that the bounds on the individual parameters are not robust. Because quasinormal mode frequencies are related to the behavior of the perturbation potential near the light ring, we propose a different strategy. Inspired by Wentzel-Kramers-Brillouin (WKB) theory, we demonstrate that the value of the potential and of its second derivative at the light ring can be robustly constrained. These constraints allow for a more direct comparison between tests based on black hole spectroscopy and observations of black hole `shadows'' by the Event Horizon Telescope and future instruments.


I. INTRODUCTION
Since the first groundbreaking direct detection of gravitational waves in 2015 [1], the LIGO-Virgo-KAGRA collaboration has reported around one hundred additional events [2][3][4][5][6]. Most of the detected signals are in agreement with the predictions of general relativity (GR) for the merger of two black holes, while some of them involve neutron stars. Besides the important implications for astrophysical formation scenarios of black hole populations, these events allow for unprecedented tests of GR in the strong field [7][8][9][10][11][12][13][14][15][16][17][18]. A key prediction from GR is that rotating astrophysical black holes should be described by the Kerr metric, and that their perturbative dynamics at late times (in the so-called "ringdown" regime) should be well described by a superposition of damped sinusoids called quasinormal modes (QNMs), with characteristic frequencies determined only by the black hole's mass and spin. If GR is the correct theory of gravity, the observed QNM frequencies should match those predicted in GR for Kerr black holes [19][20][21].
The LIGO-Virgo-KAGRA collaboration has analyzed the current catalogs, at first focusing on GW150914 and then using multiple events, with the specific aim to extract QNMs from the ringdown [14][15][16][17]. The observations are generally consistent with the presence of the * sebastian.voelkel@aei.mpg.de quadrupole ( = m = 2) fundamental mode (with "overtone number" n = 0). Reference [22] claimed evidence for the presence of overtones (higher-damping = m = 2 modes with n > 0) in GW150914, and Ref. [23] claimed evidence for the fundamental mode with = m = 3 in GW190521. Later work showed that these conclusions depend on the detector noise, data analysis methods, the choice of starting time, and nonlinear effects in the theoretical modeling of ringdown [24][25][26][27][28][29][30]. Despite the ongoing debate, the analysis of numerical relativity simulations implies that at least some overtones can in principle be extracted when the signal-to-noise ratio is large enough [31][32][33][34][35][36][37][38][39][40]. Theoretical modeling and data analysis challenges may be more subtle than anticipated, but next-generation detectors are expected to provide reliable, high-precision measurements of more than one QNM [41][42][43][44]. Therefore it is important to investigate how these measurements will inform us on possible deviations from GR.
There have been various attempts to introduce deviations from GR in the QNM spectrum in a theory-agnostic manner. These include modifications of the gravitational action around a GR black hole background [45], the addition of perturbative corrections at the level of the perturbation equations [46][47][48] or at the level of the metric [49,50] (which then require suitable assumptions for the dynamics), or the addition of free parameters in the QNM frequencies themselves [51][52][53]. While the physical interpretation of these "theory-agnostic" constraints requires a specific modified theory of gravity, it is desir-arXiv:2209.10564v2 [gr-qc] 12 May 2023 able to have a robust phenomenological framework encompassing several theories and allowing us to solve the "inverse problem" -i.e., to infer what specific GR modification caused a deviation from GR.
For non-rotating spacetimes, the inverse problem based on parametrized black hole metrics for axial (odd parity) perturbations was studied in Ref. [49]. For spinning black holes, a parametrized spectroscopy framework ("ParSpec") was introduced in Ref. [52], and applied to data in Ref. [53]. The ParSpec framework is based on a Taylor-series expansion of the QNM frequencies in the dimensionless Kerr spin parameter. In principle this can be used to stack multiple events [54], and then compare with the predictions from specific theories.
In this work we assume that deviations from the QNM spectrum in GR can be adequately captured by small modifications of the underlying perturbation equations. As a proof of principle, we focus on non-rotating or slowly rotating black holes. We adopt the parametrized formalism of Refs. [46,47], which systematically connects small deviations in the perturbation equations with the QNM spectrum. At lowest order, the idea is to consider small modifications in the perturbation potential at the linear level and to write them as a "post-Newtonian" (PN) series expansion in M/r [46]. Going to quadratic order allows one to capture (more realistically) possible couplings to additional fields [47]. While in Refs. [46,47] the focus is on the fundamental mode, Ref. [48] extended the analysis to overtones, and studied the inverse problem using a principal component analysis.
We apply a Bayesian analysis to solve the inverse problem: given a simulated set of QNM frequencies computed within the parametrized formalism, can we infer the deviation parameters in the underlying potentials? We are particularly interested in the understanding whether it is possible to constrain the more general and realistic case where many deviation parameters are varied simultaneously. Therefore we compute constraints on the individual PN-like expansion parameters twice: we first vary the parameters one at a time ("optimistic" case), and then we vary them all simultaneously ("pessimistic" case). In both cases we find that the constraints on the individual expansion parameters are not robust. The priors play an important role, especially when all parameters are varied simultaneously. In particular in the "optimistic" case the recovered bounds on the individual parameters can be biased.
However, the complex correlations in the posterior distributions in the "pessimistic" case contain valuable information. It has long been known that QNM frequencies are related to the behavior of the perturbation potential near the light ring [55][56][57]. Therefore we propose a different strategy: we map the PN series expansion to the value of the potential and its derivatives at the peak by using insights from Wentzel-Kramers-Brillouin (WKB) theory. Higher-order WKB methods relate the QNM frequencies with a Taylor expansion of the effective perturbation potential in the vicinity of its maximum, i.e., close to the light-ring [56,[58][59][60]. In fact, the (closely related) eikonal approximation was used in Refs. [61][62][63][64] to build a "post-Kerr" strategy to parametrized black hole spectroscopy. By connecting the Taylor (PN-like) expansion to the light-ring WKB expansion, we demonstrate that the value of the potential and of its second derivative at the light ring can be robustly constrained using Bayesian techniques.
Our work suggests that it may be possible to relate black hole spectroscopy tests to electromagnetic observations of black hole "shadows" by the Event Horizon Telescope and future instruments [65][66][67][68]. However, the interpretation of parametrized QNM tests of strong-field gravity will require very precise measurements and a more solid theoretical understanding of QNM excitation.
The paper is structured as follows. In Sec. II we introduce our theoretical framework and Bayesian data analysis techniques. In Sec. III we apply the method to GR and non-GR injections. In Sec. IV we discuss our results and outline possible directions for future work. Throughout the paper we adopt geometrical units (G = c = 1).

II. METHODOLOGY
In this section we give an overview of the three main ingredients of our method: the parametrized QNM framework based on a Taylor expansion of the perturbation to the GR potential (Sec. II A), the WKB expansion of the modified potential near the light right (Sec. II B), and the Bayesian inference approach used to address the inverse problem (Sec. II C).

A. Parametrized QNM framework
The perturbations of the Schwarzschild black hole were first studied in the odd-parity (axial) case by Regge and Wheeler [69], and later extended to the even-parity (polar) case by Zerilli [70]. The more general case of Kerr black holes was worked out by Teukolsky [71].
The characteristic frequencies and damping times of the QNMs for rotating and non-rotating black holes in GR have long been known (see [72][73][74][75][76] for reviews). More recently, various authors have considered black hole perturbations and QNMs in modified theories of gravity. Almost all works are limited to non-rotating or slowly rotating black holes [77][78][79][80][81][82][83], because finding exact analytical background solutions for arbitrary rotation is not always possible, and separating the perturbation equations is even harder (if at all possible). Theory-specific works using the slow-rotation approximation for gravitational perturbations include dynamical Chern-Simons gravity [84,85], Einstein-dilaton-Gauss-Bonnet gravity [86,87], and effective-field-theory extensions of GR including higher-derivative terms [88,89], although there are recent attempts at formulating generalized Teukolsky equations valid for more generic theories and arbitrary rotation [90,91].
Since Einstein's theory is very well tested in the weakand strong-field regimes, from an experimental point of view it is reasonable to treat possible modifications as small deviations from GR. Under this assumption, the parametrized framework developed in Refs. [46][47][48] allows for a convenient and efficient calculation of QNMs once the perturbation equations are cast in the form where the radial function Φ comes from a spherical harmonic decomposition of the perturbed metric or of some other perturbing field, the tortoise coordinate r * is defined in terms of the areal radius r and the function f = 1 − r H /r through dr * = dr/f , where r H is the location of the black horizon, and ω is the complex QNM frequency. For metric perturbations we can write where V GR is either the Regge-Wheeler or Zerilli potential, and If the coefficients α (k) are small, the QNMs obtained by the solution of Eq. (1) can be approximated, up to second order, by the expression where ω 0 are the GR frequencies, and the coefficients d (k) and e (ks) were first introduced in [46,47]. We compute these coefficients via a continued-fraction method, as in Ref. [48]. Note that in general the α (k) can be complex numbers and might depend on the unperturbed QNM frequency, and thus also on the overtone number itself. We are also neglecting a quadratic correction term proportional to the possible QNM dependence of the potential correction term, which can be found in Ref. [47]. By a redefinition of the field it is possible to reduce the number of parameters in the potential [92], but this was only demonstrated in the linear case, and it might not be possible in the quadratic case studied in this work (or in the more general and realistic case of rotating black holes). The α (k) coefficients must be small enough for the QNM frequencies to be adequately approximated by a quadratic expansion in α (k) . Therefore we impose the following -necessary but not sufficient -"convergence criterion" (see [46]): B. Higher-order WKB method The QNMs in this work are obtained either numerically (with a continued fraction method) or analytically (through the parametrized QNM framework). However, it will be useful to analyze our results by making use of the physical insight derived from a higher-order WKB approximation. The WKB approximation is widely used for scattering problems of the form of Eq. (1), if the potential describes a barrier with a single maximum and suitable asymptotics. Within black hole perturbation theory, the leading-order WKB approximation is known as the Schutz-Will formula [56], and it was later generalized to higher orders [58][59][60]. The higher-order WKB approximation can be written, schematically, as where n is the overtone number. The correction terms in the WKB approximationΛ i (n) are lengthy expressions [58][59][60] (for completeness, we list them in Appendix A up to third order). They involve derivatives V (p) of the effective potential with respect to the tortoise coordinate, evaluated at the maximum of the potential: Increasing the order of the WKB correction introduces even higher-order derivatives of the potential. In general the WKB approximation works rather well for the fundamental mode, and it becomes less accurate as n increases [93].
In the following we revise the calculation of the potential derivatives assuming a perturbative expansion in the coefficients α (k) . The maximum of V with respect to the tortoise coordinate is found by solving dV /dr * = 0, and it is located at r MAX =r + k α (k) δr k , wherer is the location of the peak in GR and At first order in α (k) , denoting with a p superscript the p-th derivative with respect to r * evaluated at r MAX , we have In Fig. 1 we plot δV (p) k for the = 2 axial potential (up to p = 4) as a function of k for k ≤ 40, while in Table I we list the value of these derivatives for k ≤ 10.
We can solve for ω and expand the result at first order in the α (k) . This yields the following relation for the k δV frequency at third order in the WKB approximation: where ω WKB is the frequency evaluated with the WKB method within GR, and δΛ i are the linear expansions in α of the coefficientsΛ i . By comparing Eqs. (4) and (10) order by order, we can map the coefficients d (k) to the derivatives δV (p) k of the displaced potential at the peak.

C. Bayesian approach
Bayesian analysis allows us to relate the observed data D with the underlying parameters θ of a model, or to compare different models with each other. It also allows us to quantitatively include our prior knowledge or assumptions into the analysis, and understand how they affect the posterior through Bayes' theorem which states that the posterior P θ|D is equal to the prior P (θ) times the likelihood P D|θ , divided by the evidence P (D). The evidence is often unknown or hard to compute, but it is still possible to compute the posterior via Markov-chain-Monte-Carlo (MCMC) techniques. These only require the knowledge of the prior and likelihood, and can be used to directly draw samples from the posterior. Standard MCMC techniques become very computationally expensive when the number of sampled parameters is large. The parametrized QNM framework of Eq. (4) is very beneficial in this sense, because it speeds up the calculation of the likelihood significantly by avoiding the more involved (and not necessarily always converging) calculations that are necessary in other techniques. This results in quick MCMC sampling, allowing us to study a reasonable number of parameters. We assume that our likelihood for the unknown parameters θ = [r H , α (k) ] is given by where C −1 is the inverse of the covariance matrix C of the QNM frequencies, and where r 0 H is the location of the horizon as it would be inferred in GR. We perform this rescaling since we do not know a priori the correct value for r H .
The exact form of the correlations depends on the details of the observed binary black hole system, as well as the details of the detector. However, by identifying the inverse of the covariance matrix with the Fisher matrix, it is possible to use analytic estimates derived in Ref. [21] as an approximation. For measurements of the fundamental QNM and of the first overtone (n = 0, 1) the Fisher matrix we need to estimate is a 4 × 4 matrix. The submatrix that connects the real and imaginary parts of the frequency for a fixed value of n can be approximated with the results of Ref. [21], but the correlations between different n's (i.e., the 2 × 2 off-diagonal blocks) cannot, and therefore we neglect them. However, we have analyzed how moderate random correlations between the fundamental mode and the first overtone would change the results, and we did not find substantial qualitative changes. We rescale the Fisher matrix so that the real part of the n = 0 mode is measured with a relative error of 1 %, which yields a relative error of 4.7 % for its imaginary part. Assuming that the n = 1 overtone is excited with a similar amplitude (an assumption justified by fits to numerical relativity simulations [33]) yields relative errors of 3.4 % and 8.2 % for the real and imaginary part of the overtone frequency, respectively. Ongoing work in the literature suggests that the specific model used to extract overtones can have an important impact on interpreting the constraints [29,30], but this approximate estimate of the Fisher matrix is sufficient for our purposes.
The location of the black hole horizon r H affects the QNM spectrum and depends on the underlying theory. We assume that r H is close to its GR value r 0 H and that the uncertainty on r H is relatively small (σ rH = 5 %), for example because we have a good estimate of the remnant black hole's mass from the inspiral/merger waveform. Then we multiply the likelihood of Eq. (12) by an additional factor This is equivalent to assuming a Gaussian prior for r H centered around its value r 0 H in GR. The priors P(α (k) ) cannot be chosen to be arbitrary because we are using a parametrized QNM framework. For the perturbative formalism to be valid, we must adopt bounds consistent with Eq. (5). For concreteness, in our analysis we consider two different prior realizations: α 10]. Our choice to study 11 parameters for α (k) guarantees that we can explore a large parameter space to capture modifications to GR, but at the same time we can still perform efficient MCMC sampling. To perform the MCMC analysis we use the emcee sampler [94], which is based on the affine-invariant ensemble sampler proposed by Goodman and Weare [95]. Since one would not generically expect that very large k contribute in a dominant way, neglecting higher orders is a reasonable choice. Throughout this work we assume flat (uniform) priors within these two different bounds. Considering different ranges allows us to quantify what aspects of the analysis are sensitive to the choice of priors and which are not. By using the continued-fraction numerical code discussed in Ref. [48], we have checked that the frequencies generated by taking α (k) = ±α (k) P20 , with k ∈ [0, 10], are well approximated by the quadratic expansion. However, for certain combinations of random draws from the larger prior range α (k) = ±α (k) P10 the approximation can become inaccurate. Our results below (based on the exact injection) show that this only mildly affects our reconstruction of the perturbative parameters that were used in the parametrized framework for the inference.

III. RESULTS
Using the techniques of the previous section, we now study two complementary settings. In the first setting (Sec. III A) we perform injections assuming that GR is the correct theory of gravity. In the second setting (Sec. III B) we assume a hypothetical modification with two non-zero deviation parameters -a nontrivial, but tractable case.
Within each setting we use the (complex) l = 2, n = 0, 1 QNM frequencies as hypothetical observations, with errors estimated by the Fisher matrix formalism as described in Sec. II C. To explore the impact of different priors on the posteriors, we always use two different prior ranges (α (k) P20 and α (k) P10 ). Within each of these two settings, we will address the inverse problem in two steps. We will first compare the posteriors for the parameters α (k) in the optimistic scenario (where only one parameter varies) against the pessimistic scenario (where all parameters are varied simultaneously). Then, from the reconstruction of the potential, we will evaluate V (p) and compare it with the corresponding value V (p) GR in GR by computing the relative difference

A. Results for a GR injection
First we present our results under the assumption that GR is the correct theory of gravity.
The violin plots in Fig. 2 show the optimistic (left side) vs pessimistic (right side) posteriors for the Taylor expansion coefficients α (k) in Eq. (2). Different colors correspond to posteriors obtained with the two different (flat) prior ranges. In the optimistic case, it is possible to constrain all parameters independently of the chosen prior range. While the quantitative details of each optimistic bound are different and their widths increase with k (as expected), they are qualitatively the same for all k. The situation is clearly different in the pessimistic case.
Here the posterior distributions of all α (k) have support at the prior boundaries, and they become broader when the prior range is increased. Since we allow for more parameters (12 in total: 11 α (k) 's and r H ) than observed QNMs (4, i.e., 2 complex frequencies), this is not surprising.
In Fig. 3 we show the results for the WKB deviation coefficients δ rel V (n) defined in Eq. (15). Since flat priors for α (k) do not correspond to flat priors for δ rel V (n) , it is important to compare the posteriors with the priors. We want to understand if QNM frequencies contain information on the δ rel V (n) coefficients, or if the inference is prior-dominated (in which case the QNMs would not be informative). We focus on the pessimistic posteriors, since the more general (pessimistic) assumption allows us to draw conservative conclusions. As usual, we show the two different prior ranges using different colors.
Quite remarkably, the posteriors for δ rel V (0) are very robust, and they do not depend on the prior range choice.  2. MCMC results for simulated observations of the Schwarzschild n = 0, 1 QNMs with relative errors described in the main text, and assuming a Gaussian prior for rH corresponding to 5 % at 1σ confidence. In each panel, on the left side we plot the posteriors in the optimistic case (only one α (k) is varied at a time), and on the right side we plot the posteriors in the pessimistic case (all α (k) 's are varied simultaneously). The blue (orange) color corresponds to a prior range of α Fig. 2, but here we sample the relative errors of various derivatives of the effective potential δ rel V (n) with respect to the tortoise coordinate evaluated at the maximum. The right side of each violin plot shows the sampling from the prior distribution of all α (k) simultaneously, with colors label the different α (k) prior ranges: α The posteriors for δ rel V (2) and δ rel V (3) are still more informative than the priors, but the QNM frequencies do not provide as much information as they do in the case of δ rel V (0) , and the measurement is even less informative in the case of δ rel V (4) . This is one of the main results of this paper: the "change of basis" from the α (k) 's to the δ rel V (n) 's (variables related to the light ring) is very effective, because QNMs carry physical information about the potential and its derivatives at the peak. This is in agreement with the conclusions of Ref. [48].

FIG. 3. Same MCMC results as in
Note that the derivation of these results is fully independent of the WKB approximation: we never use the WKB approximation to compute the QNMs, but rather we use Leaver's method to compute the injected QNM frequencies, and the parametrized framework of Eq. (4) to model their deviations from GR. The full MCMC analysis can capture correlations between the deviation parameters, as long as they are small enough to ensure the validity of the parametrized framework.

B. Results for a non-GR injection
As an example of a possible non-GR injection, we consider the hypothetical case in which α (2) = α (3) = 0.2, while all other α (k) 's are set to zero. We compute the QNM frequency with Leaver's method, while we use the parametrized framework for the MCMC analysis.
In Fig. 4 we report the MCMC results, following the same notation and conventions as in Fig. 2. The injected values of α (k) are marked by horizontal red lines. The plot shows that it is not possible to infer the correct values for any α (k) , except perhaps for k = 2 and (less clearly) for k = 3. The recovered values for all other α (k) 's peak away from their injected value. Some of the small-k posteriors have strong support at the smaller prior limit, but they are well captured by the larger prior limit.
Given that these results are qualitatively very different from the GR injections, it is quite remarkable that the WKB-motivated constraints on the potential near the peak are as robust as before. In Fig. 5 we show that it is indeed possible to recover the non-GR injections (shown as red horizontal lines). The posterior of the dominant term δ rel V (0) has tight support away from GR, while δ rel V (2) is broader and has some overlap with the GR hypothesis. The higher derivatives for this injection are very close to their GR value, and the corresponding violin plots are almost indistinguishable from those of Fig. 3.   Fig. 4, but here we sample the relative errors of various derivatives of the effective potential δ rel V (n) with respect to the tortoise coordinate evaluated at the maximum. The right side of each violin plot shows the sampling from the prior distribution of all α (k) simultaneously. The colors label the different α (k) prior ranges, α C. Theory-specific approach: a simple example A theory-specific analysis can, in principle, be carried out in different ways. We could avoid making use of the parametrized framework and base it on a full MCMC analysis involving all free parameters of the theory. We could also, in principle, use the parametrized framework with the theory-predicted values of α (k) , which could depend on one (or multiple) coupling constants of the modified theory of gravity of interest. Instead of repeating a full MCMC analysis, here we show that the posteriors for δ rel V (0) in the non-GR example above are already very informative when used in a "post-processing" analysis.
Let us assume, for simplicity, that the hypothetical modified theory predicts two non-zero deviation parameters, related to the only unknown parameter of the theory ζ as follows: We use this trivial example just for illustrative purposes, but no fundamental limitation prevents us from relaxing these assumptions to allow for more than two non-zero deviation parameters, or to allow for a more complicated dependence of these parameters on ζ. In fact, we would find qualitatively similar results if we assumed that α (3) is suppressed by a small factor . By using Eq. (9) we can find the linear approximation for the modified potential V (0) (ζ). The GR potential V (0) is known, so we can express the posteriors δ rel V (0) in Fig. 5 (left side of the first violin plot) in terms of V (0) . The posteriors are well approximated by a Gaussian, whose mean µ V and width σ V we can fit numerically. Then we can insert the value of µ V fitted to V (0) into Eq. (9), and invert it to find the most likely value of ζ, which we will denote by µ ζ . We can also estimate the 1σ and 2σ errors by repeating the inversion for µ V ± σ V and µ V ± 2σ V . In Fig. 6 we show the results of this procedure. The injected value of the coupling parameter ζ is well within the 1σ constraint. The values of µ ζ and σ ζ can also be used to define the theory-specific bounds on α (2) (ζ) and α (3) (ζ) in a similar way.
If we do not know the precise form of α (2) (ζ), α (3) (ζ) (i.e., in a theory-agnostic approach), we could treat them as independent variables in a post-processing analysis. From Eq. (9) we can find α (3) as a function of α (2) and V (0) , with the results shown in Fig. 7. Once again, the injected value is within the 1σ confidence levels.
Finally, we could repeat the analysis by using the pos- teriors for V (2) , which would in principle break the twoparameter degeneracy of the agnostic case. In practice, we find that the additional constraint is almost parallel to the constraint from V (0) and that it only excludes large values of the deviation parameters, for which the approximations do not hold anyway. These conclusions depend on the specific combination α (2) , α (3) that we have chosen in our example, and they may be qualitatively different for other combinations of the deviation parameters.

IV. DISCUSSION AND CONCLUSIONS
In general, one may hope that modifications to the perturbation potential proportional to α (k) r H /r k would yield smaller corrections to the QNM spectrum as k grows, as long as the coefficients α (k) are comparable in order of magnitude. In fact, while this trend is present, the contributions from terms with k ≥ 2 decay quite slowly. This lack of a strong hierarchy implies that, in general, it will be very difficult to recover the individual coefficients α (k) , at least in the absence of an underlying theory-specific ansatz. This was, indeed, the conclusion of previous work on parametrized ringdown using a principal component analysis [48].
Here we confirm this conclusion by exploring two specific cases: a GR QNM injection, and a non-GR injection. In both cases we look at two extreme scenarios: in the optimistic scenario we vary a single coefficient α (k) at a time, and in the pessimistic scenario we let multiple α (k) 's (with k = 0, . . . 10) vary simultaneously. In both scenarios we find that it is difficult to constrain the individual α (k) 's, as expected. We study a non-GR injection in which we assume that only two of the α (k) 's are non-zero. In the optimistic scenario we find that the posteriors of all of the α (k) 's disfavor GR, but only one of them is close to the correct injected value. Therefore it is difficult to identify which non-zero α (k) 's are present in the data without having more precise data, or additional (theory-specific) criteria to interpret the posteriors.
The main conclusion of this work is that the problem can be by-passed by exploiting the well-known relation between QNM frequencies and the properties of the perturbation potential at the light ring. This relation suggests that the complex correlations present when many α (k) 's are varied simultaneously can be understood by relating the α (k) 's to the value of the potential and its derivatives at the light ring. It is possible to obtain very robust constraints on the potential and its derivatives at the light ring, even in the pessimistic case in which many α (k) 's are varied simultaneously. More precisely, a measurement of the fundamental mode and of the first overtone can constrain quite precisely the value of the potential and of its second derivative evaluated at the peak. As demonstrated, this conclusion is not limited to GR injections, but it is also valid for non-GR injections.
We analyze, as a proof of principle, a hypothetical theory characterized by a single parameter ζ, such that all non-zero deviations have the form α (k) (ζ). Through a simple "post-processing" analysis based on the inferred posterior for the value of the potential at the maximum, and assuming that only two of the α (k) 's are non-zero, we can recover the injected parameter ζ and estimate confidence intervals for each of the α (k) 's. This demonstrates that MCMC constraints on the value of the potential at the maximum can be used to find theory-specific properties, without rerunning a full MCMC analysis.
The intimate link between QNM frequencies and the potential at the light ring implies that there is a close re-lation between strong-gravity tests using black hole spectroscopy and the black hole shadow observations by the Event Horizon Telescope collaboration [65,66,68]. In analogy with QNMs, the strong-field character of the shadow (which depends on the geometry near the light ring) also implies that the individual parameters in a PN expansion of the metric are hard to constrain, especially when many of them are allowed to vary simultaneously [67]. Black hole shadow measurements can still constrain linear combinations of the metric and its derivative at the light ring [67], in a way which is closely reminiscent of the WKB results reported in this work.
In closing, it is important to discuss some caveats and possible future extensions of this work. The parametrized framework [46,47] is only valid for non-rotating or slowly rotating black holes, but current gravitational-wave observations involve rotating black holes. More theorydependent studies of QNMs for rotating black holes in modified theories of gravity are necessary to overcome this nontrivial obstacle, and to understand how a generalized parametrized framework can effectively be implemented and constrained. The increase in complexity resulting from the (generally nonseparable) perturbation equations in modified gravity implies that theory-specific tests may be particularly valuable, because they typically involve a small number of free parameters. 19-ATP19-0051, 20-LPS20-0011 and 21-ATP21-0010. This research project was conducted using computational resources at the Maryland Advanced Research Computing Center (MARCC). The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas. edu [96].
Appendix A: Details for higher order WKB method In this appendix we report some of the lengthy WKB expressions used in the main text. For convenience, we report again the WKB formula (6), up to third order: Here n is the QNM overtone number, α n = n + 1/2, and the higher-order corrections read