Sky marginalization in black hole spectroscopy and tests of the area theorem

The direct detection of gravitational waves (GWs) throughout the past decade has been one of the most significant advancements in observational relativity. These observations not only confirmed the existence of GWs, but also opened the door to direct tests of General Relativity (GR) in the strong-field regime. One such prediction that can be verified is Hawking’s area theorem, which states that the remnant from a binary-black hole (BBH) merger must have an event horizon surface area greater than the sum of the progenitor horizon areas [1].

Tests of the area theorem using GWs have been carried out previously [2, 3, 4]. Generally, they involve measuring the initial two black holes’ mass and spin (from which the total initial area $A_{i}$ is derived) from the pre-merger part of the signal, during which the two black holes inspiral into each other. The area of the final black hole is independently measured using the GW that is emitted during the post-merger, or “ringdown” phase. In both cases Bayesian inference is used to produce a “posterior” probability density on the black holes’ parameters. Ideally, all parameters that describe the BBH should be allowed to vary in the analysis, in order to fully account for all statistical uncertainties. However, previous tests of the area theorem have fixed the sky position and coalescence time $t_{c}$ of the events to nominal values when doing their analysis [2, 4]. Fixing the values in this way may lead to biases in the resultant parameter estimates, obfuscating the true nature of the system [3, 5]. At the very least, it may cause an underestimate of the statistical uncertainty of measured parameters, yielding constraints on deviations from GR that are misleadingly strong.

The sky location and $t_{c}$ have been fixed in earlier studies due to technical hurdles in calculating the likelihood function. In order to independently analyze the pre- and post-merger portion of the signal it is necessary to excise the post- and pre-merger data, respectively, from the analysis. Excising the pre-merger data is also necessary when just analyzing the post-merger signal for tests of the no-hair theorem using quasi-normal modes (known as black hole spectroscopy). In either case, the excision (or “gate”) results in a modified likelihood function that cannot be solved using conventional, frequency-domain means. Existing pipelines for doing such analyses in the time domain, such as pyring [6, 7] and ringdown [8], instead calculate the likelihood by numerically inverting the noise covariance matrix for the data. This calculation can be numerically costly [9, 5]. If the sky location and $t_{c}$ are not fixed, the gate will vary in time, and the full likelihood will need to be recalculated for every unique gate time, significantly increasing computational costs. However, if the sky location and $t_{c}$ are fixed, the gate is static, and the most computationally demanding part of the calculation need only be done once.

Finch & Moore [10] devised a method to overcome this issue and vary the sky location and $t_{c}$. In their method the inspiral and merger are modeled with a wavelet series stitched to the beginning of the ringdown. This allows them to use the traditional frequency-domain likelihood, and vary the start time of the ringdown. However, this does not allow for an area theorem test, since the parameters of the initial black holes cannot be estimated from the wavelet signal model. Furthermore, the traditional frequency-domain likelihood causes the start of the ringdown to be coupled to information from just before the merger due to convolution with the whitening filter. This means that recovered pre- and post-merger parameters will not be independent measurements.

This paper presents a novel method of calculating the likelihood function using the parameter estimation code PyCBC Inference [11]. This code already utilizes “gating and in-painting” [12] to excise data from the analysis, which allows for cheaper likelihood evaluation under certain conditions as compared to the methods used by pyring and ringdown [9, 5]. However, an additional normalization factor has traditionally been omitted from this calculation, as this also required computationally expensive numerical methods to evaluate. This paper presents a method of linear interpolation that calculates this normalization factor with good approximation. Together, these methods allow for fast likelihood calculations regardless of gate position, allowing for marginalization over $t_{c}$ and sky location and a full accounting of parameter uncertainties.

This paper is structured as follows. Section II describes the modifications to the likelihood calculation in PyCBC in detail. Full waveform analyses can be conducted with these modifications to marginalize over sky location and $t_{c}$. However, the method cannot be used in partial waveform analyses; Section III documents the issues that arise in these models. Using this method, a test of the Hawking area theorem is conducted using data from GW150914 [13]. Section IV describes the configuration of the analyses used in this test, and Section V describes and analyzes the results. Section VI summarizes the findings and the implications of the marginalization method for future GW analyses.

PyCBC Inference [11] utilizes Bayesian inference to conduct parameter estimation on GW events. Bayes’ theorem is used to extract information about the parameter space $\bm{\vartheta}$ from a data set $\bm{s}$ [14]. Assuming the data set is composed of a signal $h$ and noise $n$ such that $\bm{s}=\bm{n}+\bm{h}$, the probability of observing a specific $\bm{\vartheta}$ is given by

