Poking Holes: Looking for Gaps in LIGO/Virgo's Black Hole Population

With the recent release of its second gravitational-wave transient catalog (GWTC-2), the LIGO/Virgo collaboration (LVC) has now detected 50 gravitational-wave (GW) events since the start of the advanced detector era, at least 46 of which came from binary black hole (BBH) systems (Acernese et al. 2014; Aasi et al. 2015; Abbott et al. 2019a, 2020b). GWTC-2 therefore provides a rich data set to infer properties of the astrophysical population of stellar mass black holes (Abbott et al. 2019b, 2020c). A robustly predicted feature that we can look for, specifically in the black hole (BH) mass distribution, is the theorized upper mass gap produced from effects due pair-instability supernovae (PISNe) that precludes formation of BHs with masses ∼55 M_⊙–120 M_⊙ from stellar collapse.

Stellar evolution simulations show that stars with core masses from ∼40–135 M_⊙ undergo PISNe in which the highly energetic gamma-rays produced in the core can collide with atomic nuclei and produce electron–positron pairs (Woosley et al. 2002). This production process absorbs energy that was previously counteracting the gravitational pressure causing the core to contract. Heavier core stars in the ∼65–135 M_⊙ range ignite oxygen leading to an unstable thermonuclear explosion that leaves behind no compact remnant (Heger & Woosley 2002; Heger et al. 2003). Lighter stars with core masses ∼40–65 M_⊙ can temporarily stabilize themselves after the ignition and thus go through a series of pulsations (PPISN) that shed large amounts of mass with each pulse (Woosley 2017). This continues until the mass of the star is too light to pair-produce, leaving the star to undergo a normal core-collapse supernova (CCSN) that leaves behind remnant black holes of masses ≲55 M_⊙ (Woosley 2019). Even more massive low-metallicity stars can bypass PISN completely and possibly form intermediate-mass black holes (IMBHs) with masses >120 M_⊙ (Spera & Mapelli 2017). While simulations consistently predict the existence of this proposed mass gap in the distribution of binary black hole masses, the precise locations of the lower and upper boundaries remain uncertain (Belczynski et al. 2016; Farmer et al. 2019; Stevenson et al. 2019; Belczynski 2020; Farmer et al. 2020; Vink et al. 2021). Analyses of GW data from the first and second Advanced LIGO–Virgo observing runs put constraints on the lower bound of this mass gap by modeling the black hole primary mass distribution as a power law with a sharp high-mass cutoff, and found support for this lower edge at ∼45 M_⊙ (Fishbach & Holz 2017; Talbot & Thrane 2018; Abbott et al. 2019b). While the GWTC-2 catalog is consistent with 97% of observed BBH primary masses lying below 45 M_⊙, population analyses using GWTC-2 and parameterized toy models find that there is less support for a sharp cutoff and instead preference for a shallow tail to larger masses (Abbott et al. 2020c). This shallow tail extending into the supposedly forbidden range of masses is primarily driven by the few GW detections that have posterior support in the mass gap, most notably GW190521 (Abbott et al. 2020d, 2020e). This leads to the possibility that there could be a separate population of black holes contributing in that mass range, formed through some other mechanism. One possibility is that remnants of previous black hole mergers undergo subsequent "hierarchical" mergers that, in dynamical environments such as globular clusters or active galactic nuclei, can contribute a significant fraction to the overall rate of mergers (Doctor et al. 2020; Farrell et al. 2020; Gayathri et al. 2020; Kimball et al. 2020a; McKernan et al. 2020; Romero-Shaw et al. 2020a; Secunda et al. 2020). For example, Kimball et al. (2020b) find evidence that GWTC-2 includes at least one merger with a second-generation component under certain assumptions about the first-generation black hole mass distribution.

