Inferring the Neutron Star Maximum Mass and Lower Mass Gap in Neutron Star–Black Hole Systems with Spin

It remains unclear whether NSs, specifically those in merging BNS and NSBH systems, can have significant spins. The most rapidly spinning NS in a (nonmerging) double NS system is the pulsar J1807–2500B with a period of 4.2 ms or dimensionless spin magnitude a = 0.12 (Lynch et al. 2012). Among recycled pulsars, the fastest spinning is pulsar J1748–2446ad with a period of ∼1.4 ms (Hessels et al. 2006). However, rapidly spinning NSs in which spin-down is inefficient (due to, e.g., weak magnetic fields) may have avoided electromagnetic discovery for the same reasons. In NSBH systems, it may also be possible for the NS spin to grow through accretion if the NS is born before the BH (Chattopadhyay et al. 2021) or through tidal synchronization, as has been studied in BBH systems (Qin et al. 2018).

We remain agnostic about NS spin magnitudes, modeling their distribution as a power law,

$$\begin{eqnarray}p({a}_{2}| {a}_{\max },{\beta }_{s})\propto \left\{\begin{array}{ll}{\left(1-{a}_{2}\right)}^{{\beta }_{s}} & 0\lt {a}_{2}\lt {a}_{\max }\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where a_max sets an upper limit on possible values of a₂, and β_s controls the slope. For β_s = 0, the secondary spin magnitude follows a uniform distribution; for β_s > 0, the secondary spin distribution prefers low spin. The maximum value of a_max is the breakup spin a_Kep, which is around a_Kep ≈ 0.7 for most EOSs.

We do not explicitly model NS spin tilts (the angle between the spin vector and the orbital angular momentum axis), but consider a few different assumptions and explore how they affect our inference. By default, we consider an NS spin-tilt distribution that is isotropic, or flat in $-1\lt \cos ({\mathrm{tilt}}_{2})\lt 1$. We also explore a restricted model in which NS spins are perfectly aligned with the orbit, $\cos ({\mathrm{tilt}}_{2})=1$. For the distribution of BH spins, by default, we assume that BHs are nonspinning (a₁ = 0; Fuller & Ma 2019; Mandel & Smith 2021). We alternatively assume that the BH spin distribution is uniform in spin magnitude with isotropic spin tilts.

In summary, we consider the following spin models:

• 1.  
  Zero-spin BH ("ZS," default spin model): Primary BH is nonspinning (a₁ = 0). Secondary NS spin is isotropic in spin tilt (flat in $-1\lt \cos ({\mathrm{tilt}}_{2})\lt 1$) and follows a power law in the spin magnitude a₂ (Equation (1)).
• 2.  
  Zero-spin BH + aligned-spin NS ("ZS + AS"): Same as the default, but with $\cos ({\mathrm{tilt}}_{2})=1$.
• 3.  
  Uniform and isotropic ("U+I"): Same as the default, but the primary BH spin is flat in magnitude a₁ and $\cos ({\mathrm{tilt}}_{1})$ rather than nonspinning.

Like the case with spins, we consider a few different mass models to check the robustness of our conclusions. We consider three models for NS masses, which describe the distribution of NSBH secondary masses m₂ (see Figure 1):

• 1.  
  Default: Single Gaussian distribution (panels (a), (d)–(f) of Figure 1)

  $$\begin{eqnarray*}&&\begin{array}{l}p({m}_{2}| \mu ,\sigma ,{M}_{\min },{M}_{\max })\\ =\left\{\begin{array}{ll}{{ \mathcal N }}_{T}({m}_{2}| \mu ,\sigma ) & {M}_{\min }\leqslant {m}_{2}\leqslant {M}_{\max }\\ 0 & \mathrm{otherwise},\end{array}\right.\end{array}\end{eqnarray*}$$

  where ${{ \mathcal N }}_{T}(x| \mu ,\sigma )$ denotes a truncated Gaussian distribution with mean μ and standard deviation σ.
• 2.  
  Two-component (bimodal) Gaussian distribution ("2C"), as in the Galactic NS distribution (Alsing et al. 2018, panel (b) of Figure 1),

  $$\begin{eqnarray*}&&\begin{array}{l}p({m}_{2}| { \mathcal A },{\mu }_{1},{\sigma }_{1},{\mu }_{2},{\sigma }_{2},{M}_{\min },{M}_{\max })\\ =\left\{\begin{array}{ll}{ \mathcal A }{{ \mathcal N }}_{T}({m}_{2}| {\mu }_{1},{\sigma }_{1})+ & \\ (1-{ \mathcal A }){{ \mathcal N }}_{T}({m}_{2}| {\mu }_{2},{\sigma }_{2}) & {M}_{\min }\leqslant {m}_{2}\leqslant {M}_{\max }\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray*}$$
• 3.  
  Uniform distribution ("U") with sharp cutoffs at the minimum and maximum NS mass $[{M}_{\min },{M}_{\max }]$ (panel (c) of Figure 1)

  $$\begin{eqnarray*}&&\begin{array}{l}p({m}_{2}| {M}_{\min },{M}_{\max })\\ =\left\{\begin{array}{ll}\tfrac{1}{{M}_{\max }-{M}_{\min }} & {M}_{\min }\leqslant {m}_{2}\leqslant {M}_{\max }\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray*}$$

All normal distributions (${{ \mathcal N }}_{T}$) are truncated sharply and normalized to integrate to 1 between ${M}_{\min }=1{M}_{\odot }$ and ${M}_{\max }$. In this work, we focus on inferring the maximum NS mass. While the minimum NS mass can also be inferred with GWs (Chatziioannou & Farr 2020), we fix the minimum NS mass to 1 M_⊙ in our models. If binary stellar evolution can produce NSs with extreme masses, then M_min and M_max correspond to the minimum and maximum allowable masses set by nuclear physics.

Zoom In Zoom Out Reset image size

Crucially, we allow NSs to have significant spin. Rapid uniform rotation may provide additional support to the NS, allowing it to reach masses greater than the nonspinning maximum mass M_TOV. We model the dependence of M_max on NS spin a₂ using the universal relationship from Most et al. (2020):

$$\begin{eqnarray}&&\begin{array}{l}{M}_{\max }({a}_{2},{a}_{\mathrm{Kep}},{M}_{\mathrm{TOV}})\\ \quad =\ {M}_{\mathrm{TOV}}\left(1+{A}_{2}{\left(\displaystyle \frac{{a}_{2}}{{a}_{\mathrm{Kep}}}\right)}^{2}+{A}_{4}{\left(\displaystyle \frac{{a}_{2}}{{a}_{\mathrm{Kep}}}\right)}^{4}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

with A₂ = 0.132, A₄ = 0.0071, where a_Kep corresponds to the dimensionless spin at the mass-shedding limit. For concreteness, we assume a_Kep = 0.7, which is true for most EOSs. For an NS with spin a_Kep, the maximum possible mass is around 1.2× the (nonspinning) TOV limit. To measure this relation directly from gravitational-wave data, we also optionally measure a free, linear dependence between maximum spin and critical mass (see Section 4.5):

$$\begin{eqnarray}&&{M}_{\max }({a}_{2},{a}_{\mathrm{Kep}},{M}_{\mathrm{TOV}})={M}_{\mathrm{TOV}}{A}_{1}\left(\displaystyle \frac{{a}_{2}}{{a}_{\mathrm{Kep}}}\right).\end{eqnarray} \tag{ 3 }$$

The extent to which the NS mass distribution can extend above M_TOV depends on the spin distribution. The NS mass distributions p(m₂) above include a dependence on spin and can be written as $p({m}_{2}| {M}_{\max }({a}_{2}),\theta )$, where θ includes all other parameters. Figures 1(d)–(f) shows the NS mass distribution under three variations of the spin distributions outlined in 2.1.1.

We model the primary (BH) mass distribution p(m₁) as a power law with a slope of −α and a minimum mass cutoff at M_BH:

$$\begin{eqnarray}p({m}_{1}| \alpha ,{M}_{\mathrm{BH}})\propto \left\{\begin{array}{ll}0 & x\lt {M}_{\mathrm{BH}}\\ {{m}_{1}}^{-\alpha } & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

We fix α > 0 such that the probability density decreases for increasing BH mass. The minimum BH mass represents the upper boundary of the mass gap. In order to restrict the range of m₁ to reasonable values, we optionally include a maximum BH mass of 30 M_⊙ in Equation (4). However, for most of our NSBH models, high-mass BHs are rare due to a relatively steep slope α and/or a pairing function that disfavors extreme mass-ratio pairings, and we do not explicitly model the BH maximum mass.

We assume that the pairing function between m₁ and m₂ NSBH systems follows a power law in the mass ratio m₂/m₁ = q < 1 (Fishbach & Holz 2020):

$$\begin{eqnarray}&&p(q)\propto {q}^{\beta },\end{eqnarray} \tag{ 5 }$$

where by default we assume β = 0 (Farah et al. 2021). We alternatively consider the case β = 3, which favors equal-mass pairings. Depending on the width of the mass gap, NSBHs may necessarily have unequal masses, but on a population level, a higher q may still be relatively preferred.

Putting the mass and spin distributions together, we model the distribution of NSBH masses and spins θ ≡ (m₁, m₂, a₁, a₂) given population hyperparameters Λ and model H as

$$\begin{eqnarray}&&\begin{array}{r}\pi (\theta | {\rm{\Lambda }},H)\propto p({m}_{1}| \alpha ,{M}_{\mathrm{BH}},H)p({m}_{2}| {{\rm{\Lambda }}}_{\mathrm{NS}},{a}_{2},H)\\ \ \ \ \ \ p({a}_{1}| H)p({a}_{2}| {a}_{\max },{\beta }_{s},H)p(q| \beta ,H),\end{array}\end{eqnarray} \tag{ 6 }$$

where H refers to the choice of model as described in the earlier subsections. For the extrinsic source parameters not in θ, we assume isotropic distributions in sky position, inclination and orientation, and the local-universe approximation to a uniform-in-volume distribution $p({d}_{L})\propto {d}_{L}^{2}$, where d_L is the luminosity distance.

We infer properties of the overall NSBH population with a hierarchical Bayesian approach (Loredo 2004; Mandel et al. 2019). This allows us to marginalize over the uncertainties in individual events' masses and spins (grouped together in the set θ_i for event i) in order to estimate the hyperparameters Λ describing the NS and BH mass and spin distributions. For N_det GW detections producing data d, the likelihood of the data is described by an inhomogeneous Poisson process:

$$\begin{eqnarray}&&{ \mathcal L }(d| {\rm{\Lambda }},N)={N}^{{N}_{\det }}{e}^{-N\xi ({\rm{\Lambda }})}\prod _{i=1}^{{N}_{\det }}\int { \mathcal L }({d}_{i}| {\theta }_{i})\pi ({\theta }_{i}| {\rm{\Lambda }})\ d{\theta }_{i},\end{eqnarray} \tag{ 7 }$$

where N is the total number of NSBH mergers in the universe within some observing time, ξ(Λ) is the fraction of detectable events in the population described by hyperparameters Λ (see Section 2.2.2), ${ \mathcal L }({d}_{i}| {\theta }_{i})$ is the likelihood for event i given its masses and spins θ_i , and π(θ∣Λ) describes the NSBH mass and spin distribution given population hyperparameters Λ (Equation (6). As we do not attempt to calculate event rates, we marginalize over N with a log-uniform prior and calculate the population likelihood as (Mandel et al. 2019; Fishbach et al. 2018)

$$\begin{eqnarray}&&{ \mathcal L }(d| {\rm{\Lambda }})\propto \prod _{i=1}^{{N}_{\det }}\displaystyle \frac{\int { \mathcal L }({d}_{i}| {\theta }_{i})\pi ({\theta }_{i}| {\rm{\Lambda }})\ d{\theta }_{i}}{\xi ({\rm{\Lambda }})}.\end{eqnarray} \tag{ 8 }$$

We evaluate the single-event likelihood ${ \mathcal L }(d| \theta )$ via importance sampling over N_samp parameter estimation samples θ_PE for each event:

$$\begin{eqnarray}&&\int { \mathcal L }(d| \theta )\pi (\theta | {\rm{\Lambda }})\ d\theta \simeq \displaystyle \frac{1}{{N}_{\mathrm{samp}}}\sum _{j=1}^{{N}_{\mathrm{samp}}}\displaystyle \frac{\pi ({\theta }_{\mathrm{PE},j}| {\rm{\Lambda }})}{{\pi }_{\mathrm{PE}}({\theta }_{\mathrm{PE},j})},\end{eqnarray} \tag{ 9 }$$

where π_PE(θ) is the original prior that was used in LIGO parameter estimation. We calculate the posterior on the population parameters, p(Λ∣d), from the likelihood ${ \mathcal L }(d| {\rm{\Lambda }})$, under Bayes theorem, using broad, flat priors on the parameters Λ. For prior ranges, see Table 1.

Table 1. Priors Ranges for Population Parameters

${ \mathcal A }$	[0.0, 1.0]
μ or μ1, μ2	[1.0, 3.0]
σ or σ1, σ2	[0.01, 1.5]
MTOV	[1.5, 3.5]
MBH	[1.5, 10]
α	[0, 10]
max a/aKep	[0.1, 1.0]
βs	[0.0, 5.0]
A1(optional)	[−0.5, 0.5]

Download table as:  ASCIITypeset image

While we model and measure the astrophysical source distributions, GW detectors observe only sources loud enough to be detected, i.e., sources that produce data above some threshold d > thresh. We account for this selection effect by including the term ξ(Λ), the fraction of detectable binaries from a population described by parameters Λ:

$$\begin{eqnarray}\begin{array}{rcl}\xi ({\rm{\Lambda }}) & = & {\displaystyle \int }_{d\gt \mathrm{thresh}}{ \mathcal L }(d| \theta )\pi (\theta | {\rm{\Lambda }})\ {dd}\ d\theta \\ & \equiv & \displaystyle \int {P}_{\det }(\theta )\pi (\theta | {\rm{\Lambda }})d\theta .\end{array}\end{eqnarray} \tag{ 10 }$$

To evaluate ξ(Λ), we calculate the detection probability ${P}_{\det }(\theta )$ as a function of masses and cosmological redshift following the semianalytic approach outlined in Fishbach & Holz (2017). We assume the detection threshold is a simple single-detector signal-to-noise ratio (S/N) threshold ρ_thresh = 8. We ignore the effect of spin on detectability; although systems with large aligned spins experience orbital hang-up that increases their S/N compared to small or antialigned spins, the effect is small compared to current statistical uncertainties (Ng et al. 2018).

Given the masses and redshift of a potential source, we calculate its detectability as follows. We first calculate the optimal matched-filter S/N ρ_opt using noise power spectral density (PSD) curves corresponding to aLIGO at O3 sensitivity, Design sensitivity, or A+ sensitivity (Abbott et al. 2020b); the optimal S/N corresponds to a face-on, directly overhead source. We then calculate the S/N ρ for a random sky position and orientation by generating angular factors 0 < w < 1 from a single-detector antenna pattern (Finn & Chernoff 1993) and set ρ = w ρ_opt. If ρ > ρ_thresh for a given detector noise curve, we consider the simulated source to be detected.

Finally, we estimate ξ(Λ) with a Monte Carlo integral over simulated sources. We draw simulated sources with m₁, m₂, and z according to p_draw(θ) until injection sets of 10,000 events are created. BH (m₁) is drawn from a power law with M_BH = 1.5M_⊙. NS (m₂) is drawn from a uniform distribution between 1 and 3.5M_⊙. Redshifts z are drawn uniformly in comoving volume and source-frame time. Each simulated system is labeled as detected or not based on its S/N, described above. We then approximate the integral ξ(Λ) as a sum over M_det detected simulated systems:

$$\begin{eqnarray}&&\xi ({\rm{\Lambda }})\simeq \displaystyle \frac{1}{{N}_{\mathrm{draw}}}\sum _{j=1}^{{M}_{\det }}\displaystyle \frac{\pi ({m}_{1,j},{m}_{2,j},{z}_{j}| {\rm{\Lambda }})}{{p}_{\mathrm{draw}}({m}_{1,j},{m}_{2,j},{z}_{j})}.\end{eqnarray} \tag{ 11 }$$

While the population distributions in Section 2.1 are defined in terms of m₁, m₂, a₁, and a₂, gravitational-wave detectors are most sensitive to degenerate combinations of these parameters. These include the gravitational chirp mass

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{{({m}_{1}{m}_{2})}^{3/5}}{{({m}_{1}+{m}_{2})}^{1/5}},\end{eqnarray} \tag{ 12 }$$

the symmetric mass ratio

$$\begin{eqnarray}&&\nu =\displaystyle \frac{q}{{(1+q)}^{2}},\end{eqnarray} \tag{ 13 }$$

and χ_eff, a mass-weighted sum of the component spins that is approximately conserved during the inspiral,

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{{m}_{1}{a}_{1,z}+{m}_{2}{a}_{2,z}}{{m}_{1}+{m}_{2}},\end{eqnarray} \tag{ 14 }$$

where a_1,z and a_2,z are the components of the primary and secondary spins that are aligned with the orbital angular momentum axis. If the primary is nonspinning, χ_eff reduces to $\displaystyle \frac{{m}_{2}{a}_{2,z}}{{m}_{1}+{m}_{2}}={a}_{2,z}\tfrac{q}{1+q}$.

We follow the method outlined in Chatziioannou & Farr (2020) to simulate realistic parameter estimation samples from mock GW NSBH detections. Chatziioannou & Farr (2020) use the post-Newtonian (PN) description of the GW inspiral, with PN coefficients ψ₀, ψ₂, and ψ₃ that depend on the masses and spins:

$$\begin{eqnarray}&&{\psi }_{0}({ \mathcal M })=\displaystyle \frac{3}{128{{ \mathcal M }}^{5/3}{\pi }^{5/3}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\psi }_{2}({ \mathcal M },\nu )=\displaystyle \frac{5}{96{ \mathcal M }\pi {\nu }^{2/5}}\left(\displaystyle \frac{743}{336}+\displaystyle \frac{11\nu }{4}\right)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\beta =\displaystyle \frac{1}{3}\left(\displaystyle \frac{113-76\nu }{4}{\chi }_{\mathrm{eff}}+\displaystyle \frac{76}{4}\delta m\nu {\chi }_{a}\right)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\psi }_{3}({ \mathcal M },\nu ,\beta )=\displaystyle \frac{3\left(4\beta -16\pi \right)}{128{{ \mathcal M }}^{2/3}{\pi }^{2/3}{\nu }^{3/5}},\end{eqnarray} \tag{ 18 }$$

where the mass difference δ m = (m₁ − m₂)/(m₁ + m₂) and the spin difference χ_a = (a_1,z − a_2,z )/2. The third coefficient ψ₃ encodes the spin–orbit degeneracy as β includes the spins and ν is the mass ratio. In our case, unlike in Chatziioannou & Farr (2020), the χ_a term is not negligible. For NSBH systems, especially under the assumption of a spinning secondary and nonspinning primary, both the mass difference δ m and spin difference χ_a are significant. For our mock events, we approximate the measured PN coefficients ψ_i as independent Gaussian distributions with standard deviations σ_i . As in Chatziioannou & Farr (2020), we adopt σ₀ = 0.0046ψ₀/ρ, σ₂ = 0.2341ψ₂/ρ, and σ₃ = −0.1293ψ₃/ρ, where we draw the S/N ρ according to p(ρ) ∝ ρ⁻⁴, an approximation to the S/N distribution of a uniform-in-comoving-volume distribution of sources (Chen & Holz 2014). We then sample m₁, m₂, a_1,z , and a_2,z from the ψ₀, ψ₂, and ψ₃ likelihoods, accounting for the priors induced by the change of variables by calculating the appropriate Jacobian transformations.

An example NSBH parameter estimation posterior generated according to this procedure is shown in Figure 2. We see that the masses and spins are highly correlated. In particular, the anticorrelation between the secondary mass and spin increases the uncertainty on M_TOV and the spin–maximum mass relationship.

Zoom In Zoom Out Reset image size

For each model and data set variation, we infer the minimum BH mass M_BH, the NS M_TOV, and their difference (representing the width of the mass gap), marginalizing over all other parameters of the mass and spin distribution. Results for our Default model are shown in Figures 5–10, with Figure 10 showing a corner plot over all model parameters.

Zoom In Zoom Out Reset image size

Maximum (spin-dependent) NS mass: As discussed in 2.1.2, at a given secondary spin, we model a hard cutoff in the NS mass distribution p(m₂). However, in the one-component and two-component models, some values of μ and σ taper off the mass distribution between 2 and 3 M_⊙, making it difficult to discern a sharp truncation mass M_TOV from the function's normal behavior. This results in long, flat tails to large posterior values of M_TOV (see panel (a) of Figure 5), reaching the prior bounds even if priors on M_TOV are widened. A better-measured parameter is the 99th percentile of the nonspinning NS mass distribution, M₉₉ (panel (b) of Figure 5). For models where the M_TOV cutoff is significant, the 99th percentile is essentially identical to M_TOV. For models producing a softer cutoff without significant M_TOV truncation, the 99th percentile still captures the largest NS we expect to observe, and, unlike M_TOV, the inference of M₉₉ is consistent between the three NS mass models.

For models including GW190814, we generally infer M₉₉ between 2.6 and 2.8M_⊙, with lower limits (95% credibility) of 2.4–2.5 M_⊙. Our default model (all four events, β = 0, "ZS" spin prior) measures ${M}_{99}={2.8}_{-0.2}^{+0.3}\,{M}_{\odot }$ (68% credibility); the inclusion of three additional GWTC-3 events shifts M₉₉ to ${2.9}_{-0.2}^{+0.2}{M}_{\odot }$. The cutoff mass is set by GW190814, where m₂ is extremely well constrained. Without GW190814, we estimate M_TOV between 2.0 and 2.3 M_⊙, with a lower limit (95% credibility) of 1.8–1.9 M_TOV. Without GW190814, our estimates are consistent with other estimates of M_TOV from gravitational-wave NS observations that do not consider spin.

The spin distribution affects the inferred value of M_TOV and M₉₉. For all four events, m₂ is consistent with being both nonspinning (a₂ = 0) or maximally spinning (a₂ = 0.7, a/a_Kep = 1). When the spin distribution allows or favors a maximally spinning NS, lower values of M_TOV are allowed and can still account for GW190814, the most massive secondary. When the spin distribution disfavors high spins, the spin-dependent maximum mass is lower and M_TOV must be higher in order to accommodate GW190814.

This is shown in Figure 5; the posterior on M_TOV inferred under a uniform spin distribution (β_s = 0, ${a}_{\max }/{a}_{\mathrm{Kep}}=1$), which has support at high NS spins, has a significant tail to lower values below 2.5 M_⊙ (dashed blue curve). A prior that requires GW190814 to be maximally spinning (${a}_{\min }/{a}_{\mathrm{Kep}}=0.9$) brings M_TOV estimates even lower, to ∼2.4 M_⊙, with support below 2.2 M_⊙ (green dashed curve in Figure 5). Meanwhile, requiring all NSs to be nonspinning (${a}_{\max }/{a}_{\mathrm{Kep}}=0$) means that GW190814's secondary (if it is an NS) sets the nonspinning maximum mass for the population and results in a narrower posterior preferring larger values. The difference between posteriors on M_TOV and M₉₉ modeled with GW190814 (black solid, red dotted, blue dashed) and without GW90814 (solid yellow curve) is bridged partially by models assuming GW190814's spin is near maximal. This effect is also visible in Figure 9; if GW190814 is assumed spinning, the upper end of p(m₂) visibly shifts to lower masses, and zero-spin NS mass functions truncating below GW190814's secondary's mass are allowed (see overplotted credible interval). We see that even in the absence of well-constrained a₂, modeling a spin-dependent maximum mass has significant effects on the inferred NS mass distribution.

GW190814's secondary spin: Using the posterior on M_TOV (Figure 5) inferred from the population of NSBHs excluding GW190814, we can infer the minimum secondary spin of GW190814 required for it to be consistent with the NSBH population. Results are shown in panel (a) of Figure 6. For our sample of NSBH events excluding GW190814, the results are inconclusive: Because the posterior on M_TOV is broad, GW190814 is consistent even if nonspinning (with the minimum required a₂/a_Kep = 0), but it may also be maximally spinning with a₂/a_Kep = 1. GW190814 may also be an outlier from the NSBH population, even if it is maximally spinning: For this figure, we allow min a₂/a_Kep > 1, but a₂/a_Kep > 1 would imply inconsistency with the rest of the population as ${a}_{\max }/{a}_{\mathrm{Kep}}=1$. Future GW observations of a larger population of NSBH events (see panel (b) of Figure 6) could allow a much tighter measurement of GW190814's secondary spin.

Zoom In Zoom Out Reset image size

Minimum BH mass: Across all models, the inferred BH minimum mass M_BH is between 4 and 7 M_⊙ with typical uncertainties of ± 1M_⊙. Our default model using the "ZS" spin model, all four NSBH events (GW190814, GW190426, GW200105, GW200115), and random pairing (β = 0) results in ${M}_{\mathrm{BH}}={5.4}_{-1.0}^{+0.7}\,{M}_{\odot }$ (68% credibility). At the low end, we infer ${M}_{\mathrm{BH}}={4.2}_{-1.0}^{+1.1}\,{M}_{\odot }$ using all four NSBH events, a uniform NS mass distribution, pairing function β = 3, and the "U + I" spin model. At the high end, we infer ${M}_{\mathrm{BH}}={6.7}_{-0.8}^{+0.4}\,{M}_{\odot }$ (68% credibility) using only the confident NSBH events and the spin model. The effect of the m₁, m₂ pairing function β is minimal, but assuming equal-mass pairings further reduces posterior support for low M_BH (see Figure 7).

Zoom In Zoom Out Reset image size

Mass gap: We estimate the inferred width of the lower mass gap as the difference between the minimum BH mass, M_BH, and the maximum nonspinning NS mass, M_TOV or M₉₉. The mass gap's width may range from 0 to a few M_⊙, while the mass gap's position may range from 2 to 7 M_⊙. As seen in Figures 5 and 7; the overlap between the posteriors on M_BH and M₉₉ is low, suggesting the existence of a mass gap. Similarly, panels (a) and (b) in Figure 8 show inferred (m₁, m₂) posterior predictive distributions, overplotted with the LVK m₁, m₂ posteriors. As Figure 8 illustrates, for all model variations, we find evidence for a separation between the upper end of the NS mass distribution and the lower end of the BH mass distribution.

Zoom In Zoom Out Reset image size

For our default model, we measure a mass gap of ${2.5}_{-1.0}^{+0.8}{M}_{\odot }$ (${2.3}_{-1.0}^{+0.7}{M}_{\odot }$ with 3 additional GWTC-3 events), wider than 0 M_⊙ with 97% credibility and 1 M_⊙ with 90% credibility. The inferred mass gap is widest when only using the confident NSBH events, between 3.0 and 4.5 M_⊙, and narrowest when using all four NSBH events, between 1.5 and 3.0 M_⊙. This is because the mass gap is narrowed from the NS side by the inclusion of GW190814 and from the BH side by the inclusion of GW190426 (see Figure 3. All model variations (spin prior, β, events) support the existence of a mass gap: >0 M_⊙ with 92% or higher (up to >99.9%) credibility, and >1 M_⊙ with 68% or higher (up to >99.9%) credibility.

As seen in Figure 3, additional spin assumptions (namely assuming that the BH is nonspinning and/or the NS spin is aligned) tend to prefer lower m₂ and higher m₁, which widens the inferred mass gap. When using spin priors in which the BH is assumed to be nonspinning, even when modeling all four events (including GW190814) we infer a mass gap exists with >96% credibility and that it is wider than 1 M_⊙ with >90% credibility.

In addition to the most astrophysically relevant parameters—M_BH, M_TOV, and the width of the mass gap—we also constrain other parameters of the primary and secondary mass functions. In this section, we discuss general trends in the mass distribution shape, as inferred from posterior traces (Figure 9).

Zoom In Zoom Out Reset image size

We first consider the NS mass distribution, p(m₂), which differs slightly depending on the mass model used. For the one-component model, we generally infer a broad distribution (σ ≃ 0.5) with mean μ between 1.2 and 1.6M_⊙. A broad distribution is especially necessary to explain the large secondary mass of GW190814. The two-component model generally agrees well with the one-component model, although additional substructure (see panel (a) of Figure 9), particularly a narrower peak at around 1.3 M_⊙ and a longer tail to high NS masses (above 2M_⊙), is possible. The only free parameter in the uniform model is the cutoff mass M_TOV. Though the flatness of the uniform model means we necessarily infer higher probability at masses near M_TOV, M_TOV is generally consistent with the upper limit (99th percentile M₉₉) inferred from other mass models.

The BH mass function is consistent between the three NS mass models. The most significant influence is the pairing function (β = 0 for random or β = 3 for equal-mass preference). For example, under our default model (four events), which includes random pairing (β = 0), we infer a distribution with power-law slope ${\alpha }_{\mathrm{BH}}={3.4}_{-0.9}^{+1.4}$ (${\alpha }_{\mathrm{BH}}={2.3}_{-1.0}^{+7}$ with all seven events). Under the same assumptions but preferring equal masses, β = 3, the inferred distribution shifts to significantly shallower slopes, ${\alpha }_{\mathrm{BH}}={0.9}_{-0.6}^{+1.1}$. This is because the preference for equal-mass pairing requires a shallower slope in order to account for higher-mass black holes, especially the primary of GW190814.

As seen in Figure 10, the joint posterior on β_spin and ${a}_{\max }/{a}_{\mathrm{Kep}}$ prefers low ${a}_{\max }/{a}_{\mathrm{Kep}}$ and high β_spin, but mainly recovers the flat prior, which inherently prefers steeper and smaller spin distributions. Thus, our measurement of the NS spin distribution is mostly uninformative, with a very mild preference for small spins.

Zoom In Zoom Out Reset image size