The term $p(\bm{s}|\bm{\vartheta},h)$ is the likelihood, which describes the probability of observing a signal $\bm{s}$ assuming the event has a parameter space $\bm{\vartheta}$. The term $p(\bm{\vartheta}|h)$ is the prior, a distribution that describes the a priori probability of observing $\bm{\vartheta}$ given a signal model $h$. The denominator $p(\bm{s}|h)$ is the evidence, a normalization factor used to compare analyses using different models for $h$. The resultant probability distribution on the left-hand side of the equation is known as the posterior.

The prior is chosen at the discretion of the analyst based on assumed plausible values for $\bm{\vartheta}$, whereas the likelihood is calculated directly from the data. For a system of $K$ detectors each taking $N$ time series samples of an event with a stochastic Gaussian noise background, the likelihood can be written as [15]

where $\bm{\Sigma}_{d}$ is the covariance matrix associated with $n$ in detector $d$.

To evaluate the likelihood, further assumptions must be made to calculate $\bm{\Sigma}_{d}^{-1}$. The result of these assumptions is:

where the inner product is defined by Eq. (22). (The full derivation of this likelihood function is given in Appendix A). This gives the likelihood function for a signal that is assumed to be entirely noise. To get the likelihood function with respect to the signal model $h$ evaluated for a parameter space $\bm{\vartheta}$, Eq. (3) can be rewritten by substituting $\bm{s}_{d}\rightarrow\bm{s}_{d}-\bm{h}_{d}(\bm{\vartheta})$ on the right-hand side, yielding

(Here and throughout, we use $\log$ to refer to the natural logarithm.) For a GW analysis that examines the entirety of $\bm{s}$, this expression is easily evaluated using the approximated eigenvalues of $\bm{\Sigma}_{d}$. However, many analyses only examine a portion of the time series (generally either the pre- or post-merger signal). Therefore, a region of $\bm{s}$ is omitted, or “gated”, to reduce numerical biases due to the merger. Doing this also requires excising rows and columns from $\bm{\Sigma}_{d}$, thereby breaking the Toeplitz form of the matrix and the corresponding approximations. The inverse of the truncated covariance matrix $\bm{\Sigma}_{d,tr}$ is evaluated in PyCBC using “gating and in-painting” (explained in detail in [12] and Appendix B), allowing for easy calculation of the first term of Eq. (4).

The second term, however, requires the calculation of $\log\det\bm{\Sigma}_{d,tr}$, which is a notoriously expensive numerical problem. The most widely-used matrix decomposition methods are $\sim\!O(N^{3})$ [16], which become impractical for GW time series such as GW150914, where $N\!\sim\!10^{4}$. This problem is amplified when varying sky position and $t_{c}$ in a gated analysis. The start and end times of the gate must be converted between the geocentric and individual detector frames. The time conversions directly depend on the time of merger $t_{c}$ and the sky location of the event relative to the detectors. Therefore, when varying sky location and $t_{c}$, the gate times in the detector frames will also vary. Subsequently, the elements of $\bm{\Sigma}_{d,tr}$ will vary with each unique sky location and $t_{c}$ value. Analyses that marginalize over these parameters require recalculating $\log\det\bm{\Sigma}_{d,tr}$ on multiple steps of the sampler, which would be impractical using numerical decomposition methods. This is not an issue if the sky location and $t_{c}$ are fixed, since in that case $\det\bm{\Sigma}_{d,tr}$ will not change throughout the analysis (and in fact can be ignored, as it amounts to a constant normalization term).

In order to marginalize over sky location and $t_{c}$ it is necessary to have a fast method for calculating $\det\bm{\Sigma}_{d,tr}$. We find that $\log\det\bm{\Sigma}_{d,tr}$ is strongly linearly correlated to the length of the power spectral density (PSD) of the data, which is equivalent the number of rows and columns in $\bm{\Sigma}_{d,tr}$. Fig. 1 shows this relationship using $\log\det\bm{\Sigma}_{d,tr}$ values calculated for various gate sizes applied to GW150914 data. Using a least squares fit to the points, the correlation coefficient was calculated as $0.999<R^{2}<1$, which for a sample of 8 points implies a highly significant correlation [17].

