Sky marginalization in black hole spectroscopy and tests of the area theorem
Abstract
Direct observation of gravitational waves from binary black hole (BBH) mergers has made it possible to test the laws of black-hole thermodynamics using real astrophysical sources. These tests rely on accurate and unbiased parameter estimates from the pre- and post-merger portions of a signal. Due to numerical complications, previous analyses have fixed the sky location and coalescence time when independently estimating the parameters of the pre- and post-merger signal. Here we overcome the numerical complications and present a novel method of marginalizing over sky location and coalescence time. Doing so, we find that it is not possible to model only the pre- or post-merger portions of the signal while marginalizing over timing uncertainty. We surmount this problem by simultaneously yet independently modelling the pre- and post-merger signal, with only the sky location and coalescence time being shared between the models. This allows us to marginalize over all parameters. We use our method to measure the change in area between the final and initial black holes in the BBH merger GW150914. To measure the final black hole’s area we do an analysis using quasi-normal modes (QNMs) to model the post-merger signal, and another analysis using the post-merger portion of an inspiral-merger-ringdown (IMR) template. We find excellent agreement with expectations from General Relativity. The Hawking area theorem (which states that ) is confirmed to and confidence using the QNM and IMR post-merger models, respectively. Both models yield , where is the expected change in area derived from fits to numerical relativity simulations.
I Introduction
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 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 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 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 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 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 . 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 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 . 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.
II Determinant approximation
PyCBC Inference [11] utilizes Bayesian inference to conduct parameter estimation on GW events. Bayes’ theorem is used to extract information about the parameter space from a data set [14]. Assuming the data set is composed of a signal and noise such that , the probability of observing a specific is given by
(1) |
The term is the likelihood, which describes the probability of observing a signal assuming the event has a parameter space . The term is the prior, a distribution that describes the a priori probability of observing given a signal model . The denominator is the evidence, a normalization factor used to compare analyses using different models for . 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 , whereas the likelihood is calculated directly from the data. For a system of detectors each taking time series samples of an event with a stochastic Gaussian noise background, the likelihood can be written as [15]
(2) |
where is the covariance matrix associated with in detector .
To evaluate the likelihood, further assumptions must be made to calculate . The result of these assumptions is:
(3) |
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 evaluated for a parameter space , Eq. (3) can be rewritten by substituting on the right-hand side, yielding
(4) |
(Here and throughout, we use to refer to the natural logarithm.) For a GW analysis that examines the entirety of , this expression is easily evaluated using the approximated eigenvalues of . However, many analyses only examine a portion of the time series (generally either the pre- or post-merger signal). Therefore, a region of is omitted, or “gated”, to reduce numerical biases due to the merger. Doing this also requires excising rows and columns from , thereby breaking the Toeplitz form of the matrix and the corresponding approximations. The inverse of the truncated covariance matrix 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 , which is a notoriously expensive numerical problem. The most widely-used matrix decomposition methods are [16], which become impractical for GW time series such as GW150914, where . This problem is amplified when varying sky position and 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 and the sky location of the event relative to the detectors. Therefore, when varying sky location and , the gate times in the detector frames will also vary. Subsequently, the elements of will vary with each unique sky location and value. Analyses that marginalize over these parameters require recalculating on multiple steps of the sampler, which would be impractical using numerical decomposition methods. This is not an issue if the sky location and are fixed, since in that case 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 it is necessary to have a fast method for calculating . We find that 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 . Fig. 1 shows this relationship using values calculated for various gate sizes applied to GW150914 data. Using a least squares fit to the points, the correlation coefficient was calculated as , which for a sample of 8 points implies a highly significant correlation [17].
![Refer to caption](extracted/5737897/figs/det_plots/results_H1_paper.png)
We also find that the position of the gate in the time series has no significant effect on the value of . The maximum range over which the determinant values varied was . However, since these values were generally of order , this represents a maximum fractional change of 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 . 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.
III Need for Joint Analyses
![Refer to caption](extracted/5737897/figs/inspiral_wf_comp.png)
![Refer to caption](extracted/5737897/figs/ringdown_wf_comp.png)
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 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 . When analyzing the entire time series this can only occur if the template 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 roughly 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 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 . In our case the common parameters are the right ascension , declination , and geocentric coalescence time . These parameters determine the coalescence time in each detector , 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 is:
-
1.
Generate the pre-merger (“inspiral”) template using parameters .
-
2.
Project the template into each detector’s frame using .
-
3.
Excise times after by gating and in-painting the residual using a s gate that spans .
-
4.
Calculate the pre-merger likelihood via Eq. (4).
-
5.
Repeat steps 1-4 for the post-merger (“ringdown”) template to get the post-merger likelihood . However, for the ringdown the gate spans ; i.e., it ends at whereas it starts at for the pre-merger. The data used in each domain is thereby (nearly) mutually exclusive (see Appendix B for more details).
-
6.
The total likelihood is then
(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 , both models generated independent parameter measurements. The resultant waveform templates closely match those obtained by fixing sky location and 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.
IV Methods
We performed eight analyses in total. Each analysis was differentiated by the waveform used, post-merger approximant, and whether or not sky location and 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 mode and a subdominant mode as proposed by [23]. All models apply the ringdown model starting at merger time . The priors of the QNM post-merger models were restricted such that the 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 mode erroneously matching to the mode in the signal.
Only half of the models allowed for sky location and to vary. The other half fixed these parameters to nominal values for GW150914 to replicate previous works. Specifically, the fixed parameter analyses set , , and , in accordance with the maximum likelihood values in [23].
Parameter | Description | Prior Dist. | Prior Range |
---|---|---|---|
Coalescence time | uniform | s | |
Right ascension | uniform | ||
Declination | sine angle | ||
Inclination | sine angle | ||
Source frame chirp mass | , | ||
Mass ratio | , | ||
Spin magnitude | uniform | ||
Spin polar angle | solid angle | ||
Spin azimuthal angle | solid angle | ||
Reference phase | uniform | ||
Comoving volume | uniform | Mpc3 |
Parameter | Description | Prior Dist. | Prior Range |
---|---|---|---|
Source frame final mass | uniform | ||
Final spin | uniform | ||
Initial mode amplitude | |||
mode phase | uniform | ||
Initial mode amplitude (as ratio of ) | uniform | ||
mode phase | uniform |
As a preliminary test, the sky location posteriors of the analyses with variable sky/ 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].
![Refer to caption](extracted/5737897/figs/sky_real_compare.png)
V Area theorem
The simplest area theorem test is to compare the progenitor horizon areas and to the remnant horizon area to check that
(6) |
In natural units () the area of each black hole is given by [1]:
(7) |
where the black hole’s mass and 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]
(8) |
Here, is the area of the final black hole inferred from the post-merger analysis and is the sum of the initial areas inferred from the pre-merger analysis. The expected area of the final black hole 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 is evaluated with numerical relativity using and , the denominator of Eq. (8) is positive definite. Therefore, if , then and the area theorem is violated. The confidence interval in favor of the area theorem is the fraction of the posterior for which . More robustly, indicates that signal is consistent with GR specifically, not just the more generic class of theories that satisfy the area increase law.
![Refer to caption](extracted/5737897/figs/area_thrm/real_data/chi_m_real_QNM_compare.png)
Figure 4 shows the posteriors of the final mass and final spin of the QNM post-merger models on real data. Specifically, the and posteriors from the analysis with variable sky/ 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/ run overlap almost completely with their corresponding fixed parameter posteriors.
![Refer to caption](extracted/5737897/figs/area_thrm/real_data/chi_m_real_IMR_compare.png)
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/ 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/ 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 posteriors between the four analyses of the real GW150914 data. All posteriors agree with the expected value , 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/ 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 values seen here.
![Refer to caption](extracted/5737897/figs/area_thrm/real_data/area_thrm_real_compare.png)
Table 3 summarizes the credible intervals and area theorem confidence intervals for all eight analyses. All credible intervals are in agreement with the expected value . Additionally, all analyses of the real GW150915 signal have 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 measurements approximating or exceeding 1, with very few points skewed to negative 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 values to both higher and lower values, allowing for a higher concentration of negative measurements.
GW150914 analysis | Credible interval | |
---|---|---|
Variable sky/, QNM post-merger, real waveform | 95.4% | |
Variable sky/, IMR post-merger, real waveform | 99.5% | |
Fixed sky/, QNM post-merger, real waveform | 97.4% | |
Fixed sky/, IMR post-merger, real waveform | 98.4% | |
Variable sky/, QNM post-merger, zero-noise injection | 89.1% | |
Variable sky/, IMR post-merger, zero-noise injection | 98.4% | |
Fixed sky/, QNM post-merger, zero-noise injection | 97.3% | |
Fixed sky/, IMR post-merger, zero-noise injection | 98.8% |
VI Conclusions
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 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 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 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 when doing black hole spectroscopy. In Ref. [29] we use our QNM analysis of GW150914 to investigate the evidence for the presence of the 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 . 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 . However, this is not possible when marginalizing over sky location and 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 after merger, necessitating a gap between the merger and the start of the QNM model. We plan to investigate this in a future study.
VII Acknowledgements
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.
Appendix A Likelihood Function Derivation
To calculate the posterior from Bayes’ theorem (Eq. (1)), one requires a model for both the signal and the noise . While can be determined from the linearized EFE wave solution, requires more statistical considerations.
To start, consider a network of gravitational wave detectors that each sample at a rate over a total time length . Defining the total number of samples , the full data series can be expressed as a series of -dimensional vectors such that and . To simplify calculations, assume that the signal is zero such that . Additionally, assume that noise model is a stochastic Gaussian distribution and is uncorrelated between detectors. Under these assumptions, the noise likelihood function is:
(9) |
where is the covariance matrix of the noise model for detector , defined using the ensemble average:
(10) |
(By dropping the subscripts, one obtains an equivalent expression for the full covariance matrix and data. We do so in the following steps for brevity.) This is the exact expression of the likelihood function for noise. However, this function is infeasible to calculate analytically due to the inverse covariance matrix in the numerator.
To do this, one may expand the covariance matrix definition in Eq. (10) as follows, defining :
(11) |
where in the last step the ensemble average is written out fully.
In general, this expression is dependent on time and displacement . However, one can make the assumption that the noise is wide sense stationary, where the mean and variance are both constant in time. Under this assumption, any constant can be added to the indices in Eq. (10) to obtain the same result. This makes symmetric (since the factors in Eq. (11) commute) and Toeplitz (since the elements along the diagonals are equal) [33]. Additionally, since , the elements of are even functions of .
Additionally, one can assume that the data is ergodic, meaning that new realizations of are obtained via time. Under this assumption and the properties of the elements of , the ensemble averages in Eq. (11) can be replaced with time averages:
(12) |
The last step defines the autocorrelation function , which describes the correlation between points in the time series . If goes to zero in some finite time , then all diagonals with will equal zero. This is similar to the form of a circulant matrix , a special case of a Toeplitz matrix where each row is a right-cycle permutation of the same vector (in this case, ). The eigenvectors of circulant matrices are well-known [33]:
(13) |
This generally is not true for Toeplitz matrices, but one may take advantage of Eq. (13) by recognizing that the matrix described by Eq. (12) asymptotes to a circulant matrix for large N:
(14) |
Therefore, the eigenvalues of can be evaluated using the usual eigenvalue equation as long as :
(15) |
Therefore, using the fact that is symmetric and is even:
(16) |
where is the discrete Fourier transform of the autocorrelation function:
(17) |
To simplify Eq. (16), one may impose the Wiener-Khinchin Theorem [34], which defines the power spectral density as the Fourier transform of for a wide-sense stationary stochastic process:
(18) |
One can then construct the inverse of using the eigenvalues from Eq. (18) and the eigenvectors of Eq. (13):
(19) |
using the fact that is symmetric about and defining the sample frequency . Since is real and
(20) |
one can write:
(21) |
Finally, if the inner product between two arbitrary vectors and is defined as:
(22) |
Eq. (21) can be written as an inner product:
(23) |
and the likelihood function can be expressed as:
(24) |
Appendix B Gating and In-Painting
Generally, BBH models do not require the entire waveform to be analyzed at once. For example, the analyses throughout this paper independently examine the pre- and post-merger portions of the GW150914 waveform. To maintain independence between the models, any points not corresponding to the respective model (i.e. after for the pre-merger model, or before for the post-merger model) were excised, or “gated,” from the data. Here, a gate of length applied starting at a sample will excise all samples within the range , corresponding to a gate of time length . The gated time series will therefore take the form .
By doing this, the simplifications made to derive the likelihood function are no longer valid, since is no longer Toeplitz (and, subsequently, no longer approximately circulant for large ). There are numerical methods to calculate the matrix inverse directly, but in general they can be unstable and time-intensive. Therefore, a method known as “gating and in-painting” is employed for gated waveform analyses [12].
First, the method assumes that the noise time series is the sum of , the noise series with the gated times zeroed out, and , a time series containing only the gated samples in . The goal of in-painting is to solve the following equation in the gated region:
(25) |
If the nonzero elements of are such that for all samples in the gate , then the inner product will be equal for the truncated and raw data set. Since is zero outside of the gate, will form an Toeplitz matrix containing the rows and columns of . Therefore, Eq. (25) can be rewritten within the gated region as
(26) |
and adding to will give the same result as truncating and . Unlike trying to solve for the inverse directly, this solution is readily found using a Toeplitz solver [18, 35]. Given gated data containing some gated signal , Eq. (26) can be evaluated using , and the value can be used to calculate the likelihood.
The analyses conducted in this paper utilize gating and in-painting to apply a gate starting/ending at to the waveform template. In theory, the gates should extend to the edges of the analysis segment (i.e., the pre-merger gate starts at the segment start time, and the post-merger gate ends at segment end). However, the in-painting algorithm, which is the dominant cost in our analysis, is operations for a gate of samples. Although the observable signal is only s long [13], we use an analysis segment that is s in duration in order to resolve line artifacts in the PSD. We also use a sample rate of Hz, in order to fully capture all observable signal power. The in-painting would therefore need to span samples if the gates were to extend to the beginning/end of the analysis segment. This is computationally expensive: some of our analyses would take month to complete (utilizing 64 CPU cores).
To reduce the computational cost, we instead use a 1-second gate starting (ending) at and ensuring that the remainder of the pre-merger (post-merger) template after (before) the gate is zero. This is equivalent to applying a full gate to the template when a whitening filter is applied, since the early- and late-time template will be identically zero in both cases.
We do not zero the data outside of the gates, however, since doing so would introduce additional ringing artifacts at the beginning/end of the analysis segment not accounted for by the in-painting. This means that the pre- and post-merger models will share data outside of the gates (namely, greater than 1 second before and after ). However, the data is expected to be noise-dominated in this region, yielding an insignificant contribution to the likelihood. This can be seen in Appendix C, which portrays the remnant posteriors for analyses using zero-noise waveform injections. Since the data in these analyses contain no noise, the shared data has identically zero contribution to the likelihood. These results are in strong agreement with Figs. 4, 5, and 6, indicating that any effects due to early- and late-time shared noise are indeed negligible.
Appendix C Simulation results
This section contains additional plots showing results from the analysis of the GW150914-like simulated signal in zero noise. Figure 7 shows the and posteriors for the QNM post-merger models, and Fig. 8 shows the same for the IMR post-merger models. Figure 9 shows the posteriors for the zero-noise injection models.
![Refer to caption](extracted/5737897/figs/area_thrm/injection/chi_m_inject_QNM_compare.png)
![Refer to caption](extracted/5737897/figs/area_thrm/injection/chi_m_inject_IMR_compare.png)
![Refer to caption](extracted/5737897/figs/area_thrm/injection/area_thrm_inj_compare.png)
![Refer to caption](extracted/5737897/figs/qnm_convergence.png)
![Refer to caption](extracted/5737897/figs/imr_convergence.png)
Appendix D Sampler convergence
To test the convergence of the sampler we repeated the analyses with 2000 and 4000 live points. Figure 10 shows the and posteriors for QNM post-merger models on real data with 2000 and 4000 live points, while Fig. 11 depicts the same for the IMR post-merger models. Both the QNM and IMR models were able to converge to similar distributions regardless of the number of live points used. All sky/ results are reported using 4000 live points.
References
- Bardeen et al. [1973] J. M. Bardeen, B. Carter, and S. W. Hawking, The four laws of black hole mechanics, Communications in Mathematical Physics 31, 161 (1973).
- Kastha et al. [2022] S. Kastha, C. D. Capano, J. Westerweck, M. Cabero, B. Krishnan, and A. B. Nielsen, Model systematics in time domain tests of binary black hole evolution, Physical Review D 105, 10.1103/physrevd.105.064042 (2022).
- Cabero et al. [2018] M. Cabero, C. D. Capano, O. Fischer-Birnholtz, B. Krishnan, A. B. Nielsen, A. H. Nitz, and C. M. Biwer, Observational tests of the black hole area increase law, Physical Review D 97, 10.1103/physrevd.97.124069 (2018).
- Isi et al. [2021] M. Isi, W. M. Farr, M. Giesler, M. A. Scheel, and S. A. Teukolsky, Testing the black-hole area law with GW150914, Physical Review Letters 127, 10.1103/physrevlett.127.011103 (2021).
- Wang et al. [2023] Y.-F. Wang, C. D. Capano, J. Abedi, S. Kastha, B. Krishnan, A. B. Nielsen, A. H. Nitz, and J. Westerweck, A frequency-domain perspective on GW150914 ringdown overtone, (2023), arXiv:2310.19645 [gr-qc] .
- Carullo et al. [2019] G. Carullo, W. D. Pozzo, and J. Veitch, Observational black hole spectroscopy: A time-domain multimode analysis of GW150914, Physical Review D 99, 10.1103/physrevd.99.123029 (2019).
- Carullo et al. [2023] G. Carullo, W. Del Pozzo, and J. Veitch, pyRing: a time-domain ringdown analysis python package, git.ligo.org/lscsoft/pyring (2023).
- Isi and Farr [2021a] M. Isi and W. M. Farr, maxisi/ringdown: Initial ringdown release (2021a).
- Isi and Farr [2021b] M. Isi and W. M. Farr, Analyzing black-hole ringdowns, (2021b), arXiv:2107.05609 [gr-qc] .
- Finch and Moore [2022] E. Finch and C. J. Moore, Searching for a ringdown overtone in GW150914, Physical Review D 106, 10.1103/physrevd.106.043005 (2022).
- Biwer et al. [2019] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals, Publications of the Astronomical Society of the Pacific 131, 024503 (2019).
- Zackay et al. [2021] B. Zackay, T. Venumadhav, J. Roulet, L. Dai, and M. Zaldarriaga, Detecting gravitational waves in data with non-stationary and non-gaussian noise, Physical Review D 104, 10.1103/physrevd.104.063034 (2021).
- Abbott et al. [2016] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Observation of gravitational waves from a binary black hole merger, Physical Review Letters 116, 061102 (2016).
- Sivia and Skilling [2006] D. Sivia and J. Skilling, Data analysis: a Bayesian tutorial (Oxford University Press, Oxford, 2006).
- Finn [1992] L. S. Finn, Detection, measurement, and gravitational radiation, Physical Review D 46, 5236 (1992).
- Pan et al. [1997] V. Pan, Y. Yu, and C. Stewart, Algebraic and numerical techniques for the computation of matrix determinants, Computers & Mathematics with Applications 34, 43 (1997).
- Taylor [1996] J. R. Taylor, An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements, 2nd ed. (University Science Books, 1996).
- et al. [2020a] P. V. et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17, 261 (2020a).
- Nitz et al. [2022] A. H. Nitz, S. Kumar, Y. Wang, S. Kastha, S. Wu, M. Schäfer, R. Dhurkunde, and C. D. Capano, 4-OGC: Catalog of gravitational waves from compact-binary mergers (2022), arXiv:2112.06878 [astro-ph.HE] .
- Speagle [2020] J. S. Speagle, dynesty: a dynamic nested sampling package for estimating Bayesian posteriors and evidences, Monthly Notices of the Royal Astronomical Society 493, 3132 (2020).
- et al. [2021] R. A. et al., Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo, SoftwareX 13, 10.1016/j.softx.2021.100658 (2021).
- Pratten et al. [2021] G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, and S. Husa, Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Physical Review D 103, 10.1103/physrevd.103.104056 (2021).
- Isi et al. [2019] M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, Testing the no-hair theorem with GW150914, Physical Review Letters 123, 10.1103/physrevlett.123.111102 (2019).
- Vishveshwara [1970] C. V. Vishveshwara, Scattering of Gravitational Radiation by a Schwarzschild Black-hole, Nature 227, 936 (1970).
- Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. L. Detweiler, The quasi-normal modes of the Schwarzschild black hole, Proc. Roy. Soc. Lond. A 344, 441 (1975).
- Dreyer et al. [2004] O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, Black hole spectroscopy: Testing general relativity through gravitational wave observations, Class. Quant. Grav. 21, 787 (2004), arXiv:gr-qc/0309007 .
- Cotesta et al. [2022] R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, Analysis of ringdown overtones in GW150914, Physical Review Letters 129, 10.1103/physrevlett.129.111102 (2022).
- Abbott et al. [2021] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), Tests of General Relativity with GWTC-3, (2021), arXiv:2112.06861 [gr-qc] .
- Correia et al. [2023] A. Correia, Y.-F. Wang, and C. D. Capano, Low evidence for ringdown overtone in GW150914 when marginalizing over time and sky location uncertainty, (2023), arXiv:2312.14118 [gr-qc] .
- Isi and Farr [2022] M. Isi and W. M. Farr, Revisiting the ringdown of GW150914 (2022), arXiv:2202.02941 [gr-qc] .
- Nitz et al. [2024] A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, T. Dent, L. Pekowsky, G. S. C. Davies, S. De, M. Cabero, S. Wu, A. R. Williamson, B. Machenschalk, D. Macleod, F. Pannarale, P. Kumar, S. Reyes, dfinstad, S. Kumar, M. Tápai, L. Singer, P. Kumar, veronica villa, maxtrevor, B. U. V. Gadre, S. Khan, S. Fairhurst, and A. Tolley, gwastro/pycbc: v2.3.3 release of pycbc (2024).
- Correia and Capano [2024] A. Correia and C. Capano, Sky marginalization in black hole spectroscopy (2024).
- Gray [2001] R. Gray, Toeplitz and circulant matrices: a review, Foundations and Trends in Communications and Information Theory 2 (2001).
- Chatfield [1996] C. Chatfield, The analysis of time series: an introduction, 6th ed. (CRC Press, Florida, US, 1996).
- et al. [2020b] C. R. H. et al., Array programming with NumPy, Nature 585, 357 (2020b).