The hierarchical merger scenario is not the only explanation for high-mass events like GW190521, though. Other recent work has looked more closely at the inferred source parameters of GW190521. Using gravitational waveforms with quasi-circular black hole inspirals and a standard "agnostic" prior, the LVC found the source of GW190521 to have component source-frame masses within the theorized bounds of the mass gap with ${m}_{1}^{\mathrm{src}}={85}_{-14}^{+21}\,{M}_{\odot }$ and ${m}_{2}^{\mathrm{src}}={66}_{-18}^{+17}\,{M}_{\odot }$ (Abbott et al. 2020d, 2020e). However, other waveform models and priors lead to other interpretations of this event. For example, the source of GW190521 could have had a highly eccentric orbit, been a head-on merger, or been subject to new physics allowing formation within the PISNe mass gap (Bustillo et al. 2020; Gayathri et al. 2020; Romero-Shaw et al. 2020a; Sakstein et al. 2020; Cruz-Osorio et al. 2021). Alternatively, assuming the secondary component of GW190521 comes from the same prior distribution of secondaries as other events, there is more support for the components of GW190521 to straddle the lower and upper bounds of the mass gap (Fishbach & Holz 2020). There has also been work by Nitz & Capano (2021, hereafter NC21) reanalyzing the parameter estimates of GW190521 with the recently released IMRPhenomPXHM waveform (Pratten et al. 2020), which supports mass ratios q ≡ m₁/m₂ > 4 that were not considered in the LVC analysis. NC21 found that GW190521 may be an intermediate mass ratio merger, reporting a multimodal posterior with an additional high mass ratio mode not identified in the LVC analysis. The reported source-frame component masses for the high mass ratio mode squarely puts each outside of the theorized mass gap with ${m}_{1}^{\mathrm{src}}={166}_{-35}^{+16}\,{M}_{\odot }$ and ${m}_{2}^{\mathrm{src}}={16}_{-3}^{+14}\,{M}_{\odot }$ (Nitz & Capano 2021). Figure 1 shows the 90% contours on the posterior samples from events in GWTC-2 with both the LVC GW190521 posterior samples in addition to samples from NC21 highlighted. This illustrates how differences in the analysis of GW190521 can considerably change the posterior support for its component masses to lie in the theorized PISNe mass gap. If GW190521 does "straddle" the gap, it would signal the existence of a high-mass population that could inform questions in both astrophysics and cosmology (Ezquiaga & Holz 2021).

Zoom In Zoom Out Reset image size

In this Letter we present a simple phenomenological population model parameterizing the PISNe mass gap that enforces a zero rate of BBH mergers within the gap. Our model is a complementary approach to other physically motivated models that describe the impact of PISNe on the mass spectrum (Baxter et al. 2021). Using this model we evaluate the evidence for the presence of a mass gap in LIGO/Virgo's second gravitational-wave transient catalog, and constrain its properties. We conduct each analysis twice, first using posterior samples for GW190521 released by the LVC, and alternatively using samples produced in NC21. In Section 2 we introduce our parameterized mass gap model, and the methods used to infer population properties. In Section 3 we present the results of our inference with both sets of posterior samples and two underlying mass distributions. We then discuss our interpretation of the results and astrophysical implications in Section 4 and finish with our conclusions on the support for the presence of an upper mass gap in LIGO/Virgo's BH population in Section 5.

2.1. Hierarchical Inference

We use hierarchical Bayesian inference to simultaneously infer hyperparameters of the population distribution of the primary masses (m₁), mass ratios (q), and the redshifts (z) of observed BBHs. We assume the BBH merger rate $d{ \mathcal R }$ over a given interval of masses and redshifts can be factored as

$$\begin{eqnarray}&&\displaystyle \frac{d{ \mathcal R }({m}_{1},q,z| {{ \mathcal R }}_{0},{\rm{\Lambda }})}{{{dm}}_{1}{dq}}={{ \mathcal R }}_{0}p({m}_{1}| {\rm{\Lambda }})p(q| {m}_{1},{\rm{\Lambda }})p(z| {\rm{\Lambda }}),\end{eqnarray} \tag{ 1 }$$