We also find that the position of the gate in the time series has no significant effect on the value of $\log\det\bm{\Sigma}_{d,tr}$. The maximum range over which the determinant values varied was $\sim\!10^{-6}$. However, since these values were generally of order $\sim\!10^{6}$, this represents a maximum fractional change of $\sim\!10^{-12}$ over the entire time domain. While gate position does have a minor effect on determinant value, the relative effect is so small that the value is well-approximated as a constant for a static gate length.

Using these two facts, we calculate the normalization term in Eq. (4) using a linear interpolation based solely on the size of $\bm{\Sigma}_{d,tr}$. Specifically, a least squares linear fit is generated using the determinant values for the full matrix and three differently-sized truncated matrices. The full matrix determinant is calculated using its approximated eigenvalues, while the truncated determinants are calculated numerically using SciPy [18]. The determinant can then be calculated on each step through linear interpolation based on the size of the truncated covariance matrix.

The strategy outlined in Sec. II fixes the technical hurdle of calculating the likelihood when the gate length and position are varied. Even with that, however, we find that it is not possible to only model a portion of a signal — whether it be the post-merger or the pre-merger — if $t_{c}$ or sky location are uncertain. This is because a larger likelihood can generally be obtained if the gate is shifted such that it removes as much of a signal as possible.

From Eq. (4) it is evident that the likelihood is maximized as $\bm{s}_{d}-\bm{h}_{d}(\bm{\vartheta})\rightarrow 0$. When analyzing the entire time series this can only occur if the template $\bm{h}_{d}(\bm{\vartheta})$ is similar to the signal that exists in the data (as desired). However, if a variable gate is involved, then a large likelihood can also be obtained if the gate is shifted such that it excises most (or all) of the signal in the data. The template then only needs to match the remaining noise. Since the noise manifold is larger than the signal manifold, a larger volume of prior space will generally be able to match the noise than the volume that matches the signal. For example, when using a QNM model, if the gate is shifted such that it excises the entire signal then one only needs to reduce the amplitude of the QNM template below noise level to get a relatively large likelihood. The same likelihood would then be obtained regardless of what the other parameters describing the template are set to. The end result is that the posterior probability will favor excising the entire signal even if the template that is being used can match the signal.

Specifically, in pre-merger-only analyses, the gate will be shifted to as early a time that is possible while in post-merger-only analyses, the gate will be shifted to as late as possible. These phenomena are illustrated in the middle row of Fig. 2. Shown are the results from separate pre-merger and post-merger analyses of a simulated signal in which the gate time is varied due to time and sky-location marginalization. The maximum likelihood template in the pre-merger-only analysis has a $t_{c}$ roughly $25\text{ms}$ earlier than what was injected. This significantly shortens the pre-merger signal, leading to inaccurate measurements in parameters such as total spin. Meanwhile, the maximum likelihood template for the post-merger-only analysis had a ringdown start time $\sim\!10M_{\odot}$ later than the injected value, well into the regime where the signal is expected to be noise-dominated. This leads to a “prior-in prior-out” posterior – a roughly uniform distribution across parameter space – as the signal model can fit any arbitrary model to late-time noise.

Fundamentally, this problem is due to the gated signal model not being the appropriate model for the observed data. Excising data from the analysis is mathematically equivalent to assuming that the excised data is Gaussian noise and marginalizing over all possible realizations of it. A property of multivariate Gaussian distributions is that marginalizing over a subset of dimensions yields another multivariate Gaussian distribution with the marginalized dimensions excised from the covariance matrix. This is exactly the same form as Eq. (3); in our case, each time sample is a dimension in the multivariate Gaussian. The problem is the excised times are not Gaussian noise. They contain a signal, albeit a portion of the signal that we want to ignore. Consequently, a gated signal model is not representative of the data.

This issue could be mitigated by modifying the prior to exclude the portion of parameter space that matches the noise while keeping the portion that matches the signal. With a QNM template this would mean setting a lower bound on the amplitude such that it is above noise level. However, modifying the prior in this manner is challenging to do in an unbiased way, as it would require a priori knowledge of both the noise and signal properties. It also does not solve the fundamental problem that the gated signal model is not representative of the data.

We resolve this issue by simultaneously yet independently modeling both the pre- and post-merger signals. Each domain is treated as a separate gated analysis. The signal parameters used in each domain are independent of each other, except for a set of common parameters $\bm{\vartheta}_{\rm com}$. In our case the common parameters are the right ascension $\alpha$, declination $\delta$, and geocentric coalescence time $t_{c}$. These parameters determine the coalescence time $t^{d}_{c}$ in each detector $d$, which the gate time depends on. The likelihood is calculated in each domain, then combined to get the overall likelihood.

