Bridging the Gap: Categorizing Gravitational-wave Events at the Transition between Neutron Stars and Black Holes

We search for features in the mass distribution of detected compact binary coalescences which signify the transition between neutron stars (NSs) and black holes (BHs). We analyze all gravitational-wave (GW) detections by the LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration (LVK) made through the end of the first half of the third observing run, and find clear evidence for two different populations of compact objects based solely on GW data. We confidently (99.3%) find a steepening relative to a single power law describing NSs and low-mass BHs below 2.4−0.5+0.5M⊙ , which is consistent with many predictions for the maximum NS mass. We find suggestions of the purported lower mass gap between the most massive NSs and the least massive BHs, but are unable to conclusively resolve it with current data. If it exists, we find the lower mass gap’s edges to lie at 2.2−0.5+0.7M⊙ and 6.0−1.4+2.4M⊙ . We reexamine events that have been deemed “exceptional” by the LVK collaborations in the context of these features. We analyze GW190814 self-consistently in the context of the full population of compact binaries, finding support for its secondary to be either a NS or a lower mass gap object, consistent with previous claims. Our models are the first to accommodate this event, which is an outlier with respect to the binary BH population. We find that GW200105 and GW200115 probe the edges of, and may have components within, the lower mass gap. As future data improve global population models, the classification of these events will also improve.


INTRODUCTION
The LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (LVK) continue to expand the catalog of confidently detected Gravitational-Wave (GW) transients.To date, there have been 46 unambiguous detections of binary black hole (BBH) mergers (Abbott et al. 2019a(Abbott et al. , 2021b)), two detections of binary neutron star (BNS) mergers (Abbott et al. 2017(Abbott et al. , 2020a)), and, most recently, two neutron star-black hole (NSBH) merger candidates (Abbott et al. 2021c), not including the events reported in the most recent deep exafarah@uchicago.edu* NASA Hubble Fellowship Program Einstein Postdoctoral Fellow tended catalog (Abbott et al. 2021a, GWTC-2.1).However, some of the detected sources cannot obviously be ascribed to one of these source categories.For example, GW190814 (Abbott et al. 2020b) has a secondary mass of m 2 = 2.59 +0.08 −0.09 M , making it either the most massive neutron star (NS) or the lowest mass black hole (BH) detected to date.GW190719 215514, presented for the first time in GWTC-2.1 (Abbott et al. 2021a), is in a similar position.These events lie at the edges of both the NSBH and BBH populations and therefore have the potential to probe the extremes of whichever subpopulation they belong.
However, their classification remains elusive, and in the absence of detectable tides or electromagnetic counterparts that definitively identify their secondaries as NSs, other classification schemes are necessary.Some analyses have directly used NS equation of state (EOS) constraints to compute the probability that a given component mass is below the maximum allowed NS mass (Essick & Landry 2020;Abbott et al. 2020bAbbott et al. , 2021c)).In low latency, the LVK searches used a hard cutoff in component masses of 3 M in the third observing run1 (Chatterjee et al. 2020) and 2.83 M in the second observing run (Abbott et al. 2019b), both motivated by EOS predictions for the maximum possible NS mass.
However, modern models for the EOS predict different values for the maximum allowed neutron star mass, typically ranging anywhere between 2.0-2.5 M (Margalit & Metzger 2017;Ruiz et al. 2018;Shibata et al. 2019;Legred et al. 2021).See, e.g., Chatziioannou (2020) for a review of recent observational constraints.Additionally, the maximum mass of the galactic population of NSs is currently estimated to be ∼ 2−2.6 M (Antoniadis et al. 2016;Alsing et al. 2018;Farr & Chatziioannou 2020;Fonseca et al. 2021) by electromagnetic observations of pulsars.
Meanwhile, electromagnetic observations of BHs in Xray binaries suggest an absence of objects below 5M (Bailyn et al. 1998;Özel et al. 2010;Farr et al. 2011).This mismatch between the maximum observed NS mass and minimum observed BH mass points to the possibility of a "lower mass gap": a dearth of compact objects between ∼ 2.5 and ∼ 5 M .GW detections have followed suit, with a lack of BBH detections in this range (see Fig. 1).
The low number of GW detections in the purported lower mass gap may be due to lower detector sensitivity to events with low masses (which are intrinsically quieter), but studies that take such selection effects into account find that the BBH distribution cannot be trivially extended to NS mass ranges.Fishbach et al. (2020) found that a single power law cannot fit the mass distribution of both BBHs and BNSs via an analysis of events in GWTC-1 (Abbott et al. 2019a).Similarly, Abbott et al. (2021d) find a deviation from a power law in the BBH spectrum below 4 M using BBHs in GWTC-2.
We pursue a different approach, simultaneously modeling the overall distribution of all compact binaries without first dividing them into subpopulations.In doing so, we may be less susceptible to issues arising from the placement of ad hoc boundaries.By examining the total mass distribution, we can first ask whether there is a need to subdivide the observed events into separate categories at all, such as events that fall above, below, or within the lower mass gap, instead of asserting such subpopulations exist a priori.
A full description of the lower mass gap has been previously stymied by the small number of detections available.Characterizing the nature of this feature in the overall mass spectrum of compact objects can reveal the true delineation between NSs and BHs in merging binaries, providing a useful compliment to event classification based on cutoffs that are either arbitrary, motivated by external observations of arguably different populations, or reliant on theoretical predictions.With the recent influx of GW detections, enough data is available to make such classification schemes possible.
To resolve the transition between NSs and BHs, we follow the procedure of Fishbach et al. (2020) (herein FEH20) and model the full spectrum of compact binary coalesences (CBCs).We do not separate events based on previous classifications, instead using all GW events available to determine if a global fit to component masses alone is able to find two or more distinct subpopulations.Additionally, fitting a single mass model across the entire spectrum of detected CBCs allows for the inclusion of GW events that are not definitively classified into a specific source category.
We further describe our methodology and present the population models used in Sec. 2. We show that the transition between NSs and BHs cannot be described by a single power law and explore other morphologies in Sec. 3. In Sec.4.2, we demonstrate that the models considered in this work are the first to be able to accommodate GW190814.We describe insights on the NSBH events in Sec.4.3.We then classify all low-mass events with respect to the inferred subpopulation delineations in Sec.4.1, concluding in Sec. 5