with Λ the population hyperparameters and ${{ \mathcal R }}_{0}$ the local (z = 0) merger rate. Under the condition that p(m₁∣Λ) and p(q∣m₁, Λ) are both normalized, and p(z) chosen such that p(z = 0) = 1, integrating the merger rate density across all primary masses and mass ratios at a given z, returns the total BBH merger rate density at that redshift, ${ \mathcal R }(z)$. The number density of BBH mergers can be related to the merger rate density b:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dN}({m}_{1},q,z| {{ \mathcal R }}_{0},{\rm{\Lambda }})}{{{dm}}_{1}{dqdz}}=\displaystyle \frac{{{dV}}_{c}}{{dz}}\left(\displaystyle \frac{{T}_{\mathrm{obs}}}{1+z}\right)\\ \quad \displaystyle \frac{d{ \mathcal R }({m}_{1},q,z| {{ \mathcal R }}_{0},{\rm{\Lambda }})}{{{dm}}_{1}{dq}},\end{array}\end{eqnarray} \tag{ 2 }$$

with V_c the comoving volume element and T_obs the total observing time with the factor of 1 + z converting source-frame time to detector-frame. Integrating the above number density across all primary masses, mass ratios, and redshifts out to a maximum z_max returns the expected number of BBH mergers in the universe out to z_max. Given a set of data {d_i } from N_obs gravitational-wave events, we can calculate the posterior on Λ following, e.g., Farr (2019) and Mandel et al. (2019):

$$\begin{eqnarray}&&\begin{array}{l}p\left({{ \mathcal R }}_{0},{\rm{\Lambda }}| \{{d}_{i}\}\right)\propto p({\rm{\Lambda }})p({{ \mathcal R }}_{0}){e}^{-{{ \mathcal R }}_{0}\langle {VT}{\rangle }_{{\rm{\Lambda }}}}\\ \prod _{i=1}^{{N}_{\mathrm{obs}}}\left[\displaystyle \int { \mathcal L }\left({d}_{i}| {m}_{1}^{i},{q}^{i},{z}^{i}\right)\displaystyle \frac{d{ \mathcal R }}{{{dm}}_{1}{dqdz}}({\rm{\Lambda }}){{dm}}_{1}{dqdz}\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal L }({d}_{i}| {m}_{1},q,z)$ is the single-event likelihood function used to infer each event's parameters, and 〈VT〉_Λ is the average sensitive time-volume when assuming a population corresponding to hyperparameters Λ. To estimate the 〈VT〉_Λ, we use the results of the LVC's injection campaign where the GWs from a fixed, broad population of sources were simulated, injected into real detector data, and searched for using the same analyses that produced GWTC-2. ¹ We use importance sampling over the detected simulated events to estimate 〈VT〉_Λ, marginalizing over the uncertainty in our estimate due to a finite number of simulated events, following Farr (2019). We assume that repeated sampling of 〈VT〉_Λ will follow a normal distribution (i.e., $\langle {VT}{\rangle }_{{\rm{\Lambda }}}\sim { \mathcal N }(\mu ({\rm{\Lambda }}),\sigma ({\rm{\Lambda }}))$), with μ the raw importance sampled estimate and σ standard error. Now we define N_eff, the effective number of independent draws contributing to the importance sampled estimate, as ${N}_{\mathrm{eff}}\equiv \tfrac{{\mu }^{2}}{{\sigma }^{2}}$, which we verify to be sufficiently high after reweighting to a population (i.e., ${N}_{\mathrm{eff}}\gt 4{N}_{\det }$). After marginalizing over the uncertainty estimating the sensitive time-volume, we write the marginalized posterior as