Explicitly, our algorithm to calculate the likelihood for a given set of parameter values $\bm{\vartheta}=\{\bm{\vartheta}_{\rm insp},\bm{\vartheta}_{\rm rd},\bm{\vartheta}_{\rm com}\}$ is:

1. 1.

   Generate the pre-merger (“inspiral”) template $h^{\rm insp}$ using parameters $\bm{\vartheta}_{\rm insp}$.

2. 2.

   Project the template into each detector’s frame using $\bm{\vartheta}_{\rm com}=\{\alpha,\delta,t_{c}\}$.

3. 3.

   Excise times after $t^{\rm det}_{c}$ by gating and in-painting the residual $s_{d}-h_{d}^{\rm insp}(\bm{\vartheta}_{\rm insp},\bm{\vartheta}_{\rm com})$ using a $1\,$s gate that spans $[t^{d}_{c},t^{d}_{c}+1)$.

4. 4.

   Calculate the pre-merger likelihood $p(\bm{s}_{net}|\bm{\vartheta}_{\rm insp},\bm{\vartheta}_{\rm com},h^{\rm insp})$ via Eq. (4).

5. 5.

   Repeat steps 1-4 for the post-merger (“ringdown”) template $h^{\rm rd}$ to get the post-merger likelihood $p(\bm{s}_{net}|\bm{\vartheta}_{\rm rd},\bm{\vartheta}_{\rm com},h^{\rm rd})$. However, for the ringdown the gate spans $[t^{d}_{c}-1,t^{d}_{c})$; i.e., it ends at $t^{d}_{c}$ whereas it starts at $t^{d}_{c}$ for the pre-merger. The data used in each domain is thereby (nearly) mutually exclusive (see Appendix B for more details).

6. 6.

   The total likelihood is then

   	$\displaystyle p(\bm{s}_{net}|\bm{\vartheta},h)$	$\displaystyle=$
   	$\displaystyle p(\bm{s}_{net}$	$\displaystyle|\bm{\vartheta}_{\rm insp},\bm{\vartheta}_{\rm com},h^{\rm insp})p(\bm{s}_{net}|\bm{\vartheta}_{\rm rd},\bm{\vartheta}_{\rm com},h^{\rm rd})$		(5)

We fix the issue of the gate trying to excise the signal by using this hierarchical likelihood. The two domains offset each other: in order for the post-merger gate to shift to later times the pre-merger template must match more of the signal, and vice versa. This also addresses the fundamental issue with the gated signal model highlighted above: our global model for the entire data set now contains a non-Gaussian element (the other domain’s signal model) in each domain’s excised region. Note also that no coupling occurs across the domain boundaries due to the whitening filter.

The bottom plots in Fig. 2 show the maximum likelihood waveform templates from an analysis with this configuration. Besides sky location and $t_{c}$, both models generated independent parameter measurements. The resultant waveform templates closely match those obtained by fixing sky location and $t_{c}$ to the injected values (top row of Fig.  2); the erroneous gate motion observed in the pre- and post-merger-only analyses is no longer present.

We performed eight analyses in total. Each analysis was differentiated by the waveform used, post-merger approximant, and whether or not sky location and $t_{c}$ were marginalized over. All other aspects of the analyses were kept constant. Both the pre- and post-merger signals in all models utilized the PyCBC model GatedGaussianMargPol, which inherits the normalization protocols described in Section II. The dynesty sampler [20] was used to generate posterior distributions. The samplers in these analyses used 4000 live points to ensure the convergence of each model (see Appendix C).

Half of the analyses used the original data of GW150914 with a sample rate of 2048 Hz obtained from the Gravitational Wave Open Science Center [21]. To validate these results, we repeat each run on a simulated signal in zero noise. The simulated signal was generated using the IMRPhenomXPHM waveform approximant [22] with the maximum likelihood parameters for GW150914 from Ref. [19].