METHODS
We first describe our selection criteria for which events to include within our analysis in Sec.2.1.We then describe our parametric models for the mass distribution in Sec.2.2 before describing metrics that can be used to describe features in the mass spectrum, sometimes without the need for a particular parametrization, in Sec.2.3.

Event List
We consider all published CBC detections made by the LVK to date that have a matched-filter signal-tonoise ratio (SNR) greater than 11.That is, all events in GWTC-2 (Abbott et al. 2021b) that pass this threshold, as well as the two NSBH detections from the second half of the third observing run (Abbott et al. 2021c, O3b)).None of the events newly presented in GWTC-2.1 pass this threshold (Abbott et al. 2021a).This results in 46 events, all of which are shown in Fig. 1.For the purposes of estimating sensitivity, we treat the NSBH events as if they were detected in the first half of the third observing run because data from other detections during O3b are not yet publicly available.This choice introduces a modest bias in the inferred rate of NSBHlike events, which is discussed in Appendix B

Population Models
We consider several nested mass models.That is, we describe our conclusions in the context of different sets of assumptions about the mass distribution, each of which corresponds to a specific hyperprior as described in Tab. 1.Following Fishbach &Holz (2020), FEH20, andDoctor et al. (2020) we parametrize the joint component mass distribution with separate draws for each component mass from a one-dimensional mass distribution and a pairing function.Our one-dimension mass model can be described by a two-piece broken power law adorned with several multiplicative features that chip off parts of the mass distribution to produce dips and additional features.That is, the one-dimensional distribution is described by where is a low-pass Butterworth filter with roll-off mass m 0 and the sharpness of the roll-off set by η (larger η imply sharper roll-offs), is a corresponding high-pass filter, and is a notch filter that subtracts a fraction of the signal (set by A) between γ low and γ high .Our one-dimensional mass model builds upon a basic broken power law by adding a high-pass filter at the lowest masses to allow for a fixed but smooth turn-on of the mass distribution, a low-pass filter at the highest masses to allow for a variable tapering of the mass distribution, and a notch filter in between to model a potential lower mass gap.This model has 12 parameters, which we collectively refer to as λ 1D .It was first presented in FEH20.However, we include an additional high-pass filter at low masses instead of a sharp cut-off at m min .FEH20's model is therefore nested in ours in the limit η min 1.To construct population models for both component masses, we employ a pairing function f p that relates primary and secondary masses: where Θ(•) is the Heaviside function that enforces the convention m 2 ≤ m 1 .
Note that the one-dimensional distribution's hyperparameters λ 1D are shared between the primary and secondary mass distributions.Note also that the onedimensional distribution does not, in general, correspond to the marginal distribution for either the primary or secondary mass.We use the pairing function (6) to allow for the possibility that most BBHs have different formation channels than NS-containing or lowermass-gap binaries.We note that this is a different pairing function than that used in FEH20, which employed a single power law in mass ratio for all sources.Again, their model is nested within ours in the limit For simplicity, we assume fixed, independent spin and redshift distributions.We assume sources are uniformly distributed in co-moving volume (V c ) and source frame time so that We also assume spin distributions where each component is independent and follows a distribution that is uniform in magnitude χ 1,2 and isotropic in orientation so that p(|χ|) = U(0, 1), ( 8) where U(x, y) is the uniform distribution between x and y and z 1,2 = cos θ 1,2 is the cosine of the tilt angle between component spin and a binary's angular momentum.
Although some authors have shown evidence for the evolution of the mass distribution and/or the merger rate with redshift (Fishbach et al. 2021;Abbott et al. 2021d), we do not expect these choices to affect the mass inference within current statistical uncertainties.
Eqs. 1 and 6 intentionally provide a lot of flexibility in the functional form of the mass distribution via 14 free parameters while encoding several features predicted by theory.Because the current data may not be able to constrain all the hyperparameters, we consider three nested models in order to demonstrate several of our conclusions more clearly.In order of most to least model freedom, Illustration of the primary mass distribution for the three mass models used in this work.All models are depicted with the same values for m break , α1, and α2.
• Broken Power Law + Dip (BPL+Dip): assumes the edges of low-mass features are steep but leaves the locations as free parameters, and fixes the break in the power law to occur at the lower edge of the gap.Sets η min = η high = 50 1 and assumes m break = γ low .Alternatively, at times we assume the gap is deep with a sharp lower edge but allow the upper edge to be fit as a free parameter (A = 0.98 and fits η high ).
• Broken Power Law + Drop (BPL+Drop): makes similar assumptions as BPL+Dip, except for η high = 0.This de facto models a sharp dropoff above γ low and removes the feature at γ high .
• Broken Power Law (BPL): assumes a steep turn-on at m min but removes the dip/discontinuity by setting η low = η high = 0.
An illustration of these mass models is provided in Fig. 2, and Tab. 1 summarizes the hyperprior adopted for each of these models.