$$\begin{eqnarray}&&\begin{array}{l}p\left({{ \mathcal R }}_{0},{\rm{\Lambda }}| \{{d}_{i}\}\right)\propto p({\rm{\Lambda }})p({{ \mathcal R }}_{0})\\ \prod _{i=1}^{{N}_{\det }}\left[\displaystyle \int { \mathcal L }\left({d}_{i}| {m}_{1}^{i},{q}^{i},{z}^{i}\right)p({m}_{1}^{i},{q}^{i},{z}^{i}| {\rm{\Lambda }}){{dm}}_{1}{dqdz}\right]{{ \mathcal R }}_{0}^{{N}_{\mathrm{obs}}}\\ \quad \exp \displaystyle \frac{{{ \mathcal R }}_{0}\mu ({{ \mathcal R }}_{0}\mu -2{N}_{\mathrm{eff}})}{2{N}_{\mathrm{eff}}}.\end{array}\end{eqnarray} \tag{ 4 }$$

Finally, when using the commonly chosen log-Uniform prior on ${{ \mathcal R }}_{0}$ (Abbott et al. 2020c), we can marginalize over the local merger rate, neglecting terms of ${ \mathcal O }({N}_{\mathrm{eff}}^{-2})$ or greater: (Farr 2019)

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}p\left({\rm{\Lambda }}| \{{d}_{i}\}\right)\propto \sum _{i=1}^{{N}_{\mathrm{obs}}}\mathrm{log}\left[\displaystyle \frac{1}{{K}_{i}}\sum _{j=1}^{{K}_{i}}\displaystyle \frac{p({m}_{1}^{i,j},{q}^{i,j},{z}^{i,j}| {\rm{\Lambda }})}{\pi ({m}_{1}^{i,j},{q}^{i,j},{z}^{i,j})}\right]\\ \quad -\ {N}_{\mathrm{obs}}\mathrm{log}\mu +\displaystyle \frac{3{N}_{\mathrm{obs}}+{N}_{\mathrm{obs}}^{2}}{2{N}_{\mathrm{eff}}}+{ \mathcal O }({N}_{\mathrm{eff}}^{-2})\end{array}\end{eqnarray} \tag{ 5 }$$

In the last expression we further approximated the inner integral over the individual event parameters ${m}_{1}^{i},{q}^{i},{z}^{i}$ with importance sampling over K_i single-event posterior samples generated from inference with prior $\pi ({m}_{1}^{i,j},{q}^{i,j},{z}^{i,j})$. To calculate marginal likelihoods and draw samples of the hyperparameters from the hierarchical posterior distribution shown in Equation (5), we use the Bilby (Ashton et al. 2019; Romero-Shaw et al. 2020b) and GWPopulation (Talbot et al. 2019) Bayesian inference software libraries with the Dynesty dynamic nested sampling algorithm (Speagle 2020).

2.2. Parameterized Mass Gap Model

We build on models used in Abbott et al. (2020c, 2019b) and Fishbach et al. (2018) for the mass ratio and redshift of our population, specifically $p(z| \gamma )\propto {\left(1+z\right)}^{\gamma }$ and $p(q| {m}_{1},{m}_{\min },{\beta }_{q})\propto {q}^{{\beta }_{q}}{\rm{\Theta }}({{qm}}_{1}-{m}_{\min }){\rm{\Theta }}({m}_{1}-{{qm}}_{1})$ with Θ denoting the Heaviside step function ensuring m₂ is within the range [m_min, m₁]. We choose to neglect a population model for the spins, assuming that their population follows the uniform and isotropic prior used in each event's initial analysis. For the primary mass distribution we use two different models presented in Abbott et al. (2020c), the Truncated and PowerLaw+Peak models. We choose to build upon the Truncated model for its simplicity and the Powerlaw+Peak model since Abbott et al. (2020c) found it to have the highest marginal likelihood of the models used. Additionally, it is important to include the peak (Powerlaw+Peak), as it was motivated to model the pileup of events due to PPISN mass loss that is expected from the same processes predicting the upper mass gap (Talbot & Thrane 2018). We then introduce a mass gap into both the primary and mass ratio distributions by enforcing that neither component mass can lie within the gap, which we parameterize with the location of the lower edge m_g and the width of the gap w_g :