All template models utilized IMRPhenomXPHM (abbreviated here on as IMR) to model the pre-merger signal of the waveform. Half of the models used this IMR approximant to model the post-merger signal. The other half utilized a quasinormal mode (QNM) approximant to model the post-merger section of the waveform. Table 1 lists the sampled parameters and priors used in the pre-merger models, and Table 2 lists the same for the post-merger models. The QNM approximant was configured such that the ringdown was composed of a dominant $(2,2,0)$ mode and a subdominant $(2,2,1)$ mode as proposed by [23]. All models apply the ringdown model starting at merger time $t_{c}$. The priors of the QNM post-merger models were restricted such that the $(2,2,0)$ mode contribution to the post-merger signal-to-noise ratio (SNR) was at least 2. This condition was imposed to ensure that the dominant mode was present in the model, preventing possible “label-switching” that may occur due to the $(2,2,1)$ mode erroneously matching to the $(2,2,0)$ mode in the signal.

Only half of the models allowed for sky location and $t_{c}$ to vary. The other half fixed these parameters to nominal values for GW150914 to replicate previous works. Specifically, the fixed parameter analyses set $\alpha=1.95$, $\delta=-1.27$, and $t_{c}=1126259462.408$, in accordance with the maximum likelihood values in [23].

As a preliminary test, the sky location posteriors of the analyses with variable sky/$t_{c}$ parameters on real data were plotted. As seen in Fig. 3, the analyses were able to recover most of the posterior for the full IMR analysis conducted in [19].

The simplest area theorem test is to compare the progenitor horizon areas $A_{1}$ and $A_{2}$ to the remnant horizon area $A_{f}$ to check that

In natural units ($G=c=1$) the area of each black hole is given by [1]:

where $M$ the black hole’s mass and $\chi=J/M^{2}$ is its dimensionless spin.

A more robust test can be performed by comparing the measured change in area to the expected change in area [2]

Here, $A_{f,\rm measured}$ is the area of the final black hole inferred from the post-merger analysis and $A_{i}$ is the sum of the initial areas inferred from the pre-merger analysis. The expected area of the final black hole $A_{f,\rm expected}$ is derived from the initial masses and spins measured from the pre-merger analysis. These are converted into an estimated final mass and spin using fits from numerical relativity, which is then converted into an area via Eq. (7).

Since $A_{f,expected}$ is evaluated with numerical relativity using $A_{1}$ and $A_{2}$, the denominator of Eq. (8) is positive definite. Therefore, if $H<0$, then $A_{f,measured}<A_{i}$ and the area theorem is violated. The confidence interval in favor of the area theorem is the fraction of the posterior for which $H>0$. More robustly, $H=1$ indicates that signal is consistent with GR specifically, not just the more generic class of theories that satisfy the area increase law.

Figure 4 shows the posteriors of the final mass $M_{f}$ and final spin $\chi_{f}$ of the QNM post-merger models on real data. Specifically, the $M_{f}$ and $\chi_{f}$ posteriors from the analysis with variable sky/$t_{c}$ parameters are compared with the corresponding posteriors from fixed gate analysis. All posteriors contain within them the distribution from a full IMR analysis from 4-OGC [19]. Furthermore, both posteriors from the variable sky/$t_{c}$ run overlap almost completely with their corresponding fixed parameter posteriors.

A similar plot is shown in Fig. 5, except with the analyses that used an IMR post-merger model on real data. Again, the pre-merger and post-merger posteriors from the analysis with variable sky/$t_{c}$ parameters are compared to their fixed counterparts as well as the full IMR posterior from [19]. Notably, while still in very close agreement, the variable sky/$t_{c}$ posteriors have a greater discrepancy from their fixed counterparts than seen in Fig. 4. Specifically, the post-merger posterior tends towards slightly higher masses and spins, while the pre-merger posterior includes lower spin values.

Directly comparing the pre-merger posteriors in Figs.4 and 5 shows a slight discrepancy in final mass and spin estimates. Namely, the analysis with a QNM post-merger model tends towards higher final mass and spin estimates than the IMR post-merger model. This indicates minor coupling between the pre- and post-merger models caused by their shared sky and time parameters. However, as made evident by these two plots, this is not a major effect, and both pre-merger models maintain strong consistency with the full IMR posterior.

Figure 6 compares the $H$ posteriors between the four analyses of the real GW150914 data. All posteriors agree with the expected value $H=1$, with each distribution peaking near this value. While the QNM post-merger posteriors are almost in exact agreement with each other, there is a slight discrepancy between the IMR post-merger posteriors. Specifically, the fixed sky/$t_{c}$ posterior has a much sharper peak at a lower value than the variable gate counterpart. This may be explained by Fig. 5, where the fixed IMR post-merger posterior was observed to contain lower masses and spins than the variable gate model. This could lead to the slight bias to lower $H$ values seen here.