Characterizing the Prominence of Features in the Mass Distribution
In the case of a mass gap with sharp edges (η low,high 1), the notch filter strength A is an unambiguous measure of the gap depth, and the parameters γ low and γ high are the precise locations of the beginning and end of the gap, respectively.With this in mind, we can quantify the evidence in favor of a gap with a Bayes factor between models where A = 0 and A = 0: B A =0 A=0 .However, this assumes the shape of the gap is well known a priori.This is not the case, particularly for the upper edge of the gap.In the case of smooth gap edges, the depth and location of the gap are dictated by a more complex combination of hyperparameters.We Table 1.Hyperparameters of our mass model and hyperpriors corresponding to specific nested models: Broken Power Law + Dip (BPL+Dip), Broken Power Law + Drop (BPL+Drop), and Broken Power Law (BPL).We denote the uniform distribution between x and y as U(x, y), list specific values that are fixed in some priors, and denote when hyperparameters are irrelevant to a specific nested model with "N/A".Note that there are two priors that we refer to as BPL+Dip since they explore similar phenomenology: one assumes a sharp upper edge of the gap but variable gap depth and the other assumes a deep gap but allows for smooth edges.We specify which priors were used in context within the text.† The allowed values of A within Broken Power Law + Dip are between 0 and 1, whereas the elimination of the upper edge of the notch filter in Broken Power Law + Drop allows the values of A to range between 0 and 2 in that case.
therefore define a few straightforward measures of gap prominence which can be used regardless of the exact parametrization adopted.
The first is the existence of a local minimum in the mass distribution between 1-10M .A lack of support for such a dip in the mass distribution in this region would indicate either a lack of a lower mass gap or an inability to resolve one if it does exist.Throughout this text, we report Bayes factors for the existence of a local minimum that are defined as the ratio of the percentage of draws from the hyperposterior that find a local minimum to the percentage of draws from the hyperprior that find a local minimum.
The second is the relative probability of events between γ low and γ high vs. that between m min and γ low : With this definition, r NS gap 1 indicates a prominent drop-off in the mass distribution at γ low , whereas r ∼ 1 is indicative of a smoother transition above NS masses.
The third is similar to the second: it is the relative probability of events between γ low and γ high vs. that in an equally wide mass interval starting at γ high : A value of r BH gap 1 indicates a very prominent lower mass gap, whereas r consistent with unity indicates a flat transition between NS and BH masses.

Statistical Framework
Using the parametrized distributions from Sec. 2.2, we construct a hierarchical Bayesian inference to determine the appropriate population-level parameters, Λ = {λ 1D , β low , β high } given the observed set of data {D j } for N observed events (see, e.g., Loredo 2004;Thrane & Talbot 2019;Mandel et al. 2019, for more details).We model the data as an inhomogeneous Poisson process with the rate density (expected number of events per unit time per single-event-parameter hypervolume) given by dN dm 1 dm 2 , ds 1 ds 2 dz = Rp(z)p(s 1 )p(s 2 )p(m 1 , m 2 |Λ) (12) where R acts as a normalizing constant that sets the overall magnitude of the rate.The posterior on the population hyper-parameters, assuming a prior on the overall rate of mergers p(R) ∼ 1/R and marginalizing, is where is the marginal likelihood for the j th event, is the fraction of detectable events in a population described by Λ, and P (det|m 1 , m 2 , s 1 , s 2 , z) is the probability that any individual event with parameters m 1 , m 2 , s 1 , s 2 , and z would be detected, averaged over the duration of the experiment.
If we wish to simultaneously infer the properties of individual events along with the population hyperparameters, we simply do not marginalize over those events' parameters and obtain, e.g., p(m where p(D i |m 2 ) is still marginalized over s 2 , and z (i) .We typically investigate the different levels of our hierarchical inference separately, though.That is, we examine the hyperposterior for Λ in the populationlevel inference after marginalizing over uncertainty in individual-event parameters, or we examine the eventlevel population-informed posteriors for individual event parameters after marginalizing over the uncertainty in the inferred population.
In practice, the high-dimensional integrals in Eqs. 14 and 15 are approximated via importance sampling.That is, given a set of N j event-level posterior sam-ples for the j th event drawn with a reference prior (17) Similarly, by simulating a large set of M signals drawn from an injected population p inj , we can approximate Eq. 15 with a sum over the subset of m detected signals We approximate the detectable set of signals as those that yield an optimal network signal-to-noise ratio (SNR) ≥ 11 based on networks of the LIGO Hanford and Livingston detectors during O1 and O2 as well as the LIGOs and Virgo during O3a.During each observing run, we use a single representative power spectral density (PSD) to model each detector's sensitivity separately (LVK 2015a(LVK ,b, 2017a(LVK ,b, 2020) ) and estimate the network SNR by coherently projecting simulated signals into all detectors.While this semi-analytic approach to sensitivity estimation is commonplace in the literature, we note that it makes several assumptions that may not be perfectly accurate (a single stationary PSD for each detector, the same SNR threshold for all systems, etc.).
The size of the possible systematic errors these assumptions could introduce is not known precisely, but we believe it is less than O(10%) based on comparisons of our semianalytic rate estimates with fixed populations to those previously published in the literature using real searches (Abbott et al. 2019a(Abbott et al. , 2021b)).
We sample from the posterior distribution in Eq. 13 using the approximations in Eqs. 17 and 18 to determine the shape of the mass distribution using gwpopulation (Talbot et al. 2019).Furthermore, where needed, we estimate Bayes factors via Savage-Dickey Density Ratios (Dickey & Lientz 1970;Wagenmakers et al. 2010) using the hyperposteriors and the hyperpriors described in Tab. 1.

