During the first two observing runs of Advanced LIGO and Advanced Virgo (O1 and O2), 11 gravitational wave (GW) signals were identified by the LIGO-Virgo Collaboration [1,2] including 10 from binary black hole (BBH) mergers [3]. Substantial improvements in detector sensitivity prior to the first half of the third observing run (O3a) led to the detection of 39 new GW event candidates, at least 36 of which are consistent with BBH mergers [4]. This growing population of BBH mergers allows us to probe underlying black hole formation channels [5,6].
Black holes can form directly through stellar collapse, mergers, accretion, or through collapse of dense gas in the early Universe (primordial black holes). Binary systems consisting of black holes which formed from stellar collapse are referred to as stellar-mass BBHs. In the case of stellar evolution, a mechanism known as pair-instability (PI) supernovae is expected to prevent the formation of heavier black holes [7][8][9][10][11]. In PI, the production of positron-electron pairs causes a reduction in radiation pressure. This in turn causes the core of the star to contract, raising the temperature and increasing the production of positron-electron pairs in a runaway process. Any star with a helium core mass 32 M ⊙ ≲ M HE ≲ 64 M ⊙ undergoes pulsational pair-instability (PPI). In PPI, the stellar envelope is removed and the mass of the star is reduced until the star reaches a stable state [7,[12][13][14]. Any star with a helium core mass 64 M ⊙ ≲ M HE ≲ 135 M ⊙ undergoes PI supernova which blasts the star apart, leaving no remnant compact object [7,15,16]. More massive stars are predicted to collapse directly to black holes [16][17][18].
The PI process is heavily dependent on stellar metallicity [19,20]. A metal-poor star (Z ≤ 3 × 10 -4 ) may collapse into a black hole with mass up to ∼65 M ⊙ if the star's hydrogen envelope falls onto the compact object and up to ∼45 M ⊙ otherwise [21]. Compared to a low metallicity star, a metal-rich star is expected to collapse into a lower mass compact object [8,16,[22][23][24][25][26][27]. Overall, this mechanism results in an absence of black holes in the mass spectrum known as the PI mass gap. Some theoretical uncertainties regarding the lower boundary remain [28][29][30][31][32], however, 65 M ⊙ is likely a conservative lower limit [8][9][10][11]16,33,34]. The upper boundary of the PI mass gap is more certain and is believed to be ∼120 M ⊙ [8,9,11].
Black holes can exceed the PI mass limit if they are formed from the previous merger of smaller black holes. This so-called hierarchical merger scenario can occur through dynamical encounters in dense stellar clusters such as galactic nuclei, globular clusters [35][36][37][38][39][40][41][42][43][44][45][46],o r nuclear star clusters [47,48]. Hierarchical mergers can also form through gas-capture in the disks of active galactic nuclei (AGNs) [49][50][51][52][53][54][55]. Alternatively, black holes can exceed the PI mass limit by substantial accretion, which is possible for black holes in AGN disks [56,57]. Black holes can also form within the PI mass gap when young stars within a dense young star cluster merge in a runaway collision process [58][59][60].
GW190521 is the heaviest BBH merger detected to date and marks the first direct detection of an intermediate mass black hole [61]. Reference [34] estimates the probability for at least one of the component black holes to lie within 65-120 M ⊙ to be 99.0%, assuming a binary system with quasicircular orbit. This probability increases when the analysis was repeated after taking into account the population of detected binary systems [6]. Hence, GW190521 is an appealing candidate for alternative formation scenarios-a prospect which has already been explored in various papers [62][63][64][65][66][67][68][69]. We caution, however, that analyzing the features of a potential outlier BBH event outside the context of the entire population of detected events can be misleading (a notion previously examined by Fishbach et al. in Ref. [70]).
Here, we investigate whether GW190521 or any other BBH merger is an outlier to the stellar-mass BBH population. Distinguishing between separate alternative formation channels is outside the scope of this work. We examine the BBH events detected by Coherent WaveBurst (cWB), a model waveform independent pipeline used to search for GW signals [71,72], and compare this sample of detected events against models of the stellar-mass BBH population. In the method presented below, we add simulated BBH waveforms (injections) directly into the GW detector strain data and recover them with the cWB pipeline. This allows us to directly account for detector noise when estimating the significance of potential outliers. Ultimately, we seek to answer two fundamental questions: (i) accounting for statistical fluctuations, are there any outliers to the stellar-mass BBH population? and (ii) is the distribution of detected BBH events consistent with the stellar-mass BBH population?
This paper is organized as follows. In Sec. II,w e familiarize the reader with cWB-a pipeline that uses minimal assumptions to identify GW transient signals.
Next, we specify the stellar-mass BBH population models and injections used in this study. Then, we define our test statistic and describe our statistical procedure used to classify outliers. In Sec. III, we present our results for various possible models of the stellar-mass BBH population. Finally, in Sec. IV, we state the conclusions of this study and briefly discuss our future work.
CWB is a search algorithm that identifies GW signals in the LIGO-Virgo detector data without directly using a waveform model [71,72]. The signal detection process is performed in the time-frequency domain using the Wilson Daubechies Meyer wavelet transform [73]. Here, the algorithm identifies wavelets inconsistent with detector noise and assembles them into clusters. If a cluster of wavelets is coincident in time and frequency across the detector network, then an event trigger is generated. For each event trigger, summary statistics are calculated which describe the time-frequency evolution (f 0 , ΔT s , ΔF s , M, F M , e M ), signal strength and coherence (SNR net , c c , n f ), and the likelihood fit (E c =L). Additional statistics are used to estimate the number of cycles in the reconstructed waveform (Q 0 , Q 1 ). The summary statistics used in this study are defined in Appendix C.
Through O1, O2, and O3a, cWB has contributed to the detection of 22 BBH merger candidates (displayed in Table I). Among these detections is GW190521, which cWB played a crucial role in identifying [61,74]. In this work, we compare these BBH detections against predictions made by the stellar-mass BBH population. Simulated waveforms representative of the stellar-mass BBH population are injected into the GW detector data and recovered with cWB. We restrict our analysis to events detected by cWB with a false alarm rate (FAR) less than 1 yr -1 to maintain a high purity sample set. Although the inclusion of Virgo generally improves source property reconstruction, we simplify our analysis by considering only the Hanford-Livingston detector network.
We examine two phenomenological mass models presented in Ref. [6]: the Broken Power Law model [75] and the Power Law þ Peak model [76]. Each population model is described by a set of hyperparameters whose values are estimated by applying hierarchical Bayesian inference on the full population of BBH merger candidates detected by LIGO-Virgo [77][78][79]. Below is a short description of each model:
(i) Broken Power Law: The primary mass distribution is characterized by a gentle sloping power law at low mass designed to capture the observed Salpeter initial mass function [80,81] and a second steep sloping power law at high mass which targets a mass tapering of the distribution. The mass ratio is given as a power law that favors equal mass binaries [82]. (ii) Power Law þ Peak: The primary mass distribution is characterized by a power law component at low mass and wide Gaussian component at high masses which models a build-up of black holes due to PPI supernovae [83,84]. Similar to the previous model, the mass ratio is given as a power law favoring equal mass binaries. The left plot in Fig. 1 displays the astrophysical distribution for the primary mass m 1 predicted by each model, averaged over uncertainty. These mass models describe the complete astrophysical distribution of BBHs and do not distinguish between separate formation channels. To isolate the stellar-mass BBH population, we apply a theoretically motivated cutoff on the primary mass to model the PI mass gap lower boundary. We assume that no stellar-mass BBHs form inside of the PI gap. To account for the uncertainty surrounding the PI mass limit, we consider three cases: (i) m 1 < 65 M ⊙ , (ii) m 1 < 50 M ⊙ , and (iii) m 1 < 40 M ⊙ . For brevity, we abbreviate each model and add a subscript indicating the primary mass cutoff, e.g., BPL 65 indicates the Broken Power Law model with a primary mass cutoff of m 1 < 65 M ⊙ .
An essential part of this work is that we add injections directly to GW detector strain data. This allows us to compare our detected BBH events against the population models, accounting for selection bias. We use the SEOBNRv4 waveform approximant [85] and inject simulated waveforms uniformly in time throughout the O1-O3a observing runs. In the case of O1 and O2, we use the public GW strain data available through the Gravitational Wave Open Science Center (GWOSC) [86]. At the time of this publication, the GW strain data for the complete O3a observing run was not publicly available. Instead, we used the segments of data available from the GWTC-2 GWOSC data release [86] to emulate the O3a sensitivity, and we repeated injection analysis on these data segments until realizing the full O3a observation time.
The component masses for each simulated BBH event are drawn from the tested stellar-mass model under consideration (see Sec. II B). Component black hole spins are drawn from a uniform aligned spin distribution. This spin distribution is theoretically motivated to favor isolated binaries [87][88][89][90][91][92]. The binary source distance is drawn assuming uniform density in comoving volume up to redshift z ¼ 2 [93]. Other extrinsic binary parameters including sky location and inclination angle are randomly selected.
We recover the simulated BBH events using the cWB pipeline and select events with a FAR less than 1 yr -1 . The right plot in Fig. 1 astrophysical distributions, the detected distributions favor more massive BBH mergers due to selection effects caused by the ground-based detectors and the cWB search pipeline.
In this study, we require a test statistic that captures the characteristics of an outlier to the stellar-mass BBH population. One logical statistic to consider is the source frame primary mass m 1 , since stellar-mass BBH events are not expected to have m 1 inside the PI mass gap. Ideally, one could use the source frame primary mass m 1 estimated via Bayesian inference. However, since we inject and recover Oð100;000Þ simulated BBH events in this study, this is not a feasible approach.
Instead, we utilize cWB in combination with a supervised machine learning (ML) algorithm XGBoost [94] to rapidly estimate the source frame total mass for a given event. This estimate, which we call M cWB , serves as our test statistic. The injected source frame total mass M inj is designated as the target variable in our ML regression model, and the summary statistics reconstructed by cWB are used as input features to estimate the target variable. Figure 2 shows the injected M inj and estimated M cWB for simulated BBH events drawn from the Power Law þ Peak mass model. In Appendix C.4 of Ref. [6], we use a simplified approach without ML, but we expect the method presented in this paper to be more sensitive to measure outliers of the stellar-mass BBH population.
To generate a robust estimate of the source frame total mass, we train our ML algorithm on a simulation set independent from the stellar-mass BBH population simulation sets. This way, the trained ML algorithm is conservative when applied to tested population models. We use the IMRPhenomPv2 waveform approximant [95] and inject simulated waveforms uniformly throughout O1, O2, and our emulated O3a observation run. For each training event, the primary mass m 1 is drawn from a uniform distribution m 1 ∈ ½1; 150 M ⊙ , and the mass ratio q is fixed equal to one. The component spins are drawn from an isotropic distribution. Extrinsic parameters including sky location, orientation, and distance are each drawn from uniform distributions. We again consider only events detected by cWB with FAR less than 1 yr -1 . More information regarding our ML training and tuning procedure can be found in Appendix D.
Our statistical procedure is divided into two steps. First we identify which events, if any, in our sample set are potential outliers to the stellar-mass BBH population. And second, we remove any potential outliers and determine if the remaining distribution of events is consistent with the stellar-mass BBH population. We use M cWB , defined in Sec. II D, as an outlier test statistic where a higher value of M cWB indicates a greater probability for the event to have originated inside the PI mass gap.
First, to estimate the significance of potential outliers, we consider a bootstrap resampling procedure. In this procedure, we identify the event x which has the highest value of M cWB from our sample set of N event candidates. We then randomly draw N simulated events with replacement from the model distribution and select the simulated event from the set which has the highest value of M cWB , i.e., we select the outlier from the simulated set. We repeat this routine 10,000 times, and thus estimate the distribution of outliers events O predicted by the model. So, the probability for a model to produce an outlier as significant as event x is given as:
where H 0 represents the null hypothesis that event x is consistent with the model outlier distribution O. This is the one-sided (right tail) p-value. If this probability is less than our conservatively chosen threshold α ¼ 0.01, we claim that event x is a potential outlier to the tested population model.
In the case where event x is found to be a potential outlier, we omit the event from the sample set and reanalyze the remaining N -1 events. We repeat this process until all outliers are identified.
Once we have identified the number of population outliers N out , we then study whether the remaining sample set of N -N out events are consistent with the population model. We again consider a bootstrap resampling procedure where the BBH events in our sample set are ordered according to M cWB . We then randomly draw 10,000 batches of N -N out simulated events and sort each batch according to M cWB . The nth ranked event from our sample set is compared against the nth ranked simulated event from each batch. Thus, we can identify if the nth ranked event in our sample set is consistent with predictions made by the tested population model.
To further quantify the consistency between the sample distribution and tested population model, we apply the Anderson-Darling (AD) goodness-of-fit test [96]. Compared to the Kolmogorov-Smirnov (KS) test [97,98], the AD test is expected to be more sensitive at the tails of the distribution and so is more suitable for this analysis. If the AD p-value p AD is less than β ¼ 0.01, then we find the sample distribution is inconsistent with the tested model distribution.
In this section, we present the results of our analysis under three separate assumptions of the PI mass limit. First, we examine stellar-mass BBH models under the assumption of a 65 M ⊙ PI mass limit: BPL 65 and PLP 65 , defined in Sec. II B. Our sample distribution consists of 22 BBH events detected by cWB (see Table I). We order these events according to M cWB . Here, the first rank event is GW190521 with M cWB ¼ 193.2 M ⊙ .
We compare GW190521 against outliers predicted by each model to determine if the event is a stellar-mass BBH population outlier. Specifically, we generate batches of 22 simulated events from the BPL 65 and PLP 65 models and compare GW190521 against the outlier from each batch. The rotated histograms in Fig. 3 show M cWB for GW190521 (red dotted line) compared against the distribution of M cWB for outliers predicted by BPL 65 (orange) and PLP 65 (blue). The dashed lines in each histogram denote the 98% confidence predicted by the model. Using Eq. ( 1), we find the probability for the BPL 65 and PLP 65 models to produce an outlier as significant as GW190521 to be less than 0.0001. So, assuming the PI mass limit exists as a sharp cutoff at 65 M ⊙ , we find GW190521 to be an outlier to the stellar-mass BBH population.
Next, we investigate if any other events are outliers to the stellar-mass BBH population. We omit GW190521, and reanalyze the remaining 21 detected BBH events in our sample set against batches of 21 simulated BBH events predicted by each model. The highest ranked event in this sample set is GW190706_222641 with M cWB ¼ 139.2 M ⊙ . After applying the same statistical procedure as above, we calculate a p-value of 0.109 and 0.082 for the BPL 65 and PLP 65 model, respectively, and so this event is consistent with both population models.
After identifying all potential population outliers, in this case just GW190521, we examine whether the remaining sample distribution is consistent with the tested population models. Figure 3 shows the estimated M cWB for our sample set (black) compared to predictions made by the BPL 65 (orange) and PLP 65 (blue) models. The dark orange (blue) line represents the median value of M cWB predicted by the BPL 65 (PLP 65 ) model, the shaded region indicates 90% confidence interval, and the dashed lines indicate 98% confidence interval. We find that all 21 events lie within the 98% confidence interval predicted by each model. Applying the AD goodness-of-fit test, we estimate p AD for the BPL 65 model and PLP 65 model to be greater than 0.25. So, the distribution of detected events, excluding GW190521, is consistent with both stellar-mass BBH models.
Next, we repeat our statistical analysis for stellar-mass BBH models assuming a PI mass limit of 50 M ⊙ :B P L 50 and PLP 50 . Here, we again compare the M cWB distribution of the 22 BBH events detected cWB (our sample distribution) against the M cWB distribution predicted by the two models.
As in Sec. III A, we find GW190521 to be an outlier to the tested stellar-mass BBH population models (p-value < 0.0001). So, we reanalyze the remaining 21 BBH detections. Here, the highest ranked event from this set, GW190706_222641, has an outlier p-value of 0.009 when compared to the BPL 50 model and an outlier p-value of 0.011 when compared to the PLP 50 model. By definition of our a priori p-value threshold α ¼ 0.01, we treat this event as a potential outlier to the BPL 50 model but not the PLP 50 model. No other potential outliers are identified for the BPL 50 stellar-mass model. In Fig. 4, the rotated histograms show the estimated M cWB of potential outliers compared to predictions made by each population model. Here, it is FIG. 3. Distribution of BBH events detected by cWB, shown in red/black, ranked according to M cWB and compared to predictions made by stellar-mass models BPL 65 (top) and PLP 65 (bottom), which assume a PI mass limit of 65 M ⊙ . Event names are abbreviated. Shaded region indicates 90% confidence interval, dashed lines indicate 98% confidence interval. GW190521 is identified as an outlier to both population models. The remaining BBH events are consistent with both models. evident that GW190521 is a clear outlier to both models, whereas GW190706_222641 lies just outside the 98% confidence interval for the BPL 50 model.
Next, we examine the remaining distribution of detected events. The main panels in Fig. 4 show this distribution compared to predictions made by each model. Here, we observe that the sample distribution appears to deviate from the predicted distributions, particularly the PLP 50 model distribution. That is, there are more high mass events in the sample distribution compared to the predictions made by each model. To quantify this deviation, we apply the AD test. For the BPL 50 stellar-mass model, we analyze the 20 nonoutlier BBH detections and calculate p AD to be 0.039so the detected events are consistent with the model. For the PLP 50 model, we analyze the 21 nonoutlier BBH detections, including GW190706_222641, and calculate p AD to be 0.006. This indicates that the sample distribution is inconsistent with the PLP 50 model. As a cross-check, we examine the case of evaluating GW190706_222641 as a potential outlier, since it barely passes our outlier p-value threshold. Applying the AD test for this case, we calculate p AD to be 0.036. So, if GW190706_222641 is treated as a potential outlier, then the remaining sample distribution is consistent with PLP 50 .
In the final case, we consider a PI mass limit of 40 M ⊙ . Our results for this case are summarized in Fig. 5. As in the previous two cases, we first compare GW190521 against the distribution of outliers predicted by each model. We again find GW190521 to be an outlier with a p-value less than 0.0001 for both models (leftmost rotated histogram). We then reanalyze the remaining 21 BBH events in our sample set to test if the highest ranked event in this set, GW190706_222641, is an outlier to the stellar-mass BBH models. For this event, we calculate a p-value of 0.007 and 0.006 for the BPL 40 and PLP 40 models, respectively (second rotated histogram). So, this event is also classified as a potential outlier. We then rerun our analysis for the remaining 20 detected BBH events and find no additional outliers in the sample set.
The main panels in Fig. 5 show the sample distribution, according to M cWB , compared to predictions made by each stellar-mass BBH population model. Here, we notice that more high mass events are found compared to the predictions made by either population model. Applying the AD test, we find p AD for both the BPL 40 model and the PLP 40 model to be less than 0.001. So, if the PI mass limit is at 40 M ⊙ , then we have two potential outliers-GW190521 and GW190706_222641-and out of the remaining 20 BBH events, there are more high mass events than predicted by the models.
In this paper, we introduce a method to determine if the distribution of BBH events detected by cWB is consistent with the stellar-mass BBH population. Detector noise is unavoidable and could affect the inference and FIG. 4. Distribution of BBH events detected by cWB, shown in red/black, compared to predictions made by stellar-mass models BPL 50 (top) and PLP 50 (bottom), which assume a PI mass limit of 50 M ⊙ . For the BPL 50 model, GW190521 and GW190706_222641 are evaluated as potential outliers, and the remaining events are consistent with the population model. For the PLP 50 model, only GW190521 is identified as an outlier, while the remaining distribution of events are inconsistent with the model. FIG. 5. Distribution of BBH events detected by cWB, shown in red/black, compared to predictions made by stellar-mass models BPL 40 (top) and PLP 40 (bottom). Two events-GW190521 and GW190706_222641-are found to be potential population outliers. The remaining 20 BBH events are inconsistent with both stellar-mass BBH population models.
understanding of any given event, thereby making an otherwise ordinary BBH merger appear as an outlier. Our analysis operates directly on search statistics by injecting simulated waveforms representative of the stellar-mass BBH population directly into the detector strain data and recovering them with the search pipeline. This way, any inference errors caused by the detector noise are captured within the model distribution.
The existence of the PI mass gap is widely accepted, though its bounds are not yet well defined. Following the detections of BBH mergers, and in particular GW190521, made by LIGO-Virgo, it is now possible to probe the PI mass gap. Some works have investigated stellar evolution mechanisms which may increase the PI lower limit to account for GW190521 [28,[30][31][32]. Alternative approaches have examined the mass models fit to the LIGO-Virgo BBH detections [5,6,75,76] to understand whether there is evidence for a mass gap [68] or whether outliers exist [70] (see also Appendix A). In this paper, we examine three separate theoretically motivated cutoffs in the black hole mass distribution to model the PI mass lower boundary to identify if any detected binaries merged within the mass gap.
We began this study to further investigate the origin of the GW190521 signal and to understand whether an event like GW190521 could be imitated by a stellar-mass BBH merger. We conclude that if the PI mass gap begins at or below 65 M ⊙ , then GW190521 is a clear outlier to the stellar-mass BBH population and that its outlier significance cannot be explained by statistical fluctuations. Any population study, such as the analysis presented in this paper, is dependent on the sample distribution. So, we expect the significance of stellar-mass BBH population outliers presented here to fluctuate following future observing runs. However, it is critical to note that we do not expect our outlier significance of GW190521 to change in a meaningful way (See Appendix B). Our conclusion on GW190521 is consistent with the population re-weighted posterior probability, which places the primary mass m 1 firmly above 65 M ⊙ [6].
Regarding the remaining distribution of 21 BBH mergers detected by cWB, we find it unlikely that this sample distribution is consistent with the stellar-mass BBH population, assuming the PI mass limit exists at 50 M ⊙ . In this case, there is also some evidence that GW190706_222641 could be an outlier to the stellar-mass BBH population. And if the PI mass limit is at 40 M ⊙ , differences between the sample distribution and stellar-mass BBH population models cannot be reconciled. In this case, there is overwhelming evidence that cWB has detected more high mass BBH events compared to predictions made by the stellarmass models.
In future work, we aim to expand upon this analysis by testing BBH population models corresponding to specific astrophysical channels. In this case, the inclusion of spin effects into our test statistic will become important to distinguish between different alternative scenarios.
One main conclusion of this work is that GW190521 is inconsistent with the stellar-mass BBH population, assuming a PI mass limit at or below 65 M ⊙ . This result, however, is inherently dependent on our sample size and is subject to change following additional detections in future observing runs. Nevertheless, we expect the outlier significance of GW190521 presented here will not meaningfully change following future observing runs.
To test this, we consider a hypothetical scenario where the cWB pipeline identified 100 BBH events during the O1, O2, and O3a observing runs (instead of 22 BBH events). Under this scenario, we conservatively assume that GW190521 persists as the first rank event, according to M cWB . Following the statistical procedure described in Sec. II E, we estimate the outlier p-value of GW190521 to be less than 0.0001. Though this p-value is admittedly limited by our simulated sample size, our conclusion remains that no simulated BBH event representative of the stellar-mass BBH population is as significant as GW190521.
We utilize XGBoost, a boosted decision tree ML algorithm [94], as a regression tool to estimate the source frame total mass M cWB for a given event. To estimate the target variable, ML algorithms require input features which are correlated with the target variable. In our case, the target variable is the source frame total mass M inj , and the ML features are summary statistics reconstructed by cWB. Below is the list of 13 summary statistics used in the training process, ordered according to the gain importance metric. Unsurprisingly, the most important feature is a measure of the central frequency f 0 , which is expected to be inversely proportional to the detector frame total mass. Further attempts at feature pruning reduced our evaluation score on a validation data set. Some summary statistic definitions below are equivalent to those found in Ref. [105].
(i) f 0 -Energy weighted signal central frequency.