Table 3 summarizes the $H$ credible intervals and area theorem confidence intervals for all eight analyses. All credible intervals are in agreement with the expected value $H=1$. Additionally, all analyses of the real GW150915 signal have $H>0$ at greater than 95% confidence, indicating very high agreement with the area theorem. Results from the analysis of the simulated signal are largely consistent with these results (see Appendix C for equivalent figures).

The area theorem confidence interval of any given analysis appears to be uncorrelated to whether or not the sky location and coalescence time were allowed to vary. However, IMR post-merger models tended to have slightly higher agreement with the area theorem than corresponding QNM post-merger models. This pattern may be explained by trends in the pre- and post-merger mass and spin results. The IMR post-merger models in Fig. 5 had post-merger measurements that were either in agreement with or greater than the corresponding pre-merger estimates. This corresponds to a very high concentration of $H$ measurements approximating or exceeding 1, with very few points skewed to negative $H$ values. Conversely, the post-merger posteriors shown in Fig. 4 have significant amounts of points with higher and lower mass and spin values than their pre-merger counterparts. This leads to a wider distribution of $H$ values to both higher and lower values, allowing for a higher concentration of negative measurements.

This paper presented a novel method of marginalizing over sky location and coalescence time when performing a pre- or post-merger analysis, which allows for a full accounting of uncertainties in parameter estimates. Tests of the area theorem were conducted using data from GW150914 using this method. It was found that the data for this event agrees very well with the area theorem regardless of whether sky position and $t_{c}$ were allowed to vary. The only noticeable changes in agreement in the area theorem were caused by the post-merger approximant used, though in general this was a very minor change.

While we focus on tests of the area theorem here, our method for marginalizing over sky location and $t_{c}$ applies to any analysis in which only a portion of the signal is modelled. In particular, black hole spectroscopy involves analyzing the post-merger signal using a QNM template in order to perform a test of the no-hair theorem [24, 25, 26]. As with previous tests of the area theorem, previous black hole spectroscopy studies have fixed the sky location and $t_{c}$ for the reasons discussed in Sec. II [23, 27, 28]. Our finding that a signal model is needed for the entire observable signal applies equally well to black hole spectroscopy, even though the pre-merger signal is purposely excluded in such analyses. The hierarchical method we develop in Sec. III is therefore equally relevant for marginalizing over sky location and $t_{c}$ when doing black hole spectroscopy. In Ref. [29] we use our QNM analysis of GW150914 to investigate the evidence for the presence of the $(2,2,1)$ mode, which has been hotly contested in the literature [23, 27, 30, 5].

In all our analyses here we end the pre-merger analysis when the post-merger analysis beings. This is necessitated by our finding that some model must exist for the entire observable signal when the gate time is allowed to vary. Since there was no gap between our pre- and post-merger template, we obtained agreement with the area theorem in excess of $95\%$. A more rigorous test of the area theorem is to excise the merger from the data, since any biases introduced by including the merger as part of the pre-merger model would be omitted. Refs. [4, 2] did this additional test with fixed sky location and $t_{c}$. However, this is not possible when marginalizing over sky location and $t_{c}$ for the reasons highlighted in Sec. III. Introducing a gap between the pre- and post-merger will once again favor points that excise as much of the signal as possible, even when both the pre- and post-merger were modeled simultaneously.

Introducing a gap between the pre- and post-merger models could be achieved by using three sub-domains rather than two: one each for the inspiral, merger, and ringdown. The gate for one domain would start/end at the end/start of the next. For the merger domain, an arbitrary signal model using wavelets could be used, similar to what Finch and Moore used in Ref. [10]. This would ensure that pre- and post-merger parameters are fully unbiased by the merger while ensuring that the merger itself is not arbitrarily gated out. The initial black hole areas could then be estimated using an inspiral model and the final area using a QNM model. This would also be useful in black hole spectroscopy studies involving fundamental angular QNMs. These modes are not expected to become relevant until $\sim 10M$ after merger, necessitating a gap between the merger and the start of the QNM model. We plan to investigate this in a future study.

This research was conducted using PyCBC [31]. Our data is available at [32].

A.C. was supported by funds from the Massachusetts Space Grant Consortium. C.C. acknowledges support from NSF award PHY-2309356. All computations were performed on Unity, a collaborative, multi-institutional high-performance computing cluster managed by UMass Amherst Research Computing and Data.

This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation, as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.