POPULATION-LEVEL INSIGHTS
In this section, we search for and characterize features between NS-like and BH-like masses.Following FEH20, we first determine the existence of such a feature by asking whether a single power law can describe the CBC mass spectrum between 1 − 10M in Sec.3.1.We find that it cannot.We then describe the nature of this deviation from a power law in Sec.3.3 through the features provided by the parametrized models from Sec. 2.2.

Deviation from a Single Power Law Describing BH and NS Masses
Neither FEH20 nor Abbott et al. (2021d) include GW190814, nor do they include GW190425, a BNS detected in the first half of the third LVK observing run (Abbott et al. 2020a).NSBH-like events GW200105 162426 and GW200115 042309 (abbreviated GW200105 and GW200115), which contain secondaries near the upper limits of what can be considered a NS and primaries at lower BH masses than have been previously detected in definitive BBHs, had not been published at the time of either study.These events appear to bridge the NS and BH subpopulations, so the deviations from a single power law in the range 1 − 10 M may no longer be present in the most recent data.We therefore aim to determine if the inclusion of these events removes the feature found by FEH20 and Abbott et al. (2021d).
The first mass model we fit is Broken Power Law.It has one power law spectral index α 1 at NS-like masses and another α 2 at BH-like masses, with a transition between the two at m break .Draws from the hyperposterior for Broken Power Law are shown in the top panel of Fig. 3.There is a clear transition in the mass distribution at m break = 2.4 +0.5 −0.5 M .This value is consistent with current constraints on the maximum mass of cold, non-spinning NS, which are typically between 2.0-2.5 M at 90% confidence (see, e.g., Essick et al. 2020;Legred et al. 2021;Farr & Chatziioannou 2020;Chatziioannou 2020), although larger values are still consistent with low-density nuclear physics.
We find α 2 − α 1 > 0 at 99.3% credibility, indicating that a single power law is unable to fit the mass distribution down to low masses.This is consistent with findings in FEH20 and Abbott et al. (2021d), and therefore the low mass events GW190814, GW200105, and GW200115 do not fully bridge the BH and NS populations.
As an aside, we also note that β 2 − β 1 > 0 at 99.6% credibility, indicating that the pairing functions between BBHs and CBCs with low secondary mass are distinct.However, because of the fixed location of the break between β 1 and β 2 , we do not use these parameters to locate delineations between subpopulations.

Steep Decline After NS Masses
The steepness of the NS part of the power law in Broken Power Law (i.e.large magnitude of α 1 ) may be driven primarily by a difference in the merger rates of BNS and BBH systems, rather than by the shape of the NS distribution itself.While our particular analysis has a slight bias in the relative heights of the low mass and high mass parts of the mass distribution because of the choice to treat the NSBH events as detected in the second half of the third observing run, many other studies have found the BNS merger rate to be several times higher than the BBH merger rate (Abbott et al. 2019a,c;Fishbach et al. 2020;Abbott et al. 2021b,d).
To decouple the shape of the NS distribution from the difference in merger rates between BNS and BBH systems, we we turn to Broken Power Law + Drop.This model allows for a drop-off in the mass distribution at m break .The extent of the discontinuity is parameterized by A, which is a free parameter.A = 2 corresponds to a maximal discontinuity: in this case, the rate of BHcontaining events drops to zero after m break .A = 0 corresponds to no discontinuity -in this case, Broken Power Law + Drop reduces to Broken Power Law.The resulting fit to this model is shown in the middle panel of Fig. 3.A clear drop-off is found after NS masses, which can be quantified by r NS gap , the ratio of probability mass below m break to that in a comparable interval just above it (Eq.10).We find r NS gap = 4.1 +10 −3.0 , which excludes 1 to > 90% credibility.
Fig. 4 shows the inferred posteriors for A and α 1 under the Broken Power Law + Drop framework.These confirm the need for a sharp drop-off in the mass distribution at m break , either by a large negative α 1 or by a large A. However, at the largest values of A, α 1 is unconstrained.We therefore conclude that the steep  drop-off in the merger rate after NS masses is driving the constraints on α 1 , rather than the shape of the NS distribution itself.This is consistent with the findings of Landry & Read (2021), who show that the mass distribution of NSs in binaries is not yet able to be determined from GW data alone.
Within Broken Power Law + Drop, the location of the break in the power law (and therefore the location of the discontinuity) is inferred to be at m break = 2.20 +0.51 −0.41 M , which is consistent with but shifted lower than the corresponding value inferred by Broken Power Law (2.4 +0.5 −0.5 ).This is due to the different modelling assumptions inherent in both models, and both values are equally plausible.Furthermore, given the one-dimensional posterior on A, we find that the data cannot distinguish between A = 0 and A = 0 and therefore have no clear preference for Broken Power Law or Broken Power Law + Drop.

