The ringdown signal emitted during a binary black hole coalescence can be modeled as a linear superposition of the characteristic damped modes of the remnant black hole that get excited during the merger phase. While checking the consistency of the measured frequencies and damping times against the Kerr BH spectrum predicted by general relativity GR is a cornerstone of strongfield tests of gravity the consistency of measured excitation amplitudes and phases have been largely left unexplored. For a nonprecessing quasicircular binary black hole merger we find that GR predicts a narrow region in the space of mode amplitude ratio and phase difference independently of the spin of the binary components. Using this unexpected result we develop a new null test of strongfield gravity which demands that the measured amplitudes and phases of different ringdown modes should lie within this narrow region predicted by GR. We call this the amplitudephase consistency test and introduce a procedure for performing it using information from the ringdown signal. Lastly we apply this test to the GW190521 event using the multimodal ringdown parameters inferred by Capano et al. arXiv2105.05238. While ringdown measurements errors for this event are large we show that GW190521 is consistent with the amplitudephase consistency test. Our test is particularly well suited for accommodating multiple loud ringdown detections as those expected in the near future and can be used complementarily to standard blackhole spectroscopy as a proxy for modified gravity compact objects other than black holes binary precession and eccentricity.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Theory']
Novel Ringdown AmplitudePhase Consistency Test
{'A novel ringdown amplitude-phase consistency test': "Xisco Jim'enez Forteza 1 , 2 , Swetha Bhagwat 3 , 4 , Sumit Kumar 1 , 2 , Paolo Pani 3  \n1 Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Callinstraße 38, 30167 Hannover, Germany 2 Leibniz Universitat Hannover, 30167 Hannover, Germany 3 Dipartimento di Fisica, 'Sapienza' Universit'a di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 5, 00185, Roma, Italy and 4 Institute for Gravitational Wave Astronomy & School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  \nThe ringdown signal emitted during a binary black hole coalescence can be modeled as a linear superposition of the characteristic damped modes of the remnant black hole that get excited during the merger phase. While checking the consistency of the measured frequencies and damping times against the Kerr BH spectrum predicted by General Relativity (GR) is a cornerstone of strong-field tests of gravity, the consistency of measured excitation amplitudes and phases have been largely left unexplored. For a nonprecessing, quasi-circular binary black hole merger, we find that GR predicts a narrow region in the space of mode amplitude ratio and phase difference, independently of the spin of the binary components. Using this unexpected result, we develop a new null test of strongfield gravity which demands that the measured amplitudes and phases of different ringdown modes should lie within this narrow region predicted by GR. We call this the amplitude-phase consistency test and introduce a procedure for performing it using information from the ringdown signal. Lastly, we apply this test to the GW190521 event, using the multimodal ringdown parameters inferred by Capano et al. (2021) [1]. While ringdown measurements errors for this event are large, we show that GW190521 is consistent with the amplitude-phase consistency test. Our test is particularly well suited for accommodating multiple loud ringdown detections as those expected in the near future, and can be used complementarily to standard black-hole spectroscopy as a proxy for modified gravity, compact objects other than black holes, binary precession and eccentricity.  \nIntroduction. A binary black hole (BBH) ringdown is the gravitational-wave (GW) signal emitted as the remnant black hole (BH) formed during a BBH coalescence relaxes towards a stationary configuration [2-5]. The Kerr metric [6-9] in Einstein's general theory of relativity (GR) uniquely describes this final state. The ringdown phase is modelled as the evolution of perturbations (set up during the pre-merger stage) on the Kerr metric of the remnant BH. The GW signal emitted is wellapproximated as a linear superposition of countably infinite quasi-normal modes (QNMs), i.e., exponentially damped sinusoid modes with discrete characteristic complex frequencies, which are the eigenvalues of the radial and angular Teukolsky's equations [10-13]. Each mode is characterized by its frequency f lmn , damping time τ lmn , excitation amplitude A lmn , and phase φ lmn , where the integers ( l, m ) identify the angular dependence of the mode, whereas n = 0 , 1 , 2 , .. is the overtone number (see Eq. (2) for details).  \nWhile the frequencies and damping times are solely determined by the remnant's mass and spin, the perturbation condition setup prior to ringdown phase regulate the mode excitation, namely the amplitudes and phases. For a BBH coalescence, there is an intrinsic relation between the initial binary's parameters and the perturbation condition setup for ringdown, which determine the magnitude of A lmn and φ lmn . QNMamplitudes (i.e. excitation factors) account for spacetime's geometry during merger, providing the initial data for ringdown perturbations. GW ringdown models implicitly incorporate this into A lmn and φ lmn . We can estimate A lmn and φ lmn using numerical relativity (NR) simula-  \nns [14-16] corresponding to a set of BBH masses and spins. Thus, ringdown allows us to check two key predictions of GR in the strong-field regime - a) the consistency of measured QNM spectrum in the ringdown to the expected Kerr spectrum and, b) the compatibility of measured mode excitation factors with the prediction obtained from NR BBH coalescences in GR, i.e., consistency with (pre-)merger nonlinear dynamics.  \nWhile the former is the focus of a traditional BH spectroscopy, here we concentrate on the latter prospect and devise a novel test of GR called the ringdown amplitudephase consistency (APC) test . The APC test is based on the observation that, after a suitable normalization, only a narrow region in the mode amplitude-phase space is allowed for a BBH ringdown in GR. However, the BBH ringdown amplitudes and phases in modified GR [17-21] need not be constrained to lie in this region. We also expect a similar situation when the components of binary system are not Kerr BHs [22] (e.g. in neutron-star or more exotic boson-star [23-26] coalescences) within GR. Note that in both these scenarios, the remnant can still be a Kerr BH; so while it could pass the standard BH spectroscopy tests, the QNM amplitudes and phases can provide a way to distinguish the event from a GR BBH coalescence based on the nature of the merger.  \nAkey virtue of this test is that it does not require information from the inspiral phase other than the binary extrinsic parameters. It is therefore particularly well suited for massive BBHs, where the inspiral is short and the parameter estimation of the binary intrinsic parameters (e.g. the mass ratio and spins) is uncertain, jeopardizing the accuracy of inspiral-merger-ringdown (IMR) consis-  \ntency tests [27]. GW190521 is one such event [28]. It is also the only event that has been reported to show presence of measurable subdominant angular-mode parameters [1]. Furthermore, the parameter estimation for this event yields a primary component mass lying in the pair-instability supernova BH mass gap. This could be an indication for exotic alternatives [29-31] and the data do not exclude some of these possibilities. Thus, while our main scope is to devise the APC test in general, we find GW190521 a particularly interesting test-bed because if such exotic scenarios are possible, they can affect the QNM amplitudes and phases and manifest as a violation of this null test.  \nAmplitude-phase space of a BBH ringdown. Our aim here is to show that the mode amplitude-phase space corresponding to a BBH ringdown in GR is constranined to a narrow region. For this, we must first extract the amplitudes and phases of several modes by fitting the NR simulation as a function of BBH masses and spins. While we mostly focus on non-precessing, quasi-circular BBHs for simplicity, later on we shall also discuss the effect of spin misalignment and eccentricity.  \nWe use the dominant ( l = m = 2, n = 0) mode as a baseline and work with the intrinsic (independent of sky position, distance, polarization ψ ) amplitude ratio and phase difference defined as  \nA R lmn ≡ A lmn A 220 , δφ lmn ≡ m 2 φ 220 -φ lmn , (1)  \nrespectively [32]. δφ lmn is defined such that it removes the degeneracy between φ lmn and the coalescence phase ϕ (see Supplemental Material and Ref. [33]). The ringdown waveform can be analytically written down as  \nh + + i h × = A 220 ∑ lmn ( e -i m 2 φ 220 A R lmn e iδφ lmn S lmn ( ι, ϕ ) × e i 2 πf lmn t e -t/τ lmn ) , (2)  \nwhere S lmn ( ι, ϕ ) are the spin 2-weighted spheroidal harmonics. For simplicity, we approximate them to spherical harmonic functions [34] and as discussed later, this introduces a systematic error no larger than 1% for spinning remnants with a f glyph[lessorsimilar] 0 . 9 and for the modes considered in this work. For quasi-circular BBH mergers, A R lmn and δφ lmn are functions of the binary mass ratio q = m 1 /m 2 ≥ 1 and spins. We fit for A R lmn and δφ lmn mode-wise for a set of 142 nonprecessing and noneccentric NR simulations from the SXS catalog [35] (other catalogs [36, 37] are considered in the Supplemental Material). Our simulation set is the same used to calibrate the SEOBNRv4HM model [38] and spans q ∈ [1 , 10] and the z -component of the spins χ 1 , 2 ∈ [ -0 . 9 , 0 . 9]. We fit for the ( l, m, n ) ∈ { (3 , 3 , 0) , (2 , 1 , 0) , (4 , 4 , 0) } modes and provide ready-to-use analytical fits as a function of q and a post-Newtonian motivated effective spin χ pheno ( q, χ 1 , χ 2 ) whose explicit form depends on the mode under consideration [39]. In the nonspinning limit, these fits are consistent with those obtained in [38, 40-42], while they correct  \nFIG. 1. The regions of the ( A R lmn , δφ lmn ) plane allowed for BBHs in GR for various modes. They are estimated by fitting the parameters ( A R lmn , δφ lmn ) in Eq. (2) to 142 nonprecessing waveforms, covering the parameter space q ∈ [1 , 10] and χ 1 , 2 ∈ [ -0 . 9 , 0 . 9]. Given the parity of the odd/even m modes, we plot mod( δφ lmn , π ) for the (2 , 1 , 0) and (3 , 3 , 0) modes and mod( δφ lmn , 2 π ) for the (4 , 4 , 0) mode. The dashed blue curves within each region correspond to the nonspinning limit.  \n<!-- image -->  \nthe results given in [43] for the phase difference. We refer the reader to Supplemental Material for further details.  \nThe shaded bands in Fig. 1 marks off the region in the ( A R lmn , δφ lmn ) space allowed by GR BBH ringdowns as obtained by fitting the amplitudes and phases to our NR dataset using Eq. (2). The dashed curves correspond to χ 1 , 2 = 0, wherein A R lmn and δφ lmn parametrically depend only on q . The shaded region around the dashed curve quantifies the effects of non-zero χ 1 , 2 . Given the parity and polarization conventions of the odd/even m modes used in the SXS waveforms, we plot mod( δφ lmn , π ) for the (2 , 1 , 0) and (3 , 3 , 0) modes and mod( δφ lmn , 2 π ) for the (4 , 4 , 0) mode (see Supplemental Material). For the (3 , 3 , 0) mode, we find that the GR admissible region in the ( A R lmn , δφ lmn ) space is remarkably narrow. This happens because the effects of the progenitor spins are small for (3 , 3 , 0) mode. Indeed, δφ 330 ∈ [2 . 68 , π ] and A R 330 ∈ [0 , 0 . 42] for our entire dataset. The latter range of values observed for the amplitude ratio A R 330 is about a factor two larger than what observed in the inspiral regime [40]. This shows that higher harmonics are more excited during the highly-dynamical merger and ringdown regimes. However, the effect of non-zero χ 1 , 2 on A R lmn and δφ lmn is substantial for the (4 , 4 , 0) mode and (2 , 1 , 0) mode; this leads to a less-constricted GR-permissible region for these modes. The small bump of the (2,1,0) shaded area originates from restricting the  \namplitude A R 210 to be positive at small mass ratio and positive spin (see [44] and the Supplemental Material).  \nThe ringdown APC test. All nonprecessing, quasicircular BBH ringdown governed by GR must have A R lmn and δφ lmn within in the GR permissible region on the ( A R lmn , δφ lmn ) space i.e., measurements of A R lmn and δφ lmn must lie within the shaded bands in Fig. 1. We use this to devise our APC null test of strong gravity wherein we check if the measured mode amplitudes and phases in a ringdown signal lies in the narrow GR permissible region. Practically, one can check whether the posterior distributions of the estimated A R lmn and δφ lmn have significant support in the allowed region. Since all quasi-circular BBH mergers must satisfy this constraint, the test is naturally extendable to incorporate a population of observations.  \nNote that while implementing this test, it is crucial to account for the uncertainty in the ringdown start time with respect to the global peak time t p , i.e the time at which the strain amplitude | h ( t ) = Σ lmn h lmn | maximizes [43, 45]. Estimating the time t p lmn at which each | h lmn | mode peaks from NR fits, we can shift each mode by ∆ t = t p -t p 220 , with t p 220 < t p lm 0 [43, 46]. This induces a correction -  \nA R lmn ( t p ) = A R lmn ( t p 220 ) e ∆ t τ 220 -∆ t τ lmn . (3)  \nNotice that the τ lmn are comparable for all the modes studied here. For example, for a BH with a f ≈ 0 . 85 (consistent with GW190521 [1, 28]) τ 330 /τ 220 ∼ 0 . 986, τ 440 /τ 220 ∼ 0 . 974, and τ 210 /τ 220 ∼ 0 . 956. Here a conservative choice ∆ t ≈ 10 M translates to a ≈ (1 , 2 , 3)% correction for the amplitude of the (3 , 3 , 0), (4 , 4 , 0), (2 , 1 , 0) modes, respectively; this is well below the current statistical uncertainties [1]. This conservative choice corresponds to ∆ t ≈ 2( t p 330 -t p 220 ) ≈ ( t p 210 -t p 220 ) [43] and so sets a conservative upper bound on ∆ t .  \nA similar correction needs to be accounted for the intrinsic phase, φ lmn ( t p ) = φ lmn ( t p 220 ) + ω lmn ∆ t . Using Eq. (1), this translates to  \nδφ lmn ( t p ) = δφ lmn ( t p 220 ) + ( m 2 ω 220 -ω lmn ) ∆ t. (4)  \nSince ω lmn ≈ l 2 ω 220 for l = m modes (note that this is an exact result in the eikonal l = m glyph[greatermuch] 1 limit [12, 47]), the phase correction induced by ∆ t is ≈ 10% and ≈ 20% for the (3 , 3 , 0) and (4 , 4 , 0) mode respectively, when we assume ∆ t = 10 M (see Supplemental Material). This is a conservative choice for the (3 , 3 , 0) mode as ∆ t glyph[lessorsimilar] 5 M [43]; here we expect a systematic uncertainty no larger than 4%. Therefore, for the l = m modes, the NR phase fits can be compared to the measured posteriors inference as,  \nδφ lmn ( t p ) ≈ δφ lmn ( t p 220 ) ∀ l = m. (5)  \nHowever note that for the (2 , 1 , 0) mode the last term in Eq. (4) yields a non-negligible uncertainty.  \nFIG. 2. Ringdown APC test applied to GW190521. The amplitude ratio A R 330 and the phase difference δφ 330 obtained by various parameter estimations done on the data are compared to the GR allowed region in the ( A R 330 , δφ 330 ) space (shaded gray area). The solid and dashed contours respectively represent the 67% and 95% crediblity contours obtained by Bayesian parameter estimation. The orange and blue contours correspond to the posteriors obtained by fixing the GR QNMspectrum and allowing deviations from the Kerr QNMs, respectively. The black dot-dashed contour provides the 95% credible region obtained as in [1] but fixing the polarization to the maximum likelihood value ( ψ = 0) given by the IMR analysis [48]. The red dot-dashed contour provides the 95% credible region for the NR injection SXS:0258 consistent with GW190521 and with SNR ρ = 30. The cyan cross corresponds to the ( A R 330 , δφ 330 ) fit to the maximum likelihood waveform obtained using the NRSur7dq4 [49] waveform approximant, which includes precession and binary parameters consistent to GW190521 [50].  \n<!-- image -->  \n330  \nApplication on GW190521. We exemplify our test on GW190521, the only GW event with reported subdominant angular mode in the ringdown [1]. The total signal-to-noise ratio (SNR) of this event is ρ ≈ 14, of which ρ ≈ 12 comes solely from the ringdown phase. This is a consequence of large total source mass of this binary system ( M tot = 151 +29 -17 M glyph[circledot] [28] and the system is a convenient choice to demonstrate a proof-of-concept of the APC test.  \nIn Fig. 2, we first mark the GR-permissible region obtained by the NR fits as a grey band. Note that this region also accounts for the following uncertainties: the 1 σ deviations on the best fit results, the ∼ 1% and ∼ 4% uncertainties on A R 330 and δφ 330 caused by ambiguity in ringdown start time, and the errors due to approximating the spheroidal harmonics as spherical harmonics. Next, we project the 67% and 95% credible regions of the measured posterior distribution for A R 330 and δφ 330 from [1] on to the ( A R 330 , δφ 330 ) space. The posteriors obtained by assuming that ringdown has GR predicted Kerr QNM spectrum corresponds to the orange contour whereas the blue contour corresponds to the case where the QNM frequencies and damping times are allowed to vary freely. Interestingly, this more agnostic assumption  \ndoes not deteriorate the confidence region significantly. The black dot-dashed contour provides the 95% credible region obtained by fixing the polarization angle ψ to the maximum likelihood value estimated from the full IMR analysis in [48]. This is similar to fixing right ascension and declination as done in [1] (see [51-53] for a discussion on fixing these parameters in BH spectroscopy tests). Estimating the polarization angle independently helps to break the degeneracy between ψ and δφ lmn . ψ can be estimated from the inspiral-merger regime, while the intrinsic dependence of δφ lmn on q and χ pheno arises in the ringdown phase (see the Supplemental Material).  \nLastly, we inject the NR waveform SXS:0258 into Gaussian noise at SNR ρ = 30 for a 3-detector (LIGOHanford, LIGO-Livingston and Virgo) configuration and perform a parameters estimation using the PyCBC inference library [54]. This numerical waveform has parameters compatible with GW190521 [55]) but has been injected with twice the SNR of GW190521 to estimate the quality of the test achievable in the case of higher SNR events. The red dot-dashed contour in Fig. 2 denotes the 95% credible region obtained for this case. We perform the parameter estimation at t = t p +15ms (see the Supplemental Material). We note that, as expected, the confidence region shrinks and the test is significantly more accurate with higher SNR.  \nWe presented the main result of applying APC to GW19021 in Fig. 2; We find that the 1 σ credible interval obtained in [1] has a substantial support in the gray GRpermissible region marked on the ( A R lmn , δφ lmn ) space. Therefore, we conclude that the mode amplitude and phases measured in GW19021's ringdown are compatible with the GR BBH predictions and this event passes the APC test.  \nDiscussion. The APC test provides a novel strategy for testing GR using the ringdown mode excitations. NR waveforms of BBH mergers give accurate empirical relations between A R lmn and δφ lmn as functions of the binary's mass ratio and spins. We found that only a narrow strip in the ( A R lmn , δφ lmn ) space is admissible for ringdown modes of quasi-circular BBHs within GR. We build the APC test based on this feature and present a proof-of-concept implementation of this test on GW190521. We find that the 1 σ posterior distributions of A R lmn and δφ lmn obtained in [1] for this event has substantial support in the GR-permissible region, showing that GW190521 passes the APC test. Furthermore, we verify that the combined modelling uncertainties for the (3 , 3 , 0) mode are well below the statistical uncertainties of the current GW observations. Overall, for the (3 , 3 , 0) mode the total systematic errors accumulated from the fit, ringdown start time, and spherical-harmonic approximation may reach a value ∼ 13% for both A R lmn and δφ lmn (see Supplemental Material). For GW190521, the marginalized 1 σ statistical uncertainties on these quantities are ∼ 100% - much larger than the systematic deviations accumulated from our fit uncertainties. While this holds for any l = m mode, the phase fits of l = m  \nglyph[negationslash]  \nglyph[negationslash]  \nmodes are non-negligibly affected by the shift of the peak time; this makes l = m modes not optimal for the proposed implementation of the test. However, the situation might change for louder detections as those routinely expected in the third-generation era [56, 57], in which case systematic errors of the fit might limit the accuracy of the APC test, unless the quality of NR waveforms improves. On the other hand, the measurements of the polarization angle ψ (which is degenerate with the phase difference δφ lmn ) and of other binary's intrinsic parameters are expected to improve as more interferometers are added to the network, or through an electromagnetic counterpart, and will anyway improve with third-generation detectors, therefore allowing for a more accurate test.  \nWe have focused on quasi-circular binaries with aligned spins, although we can extend a similar concept to build a more generic test. In the Supplemental Material we show that the effect of eccentricity on the fits of δφ 330 ( A R 330 ) is non-negligible only when e glyph[greaterorsimilar] 0 . 3 ( e glyph[greaterorsimilar] 0 . 6). Therefore, the current implementation of the test is robust to mild eccentricities. By comparing the posterior distribution for A R 330 shown in Fig. 2 to the eccentricity fits, we have obtained a mild bound of e glyph[lessorsimilar] 0 . 9 at the 95% level on the eccentricity of GW19052. Note also that several works using the full IMR analysis on GW190521 have reported a moderately high effective precession spin parameter [58], χ p = 0 . 68 +0 . 25 -0 . 37 (although waveform systematics and prior choices significantly affect the posterior estimates of this event [50]). Interestingly, when we fit for δφ 330 and A R 330 corresponding to the maximum likelihood waveform including precession [49], we notice that the effect of precession for GW190521 is within the measurement errors for this event. However, the best fit is marginally outside the grey shaded area (corresponding to the non-precessing scenario) in Fig. 2. Because of the large statistical error, the systematic effect of neglecting the spin precession does not affect GW190521 significantly. Therefore, we can use GW190521 as a proof of concept for the APC test, and future louder events could be used to constrain the binary precession independently from ringdown measurements only. A detailed examination of the effects of precession in the ringdown is an involved problem and requires a dedicated study.  \nInterpretation and extensions. As with any nullhypothesis consistency test, its violation suggests a departure from the adopted baseline assumptions, and so a violation of the null-test could have various origins. We spell out the viable interpretation when an event does not pass the APC test - a) Most conservatively, it might be evidence for mis-modelling the signal e.g., presence of strong spin precession or large eccentricity in the BBH; b) It could be because the observed ringdown is not BH coalescence; note this does not preclude the remnant from being a standard Kerr BH. We expect the coalescence of massive neutron stars, boson stars [23-26], and other exotic compact objects [22] to produce QNMs consistent with Kerr BHs in GR. However, the QNM amplitudes and phases can be different from GR as the merger dy-  \nnamics could be modified; these will therefore fail the APC test while being consistent with a traditional Kerr BH spectroscopy; c) Finally and most radically, it could be because the underlying coalescence dynamics is not governed by GR. Disentangling these possibilities calls for a generalization of our fits to incorporate features like precession and eccentricity, louder ringdown detections, and detecting a population of them. For instance, if the violation of the test were due to not including eccentricity/precession in our fits, out of a population of ringdowns only a subgroup would be violating it. However, if GR dynamics were under question, there could be a ubiquitous violation of the test. In this context, although measurement errors are large, it is relevant that GW190521 passes the APC test. It would be interesting to assess whether this is in tension with alternative explanations for this event, e.g. a Proca star merger [29], by fitting the ringdown amplitudes and phases for Proca star merger waveforms and performing Bayesian model selection between the two hypotheses [29].  \nWhile we focused on the ringdown signal with prior knowledge of the binary's extrinsic parameters (estimated either from the IMR analysis or from another independent sky localization), a variant of this test would be to estimate the initial binaries parameters with ringdown and check for consistency with IMR analysis. In principle one could invert the A R lmn ( q, χ pheno ) and δφ lmn ( q, χ pheno ) relations to infer an estimate of the mass ratio (and spins) from the QNM excitations. However, owing to the mild dependence of δφ lmn on the binary parameters the quality of this test is expected to be rather poor. A more promising avenue is to neglect the phases and use only the amplitude ratios of several subleading QNMs. We discuss  \n- [1] C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, A. B. Nielsen, and B. Krishnan, (2021), arXiv:2105.05238 [gr-qc].\n- [2] S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A344 , 441 (1975).\n- [3] S. A. Teukolsky, Astrophys. J. 185 , 635 (1973).\n- [4] W. H. Press and S. A. Teukolsky, Astrophys. J. 185 , 649 (1973).\n- [5] S. A. Teukolsky and W. H. Press, Astrophys. J. 193 , 443 (1974).\n- [6] J. D. Bekenstein, Phys. Rev. D7 , 2333 (1973).\n- [7] B. Carter, Phys. Rev. Lett. 26 , 331 (1971).\n- [8] S. W. Hawking, Phys. Rev. Lett. 26 , 1344 (1971).\n- [9] S. W. Hawking and J. B. Hartle, Commun. Math. Phys. 27 , 283 (1972).\n- [10] E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26 , 163001 (2009), arXiv:0905.2975 [gr-qc].\n- [11] S. Chandrasekhar, The mathematical theory of black holes (Oxford Classic Texts in the Physical Sciences, 1985).\n- [12] V. Ferrari and B. Mashhoon, Phys. Rev. D 30 , 295 (1984).\n- [13] K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2 ,  \nthis in the Supplemental Material. This is interesting for GW190521-like systems where the short signal duration and low SNR in the pre-merger part leads to controversial and model-dependent inference on the binary parameters [48, 50, 59, 60]. Also, higher sensitivity at low frequency (as expected for third-generation detectors) will improve this test significantly (see also [61] for a conceptual framework in this direction). Overall, the APC test provides an excellent arena to complement standard BH spectroscopy tests in the strong-gravity regime, especially for the nextgeneration detectors.  \nAcknowledgments. We acknowledge the Max Planck Gesellschaft for support, and we are grateful to the Atlas cluster computing team at AEI Hannover for their help. The authors are specially thankful to Lionel London, Cecilio Garc'ıa-Quiros, and Juan CalderonBustillo for the invaluable discussions and further clarifications about the NR phase alignment and phase conventions. X. Jimenez is also thankful to P. Mourier for the useful discussions about the correspondence of the fit and parameter-estimation results. S.B. is supported by the UKRI Stephen Hawking Fellowship, grant ref. EP/W005727. P.P. acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement no. DarkGRA-757480. We also acknowledge support under the MIUR PRIN (Grant 2020KR4KN2 'String Theory as a bridge between Gauge Theories and Quantum Gravity') and FARE (GWNEXT, CUP: B84I20000100001, 2020KR4KN2) programmes, and from the Amaldi Research Center funded by the MIUR program 'Dipartimento di Eccellenza' (CUP: B81I18001170001).", '2 (1999), arXiv:gr-qc/9909058 [gr-qc].': "- [14] X. Jim'enez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Purrer, Phys. Rev. D95 , 064024 (2017), arXiv:1611.00332 [gr-qc].\n- [15] C. O. Lousto and J. Healy, Phys. Rev. D93 , 124074 (2016), arXiv:1601.05086 [gr-qc].\n- [16] F. Hofmann, E. Barausse, and L. Rezzolla, Astrophys. J. 825 , L19 (2016), arXiv:1605.01938 [gr-qc].\n- [17] M. Okounkova, Phys. Rev. D 102 , 084046 (2020), arXiv:2001.03571 [gr-qc].\n- [18] M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 101 , 104016 (2020), arXiv:1911.02588 [gr-qc].\n- [19] W. E. East and J. L. Ripley, Phys. Rev. Lett. 127 , 101102 (2021), arXiv:2105.08571 [gr-qc].\n- [20] M. Elley, H. O. Silva, H. Witek, and N. Yunes, (2022), arXiv:2205.06240 [gr-qc].\n- [21] H. Lim, G. Khanna, and S. A. Hughes, (2022), arXiv:2204.06007 [gr-qc].\n- [22] V. Cardoso and P. Pani, Living Rev. Rel. 22 , 4 (2019), arXiv:1904.05363 [gr-qc].\n- [23] C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, Phys. Rev. D 96 , 104058  \n(2017), arXiv:1710.09432 [gr-qc].  \n- [24] T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, Phys. Rev. D 99 , 044046 (2019), arXiv:1802.06733 [grqc].\n- [25] M. Bezares and C. Palenzuela, Class. Quant. Grav. 35 , 234002 (2018), arXiv:1808.10732 [gr-qc].\n- [26] M. Bezares, M. Boˇskovi'c, S. Liebling, C. Palenzuela, P. Pani, and E. Barausse, Phys. Rev. D 105 , 064067 (2022), arXiv:2201.06113 [gr-qc].\n- [27] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2112.06861 [gr-qc].\n- [28] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 125 , 101102 (2020), arXiv:2009.01075 [gr-qc].\n- [29] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forn'e, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, Phys. Rev. Lett. 126 , 081101 (2021), arXiv:2009.05376 [gr-qc].\n- [30] M. Fishbach and D. E. Holz, Astrophys. J. Lett. 904 , L26 (2020), arXiv:2009.05472 [astro-ph.HE].\n- [31] J. C. Bustillo, S. H. W. Leong, K. Chandra, B. McKernan, and K. E. S. Ford, (2021), arXiv:2112.12481 [gr-qc].\n- [32] J. Calder'on Bustillo, A. Boh'e, S. Husa, A. M. Sintes, M. Hannam, and M. Purrer, (2015), arXiv:1501.00918 [gr-qc].\n- [33] V. Baibhav, E. Berti, and V. Cardoso, Phys. Rev. D 101 , 084053 (2020), arXiv:2001.10011 [gr-qc].\n- [34] E. Berti and A. Klein, Phys. Rev. D90 , 064012 (2014), arXiv:1408.1860 [gr-qc].\n- [35] The SXS Collaboration, 'SXS Gravitational Waveform Database,' (2019).\n- [36] M. Campanelli, J. Healy, C. Lousto, and Y. Zlochower, 'CCRG@RIT Catalog of Numerical Simulations,' (2022).\n- [37] K. Jani, J. Healy, J. A. Clark, L. London, P. Laguna, and D. Shoemaker, 'Georgia Tech catalog of binary black hole simulations,' (2016).\n- [38] R. Cotesta, A. Buonanno, A. Boh'e, A. Taracchini, I. Hinder, and S. Ossokine, (2018), arXiv:1803.10701 [gr-qc].\n- [39] See Supplemental Material for a discussion on the NR fits and the phase and starting time conventions.\n- [40] S. Borhanian, K. G. Arun, H. P. Pfeiffer, and B. S. Sathyaprakash, (2019), arXiv:1901.08516 [gr-qc].\n- [41] L. T. London, Phys. Rev. D 102 , 084052 (2020), arXiv:1801.08208 [gr-qc].\n- [42] I. Ota and C. Chirenti, Phys. Rev. D 101 , 104005 (2020), arXiv:1911.00440 [gr-qc].\n- [43] X. Jim'enez Forteza, S. Bhagwat, P. Pani, and V. Ferrari, Phys. Rev. D 102 , 044053 (2020), arXiv:2005.03260 [grqc].\n- [44] C. Garc'ıa-Quir'os, M. Colleoni, S. Husa, H. Estell'es, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, Phys. Rev. D 102 , 064002 (2020), arXiv:2001.10914 [gr-qc].\n- [45] S. Bhagwat, X. J. Forteza, P. Pani, and V. Ferrari, Phys. Rev. D101 , 044033 (2020), arXiv:1910.08708 [gr-qc].\n- [46] H. Estell'es, S. Husa, M. Colleoni, D. Keitel, M. MateuLucena, C. Garc'ıa-Quir'os, A. Ramos-Buades, and A. Borchers, (2020), arXiv:2012.11923 [gr-qc].\n- [47] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79 , 064016 (2009), arXiv:0812.1806 [hep-th].\n- [48] A. H. Nitz and C. D. Capano, Astrophys. J. Lett. 907 , L9 (2021), arXiv:2010.12558 [astro-ph.HE].\n- [49] V. Varma, S. E. Field, M. A. Scheel, J. Blackman,  \n- D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. Research. 1 , 033015 (2019), arXiv:1905.09300 [gr-qc].\n- [50] H. Estell'es et al. , Astrophys. J. 924 , 79 (2022), arXiv:2105.06360 [gr-qc].\n- [51] R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, (2022), arXiv:2201.00822 [gr-qc].\n- [52] M. Isi and W. M. Farr, (2022), arXiv:2202.02941 [gr-qc].\n- [53] E. Finch and C. J. Moore, (2022), arXiv:2205.07809 [grqc].\n- [54] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, Publ. Astron. Soc. Pac. 131 , 024503 (2019), arXiv:1807.10312 [astro-ph.IM].\n- [55] The progenitor parameters of the simulation are q = 2 and χ pheno = 0 . 48, with dimensionless final spin a f = 0 . 84 and we have scaled the waveform to the total detector-frame mass consistent with the ringdown. The injected waveform includes the following modes (2 , ± 2) , (2 , ± 1) , (3 , ± 3) , (3 , ± 2) , (4 , ± 4) , (4 , ± 3). We fix the polarization to the maximum likelihood value ψ = 0, consistent with [48].\n- [56] M. Maggiore et al. , JCAP 03 , 050 (2020), arXiv:1912.02622 [astro-ph.CO].\n- [57] V. Kalogera et al. , (2021), arXiv:2111.06990 [gr-qc].\n- [58] P. Schmidt, F. Ohme, and M. Hannam, Phys. Rev. D91 , 024043 (2015), arXiv:1408.1810 [gr-qc].\n- [59] S. Kastha, C. D. Capano, J. Westerweck, M. Cabero, B. Krishnan, and A. B. Nielsen, (2021), arXiv:2111.13664 [gr-qc].\n- [60] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2111.03606 [gr-qc].\n- [61] S. Bhagwat, C. Pacilio, E. Barausse, and P. Pani, (2021), arXiv:2201.00023 [gr-qc].\n- [62] X. Jim'enez Forteza and S. Bhagwat, (2022), in preparation.\n- [63] D. Keitel, X. Jimenez-Forteza, et al. , Phys. Rev. D96 , 024006 (2017), arXiv:1612.09566 [gr-qc].\n- [64] E. Barausse, A. Buonanno, S. A. Hughes, G. Khanna, S. O'Sullivan, and Y. Pan, Phys. Rev. D85 , 024046 (2012), arXiv:1110.3081 [gr-qc].\n- [65] L. London, D. Shoemaker, and J. Healy, Phys. Rev. D90 , 124032 (2014), [Erratum: Phys. Rev.D94,no.6,069902(2016)], arXiv:1404.3197 [gr-qc].\n- [66] X. J. Forteza and P. Mourier, (2021), arXiv:2107.11829 [gr-qc].\n- [67] B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. 883 , 149 (2019), arXiv:1907.09384 [astro-ph.HE].", 'Appendix A: Fits for A R lmn and δφ lmn': 'For each of the 142 NR SXS waveforms in our dataset, we fit for A R lmn and δφ lmn for ( l, m, n ) = (3 , 3 , 0), (2 , 1 , 0), and (4 , 4 , 0) modes. The waveform dataset spans q ∈ [1 , 10] and χ 1 , 2 ∈ [ -0 . 9 , 0 . 9]. We fit for the amplitudes A lmn and phases φ lmn of each h lmn mode using the following ansatz,  \nh lmn = A lmn e -iφ lmn e i 2 πf lmn t e -t/τ lmn , (A1)  \nwhere the frequencies and damping times are fixed to the values predicted by GR, while A lmn and φ lmn are amplitudes and phases with same reference starting time t = 0.  \nTABLE I. Resolution and extrapolation error estimates for NR waveforms in the SXS catalog used for our fits.  \nHowever, the ringdown modes h lmn are extracted at a reference time t p 220 corresponding to the peak of the dominant (2 , 2 , 0) mode. To have the same reference starting time, we account for a time shift ∆ t = t p -t p 220 , where t p glyph[lessorsimilar] 10 M is the time at which the strain signal peaks (global peak). As discussed in the main text, this ansatz does not fully match Eq. (2) due to the peak time ambiguity. It adds an uncertainty on the phase values δφ lmn for the l = m modes but this is typically small.  \nIn the rest of this appendix we provide a detailed discussion about the NR dataset used and the accuracy of the fits.', '1. Error estimate for NR waveforms': 'NR waveforms contain two sources of uncertainties that are of interest to our study - i) resolution uncertainties which are produced by the finiteness of the numerical grid; and ii) extrapolation errors which are produced from computing the data at a finite radii and extrapolating them to future null infinity. To quantify resolution errors, for each simulation we compute the mismatch (as defined, e.g., in Eq. (1) of [45]) between waveforms at the two highest resolutions of the simulations in the SXS catalog. To estimate the extrapolation errors, we use the highest-resolution waveform and compute the mismatch between the waveform extrapolated with second or third polynomial order. We finally compute the distribution of the mismatch across our dataset. In Table I, we give the order of magnitude of the mismatch at the median value for the distribution for each mode. We see that extrapolation errors are negligible relative to resolution errors for all modes, and the maximum mismatch is at most O (10 -3 ). A more detailed study on NR systematics in the context of ringdown amplitude and phase fits will be presented in a companion paper [62]. Next, we will also compare the results of our fits obtained using different NR BBH waveforms catalogs.', '2. Fits for A R lmn': "For a quasi-circular nonprecessing binary, the mode amplitude ratios and phase differences generically depend  \non all intrinsic binary parameters such as mass ratio q and spin amplitudes χ 1 , 2 . However, it turns out that A R lmn also depends more strongly on certain combinations of the spins, similarly to the effective spin parameters used in PN waveform modeling [40]. We obtain analytical (approximate) relations for A R lmn as a function of these BBH parameterizations using the ansatz [40]  \nA R lm 0 = a 0 δ + a 1 δ 2 + a 2 χ, ∀ odd modes (A2)  \nA R lm 0 = a 0 (1 -3 η ) + a 1 (1 -3 η ) 2 + a 2 (1 -3 η ) 3 (A3) a 3 + χ s , ∀ even modes  \nwhere the { (3 , 3 , 0) , (2 , 1 , 0) } ∈ odd modes and (4 , 4 , 0) ∈ even modes, η = q/ (1 + q ) 2 , δ = √ 1 -4 η , whereas  \nχ = χ a + χ s √ 1 -4 η 2 (A4)  \nis the favoured combination of the spin parameters and, for (anti)aligned spins, χ s,a = ( m 1 χ 1 ± m 2 χ 2 ) / ( m 1 + m 2 ) with m 1 , 2 being the progenitor BH masses. This ansatz automatically enforces A R lmn ( q ) → 0 for q → 1 for all odd modes in the χ 1 = χ 2 limit, which arises from the binary's symmetry under m 1 ↔ m 2 . Note that the above ansatz differs from the ones we have used for the nonspinning fits presented in Ref. [43].  \nWe fit the data in two hierarchical steps following [14]: we first fit the nonspinning waveforms using Eqs. (A2) with χ 1 , 2 = 0. Then, we fit for the spinning BBH waveforms, keeping the values of the coefficients obtained from the nonspinning fit to constrain the final result in the nonspinning limit. This improves the accuracy of the fit in the region of the parameter space where the NR simulations are known to be more accurate [14, 63]. We get the following analytical ready-to-use fits:  \nA R 330 = 0 . 572 √ 1 -4 η -0 . 144(1 -4 η ) + 0 . 035 χ, (A5)  \nA R 210 = ∣ ∣ ∣ 0 . 328 √ 1 -4 η +0 . 115(1 -4 η ) -0 . 414 χ ∣ ∣ ∣ , (A6)  \nA R 440 = 0 . 251 ( 1 + 59 . 773 η 3 -16 . 307 η 2 -3 η ) -0 . 011 χ s . (A7)  \nWe choose to set the amplitudes to be positive by shifting phase by a factor π i.e., δφ lmn → δφ lmn + π , for those cases where the fit provides a negative amplitude. For the (2 , 1 , 0) mode, we observe that A R 210 tends to negative values at low mass ratio q glyph[lessorsimilar] 2 and high spin χ ∼ 0 . 4. In this case we added the absolute value to the ansatz in order to keep the A R lmn > 0 convention . Note that the fits recover the test particle limit [43, 64] as q →∞ , and the nonspinning regime as χ 1 , 2 → 0. Indeed, in the nonspinning limit we verified that our fit agrees reasonably well with the fits in [41-43, 65].  \nIn the top panel of Fig. 3, we present the amplitude ratio A R lmn as a function of q for all 142 simulations. The green diamonds, red crosses, and blue dots  \nFIG. 3. Top panel: amplitude ratio A R lmn in terms of the mass ratio q obtained from fitting 142 waveforms for the (3,3,0), (2,1,0), and (4,4,0) modes. The spread of the points on the vertical direction quantifies the effects of the binary spin parameter, χ pheno = χ, χ s for the odd/even modes, respectively. The solid line shows the fit for the nonspinning case. Bottom panel: Normalised residual distributions for the three modes. Notice that effects are relatively small for the (3 , 3 , 0) and (2 , 1 , 0) compared to their typical amplitudes. Further analysis about the errors is provided in a followup paper [62].  \n<!-- image -->  \ncorrespond to A R 330 , A R 210 , and A R 440 , respectively. The solid lines denotes A R lmn for the nonspinning BBHs, i.e., A R lmn = A R lmn ( q, χ 1 , 2 = 0). For all modes considere here, A R lmn increases with the mass ratio, i.e. for more asymmetric binaries. Spins effects are small for A R 330 and A R 440 , leading to a small scatter around the solid lines. This also suggests that spin effects are generically small for these modes, even when accounting for spin misalignment.  \nIn order to show the spin and mass-ratio dependence more clearly, in the top panel of Fig. 4 we present a contour plot of A R lmn as a function of q and χ, χ s . We observe that, at variance with the other modes, A R 210 depends significantly on the binary spins. Finally, we quantify the goodness of the fit using residuals. In the bottom panel of Fig. 3, we present the normalized fit residual distributions, where residuals are the difference between the amplitude/phase obtained by the fit and that obtained from the raw NR data, i.e., Residuals =  \n( A R lmn, NR -data -A R lmn, fit ) or ( δφ lmn, NR -data -δφ lmn, fit ) for each simulation used. We see that the residuals are centered around zero with a small spread. Compared to their absolute value, these errors are small for the (3,3,0) and (2,1,0) modes and modest for the (4,4,0) mode.", 'a. Disentangling NR conventions from physical phase': 'The NR waveforms from different catalogs adopt different conventions for phases, so one must appropriately account for this to combine/compare phases across several NR waveform catalogs. These conventions arise from the choice of - a) the tetrad adopted to extract the NR waveform which adds an overall polarization angle ψ 0 ; and b) from rotations of the BH orbital plane by an angle ϕ 0 . Two waveforms (say A and B) from different catalogs with the same physical intrinsic parameters and aligned in time are related by [32, 44],  \nh A lmn ( t ) = e ι ( ψ 0 + mϕ 0 ) h B lmn ( t ) . (A8)  \nThe polarization angle across the NR codes is either ψ 0 = 0 or ψ 0 = π to preserve the rotatingcounterrotating mode symmetry, h lm = ( -1) l h ∗ l -m , for circularly-polarized, nonprecessing waveforms 1 [44]. We are interested only in the physical contribution to the phase, δ lmn ( glyph[vector] λ ), which depends only on the binary intrinsic parameters glyph[vector] λ . A generic ringdown phase φ lmn results from the sum of the three contributions [32, 41, 65]  \nφ lmn = δ lmn ( glyph[vector] λ ) + mϕ 0 + { 0 , π } . (A9)  \nUnlike the physical phase δ lmn ( glyph[vector] λ ) that depends on the BBH parameters, the extrinsic phase terms ψ 0 and ϕ 0 may vary across sets of NR simulations and codes. However, the dependence on ϕ 0 is eliminated out if we fit for the following quantity  \nδφ lmn := m 2 φ 22 n -φ lmn = m 2 δ 22 n ( glyph[vector] λ ) -δ lmn ( glyph[vector] λ ) (A10) + { 0 , ( m 2 -1 ) π } .  \nThe phase difference δφ lmn depends only on the intrinsic binary parameters glyph[vector] λ and on a global phase factor which is either zero or ( m 2 -1) π depending on the simulation. For instance, we can identify the convention used in a NR waveform by knowing that, in the low-frequency  \n<!-- image -->  \n<!-- image -->  \nχ  \n<!-- image -->  \nχ  \nϕ  \nχ  \n<!-- image -->  \n<!-- image -->  \ns  \nχ  \ns  \nχ  \nFIG. 4. Contour plots showing the amplitude ratio A R lmn (top) and phase difference δφ lmn (bottom) as functions of q and χ pheno = χ, χ s , χ φ (depending on the panels) for the ( l, m, n ) = (3 , 3 , 0) (left), (2,1,0) (middle), and (4,4,0) (right) mode. As χ pheno increases, the magnitude of A R 330 increases, while the A R 210 decreases and A R 440 mildly decreases. The small region around q ∼ 1 -2 and χ > 0 in which A R 210 increases with χ is forced by fitting | A R 210 | , which induces a sign flip on the trend of A R 210 ( χ ). In the bottom panels, δφ lmn decreases linearly with χ φ for the (3 , 3 , 0) mode while it increases as χ 2 φ for the (2 , 1 , 0) mode. The value of the phase difference δφ lmn decreases with χ pheno for the (3 , 3 , 0) and (4 , 4 , 0) modes while it increases for the (2 , 1 , 0) mode.  \n<!-- image -->  \ninspiral regime, the phase difference between the dominant (2 , 2 , 0) mode and a higher ( l, m, n ) mode satisfies mod( mφ 220 -2 φ lmn , 2 π ) = 0 , 2 π , for both the even and the odd modes (see Appendix D of [44] and [46]). For the SXS data, we have checked that this value is consistent with ψ 0 = 0. The RIT and Maya [36, 37] waveform catalogs (considered later on) adopt instead the ψ 0 = π convention [35].', 'b. Phase fits': 'Similar to the case for the amplitude ratio, we produce ready-to-use fits for δφ lmn = m 2 φ 22 n -φ lmn as a function of the BBH parameters. We use the following ansatz informed by the leading order PN expressions on δφ lmn [38]  \nδφ 330 = b 0 δ + b 1 χ φ + c 0 , (A11)  \nδφ 210 = b 0 δ + b 1 χ φ + b 2 χ 2 φ + c 0 , (A12)  \nδφ 440 = a 0 η d 0 + b 1 χ s , (A13)  \nwhere χ φ ≡ 1 2 ( χ a + √ 1 -4 η 1 -2 η χ s ) is another phenomenological fit parameter.  \nWe follow a similar hierarchical fitting procedure as the one previously described for the amplitude-ratio fits. The fits obtained for the δφ lmn read  \nδφ 330 = 2 . 759 + 0 . 406 √ 1 -4 η -0 . 055 χ φ (A14) δφ 210 = 0 . 401 + 0 . 286 √ 1 -4 η +0 . 402 χ φ +0 . 652 χ 2 φ δφ 440 = 4245 . 459 η 5 . 646 -0 . 365 χ s .  \nIn Fig. 5 we present δφ lmn (as defined in Eq. (A10)) obtained at t r -t p 22 = 10 M . Again the solid lines correspond to the nonspinning BBH case, while the spread of points around the lines quantifies the spin dependence of the result. As q increases, we observe that δφ 330 and δφ 210 mildly increase while δφ 440 rapidly vanishes. We stress that for the SXS catalog the convention is such that mod( δφ 330 , π ) = π in the low-frequency PN regime. For q ≈ 1, we see that δφ lmn develops a dependence on the BBH spin parameters, which is reduced as q increases. This is consistent with [38], where mod( δφ 330 , π ) is evaluated at t r = t p 22 . We again evaluate the fit residuals  \nfor each NR waveform; the distribution of the residuals is shown in the bottom panel of Fig. 5. Similar to Fig. 3, the largest errors are obtained for δφ 440 . Finally, in the bottom panel of Fig 4 we show δφ lmn as a function of q and χ, χ s , χ φ for the (3 , 3 , 0) (left), (2 , 1 , 0) (middle) and (4 , 4 , 0) (right) modes. As in the case for the amplitude ratio, the phase δφ 330 is only mildly affected by the spin, whereas a stronger spin dependence occurs for the (2,1,0) and (4 , 4 , 0) modes. Requiring the amplitude of the (2,1,0) mode to be non-negative makes it a nonsmooth function in the spin parameter space χ [44] and induces the sudden increasing of | A R 210 | at q ∼ 1 -2 and χ > 0. This originates the small bi-modal bump observed in the shaded area of the (2,1,0) mode of Fig. 1.  \nFIG. 5. Same as Fig. 3 but for the phase difference δφ lmn .  \n<!-- image -->', '4. Comparison of the fits using other NR catalogs': "The RIT and Maya public catalogs [36, 37] provide a large set of NR simulations that can be used for fitting and testing. In particular, the public data provided for the RIT catalog is tested to be in the convergent regime resolution errors shall dominate - and shows good global IMR agreement with the SXS data for all the modes up to l = 5. We have calibrated the amplitudes and phases  \nusing data from the SXS catalog since they provides data at different resolutions and extrapolation levels. We use these for the error estimates previously presented in this study. However, we can use NR data from the RIT and Maya catalogs to benchmark our results. This test is particularly useful since the two families of codes use significantly different numerical schemes to solve Einstein's equations for a BBH 2 . We see that both A R lmn and δφ lmn computed from different codes are affected dominantly by the finite extraction and extrapolation effects - adding up the differences between them to a few percent [63]. We observe that the systematic errors (i.e., a shift of the median value of the distribution with respect to zero), are below the fitting errors characterized by the width of the distributions. This observation holds true for all the modes considered here.  \nTo compare the data from the different catalogs, we first revisit the various conventions used in each of them - for example, the ψ 0 = 0 , π rotational factors arising from different tetrad choices in the simulations. For the RIT and Maya catalogues, we need to replace δφ lmn → -δφ lmn -m 2 π , This factor results from the different tetrad conventions used in these codes. The additional minus sign comes from the reversed definition of the imaginary component of the h lm modes between the SXS dataset and the RIT and Maya datasets, which just implies that φ SXS lmn ↔-φ RIT , Maya lmn [35].  \nIn Fig. 6, we compare amplitude ratio and phase difference of the (3 , 3 , 0) mode obtained using SXS simulation to that given by the RIT and Maya wavefroms. For the amplitude ratio (left panels of Fig. 6), we find a good match between the Maya and SXS data while the results show an offset of about A R 330 ∼ 0 . 05 when using the waveforms from the RIT catalog. This systematic offset is still smaller than the value of the uncertainty observed in GW190521, δA R 330 ∼ 0 . 1, but will become important for louder events in the future. Regarding the phase difference (right panels of Fig. 6), despite we find a resonably good agreement between all the three codes, the RIT code shows a slightly larger tails. The standard deviation obtained using the normalized residuals distributions is ∼ 0 . 1 rad for the three codes. This uncertaininty does not effect our ability to perfrom APC test for current or near-future GW observations.", '5. Identifying fit outliers in the data': 'We have observed that the BH spins have a minor impact on the values of A R lmn and δφ lmn compared to the effect of the mass ratio q . This allow us to use an ansatz linear in χ pheno for all the modes, except for δφ 210 , in  \nx  \nFIG. 6. Top panel: Amplitude ratio (left) and phase difference (right) in terms of the mass ratio q and the phenomenological spin parameter χ pheno obtained by fitting 142 waveforms for the (3,3,0) modes. Data from the SXS, RIT, and Maya catalogs are used here and the denoted by green, pink and blue respectively. The solid line joining the points represents nonspinning simulations using just the SXS catalog. Bottom panel: Normalized residual distributions obtained for the three catalogs for the amplitude ratio (left) and phase difference (right).  \n<!-- image -->  \nFIG. 7. Amplitude ratio fits A R 330 for q = 3 in terms of the effective spin χ . SXS data are denoted by green diamonds while the red dots denote the results obtained from fitting the data with the PN motivated ansatz (A2) linear in χ . One can easily identify an outlier point at A R 330 = 0 . 223 and χ = 0 . 850 (blue cross). The solid orange curve provides the results of the linear fit after removing the outlier point. In particular, this corresponds to the SXS simulation SXS:0293.  \n<!-- image -->  \nwhich the effects of χ pheno are found to be larger. While this property makes it easier to model A R lmn and δφ lmn in terms of the physical parameters, we find that presenced of outliers decrease the quality of the spin-dependent fits. The outliers seem to be dominated by numerical noise and we remove them from our calibration dataset by the following procedure. We first look for highly significant outliers in both A R lmn and δφ lmn by testing the spin fit χ pheno at each mass ratio q . More specifically, we select the data at mass ratio q = 1 , 2 , 3 , 6 , 7 , 8 to perform a bootstrap analysis over the χ pheno axis. This is achieved by computing the spinning fit at each q , that in general will contain N points, for a data set of N -1 points. We iterate the fit for the N points at each q , and we compute the value of the standard deviation σ ( q i , χ i pheno ). Then, we select the median value of ¯ σ and we discard all points beyond a conservative deviation of 4¯ σ . In Fig. 7 we show an example of this procedure applied to the amplitude ratio A R 330 . We see that the blue cross placed at the low-right corner is easily identified by this algorithm. However, it is important to keep a conservative criterion to avoid rejecting systematically false outliers.  \n330  \n440  \nFIG. 8. Difference ω lmn -l/ 2 ω 220 as a function of the remnant spin. The ω lmn ∼ l/ 2 ω 220 approximation is used in the calibration of δφ lmn for the (3,3,0) and (4,4,0) modes. For the (3,3,0) mode and a conservative choice ∆ t ∼ 10 M this adds an uncertainty of about 0 . 4 rad. However, such approximation worsens for the (4 , 4 , 0) mode.  \n<!-- image -->', '6. On the ∆ t dependence for the phase difference': 'In order to refer the NR fits to the parameter estimation to a consistent initial time we need to shift each mode by ∆ t = t p -t p 220 , where t p is the global peak time of the signal and t p 220 is an arbitrary reference time taken to be the peak time of | h 220 | . This time shift introduces a dephasing ∆ φ lmn = ω lmn ∆ t and an extra term in the phase difference δφ lmn ( t p ) (see Eq. 4)  \n( m 2 ω 220 -ω lmn ) ∆ t . (A15)  \nFortunately, as can be shown analytically in the geodesics approximation valid in the eikonal limit l = m glyph[greatermuch] 1 [12, 47], for the fundamental ( n = 0) modes the following approximation holds: ω lm 0 ∼ l/ 2 ω 220 . Thus, for the l = m modes the extra terms in Eq. (A15) is small. In Fig. 8, we show ω lmn -l/ 2 ω 220 for the (3 , 3 , 0) and (4 , 4 , 0) modes as a function of the remnant spin, showing that the difference is in the range ≈ [0 . 04 , 0 . 08] for any spin. This yields an ambiguity in δφ lmn approximately of ( ω 330 ∼ 3 / 2 ω 220 )∆ t ∼ 0 . 04∆ t/M and ( ω 440 ∼ 4 / 2 ω 220 )∆ t ∼ 0 . 07∆ t/M . Choosing a very conservative error estimate on ∆ t , namely ∆ t = 2( t p 33 -t p 22 ) = 10 M , we conclude that ignoring the dephasing introduced by ∆ t will at most introduce an overall uncertainty δφ 330 ∼ 0 . 4 rad and δφ 330 ∼ 0 . 8 rad. Note that, since ( t p 33 -t p 22 ) ∼ 5 M for all the NR simulations explored [43, 46], in practice we expect the uncertainty on δφ 330 to be typically half of this conservative estimate.  \nFIG. 9. Posteriors distributions for δφ 330 obtained from our injected signal SXS:0258 and for set of starting times t r ∈ [0 . 005 , 0 . 019] ms. The black dashed line stands for the fitting value shifted back ∆ t (3 / 2 ω 220 -ω 330 ) to account for the ∆ t = t p -t p 220 difference. For this simulation and for an inclination ι = 2 . 4, we obtain ∆ t = 0 . 013ms.  \n<!-- image -->', 'a. On the effect of the reference time for the parameter estimation': 'We require that the reference time to start the parameter estimation on our GW190521 injections satisfies t r ≥ t p 33 . This fixes a time at which the (3 , 3 , 0) mode is excited, which may source some systematic errors on δφ lmn . For GW190521 with an inclination of ι = 2 . 4, we obtain t r ≥ 0 . 005 ms. In Fig. 9 we show the posterior distributions obtained for δφ 330 for a set of reference times [0 . 005 , 0 . 019] ms. The dashed black line corresponds to the fit value, shifted to t p 22 . As expected, the posterior distributions are consistent to each other for all the times t r considered here.', '7. The effect of the eccentricity': 'The coalescence of eccentric binaries can in principle modify the initial perturbation conditions setup for the ringdown phase, and hence the QNM amplitudes and phases. Thus, the amplitude ratio A R lmn and the phase difference δφ lmn also depends on the eccentricity e . In Fig. 10 we examine the impact of the eccentricity on the values of A R 330 and δφ 330 for a set of mass ratios q = 2 , 3 , 4 - where both the standard BH spectroscopy test and our APC test are promising [1, 66]. To obtain these estimates, we use the data from the RIT catalog [36]. Notice that the values on δφ 330 are significantly modified only at relatively large values of the eccentricity with e glyph[greaterorsimilar] 0 . 3, while this value raises up to e glyph[greaterorsimilar] 0 . 6 for A R 330 . These values are still above the upper limit threshold for e ∼ 0 . 1 obtained from the search of eccentric BBHs during the first and second LIGO observation runs [67]. Thus eccentricities are not relevant for the most of the events  \nFIG. 10. The impact of the eccentricity e on the estimates of the amplitude ratio A R 330 (left) and phase difference δφ 330 (right) for a set of mass ratios q = 2 , 3 , 4 and χ pheno = 0 using the RIT catalog [36]. The shaded area represents the 95% credible intervals obtained from the posteriors shown in Fig. 2. Notice that, from the posterior distribution on A R 330 and the NR fits of eccentric data, we can constrain the eccentricity e glyph[lessorsimilar] 0 . 9 at the 95% confidence level. Moreover, the effects of the eccentricity will become relevant only when e glyph[greaterorsimilar] 0 . 3 for the δφ 330 and when e glyph[greaterorsimilar] 0 . 6 for A R 330 and for GW events louder than GW190521.  \n<!-- image -->  \nobserved by current ground based GW observatories for perfroming APC test. Finally, the shaded area provides the 95% credible intervals obtained from the posteriors shown in Fig. 2. Using these loose constraints, we obtain bounds on the eccentricity of GW190521 as e glyph[lessorsimilar] 0 . 9 at the 95% confidence level.', 'Appendix B: Consistency between mode-excitations and BBH mass ratio in GW190521': 'In this work we have proposed a new test of GR called the APC test that we demonstrate on GW190521 using the ringdown alone (and possibly prior knowledge of the binary extrinsic parameters [33]). In this section, we highlight another possible null test of GR that makes use of the amplitude ratio and its relation to the BBH mass ratio, and therefore requires the entire inspiral-mergerringdown (IMR) signal. The basic concept here is that one could estimate the mass ratio q from the ringdown by inverting the A R lmn = A R lmn ( q ) relation and then check whether the inferred value is consistent with q measured independently from the full IMR signal 3 . This is complementary to the standard IMR consistency tests performed by the LVK Collaboration [27] and is based on an idea similar to that used to design the merger-ringdown test for the BBH population presented in Ref. [61].  \nOur fits provide empirically 2-to-2 maps ( A R lmn , δφ lmn ) → ( q, χ pheno ) which can be inverted to obtain q . Note that the sensitivity of this test toward the measurement of the spin is limited since A R 330 and A R 440 have a rather mild dependence on the spins (see Fig. 3). However, even when one neglects the spin dependence, one can map A R to q to a good approximation. For instance, for q = 2 we get A R 330 = 0 . 14 , 0 . 19 for χ 1 , 2 = ± 0 . 85; the spin dependence is subleading. In principle, this kind of test could also be designed with δφ lmn . However, in practice the dependence of δφ lmn on the BBH parameters is weak, and one expects a much larger errors on the inferred BBH parameters making the test less constraining.  \nWe scrutinize GW190521 for consistency between mode-excitation and BBH mass ratio. Owing to its large total mass and short inspiral signal in the LIGO-Virgo band, the inspiral parameter estimation of GW190521 is particularly sensitive to model systematics and there is some tension among the binary parameters (including the mass ratio) inferred with different waveforms [48, 59, 60]. Thus, one cannot perform a reliable IMR-like test on this signal. In Fig. 11 we compare the different estimates for q using parameter estimation posteriors provided in these IMR studies with the expected value of q inferred from ringdown. For the ringdown estimate of q , we compare the estimate of A R 330 obtained in [1] along with our fit results. Specifically, we show  \n- (i) the parameter-estimation results on A R lmn from [1] (red) and translated to q by inverting our A R 330 mode fit. The red dot (Ours) provides its best likelihood value obtained from the marginalized distribution q -A R 330 .\n- (ii) the parameter-estimation results on q and χ 1 , 2 obtained by the LVK collaboration [60] (blue) and  \nFIG. 11. Amplitude ratio A R 330 estimated from the set of NR waveforms used in this work and restricted to q < 3. The mild dependence on the spin parameters is reflected by the relatively small spread of the points at fixed q . We add on top the different estimates of the A R 330 -q relation obtained from different mass ratio estimates: Capano+ (orange triangle), Ours (red dot) and LVK (blue square). The values of A R 330 and q predicted from NR are consistent at the 2 σ level with the distributions obtained on q and A R lmn obtained from [60] and [1].  \n<!-- image -->  \ntranslated to A R lmn by using our (3 , 3 , 0) mode fit. The blue-square (LVK) provides its best likelihood value obtained from the marginalised distribution q -A R 330 . Since the LVK does not provide an independent distribution on A R 330 , we get the elongated blue contours using our fit.  \nWe observe that Ours, Capano+, and LVK are all consistent with the parameter-estimation posterior distribution at the 1 σ and 2 σ credible level. Note that the q distribution from Capano+ is obtained from an independent set of fits [40]. We also obtained a rather flat distribution on the phenomenological spin parameter with χ ∈ [0 , 1], which is expected given the mild dependence of both A R 330 ( q, χ pheno ) and δφ 330 ( q, χ pheno ) on χ . However, this paradigm may change shortly with louder detections such as those expected from third-generation detectors [56, 57] and LISA [61].'}