Spectroscopy for asymmetric binary black hole mergers

We study Bayesian inference of black hole ringdown modes for simulated binary black hole signals. We consider to what extent different fundamental ringdown modes can be identified in the context of black hole spectroscopy. Our simulated signals are inspired by the high mass event GW190521. We find strong correlation between mass ratio and Bayes factors of the subdominant ringdown modes. The Bayes factor values and time dependency, and the peak time of the (3,3,0) mode align with those found analyzing the real event GW190521, particularly for high-mass ratio systems.


I. INTRODUCTION
The ringdown signal of a binary black hole (BBH) merger, which is related to the resonating space-time of the resulting final black hole [1][2][3], conveys extensive information about the post-merger object.Moreover, as we shall discuss in this paper, a detailed analysis of the ringdown signal allows one to infer properties of the premerger binary system.At sufficiently late stages, the ringdown waveform can be simply represented as a superposition of quasinormal modes (QNMs).The no-hair theorem in general relativity implies that the frequencies and damping times of the QNMs exclusively depend on the mass and spin of the single remnant black hole (BH).The amplitudes and phases of the QNMs depend on the initial pre-merger binary configuration.Systems with asymmetric progenitor masses tend to exhibit higher amplitudes of subdominant QNMs, while symmetric systems tend to have minimal amplitudes [4][5][6].The detection of multiple QNMs is often referred to as black hole spectroscopy.When multiple QNMs are detected simultaneously, they can be used not only to validate the no-hair theorem, but as we shall show in this paper, also to extract additional information about the mass ratio of the BBH.
In black hole spectroscopy, the detection of multiple ringdown modes enables us to assess the compatibility of a merger event with predictions from General Relativity (GR) [7,8].The detection of multiple ringdown modes obeying these predictions serves as a highly compelling test of the underlying Kerr nature of the final black hole spacetime.* jahed.abedi@uis.no The QNMs are a representation of the post-merger waveform decomposed into a spheroidal harmonic basis.The modes are enumerated by three integers (ℓ, m, n), satisfying the conditions ℓ ≥ 2, −ℓ ≤ m ≤ ℓ, and n ≥ 0. Among these modes, those with n ≥ 1 are commonly referred to as "overtones."The integers ℓ and m represent angular and azimuthal numbers, respectively, while n indicates the overtone index.
Waveforms that model the entire signal through the inspiral, merger, and ringdown (IMR) phases of a binary instead use a spherical harmonic basis to represent the signal.IMR waveforms have no overtones, and so their modes are characterized by just an ℓ and m.Consequently, whenever we use the notation (ℓ, m, n) in this paper -e.g., (2,2,0) -it corresponds to QNMs.When we use the notation ℓm -e.g., 22 -it signifies a mode of an IMR waveform model.
There are still several unresolved theoretical questions concerning black hole QNMs, which also have implications for observational studies.One of these questions relates to the starting time of the ringdown phase.After the formation of the remnant black hole, it initially possesses significant distortions from a Kerr black hole.As time progresses, these distortions dissipate; eventually, the black hole can be regarded as a linear perturbation of a Kerr black hole.The time at which this perturbative description becomes reliable, if at all, remains unclear [8][9][10].Additionally, studies exploring potential non-linear effects in the ringdown regime can be found in [11][12][13][14][15].
The distinct regimes observed in a gravitational wave signal are anticipated to have corresponding behaviors in the strong-field dynamic spacetime region near the binary system [16][17][18].Investigations into black hole horizon geometry during the post-merger phase have shown the potential identification of a ringdown regime using horizon dynamics [19][20][21][22][23]. Recent work also shows the presence of quadratic non-linearities at the horizon [24].The determination of the final black hole parameters based solely on the ringdown signal is sensitive to the assumed starting time of the ringdown [8,[25][26][27].Different choices of the starting time can lead to distinct outcomes [28,29].An analysis focused solely on the ringdown phase needs to explicitly exclude earlier portions of the signal where the perturbative description is not valid.
It has been widely anticipated that the current generation of gravitational wave detectors would predominantly detect the most prominent ringdown mode [30,31].This expectation was rooted in astrophysical assumptions concerning the mass distributions and mass ratios of binary black hole systems across the observable universe, which directly impact the amplitudes of different ringdown modes [4][5][6].Nevertheless, evidence supporting the existence of an overtone accompanying the dominant mode of GW150914 was provided in [32][33][34].These studies demonstrated the feasibility of modeling the gravitational waveform as a combination of ringdown modes, even starting in the merger phase, by incorporating overtones of the dominant mode.However, several questions related to both the data analysis methods used to detect overtones and the validity of the model at merger remain unresolved [27,35,36].The stability of overtones under small perturbations has also been widely studied [37][38][39][40][41].
GW190521 was detected by the Advanced LIGO and Advanced Virgo detectors [42].Although a multitude of scenarios have been put forward regarding the source of this event [43][44][45][46][47][48][49][50][51][52], the most conservative explanation is that it was caused by the merger of a quasi-circular BBH [42,53].Assuming GW190521 is a BBH merger, its inferred total mass would categorize it as one of the most massive binary black hole systems observed thus far [54,55]; other interpretations have even proposed higher total masses [43].The substantial mass implies that a significant portion of the inspiral phase occurs below the sensitive frequency range of the detectors, meaning the recorded signal is dominated by the post-inspiral merger and ringdown phases.Consequently, focusing the analysis on the ringdown phase becomes particularly valuable as it circumvents certain challenges associated with modeling the inspiral phase of the progenitor system.
While the nature of the progenitors of the GW190521 event remains open, the most likely outcome in most scenarios is the formation of a single black hole.The event GW190521 exhibits clear evidence of a prominent ringdown mode emitted by the final black hole following the merger [42].In the study by Capano et al. [56], an additional sub-dominant ringdown mode was identified in the signal.The dominant mode is consistent with the (ℓ, m, n) = (2, 2, 0) ringdown mode of a Kerr black hole, while the second mode is consistent with the sub-dominant fundamental (ℓ, m, n) = (3, 3, 0) mode.As detailed in [56], under the assumption of a Kerr black hole, the Bayes factor favoring the existence of both the (2, 2, 0) and (3, 3, 0) modes over just the (2, 2, 0) mode or the combination of (2, 2, 0) and (2, 2, 1) modes is esti-mated to be 56 ± 1. Subsequent work [57] based on simulated signals and real detector noise showed that the probability of falsely detecting these two modes when only one is present in a quasi-circular binary is ∼ 0.02.An additional agnostic analysis assumed no specific relationship between the ringdown modes.This indicated that only 1 in 250 simulated signals without a (3,3,0) mode produces a result at least as significant as the one observed for GW190521 [57].
A recent ringdown study by Siegel et al. (2023) [58] finds hints of multiple observable QNMs in GW190521.In their analysis, they find evidence for the presence of the (2, 1, 0) mode with an amplitude comparable to, and even exceeding, the (2, 2, 0) mode.They also find evidence for a higher-frequency subdominant mode, in agreement with Capano et al. [56].However, they identify this mode as the (3, 2, 0) QNM instead of the (3, 3, 0) mode, as in Capano et al.Their result aligns more closely with results obtained from the numerical relativity surrogate model NRSur7dq4 [59].In particular, they ascribe the large amplitude of the (2, 1, 0) to large precession in the binary.
Even under the assumption that GW190521 was caused by the merger of a quasi-circular BBH, estimates of its mass ratio vary widely depending on the waveform model used and the prior assumed.In [42], the binary was reported as being approximately equal mass.That analysis used the IMR model NRSur7dq4 (which is a surrogate model derived from numerical relativity simulations [59]) and a prior uniform in component masses.However, a subsequent analysis by Nitz & Capano [60] using a prior uniform in mass ratio m 1 /m 2 (where m 1 ≥ m 2 ) found support for larger mass ratios.Using this prior, Nitz & Capano found a bimodal posterior distribution in the mass ratio when using NRSur7dq4, with the second mode railing against m 1 /m 2 = 6, which is the upper bound of the valid range for NRSur7dq4 [59].The posterior from this analysis is replicated in the left panel of Fig. 1.Notably, the signal-to-noise ratio of the mass ratio ∼ 6 points are slightly larger than the equalmass points.This indicates that the equal-mass scenario favored in [42] is largely due to the choice of prior.
A further analysis of GW190521 by Estelles et al. [61] using the IMR model PhenomTPHM also found support for larger mass ratios, depending on the mass prior chosen.That result, which is reproduced in the right panel of Fig. 1, where a uniform prior in component masses (in detector frame) with a range m det 1,2 ∈ [10, 400]M ⊙ and default mass ratio m2/m1 constraints ∈ [0.035, 1.0] is employed, also yielded a bimodal posterior, with one mode at equal mass and another as mass ratio ∼ 5.The maximum SNR was about equal between the two.The larger mass ratio mode from the PhenomTPHM results yielded a final mass and spin estimate consistent with that found in Capano et al., whereas the NRSur7dq4 results (even at mass ratio ∼ 6) did not [57].
The discrepancy in mass ratio estimates is due to the challenging nature of GW190521.Only about one cy- A comparison of the final mass, mass ratio and SNR of samples from IMR analyses of GW190521, using the NR-Sur7dq4 [60] (Left) and IMRPhenomTPHM (Right) [61] waveform model.The solid lines in the plots represent the 50% and 90% credible contours.The dots' color represents the samples' SNR.Notably, the NRSur7dq4 results exhibit a second mode with higher SNR at larger mass ratio, which agrees with the findings from the ringdown analysis.
cle is observable above noise prior to merger.In order to match both that cycle and the ringdown, the IMR models need to go to high precession and/or larger mass ratio.However, large precession, mass ratios larger than ∼ 2, and the transition from late inspiral to ringdown is one of the hardest areas for waveform models to capture, as there are relatively few NR simulations in this space to calibrate against.For example, the NRSur7dq4 approximant is calibrated to NR simulations with mass ratios up to 4 and dimensionless spin magnitudes < 0.8, and was validated for mass ratios up to 6 by performing mismatch calculations to simulations with total masses up to 200 M ⊙ [59].For comparison, GW190521 is best matched by the NRSur7dq4 model with in-plane spin on the large object of ≳ 0.8, mass ratio ∼ 6, and total mass > 250 M ⊙ .Conversely, PhenomTPHM is calibrated to aligned-spin simulations, relying on a "twisting-up" procedure that may be less reliable late in the inspiral [61].
To better understand the discrepancy in mass ratio estimates for GW190521, in this paper we use QNM templates to analyze a range of different simulated signals.
We invert the typical question of parameter estimation: instead of asking which IMR template best matches (the largely ringdown dominated) GW190521, we ask which QNM models best match the post-merger of an IMR simulation.To that end, we analyze simulated black hole ringdown signals containing both the dominant (2, 2, 0) mode and various combinations of sub-dominant modes, including both fundamental modes and first overtones.The full set of mode combinations considered is listed in Table(I).
Our simulated signals use the IMR waveform model IMRPhenomTPHM [62].In order to better understand the underlying physics of GW190521, this study employs  zero-noise injections to mitigate complexities introduced by real detector noise.The data we analyze therefore consists of only the signal as it would be seen in the detectors without noise being present.This in effect samples over a large population of noise realisations from the assumed Gaussian distribution.While accounting for the detector response, we can then focus solely on the intrinsic characteristics of the signal.
The IMRPhenomTPHM injections are drawn from posterior samples provided in Ref. [61].We investigate the compatibility of BBH waveforms with the presence of subdominant QNM modes.Our findings indicate that within this family of quasi-circular binary systems, previously obtained QNM results favor an asymmetric system for the GW190521 event.
The key results for GW190521 in real data for various combinations of modes are shown in Fig. 2. Figs. 3 and 4 illustrate the outcome of a similar analysis as depicted in  In this paper, we analyze data within a time span of 42 milliseconds that includes both the pre-merger and postmerger regimes.Therefore, through our investigation of the quasinormal modes (QNMs), we are able to identify the point at which QNMs best model the signal.
In this paper we utilized a range of starting times for the ringdown analysis that covers a range even larger than full range of possible uncertainty that may occur in our study.This range is even larger than what we investigated in our previous studies [56,57].Further details can be found in Section III.
In Section II, we provide supplementary information regarding the data treatment in this paper.The methodology employed in this study is similar to the one utilized in Capano et al. [56,57], with the exception of an extended time span of data considered, the utilization of zero-noise data, and the inclusion of additional modes not considered in those works.Section III elucidates the procedure for generating the simulated data sets (A, B, 1 For definition of set A and set B, please refer to Section III. Signal, and Control ) for an event similar to GW190521, as well as the definition of the modified Bayes factor.Finally, in results and discussions IV, we present the results obtained in this study and provide a concise summary of these findings.

II. RINGDOWN MODEL AND SEARCH PIPELINE
In this section, we provide a summary of our data analysis procedure, focused exclusively on analyzing the ringdown component portion of the complete signal.Here, important challenges are the accurate identification of the ringdown section of the data, and restricting the analysis to this section while removing correlations with data outside the desired range.
The ringdown waveform model is described by the mathematical expression or functional form (1) The waveform model for the ringdown can be expressed in terms of the plus (h + ) and cross (h × ) polarizations of the gravitational wave.It depends on parameters such as the total mass of the remnant black hole in the detector frame M f and the source luminosity distance D L .The waveform is decomposed using the spin-2 weighted spheroidal basis −2 S ℓmn , which is a function of the remnant black hole's spin χ f , the inclination angle ι, and the azimuthal angle φ relative to the line-of-sight to the observer.The amplitude and phase of each quasi-normal mode is represented by A ℓmn and ϕ ℓmn , respectively.The complex frequency of the ringdown waveform is defined as Ω ℓmn = 2πf ℓmn +i/τ ℓmn .The characteristic frequency f ℓmn and decay time τ ℓmn are solely determined by the mass and spin of the remnant black hole, as predicted by the no-hair theorem in general relativity.
Our analysis must be restricted to the data where the QNM waveform model is valid.This is achieved using "gating and in-painting" to eliminate the influence of preringdown data [56,63].We employ the gated-Gaussian likelihood implemented in the open-source PyCBC Inference library [64,65].This likelihood function applies gating and in-painting to gravitational wave data assuming the noise is a stationary Gaussian process.
In all our analyses, a gate duration of two seconds is applied, with the gate ending at the selected start time of the ringdown phase.The power-spectral density (PSD) employed in the likelihood calculation is estimated from real detector data around the GW190521 event.This data has been publicly released by the Gravitational Wave Open Science Center [66].For all injections, the sky location (which determines the relative time, phase, and amplitude in each detector) of the ringdown template waveforms is fixed to the values obtained from the maximum likelihood result of Nitz & Capano [60].To -increasingly mirrors the pattern seen in GW190521 (illustrated in Fig. 2).Conversely, when we conduct an equivalent analysis on the Control injections (right panel), featuring only the 22 mode, we find no similarity to the observed pattern.
sample the parameter space during the inference process, we employ the dynesty nested sampler [67].This sampler efficiently explores the parameter space, provides reliable estimates of the posterior distribution and calculates the corresponding evidence.To speed convergence we numerically marginalize over the polarization of the gravitational wave using a discrete grid of 1000 points for each likelihood evaluation.This method has been shown to also improve the estimate of the evidence (and the Bayes factor, the ratio of evidences for two models) [57].
In our study, we explore various signal models that encompass a range of angular and overtone modes, as defined by Eq. 1.A complete list of all the models considered is given in Table I.In all models, the complex frequencies of the modes are determined by the Kerr hypothesis, which provides predictions for the frequencies based on the final black hole's mass and spin (here calculated using the pykerr package [68]).
The priors for all parameters utilized in this study are provided in Table II.Priors for the amplitude parameters are identical to those in Capano et al. [56].Specifically, we select the amplitude of the (2, 2, 1) mode to be within the range of [0, 5] times that of the (2, 2, 0) mode.This choice is motivated by the numerical relativity fits presented in [33].For the (3, 3, 0) mode amplitude, we employ a prior that is [0, 0.5] times the amplitude of the (2, 2, 0) mode.This choice is informed by the results of numerical simulations of binary black hole mergers described in Ref. [6].Note that none of the simulated IMR signals studied in this paper had amplitudes outside this prior boundary.For the (2, 1, 0) and (3, 2, 0) modes, we use the same prior distribution as for the (3, 3, 0) mode.We note that Siegel et al. [58] find that the (2, 1, 0) amplitude can exceed this range.Given that our objective in this study is to conduct a direct comparison of simulated signals with the real data from GW190521, specifically emphasizing the (3, 3, 0) mode, the main conclusions of our work remain unaffected in light of this finding.The maximum strain in the gravitational wave signal is expected to occur near the time of the merger.Using the numerical relativity (NR) surrogate waveform model NRSur7dq4, the LIGO Scientific and Virgo Collaborations (LVC) initially estimated the GPS time of the maximum strain in the Hanford detector to be 1242442967.4306+0.0067 −0.0106 (median ±90% credible interval) [53,59].In the ringdown analysis, we consider discrete starting times for the analysis with 3ms intervals.These starting times span [−18, 24]ms relative to the reference time t ref = 1242442967.445. 2 This extends the range of [−9, 24]ms used in previous works [56,57]

III. SELECTION OF SIMULATED SIGNALS
We produce two distinct sets of simulated signals, ("injections") denoted as set A and set B. The signals are generated using the IMRPhenomTPHM waveform model, and their parameters are drawn from an IMR posterior found analysing GW190521 using this model [61].Each set consists of two subsets Signal and Control.
Set A is constructed such that its signals feature a (3, 3, 0) ringdown mode with comparatively large amplitude.The signal parameters are drawn randomly from the posterior of [61].We then select 100 signals with an amplitude of the IMR 33 mode at least 0.2 times as large 2 Please see next section for details.
as the amplitude of the 22 mode after merger.This selection is made by comparing the modes' amplitudes at a single point after the merger.Specifically, a time interval , and equivalently for the 33 mode.Here, h p/c represents the plus/cross polarization of the geocentric IMR waveform.Due to the differences between the IMR modes and the quasi-normal modes, this method does not directly constrain the amplitude of the (3, 3, 0) QNM.However, this simple approach suffices to select samples where a significant (3, 3, 0) mode is expected.
Through the selection, the signals in this set correspond to systems with higher mass-ratios, serving the purpose to test if GW190521 is an asymmetric system.The parameter draws of set A are identical to the 100 draws used in [57].There, we implemented an additional criterion that the signal-to-noise ratio (SNR) of the IMR 33 mode be at least 4, which restricted that study to a subset of set A.
Set B is chosen to consist of signals of nearly masssymmetric systems from the GW190521 posterior.We expect this set to contain signals with comparatively small amplitudes for the subdominant ringdown modes [5,6].We construct this set by randomly selecting 9  points with mass ratios m 1 /m 2 < 2 from the posterior distribution published in Estelles et al. [61].Adding the maximum likelihood point from the m 1 /m 2 < 2 part of this posterior results in a total of 10 signals in this set.This set provides a basis for comparison with results obtained from set A to assess the significance of the asymmetric characteristics observed in GW190521.
For each set of signal parameters we generate two sets of signals.The Signal set is generated using all modes available in the approximant.The Control contains only the (IMR) 22 mode for this waveform.Both subsets of signals, Signal and Control, use the same parameter samples drawn from the GW190521 posterior.The objective of the Control study is to examine whether subdominant modes are detected in the Signal data and to investigate the start time of the late-time ringdown stage.
Similar to the findings presented in Nitz & Capano [60] with the NRSur7dq4 waveform model, Estelles et al. [61] observed a bimodal posterior distribution in the component masses for GW190521.One mode favored nearly equal masses and another mode favored mass ratios of approximately 6:1.As it is depicted in Fig. 1, these bimodal features indicate the presence of two distinct configurations or scenarios for the binary black hole system associated with GW190521, characterized by different mass distributions between the two black holes.In our analysis, the injections labeled as set A corresponds to the asymmetric configuration or mode, while the set B captures the symmetric configuration of the GW190521 event.
Note that in this context, the notation ℓm refers to the spherical harmonics, which is the basis commonly used in inspiral-merger-ringdown (IMR) waveform models.The spherical harmonics are distinct from the spheroidal harmonics employed for quasi-normal modes (QNMs).While the QNM analysis uses the spheroidal harmonics to describe the ringdown phase, the IMR models employ the spherical harmonics to represent the overall behavior of the gravitational wave signal during the inspiral, merger, and ringdown stages.Therefore, our estimations of the Bayes factor in the absence of modes other than 22 can be considered as an upper bound on the underlying (ℓ, m, 0) QNM mode.The presence of pre-cession further complicates the picture, since precession causes modes with the same ℓ but different m to mix together.However, in this paper, we make the assumption that the late-stage behavior of the 22 mode in the IMR waveform closely resembles that of the (2, 2, 0) + (2, 2, 1) QNMs.This assumption allows us to examine the behavior of the simulated signal injections specifically during the ringdown phase, while excluding the influence of QNMs other than ℓ = 2 and m = 2.
To incorporate the injections into the analysis, both the Signal and Control subsets are added to zero-noise data.This means that the data consist of only the simulated signals, while the PSD used in the analysis is estimated from data around the event GW190521.These injections are inserted at random times surrounding the estimated merger time of GW190521.The process follows the same methodology as described in Capano et al. [57].There, an offset time t offset is drawn uniformly from the range ± [4,20] s.This offset time is then added to the coalescence time t c , which is drawn from the relevant posterior distribution for each injection.The purpose of this offset is to introduce variability in the analysis and prevent any contamination from the original GW190521 data.A gap of ±4 s around GW190521 is maintained to ensure the separation between the injections and the original event.Even though the presence of this gap is not essential when using zero-noise data, we here maintain the same methodology as used in Capano et al. for consistency [57].
The ringdown analysis is then performed on a grid of times surrounding each injection.The grid, used for the validation of the Kerr Bayes factor, spans the range [−18, 24] ms.This range is wider than the previous range used in Capano et al. [56,57], which was [−9, 24] ms.
Similar to the analysis conducted in Capano et al. [56], we establish a reference time, denoted as t inj ref , for each injection to construct the grid of times used in the ringdown analyses.For each injection, the reference time is defined as t inj ref = t ref +t offset , where t ref = 1242442967.445GPS seconds represents the estimated geocentric merger time of GW190521 and t offset is a randomly chosen offset unique to each injection.The reference time t ref is ob-tained from the maximum likelihood parameters derived from the NRSurrogate analysis conducted in Nitz & Capano [60].Note that t inj ref is distinct from the injection's coalescence time t inj c .The ringdown analysis requires the identification of the time window during which the QNMs are observable.The QNM model is not valid at early times, such as before merger, where nonlinear components may be significant.However, if the analysis is conducted too long after the merger, the signal will have damped away to the extent that only the dominant mode remains observable.
In the methodology presented in [56,57], we address this challenge through the use of Bayes factors.A key objective of Bayesian inference is to derive a posterior probability distribution p(θ|d, M ) for a set of model parameters θ, given data d and a model M [69].Here, d corresponds to the data from a (gravitational-wave) measurement, and the model M describes the behaviour of noise and signal expected to produce the data.By applying Bayes' theorem, the posterior probability distribution is expressed as Here, L(d|θ, M ) is the likelihood function of the data given the parameters θ, which represents the probability of the detectors measuring data d given a model hypothesis M with source properties θ. π(θ|M ) is the prior distribution for θ, encoding the prior knowledge about the parameters before studying the data d.Z is a normalization factor called the "evidence" and quantifies how well the data aligns with the hypothesis, Z(d|M ) ≡ dθL(d|θ, M )π(θ|M ).
These provide a means to quantify the evidence that the data contain observable modes ⃗ X = (2, 2, 0), ... at a specific time t − t ref .By calculating the evidence Z ⃗ X (t) and comparing it to the evidence for the (2, 2, 0)-only model (Z 220 ) at the same time, we obtain the Bayes factor for the model.If model ⃗ X and the (2, 2, 0) model have equal prior weight (meaning that a priori the two models are considered to be equally valid) then the Bayes factor gives the odds ratio for model ⃗ X being true relative to the (2, 2, 0)-only model.
The (2, 2, 0)-only model is not a good representation of the signal at the merger [33].Therefore, if we find a large value for Z ⃗ X /Z 220 at a particular time, it is unclear whether this is due to the ⃗ X modes being a good fit for the signal or if it is simply because the (2, 2, 0)-only model is inadequate at that time.In other words, the ratio Z ⃗ X /Z 220 only indicates whether the ⃗ X modes provide a better fit than the (2, 2, 0)-only model, not necessarily whether the ⃗ X modes are truly observable.This issue becomes more pronounced as we approach the merger.
To address this concern, we leverage the observation from Ref. [33] that including overtones of the dominant mode improves the fit to the signal close to or at the merger compared to the (2, 2, 0)-only model.We modify the Bayes factor as follows: for all models ⃗ X ̸ = (2, 2, 0) + (2, 2, 1).For the (2, 2, 0) + (2, 2, 1) model, we simply use B = Z 220+221 /Z 220 .This modification allows us to identify the most likely observable modes and determine the time at which they are most observable.By considering the maximum evidence between the (2, 2, 0)-only model and the model including an additional overtone, we account for the limitations of the (2, 2, 0)-only model and obtain a more robust assessment of the observability of the different modes at various times.
When applying this methodology to GW190521, the Bayes factor B(220 + 330) peaks at t ref + 6 ms with a value of 56 ± 1 [56,57].Fig. 2 shows the result of this methodology.This indicates that the (2, 2, 0) + (3, 3, 0) model is approximately 56 times more likely to be true than the (2, 2, 0)-only model.We see that in Fig. 6 this peak time and value of bayes factor is consistent with what we obtain from GW190521.
We now apply this analysis to our injection sets.Similar to the analysis of GW190521, we conduct the analysis on a grid of times spanning t ref + [−18, 24] ms.However, to reduce computational costs given the large number of analyses, we sample the grid at intervals of 3ms instead of the 1ms intervals used in the previous work by Capano et al. [56].The results of this analysis are shown in Figs. 3, 4, and 5.

IV. RESULTS AND DISCUSSIONS
In this paper, we investigate the nature of the GW190521 event by analyzing two sets of injections labeled as A and B, which are generated from the posterior samples of an IMR analysis of GW190521.Our analysis reveals that this event does not conform to a symmetric system.This is seen in the qualitative behavior of the time-dependent Bayes factors in favor of a subdominant mode.This finding further reinforces the previous observations of a bimodal posterior distribution in the component masses [60,61] (high mass ratio preference in Fig. 1) using NRSur7dq4 and IMRPhenomTPHM waveforms.Figs. 3 and 4 demonstrate that as the mass ratio increases in the Signal subset, the pattern of multiple modes in our study becomes more similar to the pattern observed in GW190521 (shown in Fig. 2).However, when we perform the same analysis on the Control injections, which only contain the 22 mode, no resemblance to the observed pattern is found.
Remarkably, in Figs. 3 and 4 for the high mass-ratio injections of the set A Signal subset, the peak Bayes factor corresponding to the (3,3,0) mode is observed to be located near the peak of GW190521 shown in Fig. 2. The location of the peak is shown more clearly in the colorbar of Fig. 6 when considering the entire set A of injections.As the mass ratio increases, a distinct difference is observed when comparing the Signal subset to the Control subset in both sets A and B. Specifically, in Figs. 3 and  4, where high mass-ratio systems are present, there is no resemblance between the two subsets.However, the pattern in these figures closely matches that of GW190521 in Fig. 2.
In Fig. 6, it is evident that the Signal subset is shifted towards the Bayes factor value of GW190521 compared to the Control subset.This indicates that the inclusion of the (3,3,0) mode in the simulated signals for high mass ratio injections results in a similar Bayes factor value and pattern as observed for GW190521 in Fig. 2. For the low mass ratio injections in set B, both the Signal and Control subsets exhibit no Bayes factor values close to that of GW190521.This observation is further supported by Fig. 5, where the observed patterns do not resemble that for GW190521 in Fig. 2.These findings raise questions regarding the symmetry of the GW190521 event, despite the injections being drawn from a part of the posterior distribution for GW190521.In both sets A and B, as we move to low mass ratios, the Signal and Control subsets exhibit similar patterns, distinct from the observed pattern in GW190521 (Fig. 2).
In our analysis of simulated signals, we observe that the (3, 3, 0) mode exhibits a more prominent peak compared to other combinations of modes studied in this paper.The presence of a peak in Bayes factor around time-t ref ≃ 6 ms for the Signal subset, as opposed to its absence for the Control subset, suggests that at this point the linear (3,3,0) QNM is matching the full nonlinear 33 mode of the IMRPhenomTPHM model (see Fig. 3).This finding is further supported by this peak's absence for symmetric systems, in both subsets Signal and Control.Therefore, this observation highlights the specific region in the ringdown phase of the GW190521 binary black hole merger event where spectroscopy can be effectively conducted.
Furthermore, we observed that the values of the Bayes factors for high mass ratio systems in the simulated signals align with the results obtained from GW190521.
Generally, the Bayes factor serves to compare only two models, α and β.However, it's crucial to understand that a preference for model β over model α does not necessarily imply that model β is the accurate or true model.This is because both model α and model β might not be appropriate models, and hence, a preference for model β over model α does not validate model β as a suitable representation.Therefore, it is important to acknowledge that while higher-order modes can seemingly fit the premerger regime in Figs. 2 and 3, such a fit alone does not imply physical validity.The comparison of Bayes factors for the (2, 2, 0) and (2, 2, 0) + (2, 2, 1) models highlights that an improved fit is achieved with a higher number of modes, even though these models do not align well with the pre-merger time.Essentially, within the non-physical region, a greater number of modes facilitates a superior fit to the data.Additionally, noise variations will introduce uncertainty in the values of the Bayes factors, as detailed in [57].Here, we only consider injections without noise, approximating the behaviour expected when averaging over many instances of random noise.
In this analysis, since we employed injections without noise and the implications from black hole spectroscopy indicate that the late stage of black hole ringdown is accurately described by QNMs, we can confidently affirm the reliability of identifying the (3,3,0) mode using this method at later times.
In future work, we plan to study the recovery of intrinsic binary parameters from simulated signals through black hole spectroscopy.
As discussed in Section I, the NRSur7dq4 approximant is calibrated using NR simulations with mass ratios only up to 4 and valid only up to mass ratios of 6.Therefore, for achieving more precise results, more accurate waveforms are essential.
In conclusion, our findings strongly suggest that GW190521 is an asymmetric system that exhibits the presence of subdominant modes.While a more comprehensive analysis with a larger number of simulated signals is needed to obtain more precise evidence and uncertainties, our initial results indicate that GW190521 exhibits similarities to signals with mass ratios ⪆ 8.5.Moreover, our simulation results provide further support for the presence of the (3, 3, 0) quasi-normal mode in the ringdown phase of GW190521.
Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.S. K. acknowledges support from the Villum Investigator program supported by VIL-LUM FONDEN (grant no.37766) and the DNRF Chair, by the Danish Research Foundation.
FIG.1.A comparison of the final mass, mass ratio and SNR of samples from IMR analyses of GW190521, using the NR-Sur7dq4[60] (Left) and IMRPhenomTPHM (Right)[61] waveform model.The solid lines in the plots represent the 50% and 90% credible contours.The dots' color represents the samples' SNR.Notably, the NRSur7dq4 results exhibit a second mode with higher SNR at larger mass ratio, which agrees with the findings from the ringdown analysis.

Fig. 2 ,
Fig. 2, but applied to simulated signals (set A 1 ) that exhibit a relatively high amplitude of the (3,3,0) mode (left panel) and simulated signals containing only the 22 IMR mode (right panel).Fig. 5 represents the same analysis as shown in Fig. 4 while applied to simulated signals in set B. The plots in Figs. 3 are arranged in order based on different mass ratios (only for three injections), with each row representing a distinct mass ratio.Fig. 6 represents the relation between the maximum Bayes factor for the (3,3,0) mode and the initial binary's mass ratio.In this paper, we analyze data within a time span of 42 milliseconds that includes both the pre-merger and postmerger regimes.Therefore, through our investigation of the quasinormal modes (QNMs), we are able to identify the point at which QNMs best model the signal.In this paper we utilized a range of starting times for the ringdown analysis that covers a range even larger than full range of possible uncertainty that may occur in our study.This range is even larger than what we investigated in our previous studies[56,57].Further details can be found in Section III.In Section II, we provide supplementary information regarding the data treatment in this paper.The methodology employed in this study is similar to the one utilized in Capano et al.[56,57], with the exception of an extended time span of data considered, the utilization of zero-noise data, and the inclusion of additional modes not considered in those works.Section III elucidates the procedure for generating the simulated data sets (A, B,

FIG. 3 .
FIG. 3. Results of the ringdown analysis for simulated signals of set A in terms of Bayes factors as shown in Fig. 2. Set A consists of 100 injections drawn from the GW190521 IMR-posterior samples.The injected signals use the IMRPhenomTPHM waveform model.Both columns show the same injection parameters, but either using all modes implemented in IMRPhenomTPHM (left, Signal subset) or only the 22 mode (right, Control subset).Top row: Highest mass ratio injection, q = 19.2Center row: Injection with (3, 3, 0) Bayes factor closest to that for GW190521, q = 8.8.Bottom row: Lowest mass ratio injection of set A, q = 1.1.The ringdowm analysis method and selection of modes are identical for all results shown.The dashed lines show where the (3,3,0) mode Bayes factor for GW190521 peaks in Fig. 2.This figure illustrates that as the mass ratio increases within the Signal injection set (left panel), our analysis patterns of multiple modes-encompassing peak time, Bayes factor values, etc.-increasingly mirrors the pattern seen in GW190521 (illustrated in Fig.2).Conversely, when we conduct an equivalent analysis on the Control injections (right panel), featuring only the 22 mode, we find no similarity to the observed pattern.

FIG. 4 .FIG. 5 .
FIG. 4. Results of the ringdown analysis of set A, using as templates the (2, 2, 0) + (3, 3, 0) quasi-normal modes.Set A has 100 injections drawn from the GW190521 IMR-posterior samples which use the IMRPhenomTPHM waveform model.The injections' mass ratio ranges from q = 1.1 to q = 19.2.For visibility, we display only 17 randomly selected injections of the total 100, uniformly distributed in mass ratio.Injections for the left panel have all modes of IMRPhenomTPHM, while those for the right panel have only the 22 mode.The ringdowm analysis method and template modes are identical for left and right.The red curve represents the result from the (2, 2, 0) + (3, 3, 0) analysis for GW190521 and dashed lines show where its Bayes factors peak.In the left panel, this figure demonstrates that as the mass ratio increases in the Signal injections, our analysis pattern-encompassing the peak time, Bayes factor values, etc.-for (2, 2, 0) + (3, 3, 0) progressively aligns with the pattern observed in GW190521 (illustrated in red).In contrast, conducting a similar analysis on the Control injections (right panel), which feature only the 22 mode, reveals no resemblance to the observed pattern.This figure prominently demonstrates the robust support from the Bayes factor for the existence of both (2, 2, 0) + (3, 3, 0) modes over either (2, 2, 0) or (2, 2, 0) + (2, 2, 1) modes at ∼ t ref + 6 ms specifically for high mass ratio systems.
FIG. 6. Set A simulated signals: Bayes factors B 220+330 220,221 for the (2, 2, 0) + (3, 3, 0) model, maximized over time, vs. mass ratio of the corresponding injection.Bayes factors are again calculated with respect to the stronger of the two models consisting of either the (2, 2, 0) or the (2, 2, 0) + (2, 2, 1) modes.The vertical dashed line shows the Bayes factor for GW190521.The colorbar represents the time of the peak (3, 3, 0) Bayes factor relative to the reference time.Left: Injections with all IMR modes (Signal ).Right: The same injections with only the 22 IMR mode (Control ).For high mass-ratio systems (set A), the peak time and Bayes factor values for the Signal injections (left panel) are shifted towards those obtained for GW190521 (dashed vertical line), compared to those of the Control injections (right panel).
measured from the peak of the 22 IMR waveform.Here, t L1 ref = 1242442967.4243s (GPS time) denotes the estimated merger time of GW190521 at the Livingston detector (t ref = 1242442967.445s is the geocentric GPS reference time), and t L1 c stands for the coalescence time of the specific sample at the Livingston detector.In this comparison the amplitude of each IMR mode, e.g.A 22 , is defined as A 22 = (h 2 p,22 + h 2 c,22 ) 1 2

TABLE I .
Combinations of black hole ringdown modes analyzed in this work.

TABLE II .
The prior distributions of the sampling parameters for the models utilized in this study as defined by Eq. 1.