Searching for a Lower Mass Gap
Having confidently established a deviation from a single power law in the form of a steep decrease in the mass distribution after NS-like masses, we wish to identify other features.In particular, we search for signatures of the purported lower mass gap.
To do this, we employ the Broken Power Law + Dip model, described in Sec.2.2.The fit to this model is shown in the bottom panel of Fig. 3 assuming η high = 50.We infer the depth of a potential mass gap through the hyperparameter A, where A = 0 reduces to the Broken Power Law mass model and A = 1 corresponds to a maximal gap.The lower and upper edges of this gap are represented by γ low and γ high , respectively.We infer γ low = 2.2 +0.7 −0.5 M and γ high = 6.0 +2.4 −1.4 M .These are consistent with the values found in FEH20 (γ low = 2.2 +0.6 −0.5 M , γ high = 6.7 +1.0 −1.5 M ).This agreement with a GWTC-1 analysis suggests that the addition of events in GWTC-2 and the NSBH-like events are consistent with this previously inferred population and do not fully fill the lower mass gap, if one exists.These events may, however, indicate that a lower mass gap might not be completely empty.
Within the framework of Broken Power Law + Dip, we find suggestions of a lower mass gap, but are unable to conclusively detect or rule out its existence.The posterior on A is shown in Fig. 5.It peaks at large values but does not rule out the possibility that some events occupy the gap.A model with a maximal gap (A = 1) is favored over no gap (A = 0) with a Bayes factor of 55.0, and r BH gap = 1.84 +5.74 −1.14 , which is suggestive of the existence of a lower mass gap (r BH gap is defined in Eq. 11).
However, we find a Bayes factor of only 1.1 in favor of a local minimum between 1 − 10 M , indicating that the assumption of sharp gap edges may be driving the other metrics of gap prominence2 .
Due to the symmetric nature of the notch filter n(m|A, γ low , γ high , η low , η high ) used in Broken Power Law + Dip, a sharp drop-off in the mass distribution after γ low will always be accompanied by a subsequent rise at γ high when we set η low = η high} 1.Thus, if the data prefer such a sharp drop-off that is not achieved through a steep α 1 , using Broken Power Law + Dip may artificially inflate our metrics of gap significance.This systematic is due in part to the choice to set η low = η high = 50, rather than leaving the steepness of the gap edges as free parameters.We therefore seek to determine if the data actually prefer a subsequent rise in the mass distribution after the initial drop at γ low by allowing η high to be a free parameter.Positive values of η high indicate the existence of a subsequent rise, and larger positive values imply steeper rises.To lessen the number of free parameters being fit to a region with relatively few detections, we enforce a drop in the mass distribution after NS masses by fixing η low = 50 and A = 0.98.We are therefore asking, given the existence of a drop in the mass distribution at γ low , do the data necessitate a subsequent rise?Fig. 6 shows the two dimensional posteriors on γ high and η high , the parameters determining the location and steepness of the upper edge of the mass gap.We find η high > 0 at 99.8% credibility, with lower 90% credible bound at η high = 2.82, indicating that if there is a large drop in the mass distribution at γ low , the data prefer a subsequent rise at γ high rather than an additional drop (η high < 0) or a featureless transition (η high = 0).We additionally find a Bayes factor of 1.14 in favor of a local minimum between 1 − 10 M .
In this model fit, γ high = 5.98 +0.99 −0.87 M , which is completely consistent with the location of the global maximum of the BBH distribution found in Abbott et al. (2021d).They find 7.8 +1.8  −2.0 M for the Power Law + Peak model and 6.02 +0.78  −1.96 M for the Broken Power Law model, the latter of which is more morphologically similar to the high-mass behavior of the models considered in this work.
As a final test of our results, we perform a fit where all hyperparameters are allowed to vary.The results are described in Appendix A, and we find a Bayes factor of 1.63 in favor of a local minimum.This indicates slightly more evidence for a lower mass gap, but is far from conclusive.However, due to the small number of detections at low masses, we cannot meaningfully constrain more than two hyperparameters in this region.We therefore fix η high = 50 and allow A to vary for the remaining analyses, noting that a steep upper edge and variable gap depth allows for more intuitive interpretations of the event classification presented in Sec.4.1.

EVENT-LEVEL INSIGHTS
Building upon our population-level insights, we now place individual events in the context of the inferred features in the mass distribution.Sec.4.1 discusses how we do this in general, marginalizing over both the uncertainty in individual events' parameters as well as the population hyperparameters.We then discuss the implications for specific events that previous models have not been able to accommodate: GW190814 (Sec.4.2) and  the two NSBH-like events GW200105 and GW200115 (Sec.4.3).

