Accelerating modified gravity detection from gravitational-wave observations using the Parametrized ringdown spin expansion coefficients formalism
AAccelerating modified gravity detection from gravitational-wave observationsusing the Parametrized ringdown spin expansion coefficients formalism
Gregorio Carullo
1, 2 Dipartimento di Fisica “Enrico Fermi”, Universit`a di Pisa, Pisa I-56127, Italy INFN sezione di Pisa, Pisa I-56127, Italy (Dated: February 12, 2021)Harvesting the full potential of black hole spectroscopy, demands realising the importance ofcasting constraints on modified theories of gravity in a framework as general and robust as possible.Requiring more stringent – yet well-motivated – beyond General Relativity (GR) parametrizations,substantially decreases the number of signals needed to detect a deviation from GR predictionsand increases the number of GR-violating coefficients that can be constrained. To this end, weapply to LIGO-Virgo observations a high-spin version of the Parametrized ringdown spin expansioncoefficients (ParSpec) formalism, encompassing large classes of modified theories of gravity. Weconstrain the lowest-order perturbative deviation of the fundamental ringdown frequency to be δω = − . +0 . − . , when assuming adimensional beyond-GR couplings, substantially improvingupon previously published results. We also establish upper bounds (cid:96) p =2 <
23 km, (cid:96) p =4 <
35 km, (cid:96) p =6 <
42 km on the scale (cid:96) p at which the appearance of new physics is disfavoured, dependingon the mass dimension p of the ringdown coupling. These bounds exceed the ones obtained byprevious analyses or are competitive with existing ones, depending on the specific alternative theoryconsidered, and promise to quickly improve as the number of detectors, sensitivity and duty-cycle ofthe gravitational-wave network steadily increases. I. INTRODUCTION
Within the framework of GR and the current StandardModel of particle physics, black holes (BHs) are theinescapable final state of massive objects ( (cid:38) M (cid:12) )gravitational collapse. Investigating the structureof extremely compact objects, confirming their BHnature and testing the predictions of GR in the BHnear-horizon regime, is one of the main scientific goalsof gravitational-wave astronomy and BH imaging [1–5].Among the wide variety of BH characteristic signatures,a considerable amount of effort has been devoted tobound [2, 6–10] deviations in the quasi-normal modes(QNMs) dominating the ringdown process, the GWemission process driving the post-merger phase of abinary black hole (BBH) coalescence. QNMs emissionis sensitive to the space-time structure close to thelast stable photon orbit [11, 12], allowing precisiontests of gravity in the strong gravitational-field regime.Deviations from the predictions of BH dynamics inGR could either be caused by a modification to thetheory action [13] or by exotic compact objects (ECOs)mimicking BHs behaviour. The latter class of objectscould be constituted by new species of particles or byunknown stable solution of Einstein’s field equationsarising from underlying quantum effects [14].The latest catalog of observations [15] announced bythe LIGO-Virgo collaborations [16, 17] has brought thetotal number of confirmed gravitational-wave eventsto 50. The increasing amount of detections is shiftingthe focus of gravitational-wave observational sciencefrom single events properties to population properties.The signal-to-noise ratio (SNR) of these observationsis weak compared to the ones that should be achieved by both ground-based or satellite future experiments[18–20], but the number of detected events by the currentgravitational-waves detectors network will soon grow tohundreds per year [1]. Although most of these signals willhave a small SNR, thus providing weak constraints per-se,bounds on beyond-GR parametrizations will steadilyimprove with the number of additional detections [7, 21–23]. Current detector network capability of testing theBH ringdown emission could be further enhanced bycavity detuning techniques. Applied in coordination withLISA early-inspiral observations, detuning would allowground-based detectors to vastly increase their sensitivityaround the expected ringdown frequency band andboost spectroscopy measurements precision [24]. Sucha coordination strategy highlights the new and excitingpossibilities that will be opened by multi-band GWobservations, enhancing the achievable precision on testsof GR obtainable by single class of detectors alone [25–28].Given the wide range of possible modifications from GRpredictions, when performing the aforementioned tests itis important to use generic and robust parametrizations.A choice extensively employed in the literature, e.g. bythe LIGO-Virgo collaborations [2], allows to directlycast constraints on observable quantities by consideringparametrized deviations to ringdown frequencies anddamping times. Such an alternative has the advantageof capturing deviations that arise in the observedquantities of interest, without being dependent ona specific alternative theory of gravity. While beingquite general, this approach has the disadvantage ofbeing source dependent, since deviation parameterswill generally differ for different astrophysical sources.Although it is still possible to combine multiple events,construct a global figure of merit on the agreement of a r X i v : . [ g r- q c ] F e b GR with observations (Bayes Factor) and estimate thepopulation distribution of putative deviations [2, 29, 30],an approach where this additional step would not beneeded could decrease the number of signals needed to beconsidered in order to claim the detection of a deviation.Another, sometimes unappreciated, point to be consid-ered when discussing such tests is that in presence ofa GR violation not only the parametrized coefficients,but also the physical BH mass and spin would also differfrom their GR predicted value. These two additionalintrinsic parameters need to be marginalised over,greatly reducing the resolving power of the analysis. Sucha reduction stems from the possibly strong correlationsin the multi-dimensional dimensional space of intrinsicparameters. Detecting a deviation from GR, consequentlyrequires an SNR (or number of signals) larger thanexpected from uncertainty estimates solely based onconsidering GR deviation parameters.A formalism not suffering from the aforementionedresolving-power reduction and allowing to extract source-independent GR deviation parameters, has been recentlyproposed in Ref. [31]. The Parametrized ringdown spinexpansion coefficients (ParSpec) formalism, described indetail in the next sections, provides a robust and accurateparametrization covering large classes of modifiedgravity theories. In this manuscript, we will apply ahigh-spin version of the ParSpec formalism to currentLVC observations, showing how it allows to: i) placetighter constraints on BH spectra deviation parameterswith respect to the ones present in the literature, henceaccelerating the detection of possible deviations fromthe predictions of GR; ii) place competitive bounds onthe energy scale at which effects due to new physics(including corrections from Einstein-scalar-Gauss-Bonnet,dynamical Chern-Simons or Effective Field Theoriesterms in the action) can manifest in BBH systems.The paper is organised as follows: in Sec. II we reviewthe ParSpec formalism, highlighting the cases of interestfor our analysis. In Sec. III we include in the expansionhigher order terms, in order to accomodate high spins, animportant step in order to perform an analysis on LIGO-Virgo data. The main section of the paper, Sec. IV, dealswith the presentation and discussion of results obtainedby applying the ParSpec framework to merger-ringdownobservations from the LIGO-Virgo collaborations, compar-ing them to previously known bounds. Conclusions andfuture prospects are presented in Sec. V. Throughout themanuscript both c = G = 1 units, together with medianand 90% credible levels for parameters point estimatesare used, except where explicitly stated otherwise. The full set of parameters characterising the GW signal inde-pendently of the position of the observer in the sky, such as BHmasses and spins or GR deviation coefficients are referred to as intrinsic parameters.
II. PARSPECA. Basic formalism
In this section we will review the ParSpec formalism,highlighting the features especially relevant to data anal-ysis. For a complete account, see Ref. [31]. The ParSpecframework can be constructed starting from consideringthe following parametrization: ω = ω Kerr · (1 + δω ) ,τ = τ Kerr · (1 + δτ ) . (1)valid for for each QNM. Ignoring sources of systematicerrors that will be discussed in the next sections, the de-viation parameters δω, δτ can take non-zero values due toenvironmental effects, extra-charges beside the BH massand spin, an ECO mimicking the BH emission or modifica-tions to the BH GR-dynamics due to an alternative theoryof gravity. As already discussed, although the expansionpresented in Eqs. 1 is generally valid, a drawback arisesfrom the fact that the deviation parameters will dependon the specific source parameters (such as mass and spin)and thus deviation values will differ for different sources.Defining source-independent parameters would allow todirectly combine constraints from multiple observations.The parametrization efficiency could also be improved(in the sense that a smaller number of observations isneeded to confidently detect a deviation) by imposing aspecific dependence on the BH mass and spin. A way torealise these two desiderata is to consider a perturbativesetting, through the following bivariate expansion, validfor a given GW source: ω K = 1 M N max (cid:88) j =0 χ j ω ( j ) K (1 + γ δω ( j ) K ) ,τ K = M N max (cid:88) j =0 χ j τ ( j ) K (1 + γ δτ ( j ) K ) . (2)where K labels the QNM considered, ω ( j ) K , τ ( j ) K are thesource-independent numerical coefficients governing thespin expansion in GR, M, χ represent the
GR values of theBH mass and spin and N max is the maximum expansionorder considered. The specific scenario giving rise to amodification in the GR expansion is encoded into the δω ( j ) K , δτ ( j ) K coefficients and in the γ parameter. The δω ( j ) K , δτ ( j ) K coefficients set the spin-expansion structureof the modification under consideration, and are thusnumerical constants. The parameter γ is a (possiblysource-dependent) coupling constant, discussed below.Since a modification in the BH dynamics will modify notonly the spectrum, but also the BH intrinsic parameters,the expansion should in principle be written in terms ofnon-GR BH mass and spin parameters ¯ M , ¯ χ . Nonetheless,as proven in Ref. [31], at the perturbative level considered,modifications to the BH mass and spin can be absorbedin the deviation parameters already defined . Hence, theBH parameters measured by assuming that GR is thecorrect description of the BH dynamics ( M, χ ), can beused in computing the spectrum. This simple observationwill have a great impact on the analysis, as discussedin the next sections. The perturbative assumption isimposed by requiring γ δω ( j ) K (cid:28) γ δτ ( j ) K (cid:28)
1) for allvalues of { K, j } considered. The adimensional parameter γ can be expressed in terms of a (possibly dimensionful)coupling constant α . If we denote by p the mass dimensionof α , we can express γ in terms of the BH mass M s in thesource frame, which sets the energy scale of the processunder consideration: γ = αM ps = α (1 + z ) p M p := (cid:18) (cid:96) c (1 + z ) G M (cid:19) p (3)where we expressed α in terms of a parameter (cid:96) , repre-senting the typical length scale below which correctionsfrom new physics affect the system. The redshift z can be obtained from the luminosity distance assuminga standard cosmology (see Ref. [31] for a discussionon the, quite general, validity of this assumption).Considering the generally dimensionful coupling α ,two cases have to be distinguished. In the first case, α is a source-independent parameter, function of thecoupling constant(s) appearing in the action of the newtheory and it is referred to as secondary hair . Theadimensional case ( p = 0) includes certain scalar-tensortheories, Einstein-Aether and Horava gravity [32], the p = 4 case encompasses Einstein-scalar-Gauss-Bonnetand dynamical Chern-Simons gravity [33–45], whilethe p = 6 case includes the (quite general) class ofEffective Field Theories considered in Refs. [46–48]. Foran Effective Field Theory approach to scalar-tensortheories, see instead Ref. [49, 50]. In the second case, α is an additional BH charge independent from themass and spin of the BH, a primary hair . The presenceof a BH U (1) charge is included in this latter casewhen p = 2. In this scenario, the ParSpec formalismparametrizes violations of the Kerr hypothesis, due e.g.to the electric charge, in the gravitational modes whichare continuously connected to the l = m = 2 QNM inthe Schwarzschild case [51–58]. See Ref. [31, 59, 60]for an account on the classes of theories that can bemapped to the ParSpec formalism and Refs. [32, 61, 62]for extensive reviews on modified theories of gravityand their observational effects. Without listing in This implies that the deviation parameters considered will in-clude a combination of intrinsic corrections to the spectrum andcorrections to the mass and spin of the remnant BH, as discussedin Appendix A of Ref. [31]. At the perturbative level, these twotypes of corrections cannot be disentangled. detail all the implications that a measurements of theparameters considered in the ParSpec formalism entail(for a comprehensive review, see Ref. [63]), we simplyremark that constraints on these classes of theories haveimplications on the deviations from some of the theoreti-cal pillars of GR, such as the Equivalence Principle andLorentz-invariance. Furthermore, a variety of beyond-GRemission mechanisms, e.g. scalar fields emission, can leada measurable effect in the observations considered andthus be constrained with BBH observations. The outlinedapproach is perturbative in spirit and thus retains thesame merger-ringdown structure as the one present inGR. Nonetheless, deviations parameters will peak awayfrom zero even if large deviations are present or if thestructure of the merger-ringdown differs from the one ofGR (e.g. due to additional modes or non-exponentialtime-dependent terms, as in Ref. [43]). In this lattercase the values of the deviation parameter can only besubject to a phenomenological interpretation and shouldbe considered as detection indicators.
B. Comparison to alternative parametrizations
Given the extreme variety of new physics effectsthat could manifest in BBH coalescences (either inthe form of alternative compact objects or stemmingfrom a high-energy extension of GR), different choicesto parametrize possible deviation from the GR BBHhypothesis are possible.An appealing possibility would be to consider an effectivemetric in a suitable coordinate frame, encompassingboth the GR and beyond-GR case, to describe the BHspacetime. For a thorough overview of the propertiesand regions of validity of several modified Kerr metrics,see Ref. [64]. Considering deviations in the metricallows to cast constraints directly on a fundamental,although coordinate dependent, object determiningthe BH geometry. This approach was adopted inRef. [65], where observations of the M87 BH shadowobtained from the Einstein Horizon Telescope (EHT),were employed to constrain the second post-Newtonian(2 PN) terms of effective non-Kerr metrics (but seeRef. [66] for discussions about the underlying statisticalassumptions and Ref. [67] for a critical assessment of allthe underlying astrophysical assumptions entering theanalysis). An extension combining shadow measurementswith constraints from BBH inspiral signals was recentlyconsidered in Ref. [68], finding no deviations from thepredictions put forward by the Kerr metric in GR. Thesame approach could in principle be used to extractconstraints from the merger-ringdown regime of BBHremnants on an effective metric. This possibility wasexplicitly considered in the case of non-spinning BHsin Ref. [69], for a simulated dataset in the context ofnext-generation GW detectors. Unfortunately, extendingthis procedure to the rotating case, while maintaining anagnostic approach, presents non-trivial complications [69].A possibility to overcome some of these complicationsin the rotating case, would be to consider an eikonalapproximation, as described in Ref. [70] (for applicationsof the eikonal approximation to the spherically symmetriccase of alternative theories of gravity, see Refs. [71, 72]).Although coming with its own set of approximationsand thus not generally valid, this approach wouldconstitute quite a robust and generic formalism toproject the constraints obtained from observations ofringdown spectra from rotating BHs for the fundamental (cid:96) = | m | subset of modes. From the pure point of viewof data analysis, constructing templates (or even simplenumerical interpolants) of QNM frequencies as a functionof an effective metric parameters would constitute anefficient way to test gravity in this approach. Such ascenario seems challenging from the theoretical pointof view, but not unfeasible in principle once a suitableformulation of the perturbation equations is found. Theconstruction of such a template would be an extremelyinteresting avenue to be explored in the near future.Instead, solving numerically the perturbation equationsfor a given value of the effective metric parameter(s)during the analysis itself, currently seems unfeasible.In fact, the computational cost that such an operationwould add to the exploration of the parameter space, atleast within current state-of-the-art stochastic samplers,would presumably not be sustainable. Machine learningmethods would likely need to be introduced in orderto decrease computational cost. The most seriouslimitation of the aforementioned approaches focusing onparametrized metrics, stems from the far richer dynamicsthat is present in beyond-GR theories when modificationsto the field equations are consistently included [60].A set of analyses similar in spirit to the ones consideredin this work, was notably put forward by the LIGO-Virgocollaborations in the recent catalog of tests of GeneralRelativity [2]. In the latter reference, a battery of testswas considered, investigating: possible modifications inboth the generation and propagation of gravitational-waves, the nature of the coalescing objects, additionalpolarisations absent in GR and residuals in the detectorstrain after subtracting a best-fit waveform. None ofthese tests signalled a discrepancy among the predictionsput forward by GR and the observed BBH population.Of particular interest for this work are the set oftests exploring parametrized deviations in the GRpredicted waveform. The reason why presumably thesetests are not formulated as perturbative expansions interms of a ”small” parameter, is that this approachwould narrow down the set of possible theories towhich they are sensitive to. When considering a fullinspiral-plunge-merger-ringdown analysis, in the inspiralregime the deviations are applied to a set of coefficientsderived from PN theory in both the EOB and Phenomset of waveforms [73–78]; instead in the plunge-merger-ringdown regimes, the deviation coefficients are appliedto phenomenological quantities extracted from complete numerical solutions of the Einstein equations, within thePhenom family. The expansion of Eqs. 1 was insteadused both in the case of the parametrized-EOB [7]model, including in the analysis the whole signal, orin the pyRing approach [8, 9], which in order to avoidcontamination from earlier stages of the coalescence,formulates the problem completely in the time domainand restricts the analysis to the post-merger-only section.The precise set of theories, and the parameter spaceexplored within each alternative theory by such aparametrization, depends on the specific choice of priorbounds used and is difficult to compute due to the lackof accurate or complete inspiral-plunge-merger-ringdownpredictions from alternative theories. As already men-tioned, one of the disadvantages of this latter approach isthat the deviations parameters can take different valuesfor different sources and it is thus not possible to directlycombine the obtained constraints. Such a difficulty hasbeen overcome using a hierarchical approach [29, 30],which places constraints on the underlying observedpopulation distribution of such deviations (a Diracdistribution centered around the null value, if GR iscorrect and the measurement has infinite precision).This approach is the optimal – following directly fromBayes theorem – procedure of ”stacking” informationfrom multiple sources in the absence of additionalhypotheses on the underlying alternative theory ofgravity. These tests also allow the intrinsic (masses andspins) parameters of the binary to deviate from thevalues which can be extracted by assuming that GRis the correct description of the signal. This feature isrequired since without a perturbative formulation, theintrinsic parameters entering the waveform are ( ¯ M , ¯ χ ),not ( M, χ ). Once again, this has the effect of allowingfor more freedom in the reconstructed signal, but it addsa large amount of correlation among the intrinsic binaryparameters and the deviation parameters. As ubiquitousin statistics, the gain in generality that these tests havein both allowing non-perturbative deviations in thewaveform parameters and in considering deviations inthe mass and spin parameters, degrades the sensitivity ofthe tests.The ParSpec formalism instead models possible deviationsfrom GR in a perturbative spirit. Deviations in bothbinary parameters and coefficients determined by theunderlying dynamics are considered within expansionswith respect to the GR value. An implication of thisapproach is that now the expansion can be formulatedin terms of the GR binary parameters and consequentlythese parameters can be restricted within the supportof the posterior probability distribution obtained whenassuming GR, strongly reducing the amount of correlationpresent in the problem. This approach restricts thedeviations values considered, but strongly enhances theresolving power of a test, allowing to detect deviationswith a much smaller number of sources that would beneeded by the more generic tests considered by theLIGO-Virgo collaborations. Furthermore, as alreadydiscussed above, the restrictions applied in the parameterspace do not affect the capabilities of the frameworkto detect a violation from GR even for scenarios wherenon-perturbative signal deviations are present.
III. PARSPEC: HIGH-SPIN VERSIONA. Fit framework
The accuracy of the expansion contained in Eq. (2)increases with the maximum order of the expansion N max and is extendable to arbitrarily large spins. Difficulties ina high-order expansions arise from the need of simultane-ously fitting a large number of the coefficient, particularlyif the expansion is carried directly on the damping time τ ,which gives rise to a more complex structure with respectto different parametrization such as the quality factor Q = ω r ω i . These difficulties can nonetheless be circum-vented by employing sophisticated fitting algorithms ableto capture the full multi-dimensional space structure ofthe – often strongly – correlated coefficients. In this workwe make use of a Nested-Sampling algorithm, CPNest [79]able to reconstruct multi-dimensional, strongly correlatedposterior distributions (see [80] for an introduction toNested Sampling). The algorithm has been validatedwith a wide range of high dimensional analytic problems(e.g. Rosenbrock function, 50-dimensional gaussian, egg-box profile), by reproducing LIGO-Virgo measurementsof compact binary GW observations and is routinely em-ployed in the production of merger-ringdown results bythe LIGO-Virgo collaborations [2]. The algorithm has alsobeen used to extract post-Newtonian coefficients from fullnumerical simulations of the Einstein equation, describingthe coalescence of two compact objects. Once a suitableand robust parametrization is provided, analytically un-known coefficients can be extracted from numerical data,[81]. See also [82–84] for alternative approaches on similarcases. Applying the algorithm to our problem requiresto formulate the fit in terms of a simple inference prob-lem in Bayesian statistics (for a rigorous formulation ofBayesian probability theory as an extension of classicallogic, see Ref. [85]). If the coefficients to be extracted arecollectively represented by (cid:126)θ , the probability distributionon the coefficients conditioned on the available numericaldata D , the posterior distribution , is computed throughBayes theorem: p ( (cid:126)θ | d, H , I ) = p ( (cid:126)θ |H , I ) · p ( d | (cid:126)θ, H , I ) p ( d |H , I ) (4)where H represents the specific parametric template de-scribing the data (our hypothesis) and I denotes all theavailable background information. The prior distribution p ( (cid:126)θ |H , I ) encodes all the information available on the coef-ficients before the inference process (e.g. positive defined, bounded by some physical constraint, etc.). If no a-prioriinformation is available, the prior can be conservativelychosen to be uninformative. The likelihood p ( d | (cid:126)θ, H , I ) isdetermined by the error distribution of the available nu-merical data. The normalisation term Z := p ( d |H , I ), the evidence , quantifies the probability that the data d can beexplained by the model under consideration. The name evidence can be understood by the role of this term whencomparing the agreement with the data of two competingmodels H , H : p ( H | d, I ) p ( H | d, I ) = p ( H | I ) p ( d |H , I ) p ( H | I ) p ( d |H , I ) := p ( H | I ) p ( H | I ) · B (5)The first term in the above equation represents the ratioof prior probabilities of the models; the second term, the evidences ratio is called Bayes factor . When comparingGR predictions against those of a specific modified theoryof gravity, we will quote this number. The advantages ofsuch an approach with respect to standard χ methodsmainly consist in: i) providing the full multi-dimensionalprobability distributions on the parameters, encodingthe global correlation structure of the problem and thusgoing beyond simple mono-dimensional error estimateson each of the parameters; ii) avoid overfitting by usingthe Bayesian evidence as a tool to select the best numberof parameters needed to fit the dataset . B. Results
We use Eq. (2) as a template, hence defining ourhypothesis H and apply Bayes Theorem to extractposterior probability distributions on the parameters (cid:126)θ = { ω ( j ) , τ ( j ) } j =1 ,...,N max . The dataset d employed isthe one of Ref. [86], considering the modes relevant tocurrent or near-future observational analyses, namely( l, m, n ) = { (2 , , , (2 , , , (3 , , , (2 , , } . Weset uniform priors on all the coefficients with bounds[-100,100]. The likelihood function is determined underthe assumption that the error on the numerical data isdistributed as a zero-mean gaussian with a conservativestandard deviation estimate of 10 − , well below theprecision achievable by current merger-ringdown dataanalyses. Although the typical adimensional spin ofa quasi-circular BBH remnant does not reach valuesclose to extremality, the current low-SNR implies thatposterior distributions have sometimes support up to Although the latter point is important in principle, it is onlymarginally relevant for the problem at hand, since we are inter-ested in an effective low-order polynomial representation of thecomplex ringdown frequencies. Overfitting issues would be rele-vant only for very high order expansions saturating the amountof information contained in the numerical data. As shown below,current statistical uncertainties keep the number of needed coeffi-cients small enough that overfitting is not presently a concern. high spins ( a ∼ .
95) [2]. Consequently we decide toconservatively fit numerical data up to a max = 0 .
99. Theorder of the expansion N max is dictated by the requiredprecision in representing the data. Expected statisticalerrors on the δω ( δτ ) parameters are of the order of fewpercent (few tens of percent); thus, as a conservativechoice, we keep including terms in the expansion until theresiduals go below 1% (3%). The zeroth-order coefficientsare fixed from the known Schwarzschild values. Thesampler settings used are: 2048 live points, 1024 as amaximum number of Markov-Chain Montecarlo internalsteps and a pool-size composed of 256 walkers. Thestopping condition is set by requiring an estimatedprecision of 0.1 on the logarithm of the evidence.Obtained results are collected in Table I, where we showmedian and 90% credible level values for the N max = 5 (9)expansion order, which provides residuals below 1%(3%)on the ω ( τ ) fitting parameters. For high expansion or-ders, τ fitting parameters tend to rail against the priorbounds (we checked that wider bounds do not alleviatethe problem), showing large cancellations. Such behaviouris tipically observed when polynomial expansions are usedto represent highly-structured functions. This is not aconcern, because we only need an effective expansionwhich is accurate enough for our purposes; the fact thatpolynomial expansions are a poor fit of QNM complexfrequencies has long been known and indeed the mostaccurate parametric representations of these data tipicallyinvolve non-rational functions [87–89]. A polynomial ex-pansion has the advantage of being easily interpretablewhen analysing observations, but most importantly to beeasily generalisable to arbitrary alternative theories withadditional charges. As a future development it wouldbe interesting to explore other parametrizations with re-spect to the simple polynomial one considered here andin Ref. [31]. IV. OBSERVATIONAL CONSTRAINTSA. Comparison to future detector observations andcurrent expected precision
Before applying the high-spin version of ParSpec toobservational data, it is instructive to estimate the levelof precision we expect to achieve on deviation coefficients,by extrapolating the estimates obtain in Ref. [31]. In thelatter study, the ParSpec formalism was applied to a setof simulated signals expected to be detected by the ETinterferometer and the LISA mission [18, 20]. The widthof the obtained probability distributions on δω were of σ ∼ × − , at 90% credible levels, when allowing onlya single deviation parameter and combining 10 sources.In this work, we instead apply ParSpec to the only cur-rently available ringdown observations, the post-mergerphase of binary BHs coalescences observed by the LIGOand Virgo interferometers. The signal-to-noise ratio (SNR) of these sources lies in the range 5 −
15, muchsmaller compared to expected SNR values from ET/LISAobservations, O (100 − ∼ / SNR,we expect to be able to place constraints of the orderof 3 × − on δω , under the same hypothesis on thenumber of deviation parameters and number of combinedsources. In a realistic setting, as the one being consid-ered in this work, discrepancies from this simple estimatearise from: i) poor accuracy of the the aforementionedSNR scaling in the low-SNR limit; ii) imperfect wide-sense stationarity or gaussianity of the stochastic processgoverning the interferometer output, and iii) the marginal-isation of additional parameters other than parametricdeviations. The last point may also introduce, throughcorrelations, deviations in the gaussian approximation(e.g. wider tails) of the posterior probability distributionsmade in Ref. [31]. All these features will contribute towiden the posterior probability distribution and implythat the estimate presented above can only represent anoptimistic lower bound on the precision we can obtain.When comparing constraints drawn from observations ofstellar-mass BHs (GW ground-based detections) againstthose obtained from supermassive BHs (GW space-baseddetections or BH shadows observations), it should berecalled that these different mass scales actually probedifferent curvature regimes. Although ringdown signalsgenererated by mergers of BBH in LISA will be loud,the corresponding space-time curvature (quantified bythe Kretschmann invariant) associated to the emitters,massive BBHs, will be smaller compared to the one asso-ciated with stellar-mass BH mergers observed by currentand future ground-based detectors. Given that GR devia-tions are expected to be prominent in the high-curvatureregime, the high-frequency (larger than kHz) regime [90]will probe the most interesting regions of the curvature pa-rameter space. Observations of post-merger signals fromlow-mass compact objects will also unveil the formationmechanisms of remnants lying on the border of the BHlow-mass gap, extending spectroscopic measurements toECOs that may avoid gravitational collapse only in a lim-ited mass range or even to primordial BHs mergers [14].In summary, constraints placed by space and ground-based detectors should be regarded as complementary inthe curvature parameter space (for a comprehensive re-view of the possible constraints on alternative theories ofgravity when considering future ground and space-baseddetectors, see Ref. [91]). Planned ground-based detectorssuch as NEMO [90], targeting GW observations in thekHz-band, could soon open an observational window ontothis exciting frequency band. B. Analysis infrastructure: pyRing
To constrain deviation coefficients from Eq. (2), we use pyRing [8, 9], a ringdown-tailored parameter estimation python [92] package, employing a Bayesian approach and
TABLE I. Numerical results for the coefficients of the polynomial expansion considered in Eq. (2). The c n coefficients for a givenmode correspond to ω ( j ) K , τ ( j ) K coefficients in Eq. (2), for K = { (2 , , , (2 , , , (2 , , , (3 , , } , j = 0 , · · · ,
9. The lowest-ordercoefficients are obtained from perturbation theory of Schwarzschild BHs. ω (l,m,n) (2 , ,
0) (2 , ,
1) (2 , ,
0) (3 , , . . . . . +0 . − . . +0 . − . . +0 . − . . +0 . − . c − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . c . +0 . − . . +0 . − . . +0 . − . . +0 . − . c − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . c . +0 . − . . +0 . − . . +0 . − . . +0 . − . τ (l,m,n) (2 , ,
0) (2 , ,
1) (2 , ,
0) (3 , , . . . . . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . c − . +0 . − . . +0 . − . . +0 . − . − . +0 . − . c − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . c . +0 . − . . +0 . − . . +0 . − . . +0 . − . c . +0 . − . . +0 . − . . +0 . − . . +0 . − . c − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . c − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . c . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . c . +0 . − . . +0 . − . . +0 . − . . +0 . − . a fully time-domain formulation of the inference problem,both in the likelihood and in the waveform. The packageis based on the CPNest [79] sampler and relies on the
LALInference library [93] to compute projections ontodetectors. The package internally implements most ofthe state-of-the-art merger-ringdown-only analyticalwaveform templates available in the literature [94–97] andaccesses other standard inspiral-plunge-merger-ringdownwaveforms thorugh the
LALSimulation library [98]. Bothringdown-only template injections and complete inspiral-plunge-merger-ringdown injections can be performed;the latter are available both for analytical templatesand for numerical simulations accessed through the
LALSimulation library and the LIGO-Virgo NumericalRelativity Injection Infrastructure [99, 100]. To optimisecomputations, pyRing internal waveforms and alllikelihood calculations are implemented in cython [101].Computations involving the Auto-Covariance matrixare drastically optimised using some of the state-of-the-art algorithms introduced in Ref. [102], exploitingrecursive relations obtained from the Toeplitz structureof the Auto-Covariance matrix. QNM frequenciescan be pre-computed using
SciPy [103] cubic-splineinterpolation of publicly available numerical data fromRef. [87], publicly avaialable at [86]. Remnant BHmass and spin are computed using the fits presented inRef. [104–106] and accessed through the
LALSimulation library or the correspective python packages. The pyRing package was recently used to produce the first catalogof ringdown-only remnant properties measurements andconstraints on QNM complex frequencies parametricdeviations, using data from the first three observing runs of the LIGO-Virgo interferometers, see Tables VIII-IXof Ref. [2]. We adopt the same selection criterion ofRef. [2] to select the list of events on which we apply theParSpec formalism, together with the same strain data,data conditioning methods and auto-correlation-functionparameters. Additionally, when considering δτ deviations,we exclude from the analyses the events: GW170814,GW190521, GW190828 063405. The inference on thedamping time, a parameter currently only weaklymeasured, of these events has in fact been shown tosuffer from systematic errors due to behaviour of thestochastic process governing LIGO-Virgo noise [2, 107].Analyses results and settings used in Ref. [2] are publiclyavailable from the accompanying data release [108]. Thischoice sets the number of events here considered to 17(14 when considering δτ parameters) out of a total of 50BBH events detected by the LIGO-Virgo collaborations,minimising the effect of noise mimicking a GW event andcontaminating tests of GR. The waveform template usedis Kerr , a superposition of the longest-lived mode andits first overtone ( n = 0 ,
1) with free complex amplitudes,for details see Refs. [2, 9, 109]. The analysis segmentstarts at the peak of h + h × in each detector. Complexfrequencies are expressed using the fits of Eq. (2),extended to large spin values as discussed in Sec. III.Consistently with the ParSpec perturbative approach,we set uniform bounds on all deviation parameters in[ − . , .
5] range, on the (cid:96) scale in the [0 ,
75] km intervaland marginalise the remnant mass and spin withinthe 90% CLs of their GR values. To break the strongdegeneracy with the distance (which is poorly measuredfor merger-ringdown-only signals) in the dimensionful α case, we also restrict it to its 90% GR CLs interval andapply a uniform prior within this interval. All the otherpriors and analysis parameters are chosen to be identicalto those of Ref. [2]. In a realistic scenario, it must berecalled that non-zero values of deviation parameterscan be caused not only by violations of GR, but fromviolations in any of the underlying analysis’ hypotheses.One of such examples would be an incorrect descriptionof the detector noise or missing physics predicted by GR,but not taken into account into our waveform model (e.g.eccentricity, extreme precession, environmental effects).As verified in Ref. [2], all of these assumptions hold forcurrent LIGO-Virgo observations. C. Combined probability distribution functions:boundary artifacts
Being the coupling parameter positively-defined, itsposterior distribution will rail against the lower or upperbound of the prior, at least for measurements which donot exclude the extremes, as the ones considered in thiswork. Consequently, when constructing an estimate ofthe underlying (continuous) probability distribution func-tion from the discrete samples available, it is necessary totake special care in dealing with artifacts at the borders,which could bias the measurement. A continuous estimateof each event probability distribution function is neededwhen producing a combined estimate of the coupling pa-rameter using multiple events, the latter requiring themultiplication of single-event likelihoods (given the statis-tical independence of the different GW events). Since flatpriors are applied on all the parameters of interest, we canemploy posterior probability distributions in place of thelikelihood functions, renormalising appropriately. A firststep in alleviating border artifacts is to use a logarithmictransformation, for example: y ( x ) = log (cid:18) x − x m x M − x (cid:19) (6)where x m , x M represent the lower and upper border ofthe support considered. The transformation stretches thedistribution boundaries to infinity. The resulting trans-formed distribution then presents a continuous behaviourat the boundaries and producing an estimate free from bor-der artifacts is simpler. The final estimate of the originaldistribution is then obtained by multiplying the estimateof the transformed distribution with the transformationjacobian | dydx | . To compute the continuous estimate ofthe transformed distribution, we attempted using both astandard Gaussian Kernel Density Estimation (GKDE)method [103] and a Dirichlet Process Gaussian-mixture(DPGM) model, a fully Bayesian non-parametric method.For a pedagogical introduction, as well as an open-sourceimplementation of the DPGM model, see Ref. [110]. The DPGM model proved superior in alleviating border arti-facts, as directly verified on single-event measurementswhere it provided in all cases a faithful representationof the corresponding samples, and was thus chosen toproduce the combined measurements presented in theremainder of the work. D. Results
We now present constraints on the deviation coefficientsand coupling parameter appearing in Eq. (2). As GRmodel, we always retain the full GR prediction, keepingthe complete expansion up to N max = 5, for frequencyparameters ( N max = 9, for damping time parameters), asdetermined in Sec. III. Deviations from GR predictionsare instead considered up to a given order M max , whileall higher-order deviation parameters are fixed to zero.For each of the charge exponent values p = 0 , , ,
6, weconsider M max = 0 , , . In turn we consider deviations on the fol-lowing set of parameters: { ω , ω , τ , τ } . Due tocomputational constraints, we vary coefficients relativeto only one of these parameters at a time. When con-sidering M max > ω for p = 2 andM max = 2, we simultaneously sample on the coefficients: { (cid:96), δω (0)220 , δω (1)220 , δω (2)220 } . Results on the deviations param-eters of interest are reported in Table II, for differentvalues of { p, M max } , obtained by combining constraintsfrom all the events in the merger-ringdown catalog. Nostatistically significant deviation from the predictions ofGR BH dynamics is found. As expected, tighter con-straints can be placed on the δω coefficients comparedto δτ and on the the fundamental ( n = 0) mode com-pared to the first overtone ( n = 1). The first effect isdue to the larger sensitivity of GW detectors to phasedeviations with respect to amplitude variations, while thesecond to the higher SNR carried by the fundamentalmode compared to its overtones. For each mode, the un-certainty on parameters’ distributions at different ordersis consistent with the magnitude of the correspondentexpansion numerical coefficients. Statistical uncertaintiesare, as expected, larger than the previously discussed esti-mates extrapolated from the simplified setting of Ref. [31].As also recently discussed in Ref. [66] in the context ofBH shadows observations, considering multiple deviationcoefficients at the same time can lead to non-trivial cor-relations among the coefficients and, depending on the The approximation induced by the truncation of the expansionat second order has been explicitly validated in Ref. [111] inthe context of Einstein-scalar-Gauss-Bonnet gravity, where itwas shown that including terms up to the fourth order induceddifferences much smaller than the current statistical uncertaintyon the parameters. information content of the signal, possibly widen the prob-ability distributions extracted from the analysis. In mostcases, if a deviation in the BH dynamics is present, it willaffect the expansion at all orders; consequently, to per-form a complete analysis all coefficients should be variedat the same time. This, however, needs to be balancedby the limited SNR available, by computational cost ofthe analysis and by the capabilities of current samplersto handle very high-dimensional problems. Varying co-efficients one at a time, as done for example in currentparametrized tests by the LVC collaboration, is still avalid method to perform a null test: if no deviations fromGR are present, none of the coefficients should peak awayfrom zero. Increasing too much the number of free param-eters could mask small deviations by increasing statisticaluncertainties.To further highlight this point, in Fig. (1) we show re-sults on the lowest-order (non-spinning) parameters forscalar charges ( p = 0), when considering different M max .The distribution on the same coefficient may widen whenhigher-order terms are included, as expected due to multi-dimensional correlations. Such degeneracies are difficultto break due to the fact that progenitors with compara-ble masses and small spins produce remnants with spins ∼ .
7. Indeed, events with a detectable merger-ringdownsignal reported in the GWTC-2 catalog, present remnantspins in the range [0.6-0.8], as reported in Table VIII ofRef. [2]. Future observations of events with lower remnantspin (requiring a low-spin, high mass ratio merger or highcounter-aligned spins with respect to the orbital angu-lar momentum of the binary) and a detectable merger-ringdown signal will help probing the lower-spin regionand further break the aforementioned degeneracies. Toquantify the overall agreement of the observed signals withrespect to the predictions of GR, keeping into account thefull correlation structure of the signal, we compute BayesFactors among the violating hypotheses (represented byeach of the cases considered in the analysis, as discussedabove) and GR. The values are reported in Table III. TheBayes Factors generally decrease with M max , penalisingthe inclusion of more parameters, without support fromthe data. They indicate that, when considering scalar( p = 0) deviations in the fundamental frequency δω ,the data strongly disfavour a violation from GR. Follow-ing the reasoning outlined above, as expected strongerconstraints can be extracted from the frequencies withrespect to the damping times and from the fundamental( n = 0) modes with respect to the overtone ( n = 1). Theremainder of the cases generally show a mild preference forthe absence of any GR violation. In all cases no significantdeviation from the GR predictions are present within thestatistical uncertainty considered. The combined BayesFactors are computed under the hypothesis that the sameunderlying dynamics describes all the observations, withpossibly different coupling values for different events, asin the case of a BH possessing an additional U (1) charge(e.g. electric charge). As such, they do not apply to casesof mixed populations of GR BHs and alternative compact objects. To cover also these cases, we checked the cor-responding posteriors and Bayes Factor for each event,none of them showing a statistically significant preferencefor the violating-GR hypothesis. Finally, constraints onthe coupling parameter (cid:96) for different mass dimensions p and expansion orders M max are presented in Table IVand Figs. 2,3,4. Since none of the distributions excludeszero, in Table IV we only report constraints on the chargeparameter as 90% upper bounds.Similarly to what has been discussed for the scalar correc-tions, the strongest bounds come from the fundamentalfrequency, as expected. In the informative cases (mainlythe ones where deviations in δω are considered) thebound correctly shrinks for higher expansion orders, sincefor higher orders the coupling parameter enters a largernumber of terms, increasing its possible impact on thewaveform morphology and breaking the correlations witheach deviation parameter. E. Discussion
In this section we will put into context our results,comparing them with previous bounds obtained whenconsidering either agnostic parametrizations or specifictheories, using both GW and EM observations. Whencomparing the bounds presented in Tables II,IV, to thoseobtained through different analyses present in the litera-ture that considered specific alternative theories of gravity,care must be taken into ensuring that we adapt all theresults to the same conventions. In our analysis, (cid:96) is defined to be the length scale appearing in the merger-ringdown corrections as specified in Eq. (2). Instead, for agiven modified theory of gravity, (cid:96) is defined through thecoupling constant(s) appearing in the action. These twodefinitions may thus differ by multiplicative factors, whichwould also inversely affect our definition of the parametricdeviations, since the observable frequency shift is invari-ant. Also, different authors may sometimes use differentdefinitions of the couplings (again up to a multiplicativeconstant of order unity). For a review of the differentconventions used in the context of Gauss-Bonnet gravitysee Ref. [112]. The units for the quantities obtained inthis work are set by the requirement that γ correctionsenter with constant factors of O (1) and δω, δτ correctionsare comparable with the ones obtained in this work. Scalar coupling: p = – The stronger constraints im-posed by the ParSpec framework translate in a shrinkageof the uncertainties on deviation parameters. For compar-ison, in the scalar case ( p = 0), when considering a singledeviation coefficient, the uncertainty of the lowest-orderdeviation parameter on the first overtone of the funda-mental frequency δω is ∼ δω is tighterby a factor of 7 than the best known result [2], although itshould be noted that due to computational requirements,0this latter value was obtained when considering only asubset of the full dataset. Since the uncertainty on theparameters scales roughly as N − ev , this implies a reduc-tion of a factor ∼ ∼ O (1%)limit, as foreseen in Ref. [21]. Dimensionful coupling: p = – In this case we ob-tain a bound of (cid:96) (cid:46)
35 km. Two notable examples ofquadratic gravity theories with a QNM spectrum featur-ing p = 4 corrections, are Einstein-scalar-Gauss-Bonnet(EsGB) and dynamical-Chern-Simons gravity (dCS).Known bounds on the coupling constant of EsGB gravityare summarised in Table II of Ref. [91] and Table I ofRef. [114]. In the latter work, numerical simulations ofthe gravitational waveform from a BBH system consistentwith GW150914 were performed. See also Ref. [112] fora focus on the scalar-field dynamics and Ref. [115] forrecent advancements in the analytical computations of theEsGB Post-Newtonian expansion. For this theory, one ofthe best known observational bounds is (cid:96) (cid:46) −
10 km,estimated from X-ray observations of low-mass BHs inRef. [116] (although it has to be noted such systems aremuch more prone to astrophysical uncertainties with re-spect to GW observations). A comparable bound wasderived in Ref. [117, 118], by applying to this theory theconstraints obtained by the LVC from the early inspi-ral phase of a few GWTC-1 events [107]. As noted inRefs.[112, 119] a theoretical bound from the existence of aBH in EsGB gravity imposes an upper limit on the cou-pling (cid:96) c G M < √ . (cid:96) < . (cid:96) < . M = 2 . − . M (cid:12) , either the lightest black holeor the heaviest neutron star ever discovered in a doublecompact-object system . If the secondary object is indeed Where we adapted the result presented in Ref. [116] to the unitsof this paper, by comparing the ringdown corrections obtainedin this work to the results of Ref. [111], obtaining a conversionfactor of √ π · Y , where Y = 2 −
10 due to the uncertaintyin the parameters reconstruction. For Y ∼
10 our units areconsistent with those of Ref. [112]. Lacking an explicit prediction,for simplicity we also apply the same Y conversion factor to thedCS bound. This is actually a conservative upper limit, since depending onthe BH spin and on the specific coupling with the scalar-field, thebound decreases. Given the large mass ratio of the binary, tidal effects are highlysuppressed and consequently these two hypotheses could not be interpreted as a BH, this observation provides the bound (cid:96) < . tightest ever obtained on thistheory.Instead, for dCS gravity, combining NICER [122–124] andLIGO-Virgo neutron stars observations [125], Ref. [126]recently put forward the most stringent estimate to dateon the dCS length scale of (cid:96) (cid:46) −
50 km. No con-straints from the inspiral of GWTC-1 events was possiblein this case, due to the failure of the small-coupling as-sumption for the parameters of these events [117]. Seealso Ref. [44] for the prospects of constraining dCS us-ing numerical simulations of BBH mergers and LIGO-Virgo observations and Ref. [127] for an analysis usingthe non-detection of amplitude birifringence in GWTC-2 signals to obtain a bound on non-dynamical Chern-Simons (cid:96) (cid:46) km. Both EsGB and dCS gravity evadeconstraints from binary-pulsar observations [117]; solar-system measurements [128] can only place bounds ondCS several orders of magnitudes weaker than those pre-sented here, while no bound on EsGB gravity can beobtained from these latter observations [117]. It shouldbe stressed that all the aforementioned bounds rely onmeasurements obtained in the non-dynamical or low-spinregime, while measurements extracted from the merger-ringdown regimes of binary black coalescences, discussedin this manuscript, are obtained from a highly-dynamical,highly spinning remnant compact object dynamics. Dimensionful coupling: p = – For this case, theobtained bound is (cid:96) (cid:46)
42 km. The class of EFTs con-structed in Ref. [46] is one of the most important exam-ples of theories whose spectra possess a p = 6 correction.Comparing with previous analyses, this result is bothconsistent with what predicted by the QNM analysis ofRef. [47] (for a generic formalism including this classof EFTs, see Ref. [60]), and improves on the results presented in Ref. [48]. This latter work analysed theearly-inspiral phases of a few GWTC-1 observations, dis-favouring the appearance of new physics on a scale largerthan ∼
150 km.In the dimensionful cases, the measurement of the (cid:96) pa-rameter is still strongly affected by the low number ofevents considered. Louder and more numerous detec-tions from future gravitational-waves detector networkobserving runs will significantly shrink the bounds pre-sented. Improvements in the high-frequency sensitivityband of GW detectors will have the highest impact on thismeasurement, since being able to probe merger-ringdownemission from low-mass BHs will translate into sensitivityon smaller scales (higher energy) effects. The boundspresented on the Schwarzschild contribution to the spec-trum open up the possibility of constraining those the-ories for which only spherically symmetric solutions are distinguished. From the results of Ref. [47] it can be deduced that the units inwhich this results is quoted are the same as the ones used in thispaper up to factors of O (1). V. CONCLUSIONS
We presented the application of a high-spin version ofthe ParSpec formalism to observational data collectedfrom the LIGO-Virgo interferometers during the O1-O2-O3a observing runs. Theoretical frameworks as the oneconsidered in this work, allow to extract tighter con-straints on modified-gravity deviation parameters withrespect to less stringent parametrizations. We showedhow combining current detections within ParSpec allowsto start placing bounds on BH spectra at the O (5 − (cid:96) (cid:46) −
40 km, depending on the cou-pling mass dimension. The bounds obtained either exceedprevious ones present in the literature, or are comparableto some of the best-known estimates derived from otherexperiments, depending on the specific theory considered.Unlike previous results, such bounds were obtained ina highly dynamical, highly spinning regime. Upcomingdata releases and future observing runs from the LIGO-Virgo-Kagra network promise steady improvements to themeasurements presented. The measurement precision maysoon enter the regime where corrections from EffectiveField Theories of beyond-GR gravity could be detectable,offering the concrete prospect to detect modified gravitysignatures in BH spacetimes or, alternatively, constrainthe fundamental coupling constants of beyond-GR theo-ries to unprecedented precision.Possible extensions of this work include employing anunderlying GR template in which the complex amplitudes are constrained as functions of astrophysical parameters,while still allowing to fit the whole post-merger signal, asin Refs. [97, 131]. Such a template would decrease thenumber of free parameters and improve on the constraintspresented. Another possibility would be to consider differ-ent parametrizations, either using a non-polynomial formof QNM parameters or directly casting constraint on aneffective metric encompassing a wide range of beyond-GRBH solutions (possibly employing an eikonal approxima-tion). An application of this formalism to specific modifiedtheories of gravity where explicit predictions for the spec-tral deviation parameters are available is already ongoing.Finally, it would be interesting to consider an extensionof this formalism considering p as a free parameter (drop-ping the integer assumption), going beyond the currentclass of theories investigated. Such an addition will likelyintroduce non-trivial correlations in the parameter spaceand its feasibility should carefully be explored. ACKNOWLEDGMENTS
Software
LIGO-Virgo data are interfaced through
GWpy [132] and some of the post-processing steps arehandled through
PESummary [133], a sampler-agnostic python package providing interactive webpages to sim-plify results visualisation. Moreover,
PESummary meta-files are used to store the complete information (bothof the internal pyRing parameters and of the soft-ware environment) for the analysis to be completelyreproducible. The pyRing package will be pub-licly released shortly. Other open-software python packages, accessible through
PyPi , internally usedby pyRing are: corner, gwsurrogate, matplotlib,numba, numpy, seaborn, surfinBH [105, 106, 134–139].
Acknowledgments
The author would like to thankAlessandra Buonanno, Nathan Johnson-McDaniel andBangalore Sathyaprakash for useful discussions, EmanueleBerti, Walter Del Pozzo, Danny Laghi, Andrea Maselliand Paolo Pani for helpful comments on the manuscript.This work greatly benefited from discussions within the
Testing General Relativity
We would like to thank all of the essential workers whoput their health at risk during the COVID-19 pandemic,without whom we would not have been able to completethis work. TABLE II. Median and 90% credible levels of ringdown parameters deviation distributions; p indicates the coupling dimension, j the index of the expansion order, while M max the maximum expansion order considered, as defined in Eq. (2). No statisticallysignificant deviation from the predictions of BH dynamics in GR is found.( j, M max ) δω δω δτ δτ p = 0(0,0) − . +0 . − . − . +0 . − . . +0 . − . . +0 . − . (0,1) − . +0 . − . − . +0 . − . . +0 . − . . +0 . − . (0,2) − . +0 . − . − . +0 . − . . +0 . − . . +0 . − . (1,1) − . +0 . − . − . +0 . − . . +0 . − . . +0 . − . (1,2) 0 . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . (2,2) 0 . +0 . − . . +0 . − . . +0 . − . . +0 . − . p = 2(0,0) − . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . (0,1) − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . (0,2) − . +0 . − . − . +0 . − . . +0 . − . . +0 . − . (1,1) − . +0 . − . . +0 . − . − . +0 . − . . +0 . − . (1,2) − . +0 . − . − . +0 . − . . +0 . − . − . +0 . − . (2,2) 0 . +0 . − . . +0 . − . − . +0 . − . . +0 . − . p = 4(0,0) − . +0 . − . . +0 . − . − . +0 . − . . +0 . − . (0,1) − . +0 . − . . +0 . − . − . +0 . − . . +0 . − . (0,2) − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . (1,1) − . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . (1,2) − . +0 . − . . +0 . − . − . +0 . − . − . +0 . − . (2,2) 0 . +0 . − . . +0 . − . . +0 . − . . +0 . − . p = 6(0,0) − . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . (0,1) − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . (0,2) − . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . (1,1) − . +0 . − . − . +0 . − . − . +0 . − . . +0 . − . (1,2) − . +0 . − . . +0 . − . − . +0 . − . − . +0 . − . (2,2) 0 . +0 . − . . +0 . − . − . +0 . − . . +0 . − . M max = 2M max = 1 (0)221 M max = 0 (0)221 (0)220 (0)220 FIG. 1. Posterior probability distributions on the deviations from the lowest-order (non-spinning) expansion parameters ( δω, δτ )of the ( l = m = 2 , n = 0 ,
1) ringdown modes, for the p = 0 case, when combining all available LIGO-Virgo ringdown observations. TABLE III. Natural logarithm of the Bayes Factors comparing the hypotheses that all the BH spectra observed by LIGO-Virgoare described by BH dynamics in GR against an alternative dynamics; p indicates the coupling parameter dimension, M max themaximum expansion order considered, both defined in Eq. (2). The error on single-events Bayes factor is ± .
1, the error on thecombined Bayes factor is ± .
6. A large number of cases yield significantly large negative values, indicating that, with availabledata, the GR predictions are highly favoured for these parameters. In no case evidence for a statistically significant deviation isfound.M max δω δω δτ δτ p = 00 − . − . − . − . − . − . − . − . − . − . − . − . p = 20 − . − . − . − . − . − . − . − . − . − . − . − . p = 40 − . − . − . − . − . − . − . − . − . − . − . − . p = 60 − . − . − . − . − . − . − . − . − . − . − . − . (cid:96) in km units, when considering deviations in different ringdownparameters; p indicates the coupling parameter dimension, M max the maximum expansion order considered, both defined inEq. (2). For p=0 the coupling is reabsorbed in the deviation parameter. The bounds presented are either competitive orsuperior, depending on the theory considered, with those obtained through other experiments and improve upon previousbounds obtained with GW observations.M max δω δω δτ δτ p = 20 43 .
04 46 .
02 57 .
46 62 .
961 31 .
57 62 .
88 59 .
73 48 .
472 23 .
35 44 .
28 60 .
02 52 . p = 40 45 .
35 59 .
52 54 .
30 43 .
011 48 .
28 63 .
93 56 .
54 62 .
992 35 .
35 60 .
60 54 .
23 53 . p = 60 47 .
78 55 .
88 52 .
40 54 .
971 36 .
57 59 .
71 46 .
97 56 .
572 41 .
76 61 .
51 57 .
80 44 . p = 2 [km] p = 2 [km] M max = 2M max = 1 p = 2 [km] M max = 0 p = 2 [km] FIG. 2. Posterior probability distributions on the distance scale (cid:96) below which new physics effects can modify the BH emissionprocess for the p = 2 case, when combining all available LIGO-Virgo ringdown observations. p = 4 [km] p = 4 [km] M max = 2M max = 1 p = 4 [km] M max = 0 p = 4 [km] FIG. 3. Posterior probability distributions on the distance scale (cid:96) below which new physics effects can modify the BH emissionprocess for the p = 4 case, when combining all available LIGO-Virgo ringdown observations. p = 6 [km] p = 6 [km] M max = 2M max = 1 p = 6 [km] M max = 0 p = 6 [km] FIG. 4. Posterior probability distributions on the distance scale (cid:96) below which new physics effects can modify the BH emissionprocess for the p = 6 case, when combining all available LIGO-Virgo ringdown observations. [1] B. Abbott et al. (KAGRA, LIGO Scientific, VIRGO),Living Rev. Rel. , 3 (2018), 1304.0670.[2] R. Abbott et al. (LIGO Scientific, Virgo) (2020),2010.14529.[3] K. Akiyama et al. (Event Horizon Telescope), Astrophys.J. Lett. , L5 (2019), 1906.11242.[4] D. Psaltis, Gen. Rel. Grav. , 137 (2019), 1806.09740.[5] K. Akiyama et al. (Event Horizon Telescope), Astrophys.J. , L1 (2019), 1906.11238.[6] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 221101 (2016), [Erratum: Phys.Rev.Lett. 121,129902 (2018)], 1602.03841.[7] R. Brito, A. Buonanno, and V. Raymond, Phys. Rev. D98 , 084038 (2018), 1805.00293.[8] G. Carullo, W. Del Pozzo, and J. Veitch, Phys. Rev. D , 123029 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.99.123029 .[9] M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, andS. A. Teukolsky, Physical Review Letters (2019),ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.123.111102 .[10] D. Laghi, G. Carullo, J. Veitch, and W. Del Pozzo (2020),2011.03816.[11] B. Mashhoon, Phys. Rev. D , 290 (1985), URL https://link.aps.org/doi/10.1103/PhysRevD.31.290 .[12] V. Cardoso, E. Franzin, and P. Pani, Phys. Rev. Lett. , 171101 (2016), [Erratum: Phys.Rev.Lett. 117,089902 (2016)], 1602.07309.[13] D. Psaltis, Living Rev. Rel. , 9 (2008), 0806.1531.[14] V. Cardoso and P. Pani, Living Rev. Rel. , 4 (2019),1904.05363.[15] R. Abbott et al. (LIGO Scientific, Virgo) (2020),2010.14527.[16] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. ,074001 (2015), 1411.4547.[17] F. Acernese et al. (VIRGO), Class. Quant. Grav. ,024001 (2015), 1408.3978.[18] M. Punturo, M. Abernathy, F. Acernese, B. Allen,N. Andersson, K. Arun, F. Barone, B. Barr, M. Bar-suglia, M. Beker, et al., Classical and Quantum Grav-ity , 194002 (2010), URL https://doi.org/10.1088%2F0264-9381%2F27%2F19%2F194002 .[19] D. Reitze et al., Bull. Am. Astron. Soc. , 035 (2019),1907.04833.[20] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Ba-rausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bor-toluzzi, et al., Laser interferometer space antenna (2017),1702.00786.[21] G. Carullo et al., Phys. Rev. D , 104020 (2018),1805.04760.[22] S. Bhagwat, X. J. Forteza, P. Pani, and V. Ferrari, Phys.Rev. D , 044033 (2020), 1910.08708.[23] S. Bhagwat, M. Cabero, C. D. Capano, B. Krishnan,and D. A. Brown, Phys. Rev. D , 024023 (2020),1910.13203.[24] R. Tso, D. Gerosa, and Y. Chen, Phys. Rev. D ,124043 (2019), 1807.00075.[25] Z. Carson and K. Yagi, Class. Quant. Grav. , 02LT01(2020), 1905.13155.[26] Z. Carson and K. Yagi, Phys. Rev. D , 044047 (2020),1911.05258. [27] S. Datta, A. Gupta, S. Kastha, K. Arun, andB. Sathyaprakash (2020), 2006.12137.[28] A. Gupta, S. Datta, S. Kastha, S. Borhanian, K. Arun,and B. Sathyaprakash, Phys. Rev. Lett. , 201101(2020), 2005.09607.[29] A. Zimmerman, C.-J. Haster, and K. Chatziioannou,Phys. Rev. D , 124044 (2019), 1903.11008.[30] M. Isi, K. Chatziioannou, and W. M. Farr, Phys. Rev.Lett. , 121101 (2019), 1904.08011.[31] A. Maselli, P. Pani, L. Gualtieri, and E. Berti, Phys.Rev. D , 024043 (2020), 1910.12893.[32] E. Berti et al., Class. Quant. Grav. , 243001 (2015),1501.07274.[33] L. G. Collodel, B. Kleihaus, J. Kunz, and E. Berti, Class.Quant. Grav. , 075018 (2020), 1912.05382.[34] S. Alexander and N. Yunes, Phys. Rept. , 1 (2009),0907.2562.[35] D. Ayzenberg and N. Yunes, Phys. Rev. D ,069905 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.069905 .[36] N. Yunes and F. Pretorius, Phys. Rev. D ,084043 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.79.084043 .[37] A. Maselli, P. Pani, R. Cotesta, L. Gualtieri, V. Fer-rari, and L. Stella, The Astrophysical Journal , 25(2017), URL https://doi.org/10.3847%2F1538-4357%2Faa72e2 .[38] P. Pani, C. F. B. Macedo, L. C. B. Crispino, and V. Car-doso, Phys. Rev. D , 087501 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.84.087501 .[39] A. Maselli, P. Pani, L. Gualtieri, and V. Ferrari, Phys.Rev. D , 083014 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.92.083014 .[40] N. Yunes and L. C. Stein, Phys. Rev. D ,104002 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.104002 .[41] F.-L. Juli´e and E. Berti, Phys. Rev. D ,104061 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.104061 .[42] M. Okounkova, L. C. Stein, M. A. Scheel, and D. A.Hemberger, Phys. Rev. D , 044020 (2017), 1705.07924.[43] M. Okounkova, L. C. Stein, M. A. Scheel, and S. A.Teukolsky, Phys. Rev. D , 104026 (2019), 1906.08789.[44] M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel,and S. A. Teukolsky, Phys. Rev. D , 104016 (2020),1911.02588.[45] P. A. Cano, K. Fransen, and T. Hertog, Phys. Rev. D , 044047 (2020), 2005.03671.[46] S. Endlich, V. Gorbenko, J. Huang, and L. Senatore,JHEP , 122 (2017), 1704.01590.[47] V. Cardoso, M. Kimura, A. Maselli, and L. Sen-atore, Phys. Rev. Lett. , 251105 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.251105 .[48] N. Sennett, R. Brito, A. Buonanno, V. Gorbenko,and L. Senatore, Phys. Rev. D , 044056 (2020),1912.09917.[49] G. Franciolini, L. Hui, R. Penco, L. Santoni, andE. Trincherini, JHEP , 127 (2019), 1810.07706.[50] J. Noller, L. Santoni, E. Trincherini, and L. G. Trom-betta, Phys. Rev. D , 084049 (2020), 1911.11671. [51] K. Kokkotas, Nuov Cim B , 991–998 (1993).[52] E. Berti and K. D. Kokkotas, Phys. Rev. D , 124008(2005), gr-qc/0502065.[53] P. Pani, E. Berti, and L. Gualtieri, Phys. Rev. D , 064048 (2013), URL https://link.aps.org/doi/10.1103/PhysRevD.88.064048 .[54] Z. Mark, H. Yang, A. Zimmerman, and Y. Chen, Phys.Rev. D , 044025 (2015), 1409.5800.[55] M. Zilhao, V. Cardoso, C. Herdeiro, L. Lehner, andU. Sperhake, Phys. Rev. D , 124062 (2012), 1205.1063.[56] O. J. C. Dias, M. Godazgar, and J. E. Santos, Phys.Rev. Lett. , 151101 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.151101 .[57] V. Cardoso, C. F. B. Macedo, P. Pani, and V. Ferrari,jcap , 054 (2016), 1604.07845.[58] G. Bozzola and V. Paschalidis (2020), 2006.15764.[59] V. Cardoso, M. Kimura, A. Maselli, E. Berti, C. F. B.Macedo, and R. McManus, Phys. Rev. D , 104077(2019), 1901.01265.[60] R. McManus, E. Berti, C. F. Macedo, M. Kimura,A. Maselli, and V. Cardoso, Phys. Rev. D , 044061(2019), 1906.05155.[61] K. Yagi and L. C. Stein, Class. Quant. Grav. , 054001(2016), 1602.02413.[62] N. Yunes and X. Siemens, Living Rev. Rel. , 9 (2013),1304.3473.[63] N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D ,084002 (2016), 1603.08955.[64] T. Johannsen, Phys. Rev. D , 124017 (2013),1304.7786.[65] D. Psaltis, L. Medeiros, P. Christian, F. ¨Ozel,K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay,D. Ball, et al. (EHT Collaboration), Phys. Rev. Lett. , 141104 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.141104 .[66] S. H. V¨olkel, E. Barausse, N. Franchini, and A. E. Brod-erick (2020), 2011.06812.[67] S. E. Gralla, Phys. Rev. D , 024023 (2021),URL https://link.aps.org/doi/10.1103/PhysRevD.103.024023 .[68] D. Psaltis, C. Talbot, E. Payne, and I. Mandel (2020),2012.02117.[69] S. H. V¨olkel and E. Barausse, Phys. Rev. D , 084025(2020), 2007.02986.[70] K. Glampedakis, G. Pappas, H. O. Silva, and E. Berti,Phys. Rev. D , 064054 (2017), 1706.07658.[71] K. Glampedakis and H. O. Silva, Phys. Rev. D ,044040 (2019), 1906.05455.[72] H. O. Silva and K. Glampedakis, Phys. Rev. D ,044051 (2020), 1912.09286.[73] A. Boh´e et al., Phys. Rev. D , 044028 (2017),1611.03703.[74] P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D ,104063 (2012), 1207.3088.[75] M. Hannam, P. Schmidt, A. Boh´e, L. Haegel, S. Husa,F. Ohme, G. Pratten, and M. P¨urrer, Phys. Rev. Lett. , 151101 (2014), 1308.3271.[76] S. Khan, K. Chatziioannou, M. Hannam, and F. Ohme,Phys. Rev. D , 024059 (2019), 1809.10113.[77] S. Khan, S. Husa, M. Hannam, F. Ohme, M. P¨urrer,X. Jim´enez Forteza, and A. Boh´e, Phys. Rev. D ,044007 (2016), 1508.07253.[78] S. Khan, F. Ohme, K. Chatziioannou, and M. Hannam,Phys. Rev. D , 024056 (2020), 1911.06050. [79] W. Del Pozzo and J. Veitch, CPNest : an efficientpython parallelizable nested sampling algorithm , https://github.com/johnveitch/cpnest .[80] D. Sivia and J. Skilling, Data analysis: a Bayesiantutorial (Oxford Science Publications, 2006).[81] M. Breschi,
Hierarchical Bayesian Analysis of Post-Newtonian Gravitational Waveforms , https://etd.adm.unipi.it/theses/available/etd-09182018-183813/ .[82] F. Messina and A. Nagar, Phys. Rev. D , 124001(2017), [Erratum: Phys.Rev.D 96, 049907 (2017)],1703.08107.[83] F. Messina, R. Dudi, A. Nagar, and S. Bernuzzi, Phys.Rev. D , 124051 (2019), 1904.09558.[84] M. van de Meent and H. P. Pfeiffer, Phys. Rev. Lett. , 181101 (2020), 2006.12036.[85] E. T. Jaynes, Probability Theory: The Logic of Science (Cambridge University Press, 2003).[86] E. Berti, URL https://pages.jh.edu/~eberti2/ringdown/ .[87] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev.
D73 ,064030 (2006), gr-qc/0512160.[88] A. Nagar et al., Phys. Rev. D , 104052 (2018),1806.01772.[89] L. London and E. Fauchon-Jones, Class. Quant. Grav. , 235015 (2019), 1810.03550.[90] K. Ackley et al., Publ. Astron. Soc. Austral. , e047(2020), 2007.03128.[91] S. E. Perkins, N. Yunes, and E. Berti (2020), 2010.09010.[92] G. Van Rossum and F. L. Drake, Python 3 ReferenceManual (CreateSpace, Scotts Valley, CA, 2009), ISBN1441412697.[93] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vi-tale, B. Aylott, K. Blackburn, N. Christensen, M. Cough-lin, et al., Phys. Rev. D , 042003 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.042003 .[94] I. Kamaretsos, M. Hannam, and B. Sathyaprakash, Phys.Rev. Lett. , 141102 (2012), 1207.0399.[95] L. London, D. Shoemaker, and J. Healy, Phys. Rev. D , 124032 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.90.124032 .[96] L. T. London (2018), 1801.08208.[97] T. Damour and A. Nagar, Phys. Rev. D , 024054(2014), 1406.0401.[98] LIGO Scientific Collaboration, LIGO Algorithm Library- LALSuite , free software (GPL) (2018).[99] C. R. Galley and P. Schmidt (2016), 1611.07529.[100] P. Schmidt, I. W. Harry, and H. P. Pfeiffer (2017),1703.01076.[101] S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D. Selje-botn, and K. Smith, Computing in Science Engineering , 31 (2011), ISSN 1521-9615.[102] S. Maran`o, B. Edwards, G. Ferrari, andD. F¨ah, Bulletin of the Seismological Societyof America , 276 (2017), ISSN 0037-1106,https://pubs.geoscienceworld.org/bssa/article-pdf/107/1/276/2644797/276.pdf, URL https://doi.org/10.1785/0120160030 .[103] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland,T. Reddy, D. Cournapeau, E. Burovski, P. Peterson,W. Weckesser, J. Bright, et al., Nature Methods (2020).[104] X. Jim´enez-Forteza, D. Keitel, S. Husa, M. Han-nam, S. Khan, and M. P¨urrer, Phys. Rev. D ,064024 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.064024 . [105] V. Varma, D. Gerosa, L. C. Stein, F. H´ebert, andH. Zhang, Phys. Rev. Lett. , 011101 (2019),1809.09125.[106] V. Varma, S. E. Field, M. A. Scheel, J. Blackman,D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer,Phys. Rev. Research. , 033015 (2019), 1905.09300.[107] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. D , 104036 (2019), 1903.04467.[108] LIGO Scientific Collaboration and Virgo Collaboration, Data release for tests of general relativity with gwtc-2 (2020).[109] M. Giesler, M. Isi, M. A. Scheel, and S. Teukolsky, Phys.Rev. X9 , 041060 (2019), 1903.08284.[110] W. Del Pozzo, C. P. L. Berry, A. Ghosh, T. S. F. Haines,L. P. Singer, and A. Vecchio, Monthly Notices of theRoyal Astronomical Society (2018), ISSN 1365-2966,URL http://dx.doi.org/10.1093/mnras/sty1485 .[111] Z. Carson and K. Yagi, Class. Quant. Grav. , 215007(2020), 2002.08559.[112] H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, Phys.Rev. D , 064035 (2019), 1810.05177.[113] T. Akutsu et al. (KAGRA) (2020), 2005.05574.[114] M. Okounkova, Phys. Rev. D , 084046 (2020),2001.03571.[115] B. Shiralilou, T. Hinderer, S. Nissanke, N. Ortiz, andH. Witek (2020), 2012.09162.[116] K. Yagi, Phys. Rev. D , 081504 (2012), 1204.4524.[117] R. Nair, S. Perkins, H. O. Silva, and N. Yunes, Phys.Rev. Lett. , 191101 (2019), 1905.00870.[118] S. Tahura, K. Yagi, and Z. Carson, Phys. Rev. D ,104001 (2019), 1907.10059.[119] A. Maselli, L. Gualtieri, P. Pani, L. Stella, and V. Ferrari,Astrophys. J. , 115 (2015), 1412.3473.[120] P. Pani, E. Berti, V. Cardoso, and J. Read, Phys. Rev.D , 104035 (2011), 1109.0928.[121] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J.Lett. , L44 (2020), 2006.12611.[122] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M.Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V.Bilous, D. Chakrabarty, et al., The Astrophysical Jour-nal , L21 (2019), URL https://doi.org/10.3847/2041-8213/ab481c .[123] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov,Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K.Harding, W. C. G. Ho, J. M. Lattimer, et al., TheAstrophysical Journal , L24 (2019), URL https: //doi.org/10.3847/2041-8213/ab50c5 .[124] Z. Arzoumanian, K. Gendreau, C. Baker, T. Cazeau,P. Hestnes, J. Kellogg, S. Kenyon, R. Kozon, K.-C. Liu,S. Manthripragada, et al. (2014), vol. 9144, p. 914420.[125] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 161101 (2017), 1710.05832.[126] H. O. Silva, A. M. Holgado, A. C´ardenas-Avenda˜no, andN. Yunes (2020), 2004.01253.[127] M. Okounkova, W. M. Farr, M. Isi, and L. C. Stein(2021), 2101.11153.[128] C. M. Will, Living Rev. Rel. , 4 (2014), 1403.7377.[129] A. G. Suvorov and S. H. V¨olkel (2021), 2101.09697.[130] K. Chamberlain and N. Yunes, Phys. Rev. D , 084039(2017), 1704.08268.[131] X. J. Forteza, S. Bhagwat, P. Pani, and V. Ferrari (2020),2005.03260.[132] D. Macleod, A. L. Urban, S. Coughlin, T. Massinger,M. Pitkin, paulaltin, J. Areeda, E. Quintero, T. G. Bad-ger, L. Singer, et al., gwpy/gwpy: 1.0.1 (2020), URL https://doi.org/10.5281/zenodo.3598469 .[133] C. Hoy and V. Raymond (2020), 2006.06639.[134] D. Foreman-Mackey, A. Price-Whelan, W. Vousden,G. Ryan, M. Pitkin, V. Zabalza, jsheyl, A. Smith,G. Ashton, M. Smith, et al., dfm/corner.py: cor-ner.py v2.1.0.rc1 (2020), URL https://doi.org/10.5281/zenodo.3937526 .[135] S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, andM. Tiglio, Phys. Rev. X , 031006 (2014), 1308.3565.[136] J. D. Hunter, Computing in Science & Engineering ,90 (2007).[137] S. K. Lam, A. Pitrou, and S. Seibert, in Proceedings of theSecond Workshop on the LLVM Compiler Infrastructurein HPC (Association for Computing Machinery, NewYork, NY, USA, 2015), LLVM ’15, ISBN 9781450340052,URL https://doi.org/10.1145/2833157.2833162 .[138] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gom-mers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor,S. Berg, N. J. Smith, et al., Nature , 357 (2020),URL https://doi.org/10.1038/s41586-020-2649-2 .[139] M. Waskom, O. Botvinnik, D. O’Kane, P. Hob-son, S. Lukauskas, D. C. Gemperline, T. Augspurger,Y. Halchenko, J. B. Cole, J. Warmenhoven, et al., mwaskom/seaborn: v0.8.1 (september 2017) (2017), URL https://doi.org/10.5281/zenodo.883859https://doi.org/10.5281/zenodo.883859