$$\begin{eqnarray*}p({m}_{1}| {m}_{g},{w}_{g},{\rm{\Lambda }})\propto \left\{\begin{array}{ll}0 & {m}_{g}\leqslant {m}_{1}\leqslant {m}_{g}+{w}_{g}\\ p({m}_{1}| {\rm{\Lambda }}) & \mathrm{otherwise}\end{array}\right.,\end{eqnarray*}$$

$$\begin{eqnarray*}p(q| {m}_{g},{w}_{g},{m}_{1},{\rm{\Lambda }})\propto \left\{\begin{array}{ll}0 & {m}_{g}\leqslant {m}_{1}q\leqslant {m}_{g}+{w}_{g}\\ p(q| {m}_{1},{\rm{\Lambda }}) & \mathrm{otherwise}\end{array}\right..\end{eqnarray*}$$

Our model prescribes a zero rate within the mass gap, which might be expected if the entire population of sources is formed through stellar collapse. We enforce an overall maximum BH mass of 200 M_⊙ so that if the upper edge of the gap is not constrained (i.e., m_g + w_g ≥ 200 M_⊙), it is equivalent to the underlying primary mass model with a maximum mass at m_g . The hyperparameters' descriptions according to each of the models used along with their corresponding priors can be found in Table 1.

Table 1. Prior Choices and Description of Hyperparameters for Used Population Models

Model	Parameter	Description	Prior
Primary Mass Model Parameters
Truncated	α	slope of the powerlaw	U(−4, 12)
	mmin	minimum mass cutoff (M⊙)	U(2 M⊙, 10 M⊙)
	mmax	maximum mass cutoff (M⊙)	200 M⊙
PowerLaw+Peak	α	slope of the powerlaw	U(−4, 12)
	mmin	minimum mass cutoff (M⊙)	U(2 M⊙, 10 M⊙)
	mmax	maximum mass cutoff (M⊙)	200 M⊙
	μp	mean of Gaussian peak (M⊙)	U(20 M⊙, 70 M⊙)
	σp	width of the Gaussian peak (M⊙)	U(0.4 M⊙, 10 M⊙)
	λp	fraction of BBH in the Gaussian component	U(0, 1)
Mass Ratio Model Parameters
PowerLaw MassRatio	βq	slope of the mass ratio powerlaw	U(−4, 12)
Redshift Evolution Model Parameters
PowerLaw Redshift	γ	slope of redshift evolution powerlaw ${\left(1+z\right)}^{\gamma }$	U(−6, 6)
Mass Gap Parameters
Agnostic MassGap	mg	lower bound of PISNe mass gap (M⊙)	U(40 M⊙, 100 M⊙)
	wg	width of the PISNe mass gap (M⊙)	U(0 M⊙, 160 M⊙)
Informed MassGap	${m}_{{\rm{g}},\min }$	lower bound of the PISNe mass gap (M⊙)	${ \mathcal N }(\mu =55\,{M}_{\odot },\sigma =10\,{M}_{\odot })$
	${m}_{{\rm{g}},\max }$	upper bound of the PISNe mass gap (M⊙)	${ \mathcal N }(\mu =120\,{M}_{\odot },\sigma =20\,{M}_{\odot })$

Download table as:  ASCIITypeset image

2.3. Injections and Sensitivity Estimates

We fit our population models to the 46 definitive BBH mergers in GWTC-2 (i.e., excluding GW170817, GW190425, GW190814, and GW190426_15215; Abbott et al. 2017, 2020a, 2020f), and since we are focused on the details of the high-mass population we neglect the low-mass smoothing feature of the models used in Abbott et al. (2020c) for simplicity. Figure 2 (left) shows the inferred merger rate density as a function of primary mass, when using the gap model on top of the underlying Truncated model, in which there is clear inference of a mass gap. When using the Truncated model, the choice of GW190521 posterior samples (either from the LVC or NC21) does not significantly affect the outcome, but the gap is inferred to be narrower and with an upper edge at lower masses when using the LVC samples. Figure 2 (right) shows the 1-d and 2-d marginal posterior distributions for the two gap parameters. Here we can see that the posterior distribution for the width of the gap is less constrained using NC21 samples, but both cases show little support for a zero-width gap. Using the LVC GW190521 samples we find support for the Truncated+Gap over Truncated model with Bayes factor $\mathrm{ln}{ \mathcal B }=2.79$, with lower and upper bounds at ${55.12}_{-4.38}^{+7.54}\,{M}_{\odot }$ and ${103.74}_{-6.32}^{+17.01}\,{M}_{\odot }$, respectively. While the gap model is clearly favored in that comparison, when using samples from NC21's for GW190521 the gap is more clearly favored with $\mathrm{ln}{ \mathcal B }=6.5$, with lower and upper bounds at ${55.33}_{-4.21}^{+5.21}\,{M}_{\odot }$ and ${126.03}_{-22.65}^{+30.25}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Figure 3 (left) shows the inferred merger rate density as a function of primary mass when imposing a gap on the most favored mass model in Abbott et al. (2020c), Powerlaw+Peak. When including the Gaussian peak in our primary mass distribution, the support for the upper mass gap significantly reduces, regardless of which GW190521 samples are used. We find log Bayes factors for inclusion of the gap to be $\mathrm{ln}{ \mathcal B }=-0.5$ and $\mathrm{ln}{ \mathcal B }=0.5$ when using the LVC and Nitz & Capano (2021) GW190521 posterior samples, respectively. Figure 3 (right) shows the 1-d and 2-d marginal posterior distributions for the gap parameters, which show poorer constraints on the gap in the Powerlaw+Peak+Gap model relative to Powerlaw+Gap. In this case, both choices of posterior samples show support for a zero-width gap, as reflected in the Bayes factors.

Zoom In Zoom Out Reset image size

Our results are inconclusive about the existence of a high-mass mass gap. While a gap is clearly inferred when using a pure power-law model of the population, adding a Gaussian peak to the mass distribution washes away the need for the gap. Furthermore, differences in parameter estimates with different priors and waveforms give rise to different inferences on the gap parameters if a gap indeed exists. These results are summarized in Table 2 through Bayes factors comparing the marginal likelihood of each model to that of the model with highest marginal likelihood (which therefore has $\mathrm{ln}{ \mathcal B }=0$). The Bayes factors for LVC and NC21 parameter estimates are treated separately in the table. We also include Bayes factors for analyses with "informed" priors on the gap boundary, where we place Gaussian priors for ${m}_{{\rm{g}},\min }$ and ${m}_{{\rm{g}},\max }$ centered on the approximate expected gap bounds (i.e., $p({m}_{{\rm{g}},\min })\sim { \mathcal N }(\mu =55\,{M}_{\odot },\sigma =10\,{M}_{\odot })$ and $p({m}_{{\rm{g}},\max })\sim { \mathcal N }(\mu =120\,{M}_{\odot },\sigma =20\,{M}_{\odot })$). With the smaller prior volume in these runs, the Bayes factors are higher than with the uninformed gap priors. Nevertheless, these Bayes factors do not increase enough to change the general finding of this work that the gap is favored with the Truncated model, but its existence is unclear when considering the Powerlaw+Peak model.

While the gap (if it exists) is difficult to resolve at present due to low number statistics, future detections will enable a finer look at the gap. Three extensions should be made when more high-mass detections are available, which we have eschewed for now due to the single event GW190521 driving the inference:

• 1.  
  Allow a different mass ratio distribution for high-total-mass events than for low-total-mass.
• 2.  
  Allow a nonzero rate in the gap, possibly with spins enforced to be near χ ∼ 0.7 to account for hierarchical mergers (Fishbach et al. 2017; Gerosa & Berti 2017; Kimball et al. 2020a; Doctor et al. 2021).
• 3.  
  Allow the merger rate normalization above the gap to be a free parameter. ²

The location of the lower edge of the PISNe mass gap has been found to be insensitive to many variations in stellar physics, especially metallicities (Farmer et al. 2019). A metallicity-independent feature in the BH mass spectrum can provide a "standard siren" that allows for independent measurements of redshift and luminosity distances to GW sources to directly measure the Hubble constant (Farr et al. 2019). The lower edge of the PISNe mass gap has been found to be very sensitive to variations in the ${}^{12}{\rm{C}}{\left(\alpha ,\gamma \right)}^{16}{\rm{O}}$ reaction rate, with some choices of rate pushing the lower bound up to ∼90 M_⊙, illustrating that constraints on the lower bound can also be used to put constraints on nuclear physics going on inside stars' cores. (Farmer et al. 2019, 2020). These astrophysical implications rely on the CO–BH mass relation from Farmer et al. (2019) that predicts a pileup of BHs below the onset of the PISNe mass gap, implying constraints on the gap lower bound that neglect a pileup could not be reliably used to constrain nuclear physics. The upper edge of the mass gap is currently not well constrained, but Ezquiaga & Holz (2021) argue LIGO/Virgo (at A+ sensitivity) will be sensitive to BBHs with component masses that could lie above the PISNe gap. Future constraints on the upper edge may also provide a novel probe of physics beyond the standard model (Croon et al. 2020).

Black holes formed through stellar collapse are expected to have a gap in their mass spectrum from ∼55 M_⊙–120 M_⊙. We assess the support, or lack thereof, for the existence of such a gap in the GWTC-2 catalog, using two parameterized black hole binary merger population models. Our population models build on the Truncated and Powerlaw+Peak models previously fit to these catalogs, and explicitly allow for a zero-rate mass gap with a population of black holes above the gap. Our analyses also consider two separate inferences of GW190521 parameters, one from the LVC and the other from Nitz & Capano (2021). We find that the results of our inference regarding the existence of a gap are contingent in part on the choice of population model and GW190521 parameter estimation results.

If a pure power law is used to describe the distribution of primary masses, we infer a mass gap from ${56.12}_{-4.38}^{+7.54}\,{M}_{\odot }$ to ${103.74}_{-6.32}^{+17.01}\,{M}_{\odot }$; however, if the data support more unequal masses for GW190521 as suggested in Nitz & Capano (2021), we infer a mass gap from ${55.33}_{-4.21}^{+5.21}\,{M}_{\odot }$ to ${126.03}_{-22.65}^{+30.25}\,{M}_{\odot }$. When using a power law with an additional Gaussian component, we no longer find significant support for a zero-rate mass gap. This does not, however, imply the nonexistence of a mass gap due to PISNe but points toward there being a secondary population of BHs that LIGO/Virgo is observing not formed through isolated stellar evolution. Future studies may be able to distinguish between these multiple formation channels in part by looking for a zero-rate gap in BH subpopulations while additionally using informed constraints on expected properties that a hierarchically formed population of BHs would have.

This research has made use of data, software, and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. This work benefited from access to the University of Oregon high performance computer, Talapas. This material is based upon work supported in part by the National Science Foundation under grant PHY-1807046 and work supported by NSF's LIGO Laboratory, which is a major facility fully funded by the National Science Foundation.

Software: Astropy (Astropy Collaboration et al. 2018), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), seaborn (Waskom 2021), bilby (Ashton et al. 2019), GWPopulation (Talbot et al. 2019).