Event Classification
Fits to the full spectrum of CBCs allow us to classify events with respect to features inferred from the detected population.
Tab. 2 uses the Broken Power Law model fit to characterize events in relation to m break .It lists the probabilities that each event's component masses lie in regions of interest.For example, if one interprets m break to be the maximum NS mass, the probability that an event contains a NS is P (m 2 ≤ m break ) and the probability that an event is a BNS is P (m 1 ≤ m break ), since m 2 ≤ m 1 by definition.Similarly, if we consider objects with masses larger than m break to be BHs, the probability that an event is a NSBH is P (m 2 ≤ m break < m 1 ).Using these assumptions, we find that all probabilities presented in Tab. 2 are qualitatively consistent with the classifications given to each event by the LVK.
GW170817 is found here to have P (m 1 ≤ m break ) = 0.99, which is consistent with the LVK's classification of this event as a BNS.GW190425 has P (m 1 ≤ m break ) = 0.9, with < 10% support for a NSBH classification.This is consistent with this event's classification as a BNS.NSBH events GW200105 and GW200115 both have larger P (m 2 ≤ m break < m 1 ) than their respective P (m 1 < m break ) or P (m 2 > m break ), so under this framework they are more likely NSBHs than BNSs or BBHs.The secondary mass of GW190814 is neither definitively a NS nor a BH under this classification scheme: the posterior on m break overlaps almost completely with that of GW190814's m 2 .
Tab. 3 uses the Broken Power Law + Dip model fit to characterize events in relation to γ low and γ high .Events that have been previously identified by the LVK as containing a NS (GW170817, GW190425, GW200105, GW200115) all have P (m 2 ≤ γ low ) > 0.95.BNS event GW170817 has P (m 1 ≤ γ low ) > 0.95 as well, but BNS event GW190425 has ∼ 10% support of its primary being in the mass gap rather than a NS.This is likely due to the relatively large uncertainties on both the component masses of GW and on γ low .Events previously identified by the LVK as being consistent with NSBHs (Abbott et al. 2021c, GW200105, GW200115) receive a similar classification in this analysis, with large P (m 2 ≤ γ low ), and near-vanishing P (m 1 ≤ γ low ).
Additionally, the component masses of three events have considerable support in the mass gap: the secondary of GW190814, the secondary of GW190924 021846, and the primary of GW200115.While a totally empty gap is not ruled out -none of the three have more than 75% support in the gap -the combined effect of these three events is to decrease the inferred depth of the lower mass gap.The secondary mass of GW190814 remains ambiguously classified: it has non-negligible support of being either in the lower mass gap or being a NS.As the posterior on m 2 for GW190814 is well constrained, the ambiguous classification is instead due to the inference's inability to determine whether the event should be positioned in the first power law or in the notch filer of Broken Power Law + Dip.This results in bimodality in γ low 's posterior distribution, which we discuss further in Sec.4.2.As more detections of GW events constrain the global mass distribution, uncertainties in the locations of γ low and γ high will decrease, allowing for more a more definitive classification of this event.
The probabilities in Tabs. 2 and 3 are marginalized over the uncertainties in each event's component masses and the uncertainties in the relevant model hyperparameters.
We further examine three events that were published as "exceptional" by the LVK.

GW190814
GW190814 was observed in the first half of the third observing run (Abbott et al. 2020b) and has the most asymmetric mass ratio of any gravitational wave event to date.Its secondary has a mass of m 2 = 2.59 +0.08 −0.09 M , leaving the nature of this component ambiguous: depending on the choice of maximum allowed NS mass, it can either be classified as a low-mass BH or the most massive NS observed to date.Previous "leave-one-out" studies (Abbott et al. 2021d;Essick et al. 2021) have definitively identified GW190814 as an outlier with respect to the detected BBH population.Additionally, Landry & Read (2021) find that including the secondary of GW190814 in the population of neutron stars in merging binaries considerably changes the inferred mass distribution when compared to a fit that excludes this event, though uncertainties in both fits are large.It is therefore desirable to find a model that can account for such an extreme event.
We perform a similar study with GW190814 by comparing the mass distributions inferred with the full event list to those inferred by excluding GW190814.We find that in the case of both Broken Power Law and Broken Power Law + Dip, the inferred distributions without GW190814 are consistent with those with GW190814.In the framework of Broken Power Law, the posterior of m break is largely unchanged with the exclusion of this event, but the posterior distribution of β 1 shifts to higher values when excluding GW190814.In the framework of Broken Power Law + Dip, the posterior of γ low when GW190814 is included in the inference is in some disagreement with that inferred without GW190814.When GW190814 is included, the posterior on γ low is bimodal, with the lower mode corresponding to putting GW190814's secondary in the gap and the higher mode corresponding to locating the gap at higher masses than GW190814's secondary.Excluding GW190814 keeps only the lower mode.Additionally, the posterior on A peaks at a slightly higher value when GW190814 is excluded.However, the inferred values with and without GW190814 are still consistent with one another, as can be seen in Fig. 7. Therefore, both Broken Power Law and Broken Power Law + Dip successfully account for this event.
Note that such leave-one-out analyses can introduce biases that may lead to the false classification of the left-out event as an outlier unless modifications to the likelihood are used to mitigate this effect (Essick et al. 2021).However, since these modifications tend to bring the hyperposterior inferred without the event in question closer to that inferred with all events, we assume the conclusions derived from Fig. 7 will remain the same with the use of the methods presented in Essick et al. (2021).

GW200105 and GW200115
GW200105 and GW200115 are the most definitive detections of NSBH-like systems via gravitational waves to date, observed by the LVK during the second half of the third observing run (Abbott et al. 2021c).GW200105 has component masses m 1 = 8.9 +1.2 −1.5 M and m 2 = 1.9 +0.3 −0.2 M , and GW200115 has component masses m 1 = 5.7 +1.8 −2.1 M and m 2 = 1.5 +0.7 −0.3 M .Both events therefore have components near the edges of the purported lower mass gap, and analyzing them in the context of the Broken Power Law + Dip model is helpful in understanding their classification.Likewise, these events play a large role in probing the values of γ low and γ high inferred in this analysis.
We first determine if GW200105 and GW200115 are consistent with the rest of the population of CBCs by performing a model fit to Broken Power Law + Dip without these events and comparing it to the population inferred with these events.Appendix B shows these two fits to the data.We find that the two inferred populations are consistent with one another, and that the NSBH events have the effect of constraining γ low and γ high towards a narrower gap.This effect, as well as the fact that both events have P (γ low ≤ m i ≤ γ high ) < 0.5, indicate that GW200105 and GW200115 are both consistent with being "gap-straddling binaries," though to the inference when GW190814 is included, and grey contours correspond to the inference when it is excluded.Almost all hyperparameters remain unchanged with the inclusion vs exclusion of GW190814, indicating that this event is consistent with the inferred population using this model.The only parameter that exhibits a noticeable change is γ low , where the bimodality of the posterior distribution disappears with the exclusion of GW190814.there is still a definite possibility that either or both events have one component in the lower mass gap.

DISCUSSION
A lower mass gap between NSs and BHs has been suggested by observations of X-ray binaries that find a lack of black holes below 5 M (Bailyn et al. 1998;Özel et al. 2010;Farr et al. 2011) despite predictions for the maximum allowed neutron star mass by the neutron star equation of state (Margalit & Metzger 2017;Ruiz et al. 2018;Shibata et al. 2019;Chatziioannou 2020;Legred et al. 2021) of ∼ 2 -2.5 M .Such a dearth of systems just above the maximum neutron star mass is in contrast with theoretical studies of compact object formation scenarios which predict a continuous distribution of supernova remnant masses that have a smooth transition from NSs to BHs (e.g., Fryer & Kalogera 2001;Fryer et al. 2012;Zevin et al. 2020;Drozda et al. 2020).
In this work, we explore whether there is evidence for a lower mass gap in the existing gravitational wave data.Inferences on the edges of such a gap have previously been carried out by examining the extrema of separate black hole and NS distributions.However, this approach lessens the amount of information available to both analyses, and can lead to biased inferences.It is insufficient to simply look at the largest neutron star mass and smallest black hole mass in the catalog, and make statements about a mass gap.The best estimates for individual compact object masses can only be achieved through an analysis of the full population, including systems on both sides of the gap as well as those in the gap (if any).Fitting a single model across the entire spectrum of compact binary coalescences allows us to find features that are difficult to infer from fits to the BBH, BNS, or NSBH mass distributions alone.
We present global fits to the mass distribution of NSs and BHs in merging binary systems in Fig. 3.All models considered find a sharp drop in the mass distribution at ∼ 2-3 M , either via a steep power law (top panel) or a sharper "discontinuity" (middle and bottom panels).Although this feature has been determined solely from the mass distribution of GW systems, it presumably relates to the maximum mass of neutron stars in GW binaries.It is suggestive that this break occurs precisely in the mass range predicted from equation-of-state calculations Legred et al. (2021).More detections at NS masses will inform the shape of the mass spectrum below 3 M , allowing us to better resolve the nature of this drop-off.
For models that allow for a feature in the mass distribution after the initial "neutron star" drop, we find a subsequent rise at ∼ 4.5-8.5 M (bottom panel of Fig. 3).While this hints at a potential mass gap between neutron stars and black holes, we cannot confidently detect nor rule out the existence of such a feature with the limited detections and low sensitivity in that region.
As we emphasize above, understanding a single GW event can only be done in the context of the full underlying population.Given our population model which bridges the gap between NSs and BHs, we are able to analyze individual events in the context of this model, and thereby classify them as BNS, NSBH, BBH, or binaries with one component in the mass gap.
Using this population-based classification scheme, we find that GW200105 likely straddles the mass gap, with its primary component confidently higher than the inferred upper edge, and the secondary likely falling below the edge, though uncertainties on the lower edge location make the exact classification uncertain.We similarly reaffirm that GW200115 is an NSBH, with its secondary component below the lower edge of the mass gap and its primary having roughly equal support below and above the upper edge.Our inference is unable to determine if GW190814's secondary mass lies below or above the lower edge of the mass gap, but slightly prefers it to lie within the mass gap.This ambiguity is driven by uncertainty in the location of the lower edge of the gap.Thus, the nature of GW190814 may emerge as future data constrain global population models.We re-emphasize that these results are based solely on GW data; the existence of baryons has not been assumed in this analysis.Nonetheless, the conclusions are broadly consistent with expectations from NS EOS studies as well as analyses of galactic X-ray binaries.
All of the models considered in this work can accommodate the full list of observed events without any obvious outliers.In particular, we find that GW190814 is consistent with the full set of compact binary coalescences, even if we are not yet able to identify to which subpopulation it belongs.In this work we focused on the component mass distribution to identify the component objects of binaries.However, future investigations of the mass ratio distribution may help us distinguish between NSBH and BBH population with NSBH systems in recent population synthesis work tending to prefer mass ratios < 0.5 (Broekgaarden et al. 2021).
The delineation between neutron stars and black holes is a fundamental observable of GW binary populations.Observable constraints will inform the astrophysics of compact object formation, neutron star equation-ofstate studies, and the formation channels of GW binaries (e.g.Fryer et al. 2002;Belczynski et al. 2012;Fryer et al. 2012;Müller et al. 2016;Breivik et al. 2019;Shao & Li 2021;Liu et al. 2021).For example, the confident identification of a feature in the overall mass distribution that is distinct from the maximum neutron star mass inferred from nuclear theory would have important implications for specific supernova explosion mechanisms (Fryer et al. 2012;Kiziltan et al. 2013).From GW observations alone, we have identified a sharp drop in the mass distribution of compact objects at ∼ 2-3 M , as might be expected from the maximum allowed mass of neutron stars.There is insufficient data to confidently constrain the existence of a subsequent mass gap between neutron stars and black holes.However, with the second half of LVK's third observing run completed, and a fourth observing run planned, the GW community is positioned to definitively resolve or disprove this mass gap. the events GW200105 and GW200115.The mass distribution does not change significantly with the exclusion of these events, indicating that the NSBHs are consistent with the rest of the detected population.When the NSBHs are included, the upper edge of the mass gap, γ high , is slightly better constrained and support for a wider gap (low γ low and high γ high ) diminishes.This suggests that these events probe the gap location.Hyperparameters on which the exclusion of the NSBH events has no effect are not included in this plot.

Figure 1 .
Figure 1.90% posterior credible intervals for the component masses for all CBCs included in this population study assuming uniform priors in detector-frame masses.Events classified by the LVK as BBHs, BNSs and NSBHs are shown in black, green, and blue, respectively.The ambiguously classified event GW190814 is shown in orange.The grey band indicates the approximate location of the purported "lower mass gap."GW190814 is the only event within this region at more than 90% credibility.
Figure 2.Illustration of the primary mass distribution for the three mass models used in this work.All models are depicted with the same values for m break , α1, and α2.
mmin roll-off mass for high-pass filter at the lowest masses allowed U(1.0 M , 1.4 M ) filter ηmin sharpness of high-pass filter at the lowest masses allowed ηmin =

Figure 3 .
Figure 3. Draws from the inferred mass distributions between 1−10M using the Broken Power Law (top), Broken Power Law + Drop (middle), and Broken Power Law + Dip (bottom) models.Each panel has 200 random draws from the hyperposterior (thin black lines).As an illustration, the draw that maximizes the posterior is highlighted in color.90% credible intervals on the locations of features in the mass distribution are indicated by colored bands.All three models find a steep falloff in the mass distribution between ∼ 2.2 − 3M and Broken Power Law + Dip finds a subsequent rise at ∼ 6M .

Figure 4 .
Figure 4. Values of the spectral index governing the low mass part of the spectrum α1 and drop depth A inferred by the Broken Power Law + Drop model.Contours indicate 1-σ, 2-σ, and 3-σ credible regions.The two parameters are correlated and the lower right region is excluded: mild power law slopes are only allowed when A is large.Thus, the inference finds a dramatic drop-off in the merger rate, either through a steep negative initial power-law slope (α1), or a discontinuity (large A).The two scenarios are both allowed by the data due to the inability to resolve the shape of the NS distribution.

Figure 5 .
Figure 5. Posterior on the gap depth parameter A under the Broken Power Law + Dip model including all events.A = 1 indicates a maximal gap, and A = 0 indicates no mass gap.A peaks at 0.90, with 90% highest density interval between 0.55 and 1.0 (vertical dashed lines), indicating support for both a partially-filled and totally empty mass gap.

Figure 6 .
Figure 6.Corner plot of the parameters governing the location (γ high ) and steepness (η high ) of the upper edge of the mass gap in the case where a steep decline in the mass distribution is enforced after NS masses (A = 0.98).The 90% highest probability density interval for η high is indicated by vertical dashed lines.η high is found to be positive, indicating a rise in the mass distribution near γ high after the drop enforced at γ low .

Figure 7 .
Figure7.Effects of leaving out event GW190814 under the Broken Power Law + Dip model.Green contours correspond to the inference when GW190814 is included, and grey contours correspond to the inference when it is excluded.Almost all hyperparameters remain unchanged with the inclusion vs exclusion of GW190814, indicating that this event is consistent with the inferred population using this model.The only parameter that exhibits a noticeable change is γ low , where the bimodality of the posterior distribution disappears with the exclusion of GW190814.

Figure 8 .
Figure8.A subset of hyperposteriors for the Broken Power Law + Dip model inferred with (green) and without (brown) the events GW200105 and GW200115.The mass distribution does not change significantly with the exclusion of these events, indicating that the NSBHs are consistent with the rest of the detected population.When the NSBHs are included, the upper edge of the mass gap, γ high , is slightly better constrained and support for a wider gap (low γ low and high γ high ) diminishes.This suggests that these events probe the gap location.Hyperparameters on which the exclusion of the NSBH events has no effect are not included in this plot.

Table 2 .
P (m2 ≤ m break ) P (m1 ≤ m break ) P (m2 ≤ m break < m1) Events classified in the Broken Power Law formalism.Values are the probabilities that each event's component masses lie in a given region of parameter space, considering uncertainty in the component masses as well as in the inferred values of m break .If m break is interpreted as the delimiter between NS and BH masses, the probabilities can be understood as labeled in the second row.The uncertainties in these probabilities are estimated to be less than 1%.Only events with > 5% probability of containing a NS according to this classification are included in this table.

Table 3 .
Events classified in the Broken Power Law + Dip formalism.Values are the probabilities that each event's component masses lie in a given region of parameter space, considering uncertainty in the component masses as well as in the inferred values of γ low and γ high .If γ low and γ high are interpreted as the lower and upper edges of the mass gap and γ low is therefore interpreted as the maximum observed NS mass in merging binaries, the probabilities can be understood as labeled in the second row.The uncertainties in these probabilities are estimated to be less than 1%.Only events with > 5% probability of being in one of these categories are included in this table.