Abstract
We report the observation of a coalescing compact binary with component masses 2.5–4.5 M⊙ and 1.2–2.0 M⊙ (all measurements quoted at the 90% credible level). The gravitational-wave signal GW230529_181500 was observed during the fourth observing run of the LIGO–Virgo–KAGRA detector network on 2023 May 29 by the LIGO Livingston observatory. The primary component of the source has a mass less than 5 M⊙ at 99% credibility. We cannot definitively determine from gravitational-wave data alone whether either component of the source is a neutron star or a black hole. However, given existing estimates of the maximum neutron star mass, we find the most probable interpretation of the source to be the coalescence of a neutron star with a black hole that has a mass between the most massive neutron stars and the least massive black holes observed in the Galaxy. We provisionally estimate a merger rate density of for compact binary coalescences with properties similar to the source of GW230529_181500; assuming that the source is a neutron star–black hole merger, GW230529_181500-like sources may make up the majority of neutron star–black hole coalescences. The discovery of this system implies an increase in the expected rate of neutron star–black hole mergers with electromagnetic counterparts and provides further evidence for compact objects existing within the purported lower mass gap.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
In 2023 May, the fourth observing run (O4) of the Advanced LIGO (Aasi et al. 2015), Advanced Virgo (Acernese et al. 2015), and KAGRA (Somiya 2012; Aso et al. 2013) observatory network commenced following a series of upgrades to increase the sensitivity of the network. The prior three observing runs opened the field of observational gravitational-wave (GW) astronomy, with 90 probable compact binary coalescence (CBC) candidates reported by the LIGO Scientific, Virgo, and KAGRA Collaboration (LVK) at the conclusion of the third observing run (O3; Abbott et al. 2023a) and further candidates found by external analyses (e.g., Mehta et al. 2023b; Nitz et al. 2023; Wadekar et al. 2023). These included the first observation of merging stellar-mass black holes (Abbott et al. 2016a), the first observation of a stellar-mass black hole merging with a neutron star (Abbott et al. 2021a), and the first observation of two merging neutron stars (Abbott et al. 2017a). The first observation of two merging neutron stars was also a multimessenger event accompanied by emission across the electromagnetic (EM) spectrum (Abbott et al. 2017b; Margutti & Chornock 2021). The continued discovery of CBCs by the international GW detector network in O4 and beyond promises to reveal new information about the formation pathways of compact binaries and the physics of their evolution.
While GW observations have enabled detailed characterization of the population of stellar-mass compact-object binary mergers overall (e.g., Abbott et al. 2023b), the small number of observed mergers at the low-mass end of the black hole mass spectrum leads to considerable uncertainty in this region of the mass distribution. Dynamical mass measurements of X-ray binary systems within the Milky Way suggest a paucity of compact objects with masses between ∼3 and 5 M⊙, and hence a lower mass gap that divides the population of neutron stars (with masses less than ∼3 M⊙; Rhoades & Ruffini 1974; Kalogera & Baym 1996) and stellar-mass black holes (observed to have masses above 5 M⊙; Bailyn et al. 1998; Ozel et al. 2010; Farr et al. 2011b). Although a number of recent observations of noninteracting binary systems (Thompson et al. 2019; Jayasinghe et al. 2021), radio pulsar surveys (Barr et al. 2024), and GW observations (R. Abbott et al. 2020c; Abbott et al. 2023b) have found hints of compact objects residing in the lower mass gap, GW observations have yet to quantify the extent and potential occupation of this gap (Fishbach et al. 2020; Farah et al. 2022; Abbott et al. 2023b).
Here we report on the compact binary merger signal GW230529_181500, henceforth abbreviated as GW230529, which was detected by the LIGO Livingston observatory on 2023 May 29 at 18:15:00 UTC; all other observatories either were offline or did not have the sensitivity required to observe this signal. We find that the compact binary source of GW230529 had component masses of and (all measurements are reported as symmetric 90% credible intervals around the median of the marginalized posterior distribution with default uniform priors, unless otherwise specified). Although we cannot definitively determine the nature of the higher-mass (primary) compact object in the binary system, if we assume that all compact objects with masses below current constraints on the maximum neutron star mass are indeed neutron stars, the most probable interpretation for the source of GW230529 is the coalescence between a 2.5 to 4.5 M⊙ black hole and a neutron star. GW230529 provides further evidence that a population of compact objects exists with masses between the heaviest neutron stars and lightest black holes observed in the Milky Way. Furthermore, if the source of GW230529 was a merger between a neutron star and a black hole, its masses are significantly more symmetric than the neutron star–black holes (NSBHs) previously observed via GWs (Abbott et al. 2023a), which increases the expected rate of NSBHs that may be accompanied by an EM counterpart.
We report on the status of the detector network at the time of GW230529 in Section 2 and provide details of the detection in Section 3. In Section 4 we present estimates of the source properties, along with a discussion of the inferred masses, spins, tidal effects, and consistency of the signal with general relativity. We provide updated constraints on merger rates and the inferred properties of the compact binary and NSBH populations, as well as population-informed posteriors on source properties, in Section 5. Section 6 provides analysis and interpretation of the physical nature of the source components. Implications for multimessenger astrophysics and the formation of low-mass black holes are discussed in Sections 7 and 8, respectively. Section 9 summarizes our findings. Data from the analyses in this work are available on Zenodo (Abbott et al. 2024).
2. Observatory Status and Data Quality
At the time of GW230529 (2023 May 29 18:15:00.7 UTC), LIGO Livingston was in observing mode; LIGO Hanford was offline, having gone out of observing mode 1.5 hr prior to the detection. The Virgo observatory was undergoing upgrades and was not operational at the time of the detection. The KAGRA observatory was in observing mode, but its sensitivity was insufficient to impact the analysis of GW230529. Hence, only the data from the LIGO Livingston observatory are used in the analysis of GW230529.
LIGO Livingston was observing with stable sensitivity for ≃66 hr up to and including the time of the GW signal. At the time of the detection, the sky-averaged binary neutron star (BNS) inspiral range was 150 Mpc. From the start of O4 until this event, the LIGO Livingston BNS inspiral range (Chen et al. 2021a) varied between ≃140–160 Mpc, a 4.5%–19.4% increase compared to the median BNS range in O3 (Buikema et al. 2020; Abbott et al. 2023a). Additional details on upgrades to the LIGO observatories for O4 can be found in Appendix A.
The Advanced LIGO observatories are laser interferometers that measure strain (Aasi et al. 2015). The observatories are calibrated via photon radiation pressure actuation. An amplitude-modulated laser beam is directed onto the end test masses, inducing a known change in the arm length from the equilibrium position (Karki et al. 2016; Viets et al. 2018). For the strain data used in the GW230529 analysis, the maximum 1σ bounds on calibration uncertainties at LIGO Livingston were 6% in amplitude and in phase for the frequency range 20–2048 Hz.
Detection vetting procedures similar to those of past GW candidates (Abbott et al. 2016b; Davis et al. 2021), when applied to GW230529, find no evidence that instrumental or environmental artifacts (Effler et al. 2015; Nguyen et al. 2021; Helmling-Cornell et al. 2023) could have caused GW230529. We find no evidence of transient noise that is likely to impact the recovery of the GW signal in the 256 s LIGO Livingston data segment containing the signal.
3. Detection of GW230529
GW230529 was initially detected in low latency in data from the LIGO Livingston observatory. The signal was detected independently by three matched-filter search pipelines: GstLAL (Messick et al. 2017; Sachdev et al. 2019; Hanna et al. 2020; Cannon et al. 2021; Ewing et al. 2024; Tsukada et al. 2023), MBTA (Adams et al. 2016; Aubin et al. 2021), and PyCBC (Allen 2005; Allen et al. 2012; Usman et al. 2016; Nitz et al. 2017; Davies et al. 2020; Dal Canton et al. 2021). Although GW230529 occurred when only a single detector was observing, it was detected by all three pipelines with high significance, and stands out from the background distribution of noise triggers. More details are given in Appendix B.
All three pipelines have a similar matched-filter-based approach to identify GW candidates but differ in the details of implementation. Each pipeline begins by performing a matched-filtering analysis on the data from the observatory with a bank of GW templates (Roy et al. 2019; Florian et al. 2024; Sakon et al. 2024). Times of high signal-to-noise ratio (S/N) are identified, after which each pipeline performs its own set of signal consistency tests, combines them with the S/N to make a single ranking statistic for each candidate GW event, and calculates a false alarm rate (FAR) by comparing the ranking statistic of the candidate with the background distribution. The FAR is the expected rate of triggers caused by noise with a ranking statistic greater than or equal to that of the candidate.
For each pipeline, this procedure can be done in two modes: a low-latency or online mode, and an offline mode. In the online mode, the pipelines only use the background information available at the time of a candidate to estimate its significance, along with low-latency data-quality information. In the offline mode, pipelines use more background information gathered from subsequent observing times to estimate the significance of candidates and use more refined data-quality information. The offline mode enables more robust and reproducible results at the cost of greater latency.
In both cases, the background distribution is extrapolated to higher significances. This enables the inverse FAR to be greater than the time for which the background was collected. Pipelines have differing extrapolation methods and differing methods for calculating the FAR of single-detector GW candidates, details of which can be found in Appendix C. These differing methods can cause the FARs reported by different pipelines for the same candidate to vary to a significant degree, as seen in the case of GW230529. However, all three pipelines recover a nearly identical S/N for GW230529, as expected for an astrophysical signal with a high match to the search templates. This, in concert with the fact that the inverse FARs from the offline analyses of all pipelines are much greater than the duration of the first half of the fourth observing run (O4a), makes it highly likely that GW230529 is of astrophysical origin. Table 1 shows the online S/N, online inverse FAR, and offline inverse FAR for each search pipeline.
Table 1. Properties of the Detection of GW230529 from Each Search Pipeline
| GstLAL | MBTA | PyCBC | |
|---|---|---|---|
| Online S/N | 11.3 | 11.4 | 11.6 |
| Online inverse FAR (yr) | 1.1 | 1.1 | 160.4 |
| Offline inverse FAR (yr) | 60.3 | >1000 | >1000 |
Note. Significance is measured by the inverse of the FAR.
Download table as: ASCIITypeset image
4. Source Properties
The source properties of GW230529 are inferred using a Bayesian analysis of the data from the LIGO Livingston observatory. We analyze 128 s of data including frequencies in the range of 20–1792 Hz in the calculation of the likelihood. To describe the detector noise, we assume a noise power spectral density (PSD) given by the median estimate provided by BayesWave (Cornish & Littenberg 2015; Littenberg & Cornish 2015; Cornish et al. 2021). We use the Bilby (Ashton et al. 2019; Romero-Shaw et al. 2020) or Parallel Bilby (Smith et al. 2020) inference libraries to generate samples from the posterior distribution of the source parameters using a nested sampling (Skilling 2006) algorithm, as implemented in the Dynesty software package (Speagle 2020).
Given the uncertain nature of the compact objects of the source, we analyze GW230529 using a range of waveform models that incorporate a number of key physical effects. For our primary analysis, we employ binary black hole (BBH) waveform models that include higher-order multipole moments, the effects of spin-induced orbital precession, and allow for spin magnitudes on both components up to the Kerr limit but do not include tidal effects on either component.
Systematic errors in inferred source properties due to waveform modeling may be significant for NSBH systems (Huang et al. 2021). We mitigate these effects by combining the posteriors inferred using two different signal models: the phenomenological frequency-domain model IMRPhenomXPHM (García-Quirós et al. 2020; Pratten et al. 2020, 2021), and a time-domain effective-one-body (EOB) model SEOBNRv5PHM (Khalil et al. 2023; Pompili et al. 2023; Ramos-Buades et al. 2023; van de Meent et al. 2023). The posterior samples obtained independently using the two signal models are broadly consistent. Unless otherwise noted, we present results throughout this work obtained by an equal-weight combination of the posterior samples from both models under the default priors described in Appendix D. Our measurements of key source parameters for GW230529 are presented in Table 2.
Table 2. Source Properties of GW230529 from the Primary Combined Analysis (BBH Waveforms, High-spin, Default Priors)
| Parameter | Value |
|---|---|
| Primary mass m1/M⊙ | |
| Secondary mass m2/M⊙ | |
| Mass ratio q = m2/m1 | |
| Total mass M/M⊙ | |
| Chirp mass | |
| Detector-frame chirp mass | |
| Primary spin magnitude χ1 | |
| Effective inspiral-spin parameter χeff | |
| Effective precessing-spin parameter χp | |
| Luminosity distance DL/Mpc | |
| Source redshift z |
Note. We report the median values together with the 90% symmetric credible intervals at a reference frequency of 20 Hz.
Download table as: ASCIITypeset image
Analysis details and results from other waveform models we consider are reported in Appendix D; we find that the key conclusions of the analyses presented here are not sensitive to the choice of signal model. In particular, the use of BBH models is validated by comparison to waveform models that include tidal effects, finding no evidence that the BNS or NSBH models are preferred, consistent with previous observations (Abbott et al. 2021a). This is expected given the moderate S/N with which GW230529 was detected (Huang et al. 2021).
The analysis of GW230529 indicates that it is an asymmetric compact binary with a mass ratio and source component masses and . The primary is consistent with a black hole that resides in the lower mass gap (3 M⊙ ≲ m1 ≲ 5 M⊙; Ozel et al. 2010; Farr et al. 2011b), with a mass <5 M⊙ at the 99% credible level. The posterior distribution on the mass of the secondary is peaked around ∼1.4 M⊙ with an extended tail beyond 2 M⊙, such that P(m2 > 2 M⊙) = 5%. The mass of the secondary is consistent with the distribution of known neutron star masses, including Galactic pulsars (Antoniadis et al. 2016; Özel & Freire 2016; Alsing et al. 2018; Farrow et al. 2019) and extragalactic GW observations (Landry & Read 2021; Abbott et al. 2023b). Figure 1 shows the component mass posteriors of GW230529 relative to other BNSs (GW170817 and GW190425) and NSBHs (GW200105_162426 and GW200115_042309, henceforth abbreviated as GW200105 and GW200115) observed by the LVK, as well as GW190814 (Abbott et al. 2017a, 2020a; R. Abbott et al. 2020c; Abbott et al. 2021a).
Figure 1. The one- and two-dimensional posterior probability distributions for the component masses of the source binary of GW230529 (teal). The contours in the main panel denote the 90% credible regions, with vertical and horizontal lines in the side panels denoting the 90% credible interval for the marginalized one-dimensional posterior distributions. Also shown are the two O3 NSBH events GW200105_162426 and GW200115_042309 (orange and blue, respectively; Abbott et al. 2021a) with FAR < 0.25 yr−1 (Abbott et al. 2023a), the two confident BNS events GW170817 and GW190425 (pink and green, respectively; Abbott et al. 2017a, 2019a, 2020a, 2024b), and GW190814 (red; R. Abbott et al. 2020c; Abbott et al. 2024b), where the secondary component may be a black hole or a neutron star. Lines of constant mass ratio are indicated by dotted gray lines. The gray shaded region marks the 3–5 M⊙ range of primary masses. The NSBH events and GW190814 use combined posterior samples assuming a high-spin prior analogous to those presented in this work. The BNS events use high-spin IMRPhenomPv2_NRTidal (Dietrich et al. 2019a) samples.
Download figure:
Standard image High-resolution imageTo capture dominant spin effects on the GW signal, we present constraints on the effective inspiral spin χeff, which is defined as a mass-weighted projection of the spins along the unit Newtonian orbital angular momentum vector (Damour 2001; Racine 2008; Ajith et al. 2014),
where the dimensionless spin vector χ i of each component is related to the spin angular momentum S i by . If χeff is negative, it indicates that at least one of the spin component projections must be antialigned with respect to the orbital angular momentum, i.e., . We measure an effective inspiral spin of , which is consistent with a binary in which one of the spin components is antialigned or a binary with negligible spins. The measurement is primarily driven by the spin component of the primary compact object , with a probability that χ1,z < 0 of 83%. However, there is a degeneracy between the measured masses and spins of the binary components such that more comparable mass ratios correlate to more negative values of χeff (Cutler & Flanagan 1994) for this system. We show the correlation between χ1,z and the mass ratio and primary mass in Figure 2, with more negative values of χ1,z corresponding to more symmetric mass ratios and smaller primary masses. The secondary spin is only weakly constrained and broadly symmetric about 0, . We find no evidence for precession, with the posteriors on the effective precessing spin χp (Schmidt et al. 2015) being uninformative.
Figure 2. Selected source properties of GW230529. The plot shows the one-dimensional (diagonal) and two-dimensional (off-diagonal) marginal posterior distributions for the primary mass m1, the mass ratio q, and the spin component parallel to the orbital angular momentum . The shaded regions denote the posterior probability, with the solid (dashed) curves marking the 50% and 90% credible regions for the posteriors determined using a high-spin (low-spin) prior on the secondary of χ2 < 0.99 (χ2 < 0.05). The vertical lines in the one-dimensional marginal posteriors mark the 90% credible intervals.
Download figure:
Standard image High-resolution imageThe presence of a neutron star in a compact binary imprints tidal effects onto the emitted GW signal (Flanagan & Hinderer 2008). The strength of this interaction is governed by the tidal deformability of the neutron star, which quantifies how easily the star will be deformed in the presence of an external tidal field. In contrast, the tidal deformability of a black hole is zero (Binnington & Poisson 2009; Damour & Nagar 2009; Chia 2021), offering a potential avenue for distinguishing between a black hole and a neutron star. We investigate the tidal constraints for both the primary and secondary components using waveform models that account for tidal effects (Dietrich et al. 2019a; Matas et al. 2020; Thompson et al. 2020a), which do not qualitatively change the mass and spin conclusions discussed above. Irrespective of whether we analyze GW230529 with an NSBH model that assumes only the tidal deformability of the primary compact object to be zero or a BNS model that includes the tidal deformability of both objects, we find the tidal deformability of the secondary object to be unconstrained. The dimensionless tidal deformability of the primary peaks at zero, consistent with a black hole. The constraints on this parameter are also consistent with dense matter equation of state (EOS) predictions for neutron stars in this mass range.
We also perform parameterized tests of the GW phase evolution to verify whether GW230529 is consistent with general relativity and find no evidence of inconsistencies. More detailed information on tidal deformability analyses and testing general relativity can be found in Appendices E.3 and F, respectively.
5. Impact of GW230529 on Merger Rates and Populations
We provide a provisional update to the NSBH merger rate reported in our earlier studies (Abbott et al. 2021a, 2023b) by incorporating data from the first 2 weeks of O4a using two different methods. In the first, event-based approach, we consider GW230529 to be representative of a new class of CBCs and assume its contribution to the total number of NSBH detections to be a single Poisson-distributed count (Kim et al. 2003; Abbott et al. 2021a) over the span of time from the beginning of the first observing run (O1) through the first 2 weeks of O4a. We find the rate of GW230529-like mergers to be . When computing the rates of the significant NSBH events in O3 detected with FAR < 0.25 yr−1 (Abbott et al. 2021a, 2023a) using the same method, we find a total event-based NSBH merger rate of . Using this selection criterion, the same as used in Abbott et al. (2023b), the population of NSBH mergers includes GW200105 and GW200115; we do not include GW190814, as its source binary is most probably a BBH (R. Abbott et al. 2020c; Abbott et al. 2023b; Essick & Landry 2020).
For the second, broad population-based approach, we consider GW200105, GW200115, and GW230529 to be members of a single CBC class, together with an ensemble of less significant candidates. The definition of this class is determined by a simple cut on masses and spins resulting in an NSBH-like region of parameter space. We aggregate data from all the triggers found by GstLAL from the beginning of O1 through the first 2 weeks of O4a, assess their impact on the estimated merger rates of NSBHs while also accounting for the possibility that some of them are of terrestrial origin (Farr et al. 2015; Kapadia et al. 2020), and update the NSBH merger rate obtained at the end of O3 (Abbott et al. 2023a) to . Further details regarding the classification of triggers in the population-based rate method can be found in Appendix G.
We show updates to both the event-based and population-based NSBH rate constraints in Figure 3. The event-based rate estimates highlight that GW230529-like systems merge at a similar (or potentially higher) rate to the more asymmetric NSBHs identified in GWTC-3. The population-based rate estimate is more representative of the full NSBH merger rate, as it includes less significant GW events in the data. Both of our updated estimates are consistent with the findings of Abbott et al. (2021a) within the measurement uncertainties. Further details about rate estimates can be found in Appendix G.
Figure 3. Posterior on the merger rates of NSBH systems. The solid and dashed lines represent the broad population-based rate calculation and the event-based rate calculation, respectively.
Download figure:
Standard image High-resolution imageIn addition to updates to the overall merger rate of NSBHs, we also study the impact of GW230529 on the mass and spin distributions of the compact binary population as inferred from GWTC-3 (Abbott et al. 2023b). We employ hierarchical Bayesian inference to marginalize over the properties of individual events and infer the parameters of a given population model (e.g., Mandel et al. 2019; Thrane & Talbot 2019; Vitale et al. 2020). The updates to population model results in this work are provisional because they only include one GW signal from O4a, although the biases resulting from this selection will not be severe since GW230529 occurred near the start of the observing run. We quantify this by comparing the event-based rate estimate (which accounts for O4a sensitivity to compute its detectable time–volume) with the NSBH rates attained by the various population analyses considered in this work, finding that they are consistent with each other.
We use three different population models in our analysis. The first two models consider the population of compact-object binaries as a whole without distinguishing by source classification, using either the parameterized Power law + Dip + Break model (Fishbach et al. 2020; Farah et al. 2022; Abbott et al. 2023b) or the nonparametric Binned Gaussian Process model (Mohite 2022; Abbott et al. 2023b; Ray et al. 2023b). The Power law + Dip + Break model is designed to search for a separation in masses between neutron stars and black holes by explicitly allowing for, but not enforcing, a dip in the component mass distribution. The Binned Gaussian Process model is designed to capture the structure of the mass distribution with minimal assumptions about the population. The broad CBC population analyses include all candidates reported in GWTC-3 with FAR < 0.25 yr−1, the same selection criterion used in Abbott et al. (2023b, see Table 1 therein). The third model we investigate (NSBH-pop; Biscoveanu et al. 2022) considers only the population of NSBH mergers with FAR < 0.25 yr−1. NSBH-pop is a parametric model designed to constrain the population distributions of NSBH masses and black hole spin magnitudes. This model assumes that all analyzed events have a black hole primary and neutron star secondary; we do not include GW190814 in this analysis. Further details regarding the population model parameterizations and priors can be found in Appendix H.
The inclusion of GW230529 in our population analyses has several effects on the inferred properties of NSBH systems and the CBC population as a whole.
The inferred minimum mass of black holes in the NSBH population decreases with the inclusion of GW230529. The mass spectrum of black holes in the NSBH population with and without GW230529 inferred using the NSBH-pop model is shown in Figure 4. As the source of GW230529 is the NSBH with the smallest black hole mass observed to date, the minimum mass of black holes in NSBH mergers shows a significant decrease with the inclusion of this candidate: with GW230529 compared to without. In contrast, the parameter that governs the lower edge of the dip feature in the Power law + Dip + Break model, which represents the minimum black hole mass in the CBC population, does not shift significantly with the inclusion of GW230529. This difference is because the Power law + Dip + Break model makes no assumptions about the classification of the components and therefore does not enforce sharp features at the edges of each subpopulation, meaning that the source of GW230529 is not necessarily an NSBH in this model. However, the Binned Gaussian Process and power law + dip + break models are designed to capture the structure of the full compact binary mass spectrum; because they do not assign a source classification to either of the binary components, features present in these population models can have differing astrophysical interpretations from the NSBH-pop model.
Figure 4. The differential merger rate of NSBH systems as a function of black hole masses (solid curves: mean; shaded regions: 90% credible interval) and minimum black hole mass (dashed lines: 90% credible interval) using the NSBH-pop model with (teal) and without (gray) GW230529. The minimum black hole mass is a parameter of the NSBH-pop population model; see Table 7. While the median rate is always inside the credible region, the mean can be outside the credible region.
Download figure:
Standard image High-resolution imageGW230529 increases the inferred rate of compact binary mergers with a component in the 3–5 M⊙ range. A region of interest in the mass distribution is the border between the masses of neutron stars and black holes. We choose 3–5 M⊙ to represent the nominal gap between these populations. In Figure 5 we show the posterior on the rate of mergers with one or both component masses in the 3–5 M⊙ range, with and without GW230529. For the Power law + Dip + Break model there is a small increase in the merger rate within this mass range, with GW230529 versus without. For the Binned Gaussian Process model there is a larger increase in the merger rate, with GW230529 versus without. The differing degree of change between the two models is due to different assumptions for the mass ratio distribution of merging compact binaries, as well as the potential dip in the merger rate at low black hole masses built into the Power law + Dip + Break model. Binned Gaussian Process does not fit for a specific pairing function, while Power law + Dip + Break assumes the same mass ratio distribution throughout the whole population of CBCs. As most observed BBHs favor equal masses, any population model conditioned on those observations should also favor equal masses in the BBH mass range. The implicit assumptions within the Power law + Dip + Break model require this preference to be imposed for all masses, including relatively low mass systems like GW230529. However, the assumptions in Binned Gaussian Process do not broadcast this preference as strongly across different mass scales and therefore may support more asymmetric mass ratios at lower masses. Regardless of the mass ratio distribution assumptions made by each model, we find that the inclusion of GW230529 provides further evidence that the ∼3–5 M⊙ region is not completely empty.
Figure 5. Posterior on the merger rate of binaries with one or both components between 3 and 5 M⊙. The solid curves show the results from the Binned Gaussian Process analysis, and the dashed curves show the results from the Power law + Dip + Break analysis. Both models analyze the full black hole and neutron star mass distribution. The teal and gray curves show the analysis results with and without GW230529, respectively.
Download figure:
Standard image High-resolution imageGW230529 is consistent with the population inferred from previously observed CBC candidates. In Figure 6, we show the population distributions of the full black hole and neutron star mass spectrum for the primary component of compact binary mergers using the Binned Gaussian Process (top panel) and Power law + Dip + Break (bottom panel) population models. We qualitatively see that GW230529 is not an outlier with respect to the masses of previously observed compact-object binaries because the inclusion of GW230529 in the population does not significantly alter the full-population posterior constraints. This differs from the detection of GW190814, which was an outlier with respect to the rest of the observed BBH population at the time due to its small secondary mass (Essick et al. 2022). The observation of GW190814 strongly suggested that the region between the masses of neutron stars and black holes was populated (Abbott et al. 2023b), a conclusion that is strengthened with the detection of GW230529 (even though it is not an outlier).
Figure 6. The differential binary merger rate as a function of the mass of the primary component using the Binned Gaussian Process model (top panel) and the Power law + Dip + Break model (bottom panel) for the full compact binary population (solid curves: mean; shaded regions: 90% credible interval) with (teal) and without (gray) GW230529.
Download figure:
Standard image High-resolution imageThe component masses inferred for GW230529 differ across differing population-informed priors. In Figure 7, we show the mass and spin posteriors obtained using priors informed by each of the population models considered in this work. The priors informed by the NSBH-pop model prefer unequal mass ratios and small black hole spins; they suppress the extended posterior tail out to equal masses and antialigned spins. The Binned Gaussian Process model also pulls the posteriors to more asymmetric mass ratios and has less support for antialigned spins. As shown in the top panel of Figure 6, the merger rate density inferred using the Binned Gaussian Process model is nearly flat across the region of parameter space covered by GW230529, meaning that the priors informed by this population model have less of an impact on the shape of the posterior compared to the parameterized models considered. Unlike the NSBH-pop model, the Power law + Dip + Break model has a sharp drop in the merger rate above 3 M⊙. This sharp feature results in the posterior on the primary component being pulled below 3 M⊙ and the posterior on the secondary being pulled to higher masses. The Binned Gaussian Process model has a similar feature, although it is closer to 2 M⊙ and therefore too low to significantly affect the mass estimates for this system. The Power law + Dip + Break model also assumes the same preference for equal-mass systems across the entire mass spectrum and thus, given that the majority of the CBC population is consistent with symmetric mass ratios, infers that the GW230529 binary components are more similar in mass. The combination of the drop in the differential merger rate and the preference for symmetric mass ratios increases the support for more extreme spins oriented in the hemisphere opposite the orbital angular momentum. From these models we find that the binary source of GW230529 either has asymmetric components with low values of χeff or has components that are similar in mass with χeff ∼ −0.3. These results highlight that the choice of prior has a significant impact on the inferred masses and spins of the GW230529 source and hence the inferred nature of the binary components. We assess the nature of the components further in Section 6. More details on the population prior reweighting are given in Appendix H.4.
Figure 7. Posterior distributions of GW230529 under various prior assumptions. The solid teal curve shows the posterior distributions using the default high-spin priors described in Appendix D. The various dashed and dotted curves show posteriors obtained using population priors based on the NSBH-pop model (orange), Binned Gaussian Process (blue), and Power law + Dip + Break (red) models.
Download figure:
Standard image High-resolution imageWe assess whether any of our population models are favored as priors for GW230529 by calculating Bayes factors. Given that the population models consider different sets of GW events, Bayesian evidences for each population model are not directly comparable. However, we may compare the evidence for the single event GW230529 under different models. In this calculation we only consider the effect of the different shapes of the population models, not their overall normalization over the broader mass parameter space considered by each model. Thus, we normalize each of the population models over the same range of component masses as the original parameter estimation priors. We find no significant preference between the population models, with Bayes factors between all three models of .
6. Nature of the Compact Objects in the GW230529 Binary
Without clear evidence for or against tidal effects in the signal, the physical nature of the compact objects in the GW230529 source binary can be assessed by comparing the measured masses and spins of each component with the maximum masses and spins of neutron stars allowed by previous observational data. However, statistical uncertainties in component masses make it especially difficult to determine whether compact objects with masses between ∼2.5 and 5 M⊙ are consistent with being black holes or neutron stars (Hannam et al. 2013; Littenberg et al. 2015). Nevertheless, we assess the nature of the source components by marginalizing over our uncertainties in the masses and spins of GW230529, in the population of merging binaries, and in the supranuclear EOS, to compute the posterior probability that each component had a mass and spin less than the maximum mass and spin supported by the EOS. We follow the procedure introduced in the context of GW190814 (R. Abbott et al. 2020c; Essick & Landry 2020), which relied on only the maximum neutron star mass and has subsequently been extended to include the effects of spin (R. Abbott et al. 2020c; Abbott et al. 2021a, 2023b). Our analysis provides only an upper limit on the probability that an object is a neutron star, as it assumes that all objects consistent with the maximum mass and spin of a neutron star are indeed neutron stars (Essick & Landry 2020).
To assess the nature of the component masses of the GW230529 source, we consider two versions of an astrophysically agnostic population model that is uniform in source-frame component masses and spin magnitudes and isotropic in spin orientations: one with component spin magnitudes χi ≡ ∣ χ i ∣ that are allowed to be large (χ1, χ2 ≤ 0.99), and one where they are restricted a priori to be small (χ1, χ2 ≤ 0.05). We also consider the Power law + Dip + Break population model (Appendix H.2) fit to the GW candidates from GWTC-3 (Abbott et al. 2023a); including GW230529 in the population model does not significantly affect the results. For each population, we marginalize over the uncertainty in the maximum neutron star mass conditioned on the existence of massive pulsars and previous GW observations (Landry et al. 2020) using a flexible Gaussian process (GP) representation of the EOS (Landry & Essick 2019; Essick et al. 2020b). More information on the EOS choices can be found in Appendix I.
Table 3 reports the probability that each component of the merger is consistent with a neutron star. In general, we find that the secondary is almost certainly consistent with a neutron star, and the primary is most probably a black hole. However, when incorporating information from the Power law + Dip + Break population model, we find that there is a ∼1 in 10 chance that the primary is consistent with a neutron star. If we further relax the fixed spin assumptions implicit within the Power law + Dip + Break model for objects with masses ≤2.5 M⊙ from χi ≤ 0.4 to χi ≤ 0.99 (see Appendix H.2), we can find probabilities as high as P(m1 is NS) = 27.3% ± 3.8%. This ambiguity is similar to the secondary component of GW190814 (2.50 ≤ m2/M⊙ ≤ 2.67), which is consistent with a neutron star if it was rapidly spinning (e.g., R. Abbott et al. 2020c; Essick & Landry 2020; Most et al. 2020a).
Table 3. Probabilistic Source Classification Based on Consistency of Component Masses with the Maximum Neutron Star Mass and Spin
| χ1, χ2 ≤ 0.99 | χ1, χ2 ≤ 0.05 | Power law + Dip + Break | |
|---|---|---|---|
| P(m1 is NS) | 2.9% ± 0.4% | <0.1% | 8.8% ± 2.8% |
| P(m2 is NS | 96.1% ± 0.4% | >99.9% | 98.4% ± 1.3% |
Note. All estimates marginalize over uncertainty in the masses, spins, and redshift of the source as well as uncertainty in the astrophysical population and the EOS. We consider three population models: two distributions that use astrophysically agnostic priors and consider either large spins (χ1, χ2 ≤ 0.99) or small spins (χ1, χ2 ≤ 0.05), and a population prior using the Power law + Dip + Break model fit with only the events from GWTC-3 (Abbott et al. 2023b). We use an EOS posterior conditioned on massive pulsars and GW observations (Landry et al. 2020). All errors approximate 90% uncertainty from the finite number of Monte Carlo samples used with the exception of the low-spin results, for which we only place an upper or lower bound.
Download table as: ASCIITypeset image
The differences observed between population models primarily reflect the uncertainty in the mass ratio and spins of the GW230529 source. For example, incorporating the Power law + Dip + Break population model as a prior updates the posterior for m1 from to . Additional observations of compact objects in or near the lower mass gap may clarify the composition of GW230529 by further constraining the exact shape of the distribution of compact objects below ∼5 M⊙ or the supranuclear EOS.
7. Implications for Multimessenger Astrophysics
In NSBH mergers, the neutron star can either plunge directly into the black hole or be tidally disrupted by its gravitational field. Tidal disruption would leave some remnant baryonic material outside the black hole that could potentially power a range of EM counterparts, including a kilonova (Lattimer & Schramm 1974; Li & Paczynski 1998; Tanaka & Hotokezaka 2013; Tanaka et al. 2014; Kawaguchi et al. 2016; Fernández et al. 2017) or a gamma-ray burst (Mochkovitch et al. 1993; Janka et al. 1999; Paschalidis et al. 2015; Shapiro 2017; Ruiz et al. 2018). The conditions for tidal disruption are determined by the mass ratio of the binary, the component of the black hole spin aligned with the orbital angular momentum, and the compactness of the neutron star (Pannarale et al. 2011; Foucart 2012; Foucart et al. 2018; Krüger & Foucart 2020). While the disruption probability of the neutron star in GW230529 can be inferred based on the binary parameters, we are unlikely to directly observe the disruption in the GW signal without next-generation observatories (Clarke et al. 2023).
We use the ensemble of fitting formulae collected in Biscoveanu et al. (2022), including the spin-dependent properties of neutron stars (Cipolletta et al. 2015; Breu & Rezzolla 2016; Foucart et al. 2018; Most et al. 2020b), to constrain the remnant baryon mass outside the final black hole following GW230529, assuming that it was produced by an NSBH merger. We additionally marginalize over the uncertainty in the GP-EOS results obtained using the method introduced in Section 6 (Legred et al. 2021, 2022). Using the high-spin combined posterior samples obtained with default priors, we find a probability of neutron star tidal disruption of 0.1, corresponding to an upper limit on the remnant baryon mass produced in the merger of 0.052 M⊙ at 99% credibility. The low secondary spin priors (χ2 < 0.05) yield a tidal disruption probability and remnant baryon mass upper limit of 0.042 and 0.011 M⊙, respectively. A rapidly spinning neutron star is less compact than a slowly spinning neutron star of the same gravitational mass under the same EOS. This decrease in compactness leads to a larger disruption probability and a larger remnant baryon mass following the merger, explaining the trend we see when comparing the results obtained under the low secondary spin and high-spin priors. The source binary of GW230529 is the most probable of the confident NSBHs reported by the LVK to have undergone tidal disruption because of the increased symmetry in its component masses. However, the exact value of the tidal disruption probability and the remnant baryon mass for this system are prior dependent.
We can also gauge how the inclusion of GW230529 impacts estimates of the fraction of NSBH systems detected in GWs that may be accompanied by an EM counterpart, fEM-bright. Using the mass and spin distributions inferred under the NSBH-pop population model described in Section 5, we find a 90% credible upper limit on the fraction of NSBH mergers that may be EM bright of fEM-bright ≤ 0.18 when including GW230529, an increase relative to fEM-bright ≤ 0.06 obtained when excluding GW230529 from the analysis. When including GP-EOS constraints additionally conditioned on NICER observations (Legred et al. 2021, 2022), the posterior on the EM-bright fraction further increases to peak away from zero, . Additional details on analyses with this alternative choice of EOS constraint are given in Appendix I. These estimates assume that any remnant baryon mass could power a counterpart, although the actual threshold value is astrophysically uncertain. Figure 8 shows the posterior for fEM-bright with and without the inclusion of GW230529. While the exact value of fEM-bright depends on the assumed population model, the increase in fEM-bright for NSBHs upon the inclusion of GW230529 in the population is robust against modeling assumptions.
Figure 8. Posterior on the fraction of NSBH systems detected with GWs that may be EM bright, fEM-bright, depending on the threshold remnant mass required to power a counterpart, . The solid and dashed curves represent different values of the minimum remnant mass . The teal and gray curves show the analysis results with and without GW230529, respectively.
Download figure:
Standard image High-resolution imageUsing these updated multimessenger prospects, we can infer the contribution of NSBH mergers to the production of heavy elements (Biscoveanu et al. 2022) and the generation of gamma-ray bursts (Biscoveanu et al. 2023). Assuming that all remnant baryon mass produced in NSBH mergers is enriched in heavy elements via r-process nucleosynthesis (Lattimer & Schramm 1974, 1976), we infer that NSBH mergers contribute at most 1.1 M⊙ Gpc−3 yr−1 to the production of heavy elements and that the rate of gamma-ray bursts with NSBH progenitors is at most 23 Gpc−3 yr−1 at 90% credibility. This likely represents a small fraction of all short gamma-ray bursts, for which astrophysical beaming-corrected rate estimates are in the range of (Mandel & Broekgaarden 2022), and of the total r-process material produced by compact-object mergers (Côté et al. 2018; Chen et al. 2021b).
No significant counterpart candidates have been reported for GW230529 (IceCube Collaboration 2023; Karambelkar et al. 2023; Lesage et al. 2023; Lipunov et al. 2023; Longo et al. 2023; Savchenko et al. 2023; Sugita et al. 2023a, 2023b; Waratkar et al. 2023). This is unsurprising given that it was only observed by a single detector and hence was poorly localized on the sky.
8. Astrophysical Implications
Since the late 1990s, there have been claims about the existence of a mass gap between the maximum neutron star mass and the minimum black hole mass (∼3–5 M⊙) based on dynamical mass measurements of Galactic X-ray binaries (Bailyn et al. 1998; Ozel et al. 2010; Farr et al. 2011b). The lower edge of this purported mass gap depends on the maximum possible mass with which a neutron star can form in a supernova explosion, which cannot exceed the maximum allowed neutron star mass given by the EOS. Some EOSs support masses up to ∼3 M⊙ for nonrotating neutron stars (Mueller & Serot 1996; Godzieba et al. 2021) and even larger masses for rotating neutron stars (Friedman & Ipser 1987; Cook et al. 1994), although such large values are disfavored by the tidal deformability inferred for GW170817 (Abbott et al. 2019b) and Galactic observations (Alsing et al. 2018; Farr & Chatziioannou 2020). The upper edge and extent of the mass gap depend on the minimum black hole mass that can form from stellar core collapse. However, it remains an open question whether observational or evolutionary selection effects inherent to the detection of Galactic X-ray binaries can lead to the observed gap in compact-object masses (Fryer & Kalogera 2001; Kreidberg et al. 2012; Siegel et al. 2023).
In recent years, new EM observations have unveiled a few candidates in the ∼3 M⊙ region, mostly from noninteracting binary systems (Thompson et al. 2019, 2020b; van den Heuvel & Tauris 2020) and radio surveys for pulsar binary systems (Barr et al. 2024). Microlensing surveys do not support the existence of a mass gap but cannot exclude it either (Wyrzykowski et al. 2016; Wyrzykowski & Mandel 2020). The LVK has already observed one component of a merger whose mass falls between the most massive neutron stars and least massive black holes observed in the Galaxy: the secondary component of GW190814 (2.5–2.7 M⊙ at the 90% credible level; R. Abbott et al. 2020c; Abbott et al. 2024b). The secondary components of the GW200210_092254 and GW190917_114630 source binaries also have support in the lower mass gap (Abbott et al. 2023a, 2024b), although their mass estimates are also consistent with a high-mass neutron star. Unlike GW230529, the primary components of all three of these binaries can be confidently identified as black holes. Overall, the existence of a mass gap between the most massive neutron stars and least massive black holes still stands as an open question in astrophysics.
GW230529 is the first compact binary with a primary component that has a high probability of residing in the lower mass gap. Hence, GW230529 reinforces the conclusion that the 3–5 M⊙ range is not completely empty (see Figure 5 in Section 5). However, the 3–5 M⊙ range may still be less populated than the surrounding regions of the mass spectrum. This conclusion is consistent with previous population analyses and rate estimates based on GWTC-3 (Abbott et al. 2023b).
The formation of GW230529 raises a number of questions. Given our current understanding of core collapse in massive stars (O'Connor & Ott 2011; Janka 2012; Ertl et al. 2020), it is unlikely that the primary component formed via direct collapse because of its low mass, although stochasticity in the physical mechanisms that determine remnant mass may populate the low-mass end of the black hole mass spectrum (Mandel & Müller 2020; Antoniadis et al. 2022). Formation by fallback is a viable scenario: recent numerical models of core-collapse supernovae suggest that the formation of 3–6 M⊙ black holes via substantial fallback is rare but still possible (e.g., Sukhbold et al. 2016; Ertl et al. 2020; Vigna-Gómez et al. 2021). One-dimensional hydrodynamical simulations of core collapse adopting pure helium star models predict that there is no empty gap, only a less populated region between 3 and 5 M⊙, with the lowest mass of black holes produced by prompt implosions starting at ≈6 M⊙ (Ertl et al. 2020). Another relevant parameter is the timescale for instability growth and launch of the core-collapse supernova; if this is long enough (≳200 ms), the proto−neutron star might accrete enough mass before the explosion to become a mass-gap object (Fryer & Kalogera 2001; Fryer et al. 2012). Compact binary population synthesis models that allow for longer instability growth timescales naturally produce merging NSBH systems with the primary component in the lower mass gap (Belczynski et al. 2012; Zevin et al. 2020; Broekgaarden et al. 2021; Chattopadhyay et al. 2021; Olejak et al. 2022). It may also be possible to form mass-gap objects through accretion onto a neutron star. If the first-born neutron star in the binary accretes enough material prior to the formation of the second-born compact object, it can trigger an accretion-induced collapse into a black hole, yielding a mass-gap object (Siegel et al. 2023). This scenario may be aided by super-Eddington accretion onto the first-born neutron star, which has been proposed to explain NSBHs with mass-gap black holes like the source of GW230529 (Zhu et al. 2024). Considering the large number of uncertainties about the outcome of core-collapse supernovae (e.g., Burrows & Vartanyan 2021), the primary mass of GW230529 provides a piece of evidence to inform and constrain future models. Overall, astrophysical models in the past have preferentially adopted prescriptions that enforce the presence of a lower mass gap (e.g., Fryer et al. 2012), and the inferred rate of GW230529-like systems urges a change of paradigm in such model assumptions.
Alternatively, the primary component of GW230529 might be the result of the merger of two neutron stars. For instance, the 90% credible intervals for the remnant mass of GW190425 and the primary mass of GW230529 overlap (Abbott et al. 2020a). This may hint at a scenario where the primary component could be either the member of a former triple or quadruple system (Fragione et al. 2020; Lu et al. 2021; Vynatheya & Hamers 2022; Gayathri et al. 2023) or the result of a dynamical capture in a star cluster (Clausen et al. 2013; Arca Sedda 2020, 2021; Gupta et al. 2020; Rastello et al. 2020) or an active galactic nucleus disk (Yang et al. 2020; Tagawa et al. 2021). This scenario was proposed for the formation of a compact object discovered in a binary with a pulsar in the globular cluster NGC 1851, which is measured to have a mass of 2.09–2.71 M⊙ at 95% confidence (Barr et al. 2024). However, dynamically induced BNS mergers are predicted to be rare in dense stellar environments ( e.g., Ye et al. 2020; Samsing & Hotokezaka 2021). Thus, the rate of mergers between a BNS merger remnant and a neutron star may be at least several orders of magnitude lower than the rate we infer for GW230529-like mergers, making this scenario improbable.
Non-stellar-origin black hole formation scenarios such as primordial black holes (e.g., Clesse & Garcia-Bellido 2022) remain a possibility. However, there are significant uncertainties in the predicted mass spectrum and merger rate of primordial black hole binaries, thus making it difficult to attribute a primordial origin to compact objects that have masses consistent with predictions from massive-star core collapse. Furthermore, results from microlensing surveys indicate that primordial black hole mergers cannot be a dominant source of GWs in the local universe (Mroz et al. 2024).
It has also been suggested that mergers apparently involving mass-gap objects could instead be gravitationally lensed BNSs (Bianconi et al. 2023; Magare et al. 2023), with the lensing magnification making them appear heavier and closer (Wang et al. 1996; Dai et al. 2017; Hannuksela et al. 2019). This scenario is difficult to explicitly test in the absence of tidal information (Pang et al. 2020) or EM counterparts (Bianconi et al. 2023), but the expected relative rate of strong lensing is low at current detector sensitivities (Abbott et al. 2023c; Magare et al. 2023; Smith et al. 2023).
Finally, we also find mild support for the possibility that the primary component is a neutron star rather than a black hole when considering the population-based Power law + Dip + Break prior that incorporates the potential presence of a gap at low black hole masses. In this case, GW230529 would be the most massive neutron star binary yet observed, with both components ≳2.0 M⊙, and have nonzero spins that are antialigned with the orbital angular momentum. The effective inspiral spin of the BNSs in this scenario would differ significantly from BNS sources previously observed in GWs (Abbott et al. 2017a, 2020a), as well as those inferred for Galactic BNSs if they were to merge (Zhu et al. 2018). Spins oriented in the hemisphere opposite the binary orbital angular momentum could be the result of supernova natal kicks (Kalogera 2000; Farr et al. 2011a; O'Shaughnessy et al. 2017; Chan et al. 2020) or spin-axis tossing (Tauris 2022) at birth. For example, one of the neutron stars in the binary pulsar system J0737–3039 is significantly misaligned with respect to both the spin of its companion and the orbital angular momentum of the binary system, which may require off-axis kicks (Farr et al. 2011a). Alternatively, isotropically oriented component spins that result from random pairing in dynamical environments could lead to the significant spin misalignment (e.g., Rodriguez et al. 2016).
9. Summary
GW230529 is a GW signal from the coalescence of a 2.5–4.5 M⊙ compact object and a compact object consistent with neutron star masses. The more massive component in the merger provides evidence that compact objects in the hypothesized lower mass gap exist in merging binaries. Based on mass estimates of the two components in the merger and current constraints on the supranuclear EOS, we find the most probable interpretation of the GW230529 source to be an NSBH coalescence. In this scenario, the source binary of GW230529 is the most symmetric-mass NSBH merger yet observed, and the primary component of the merger is the lowest-mass primary black hole observed in GWs to date. Because NSBHs with more symmetric masses are more susceptible to tidal disruption, the observation of GW230529 implies that more NSBHs than previously inferred may produce EM counterparts. However, we cannot rule out the contrasting scenario that the source of GW230529 consisted of two neutron stars rather than a neutron star and a black hole. In this case, the source of GW230529 would be the only BNS coalescence observed to have strong support for nonzero and antialigned spins, as well as the highest-mass BNS system observed to date. Regardless of the true nature of the GW230529 source, it is a novel addition to the growing population of CBCs observed via their GW emission and highlights the importance of continued exploration of the CBC parameter space in O4 and beyond.
Strain data containing the signal from the LIGO Livingston observatory are available from the Gravitational Wave Open Science Center. 311 Specifically, we release the L1:GDS-CALIB_STRAIN_CLEAN_AR channel, where the AR designation means that the strain data are analysis ready; this strain channel was also released at the end of O3. Samples from the posterior distributions of the source parameters, hyperposterior distributions from population analyses, and notebooks for reproducing all results and figures in this Letter are available on Zenodo (Abbott et al. 2024). The software packages used in our analyses are open-source.
Acknowledgments
This material is based on work supported by NSF's LIGO Laboratory, which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO 600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS), and the Netherlands Organization for Scientific Research (NWO) for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies, as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación (AEI), the Spanish Ministerio de Ciencia, Innovación y Universidades, the European Union NextGenerationEU/PRTR (PRTR-C17.I1), the ICSC - CentroNazionale di Ricerca in High Performance Computing, Big Dataand Quantum Computing, funded by the European Union NextGenerationEU, the Comunitat Autonòma de les Illes Balears through the Direcció General de Recerca, Innovació i Transformació Digital with funds from the Tourist Stay Tax Law ITS 2017-006, the Conselleria d'Economia, Hisenda i Innovació the FEDER Operational Program 2021–2027 of the Balearic Islands, the Conselleria d'Innovació Universitats, Ciència i Societat Digital de la Generalitat Valenciana and the CERCA Programme Generalitat de Catalunya, Spain, the National Science Centre of Poland and the European Union—European Regional Development Fund; Foundation for Polish Science (FNP), the Polish Ministry of Science and Higher Education, the Swiss National Science Foundation (SNSF), the Russian Science Foundation, the European Commission, the European Social Funds (ESF), the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the French Lyon Institute of Origins (LIO), the Belgian Fonds de la Recherche Scientifique (FRS-FNRS), Actions de Recherche Concertées (ARC) and Fonds Wetenschappelijk Onderzoek—Vlaanderen (FWO), Belgium, the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, the Natural Science and Engineering Research Council Canada, Canadian Foundation for Innovation (CFI), the Brazilian Ministry of Science, Technology, and Innovations, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the National Science and Technology Council (NSTC), Taiwan, the United States Department of Energy, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN, and CNRS for provision of computational resources. This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133 and 20H05639, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, National Research Foundation (NRF), Computing Infrastructure Project of Global Science experimental Data hub Center (GSDC) at KISTI, Korea Astronomy and Space Science Institute (KASI), and Ministry of Science and ICT (MSIT) in Korea, Academia Sinica (AS), AS Grid Center (ASGC) and the National Science and Technology Council (NSTC) in Taiwan under grants including the Rising Star Program and Science Vanguard Research Program, Advanced Technology Center (ATC) of NAOJ, and Mechanical Engineering Center of KEK. We thank the anonymous journal referee for helpful comments.
Software: Calibration of the LIGO strain data was performed with a GstLAL-based calibration software pipeline (Viets et al. 2018). Data-quality products and event-validation results were computed using the DMT (Zweizig 2006), DQR (LIGO Scientific Collaboration & Virgo Collaboration 2018), DQSEGDB (Fisher et al. 2021), gwdetchar (Urban et al. 2021), hveto (Smith et al. 2011), iDQ (Essick et al. 2020a), Omicron (Robinet et al. 2020), and PythonVirgoTools (Virgo Collaboration 2021) software packages and contributing software tools. Analyses in this catalog relied on software from the LVK Algorithm Library Suite (LIGO Scientific, Virgo, & KAGRA Collaboration 2018). The detection of the signals and subsequent significance evaluations were performed with the GstLAL-based inspiral software pipeline (Messick et al. 2017; Sachdev et al. 2019; Hanna et al. 2020; Cannon et al. 2021), with the MBTA pipeline (Adams et al. 2016; Aubin et al. 2021), and with the PyCBC (Usman et al. 2016; Nitz et al. 2017; Davies et al. 2020) packages. Estimates of the noise spectra and glitch models were obtained using BayesWave (Cornish & Littenberg 2015; Littenberg et al. 2016; Cornish et al. 2021). Low-latency source localization was performed using BAYESTAR (Singer & Price 2016). Source-parameter estimation was primarily performed with the Bilby and Parallel Bilby libraries (Ashton et al. 2019; Romero-Shaw et al. 2020; Smith et al. 2020) using the Dynesty nested sampling package (Speagle 2020). SEOBNRv5PHM waveforms used in parameter estimation were generated using pySEOBNR (Mihaylov et al. 2023). FTI and TIGER waveforms used for testing general relativity were generated using Bilby TGR (Ashton et al. 2024). PESummary was used to post-process and collate parameter estimation results (Hoy & Raymond 2021). The various stages of the parameter estimation analysis were managed with the Asimov library (Williams et al. 2023). Plots were prepared with Matplotlib (Hunter 2007), seaborn (Waskom 2021) and GWpy (Macleod et al. 2021). NumPy (Harris et al. 2020) and SciPy (Virtanen et al. 2020) were used for analyses in the manuscript.
Appendix A: Upgrades to the Detector Network for O4
The Advanced LIGO, Advanced Virgo, and KAGRA detectors have all undergone a series of upgrades to improve the network sensitivity in preparation for O4. In this appendix we focus on upgrades made to the LIGO Livingston observatory, as it was the only detector to observe GW230529; similar upgrades were made at LIGO Hanford. Some of these improvements are a part of the Advanced LIGO Plus (A+) detector upgrades (Abbott et al. 2020b).
The principal commissioning work done in preparation for O4 was to enhance the optics systems in the detector. The pre-stabilized laser was redesigned for O4 with amplifiers allowing input power up to 125 W (Bode et al. 2020; Cahillane & Mansell 2022). LIGO Livingston operated at 63 W in O4a. Both end test masses at LIGO Livingston were replaced because their mirror coatings contained small, pointlike absorbers that limited detector sensitivity by degrading the power-recycling gain (Brooks et al. 2021). For O4, frequency-dependent squeezing was implemented with the addition of a 300 m filter cavity to rotate the squeezed vacuum quadrature across the bandwidth of the detector (McCuller et al. 2020; Dwyer et al. 2022). With frequency-dependent squeezing, uncertainty is reduced in both the amplitude and phase quadratures, allowing for a broadband reduction in quantum noise (Ganapathy et al. 2023). Squeezing was implemented independently of frequency in O3, which only reduced the quantum noise at high frequencies (Tse et al. 2019).
Other commissioning work was done to improve data quality and reduce noise at low frequencies. For O4, scattered-light noise was mitigated by removing some problematic scattering surfaces and improving the resonant damping of other scatterers (Soni et al. 2020; Davis et al. 2021). Low-frequency technical noise (i.e., noise that is not intrinsic to the detector design but can limit the detector sensitivity and performance) was reduced with the commissioning of feedback control loops, noise subtraction, and better electronics. Overall, the upgrades made during the commissioning period for O4 led to a broadband reduction in noise and improvement in detector sensitivity and performance.
Appendix B: Detection Time Line and Circulars
The LVK issues low-latency alerts to facilitate prompt follow-up of GW candidates (Chaudhary et al. 2023). GW230529 was initially identified in a low-latency search using data from the LIGO Livingston observatory. The candidate was given the name S230529ay in the Gravitational-Wave Candidate Event Database (GraceDB). 312
After the detection of GW230529, an initial notice was sent out to astronomers through NASA's General Coordinates Network (GCN; LIGO Scientific, Virgo, & KAGRA Collaboration 2023a). The sky localization computed in low latency by BAYESTAR (Singer & Price 2016; Singer et al. 2016) had a 90% credible area of ≈24,200 deg2; the large credible area was due to GW230529 only being observed by a single detector. The initial alert also included low-latency estimates of the mass-based source classification (BNS, BBH, and NSBH) produced by the PyCBC pipeline (Villa-Ortega et al. 2022). Additional estimates that the source binary of GW230529 contained at least one neutron star (> 99.9%), contained at least one compact object in the lower mass gap 3–5 M⊙ (98.5%), and had neutron star matter ejected outside the final compact object (12.1%) were produced by a machine-learning method trained to infer source properties from the search pipeline results (Chatterjee et al. 2020).
Low-latency parameter estimation was performed with Bilby (Ashton et al. 2019; Romero-Shaw et al. 2020) and used to produce an updated sky localization and estimates of source properties sent out in another alert (LIGO Scientific, Virgo, & KAGRA Collaboration 2023c). The updated sky localization had a 90% credible area of ≈25,600 deg2. The updated estimate of source properties led to decreases in the probabilities that the source binary included at least one neutron star (98.5%), one compact object in the mass gap (69.4%), and ejected matter outside the remnant object (1.1%).
Appendix C: Quantifying Significance for a Single-detector Candidate Event
Each search pipeline has its own method for calculating the FAR of a single-detector GW candidate. In this appendix we will summarize these methods.
C.1. GstLAL
The GstLAL pipeline assigns a likelihood ratio to all of its candidate signals as the ranking statistic to estimate their significance. The likelihood ratio is the ratio of the probability of obtaining the candidate under the signal hypothesis to the probability of obtaining the candidate under the noise hypothesis. The latter probability is largely calculated by collecting S/N and ξ2 statistics of noise triggers during the analysis, where ξ2 is a signal consistency test statistic (Messick et al. 2017). Similar templates in the template bank are grouped together on the basis of the post-Newtonian (PN) expansion of their waveforms (Morisaki & Raymond 2020; Sakon et al. 2024). Each template bin collects its own S/N and ξ2 statistics, which get used to rank candidates recovered by templates in that bin. Further details about the calculation of the probability of obtaining the candidate under the noise hypothesis, as well as the GstLAL likelihood ratio in general, can be found in Tsukada et al. (2023). The S/N–ξ2 statistics collected from noise triggers from the same template bin as GW230529 with GW230529 superimposed are shown in Figure 9. GW230529 clearly stands out from the background, and there are no noise triggers at its position in S/N–ξ2 space.
Figure 9. The noise background for the LIGO Livingston observatory collected by the GstLAL pipeline in offline mode, with the S/N–ξ2 values of GW230529 overlaid. The test statistic ξ2 measures signal consistency. The background is collected from templates that are part of the same bin as the template that recovered GW230529 and is consequently the background used to rank GW230529 and calculate its significance. The background distribution has effectively no support at the position of GW230529. The color bar represents the probability of producing a certain S/N and ξ2 from noise triggers.
Download figure:
Standard image High-resolution imageSince nonstationary noise transients known as glitches (Nuttall 2018) are not expected to be correlated across detectors, single-detector GW candidates, whether during times when a single or multiple detectors are observing, are more likely to be glitches than coincident GW candidates. To account for this, the GstLAL pipeline downweights the logarithm of the likelihood ratio of single-detector candidates by subtracting an empirically tuned factor of 13, which is optimized to maximize the recovery of candidates with an astrophysical origin while minimizing the recovery of glitches (Abbott et al. 2020a). The FAR of the GstLAL pipeline for GW230529 quoted in Table 1 is calculated after the application of this penalty.
Finally, the FAR for a candidate in the search is computed by comparing its likelihood ratio with the likelihood ratios of noise triggers not found in coincidence, after accounting for the live time of the analysis. These noise triggers not found in coincidence allow the pipeline to extrapolate the background distribution of likelihood ratios to large values, enabling the inverse FAR of a candidate to be greater than the duration of the analysis. While in the low-latency mode, the GstLAL FAR calculation uses the background from the start of the analysis period until the time of the candidate. In contrast, for its offline analysis, the GstLAL pipeline uses background collected from the full analysis period, including the entirety of O4a, to rank candidates. This differs from the other search pipelines presented in this work, which only use background from the first 2 weeks of O4a for their offline search results presented here. The GstLAL offline analysis was performed using the same template bank as the online analysis (Ewing et al. 2024; Sakon et al. 2024). Details are liable to change for future offline GstLAL analyses in O4a. However, for such a significant candidate as GW230529, we do not expect these changes to impact its interpretation as highly likely being of astrophysical origin.
C.2. MBTA
The MBTA single-detector candidate search for O4a focuses on a subset of the full MBTA parameter space. The goal is to focus on BNS and NSBH signals, which are most interesting because of their rarity and possible EM counterparts. Their long signal duration also makes them easier to identify and allows for a better sensitivity to be reached when compared to a single-detector candidate search using the whole MBTA parameter space. The parameter space for the O4a online single-detector search was defined based on the detector-frame (redshifted) primary mass and chirp mass of the templates and motivated by the computation of the probability for whether a nonzero amount of neutron star material remained outside the final remnant compact object (Chatterjee et al. 2020) introduced in Appendix B. The considered space is (1 + z)m1 < 50 M⊙ and (Juste 2023). For the offline analysis, this space was adapted and we consider only . MBTA online also applied selection criteria based on data-quality tests to its single-detector candidates (Juste 2023). This was changed to a reweighting of the ranking statistic for the offline analysis.
The ranking statistic for MBTA single-detector triggers is a reweighted S/N. The reweighting is based on the computation of a quantity called auto χ2 (Aubin et al. 2021), which tests the consistency of the time evolution of the S/N time series. The additional reweighting used in the offline analysis relies on the computation of a quantity that identifies an excess of S/N for single-detector triggers. The excess of S/N is larger for triggers that have a noise origin compared to astrophysical signals or injections and therefore allows for discrimination between astrophysical signals and noise. The ranking statistic distribution of the MBTA offline analysis for single-detector triggers during the first 2 weeks of O4a is shown in Figure 10. Other triggers that produced significant public alerts have been removed from the plotted distribution. GW230529 stands out from the background with a high ranking statistic that reflects the auto χ2 value being completely consistent with a signal origin.
Figure 10. Ranking statistic distribution of MBTA single-detector triggers in the LIGO Livingston observatory during the first 2 weeks of O4, excluding significant public alerts aside from GW230529.
Download figure:
Standard image High-resolution imageThe FAR in MBTA is a function of pastro (Andres et al. 2022), the probability of astrophysical origin of the candidate. It is derived from the combined parameterizations of the FAR as a function of the ranking statistic and of pastro as a function of the ranking statistic. Inverting the latter gives the ranking statistic as a function of pastro and eventually the FAR as a function of pastro. The background estimation for MBTA single-detector triggers was computed differently during the online and offline analyses. The online analysis relies on the computation of simulated single-detector triggers obtained through random combinations of single frequency band data (Juste 2023). O4a online analysis and the method of randomly combining single frequency band data showed that the background for MBTA single-detector triggers follows an exponential distribution and is stable in time beyond statistical fluctuations. This prompted us to update the model we use for the offline (and future) analyses to reach greater sensitivity. This new model involves extrapolating the observed distribution of single-detector triggers and removing the significant triggers. A safety margin is used in the extrapolation such that the FAR is overestimated relative to the best-fit extrapolated distribution. This change in methods, in addition to the difference in handling of single-detector triggers online and offline, explains the difference in inverse FAR for the online (1.1 yr) and offline (>1000 yr) analyses.
C.3. PyCBC
In PyCBC, each GW candidate is assigned a ranking statistic, and then a FAR is calculated by comparing the ranking statistic of the candidate with the background distribution. The details of this procedure differ between the online and offline versions of the PyCBC search.
In the online pipeline, we consider single-detector candidates only from templates with duration greater than 7 s above a starting frequency of 17 Hz, corresponding to a range of masses and spins for which EM emission due to neutron star ejecta might be expected. Low-latency data-quality time series produced by iDQ, a machine-learning framework for autonomous detection of noise artifacts using only auxiliary data channels insensitive to GWs (Essick et al. 2020a), are used to veto candidates at times when a glitch is likely to be present in the data. The remaining single-detector candidates are ranked by the reweighted S/N (Nitz 2018). Only candidates with a χ2 statistic (Allen 2005) <2.0 and reweighted S/N > 6.75 are kept. The rate density of single-detector triggers above the reweighted S/N threshold is fit with a decreasing exponential.
In offline PyCBC, we only consider single-detector candidates from templates with duration greater than 0.3 s above a starting frequency of 15 Hz. We also require that candidates have a χ2 statistic <10.0, a reweighted S/N > 5.5, and a PSD variation statistic (Mozzon et al. 2020) <10.0, in order to exclude high-amplitude noise transients. Each candidate is assigned a ranking statistic equal to the logarithm of the ratio of astrophysical signal likelihood to detector noise likelihood. The PyCBC O4 offline search introduced two changes to the ranking statistic relative to the O3 calculation. First, an explicit model of the signal distribution covering the space of binary masses and spins is included (Kumar & Dent 2024). The model is designed to maximize the number of detected signals from the known compact binary distribution, while also maintaining sensitivity to signals in previously unpopulated regions. Second, the model of the rate density of events caused by detector noise now includes a term describing the variation in rate during times of heightened detector noise (Davis et al. 2022). The rate density of single-detector candidates above a ranking statistic threshold of 0 is fit with a decreasing exponential.
In both the online and offline versions of PyCBC, the exponential fit of trigger rate density above the relevant ranking statistic threshold is used to extrapolate the FAR of single-detector candidates beyond the observing time of the search (Cabourn Davies & Harry 2022). The FAR calculations for GW230529 used triggers from the first 2 weeks of O4 at the LIGO Livingston observatory. The distribution of ranking statistics for offline single-detector candidates at the LIGO Livingston observatory is shown in Figure 11. Similar to the other searches, GW230529 clearly stands out from the background distribution.
Figure 11. Ranking statistic distribution of PyCBC offline single-detector triggers in the LIGO Livingston observatory during the first 2 weeks of O4. The plotted distribution may include triggers from other less significant signals arriving during the period analyzed.
Download figure:
Standard image High-resolution imageAppendix D: Priors, Waveform Systematics, and Bayes Factors
Given the uncertain nature of the compact objects, we analyze GW230529 with a suite of models that incorporate a number of key physical effects. The choices of data duration (128 s) and frequency bandwidth (20–1792 Hz) analyzed were informed by comparing waveforms spanning the mass range recovered in preliminary analyses with the detector noise PSD at the time of the event. The high-frequency cutoff is chosen to avoid loss of power at high frequencies due to low-pass filtering of the data. For a subset of the signal models below, we employ a range of techniques to speed up the evaluation of the likelihood, including heterodyning (also known as relative binning; Cornish 2010, 2021; Zackay et al. 2018; Krishna et al. 2023), multibanding (García-Quirós et al. 2021; Morisaki 2021), and reduced-order quadratures (Canizares et al. 2015; Smith et al. 2016; Morisaki et al. 2023).
The physical effects included and the spin prior ranges for each waveform model considered in this work are shown in Table 4. All analyses use mass priors that are flat in the redshifted component masses with chirp masses and mass ratios q ∈ [0.125, 1]. The luminosity distance prior is uniform in comoving volume and source-frame time in the range DL ∈ [1, 500] Mpc. The priors on the tidal deformability parameters are chosen to be uniform in the component tidal deformabilities over Λ ∈ [0, 5000]. Standard priors (e.g., Romero-Shaw et al. 2020; Abbott et al. 2024b) are used for all the other extrinsic binary parameters. Below we provide further motivation for considering this suite of waveform models.
![${ \mathcal M }\,={({m}_{1}{m}_{2})}^{3/5}/{({m}_{1}+{m}_{2})}^{1/5}\in [2.0214,2.0331]\,{M}_{\odot }$](https://content.cld.iop.org/journals/2041-8205/970/2/L34/revision2/apjlad5bebieqn50.gif)
Table 4. Summary of Parameter Estimation Analysis Choices for GW230529, and the Physical Content of Each Waveform Model Used
| Waveform Model | Precession | Higher Multipoles | Tides | Disruption | Spin Prior |
|---|---|---|---|---|---|
| IMRPhenomNSBH | ⋯ | ⋯ | ✓ | ✓ | χ1 < 0.50, χ2 < 0.05 |
| IMRPhenomPv2_NRTidalv2 | ✓ | ⋯ | ✓ | ⋯ | χ1 < 0.99, χ2 < 0.05 |
| IMRPhenomXPHM | ✓ | ✓ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.99 |
| SEOBNRv5PHM | ✓ | ✓ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.99 |
| SEOBNRv4_ROM_NRTidalv2_NSBH | ⋯ | ⋯ | ✓ | ✓ | χ1 < 0.90, χ2 < 0.05 |
| IMRPhenomXPHM | ✓ | ✓ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.05 |
| IMRPhenomXP | ✓ | ⋯ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.99 |
| IMRPhenomXHM | ⋯ | ✓ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.99 |
| IMRPhenomXAS | ⋯ | ⋯ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.99 |
| IMRPhenomXAS | ⋯ | ⋯ | ⋯ | ⋯ | χ1 < 0.50, χ2 < 0.05 |
| IMRPhenomPv2_NRTidalv2 | ✓ | ⋯ | ✓ | ⋯ | χ1 < 0.05, χ2 < 0.05 |
| IMRPhenomXPHM | ✓ | ✓ | ⋯ | ⋯ | χ1 < 0.05, χ2 < 0.05 |
| SEOBNRv5PHM | ✓ | ✓ | ⋯ | ⋯ | χ1 < 0.99, χ2 < 0.05 |
Note. Here tides refers to modeling of neutron star tidal deformability and disruption when tidal forces overcome the self-gravity of the neutron star. The spin prior denotes any restrictions on the spin magnitude of the binary components.
Download table as: ASCIITypeset image
The primary analysis in this work uses the IMRPhenomXPHM and SEOBNRv5PHM BBH waveform models (corresponding to the third and fourth rows in Table 4), which provide an accurate description of BBHs but do not incorporate tidal effects or the potential tidal disruption of a companion neutron star. We did not consider an extension of IMRPhenomXPHM that incorporates additional physics beyond the two precessing higher-order multipole moment BBH models used here, as these effects only enter at high frequencies and are not relevant for this event (Thompson et al. 2024). In order to quantify the impact of neglecting tidal physics, we analyze GW230529 with NSBH and BNS waveform models that incorporate different tidal information. The NSBH models only incorporate tidal information from the less massive compact object, are restricted to spins aligned with the orbital angular momentum, and only model the dominant ℓ = m = 2 multipole moment. However, they do model the possible tidal disruption of the neutron star, which occurs when tidal forces of the black hole dominate over the self-gravity of the neutron star.
We use two NSBH models, a frequency-domain phenomenological model IMRPhenomNSBH (Thompson et al. 2020a) and a frequency-domain EOB surrogate model SEOBNRv4_ROM_NRTidalv2_NSBH (Matas et al. 2020). IMRPhenomNSBH uses the BBH models IMRPhenomC (Santamaria et al. 2010) and IMRPhenomD (Husa et al. 2016; Khan et al. 2016) as baselines for the amplitude and phase, respectively, incorporating corrections to the phase due to tidal effects following the NRTidal model (Dietrich et al. 2017, 2019a, 2019b) and to the amplitude following Pannarale et al. (2013, 2015). SEOBNRv4_ROM_NRTidalv2_NSBH uses the SEOBNRv4 BBH model as a baseline (Bohé et al. 2017) and applies the same corrections to account for tidal deformability and disruption as IMRPhenomNSBH but is additionally calibrated against numerical relativity (NR) simulations of NSBH mergers (Kyutoku et al. 2010, 2011; Foucart et al. 2013, 2014, 2019). As the NSBH waveform models do not capture spin-induced orbital precession, we analyze GW230529 with the aligned-spin BBH waveform models IMRPhenomXAS (Pratten et al. 2020), which only contains the dominant harmonic, and IMRPhenomXHM (García-Quirós et al. 2020), which includes higher-order multipole moments. Models for NSBH binaries that contain both higher-order multipole moments and precession are under active development (Gonzalez et al. 2023) and could potentially allow for a more consistent model selection between the different source categories.
Given that the mass of the primary overlaps with current estimates of the maximum known neutron star mass (Landry & Read 2021; Romani et al. 2022), we also analyze GW230529 with a BNS waveform model that allows for tidal interactions on both the primary and the secondary components of the binary. We use IMRPhenomPv2_NRTidalv2 (Dietrich et al. 2019a), which adds a model for the tidal phase that is calibrated against both EOB and NR simulations to the underlying precessing-spin point-particle waveform model IMRPhenomPv2 (Hannam et al. 2014). A limitation is that IMRPhenomPv2_NRTidalv2 is calibrated against a suite of equal-mass, nonspinning BNS NR simulations, where the maximum neutron star mass is only 1.372 M⊙. Nevertheless, the model has been validated against a catalog of EOB–NR hybrid waveforms that includes asymmetric configurations and heavier neutron star masses (Dietrich et al. 2019a; Abac et al. 2024). Because of large uncertainties in BNS postmerger waveform modeling, IMRPhenomPv2_NRTidalv2 includes an amplitude taper from 1 to 1.2 times the estimated merger frequency (Dietrich et al. 2019b). The resulting suppression of S/N at these frequencies may be responsible for the posterior differences between IMRPhenomPv2_NRTidalv2 and the BBH and NSBH waveform models we consider (Figure 12).
Figure 12. The two-dimensional q–χeff posterior probability distributions for GW230529 using various BBH, NSBH, and BNS signal models. The vertical dotted lines indicate primary masses that have been mapped to the mass ratio given the median chirp mass estimated from the IMRPhenomXPHM high-spin analysis.
Download figure:
Standard image High-resolution imageFrom the mass estimates alone, the secondary object is consistent with being a neutron star. As such, we follow earlier analyses and analyze GW230529 with two spin priors (Abbott et al. 2021a): an agnostic high-spin prior and an astrophysically motivated low-spin prior. The high-spin prior assumes that the spins on both objects are isotropically oriented and have dimensionless spin magnitudes that are uniformly distributed up to χ1, χ2 ≤ 0.99. The low-spin prior, on the other hand, is inspired by the extrapolated maximum spin observed in Galactic BNSs that will merge within a Hubble time (Burgay et al. 2003; Stovall et al. 2018) and restricts the spin magnitude of the secondary to be χ2 ≤ 0.05 while keeping χ1 ≤ 0.99. As the nature of the primary as a black hole or neutron star is uncertain, we also use a third choice of prior where both spins are restricted to χ1, χ2 ≤ 0.05 for the IMRPhenomPv2_NRTidalv2 and IMRPhenomXPHM waveform models (rows 11 and 12 of Table 4, respectively). The purpose of using multiple priors is to help gauge whether astrophysical assumptions on the neutron star spin impact any of the statements on the probability that the primary object lies in the lower mass gap. In Figure 12, we show the impact of waveform systematics and prior assumptions on the strongly correlated mass ratio effective inspiral spin distribution.
The degree to which a given waveform model under certain assumptions matches the data can be gauged using the Bayes factor. However, Bayes factors penalize extra degrees of freedom in models that do not improve the fit when those extra degrees of freedom are constrained by the data (Mackay 2003). We find no significant preference () between the high-spin and low-spin prior on the secondary; we similarly do not find a statistical preference for or against the effects of precession, higher-order multipole moments, or tidal deformation or disruption of the secondary object.
Appendix E: Additional Source Properties
E.1. Component Spins and Precession
Assuming a high-spin prior, the posteriors for the primary spin magnitude are only weakly informative, disfavoring zero and extremal spins with χ1 = 0.07–0.85. For the secondary, the posteriors are even less informative. Under the low-spin prior, the primary spin magnitude posterior is peaked at slightly higher values, with zero and extremal spins being disfavored to a larger degree, χ1 = 0.10–0.84. The secondary spin is completely uninformative over the restricted range of spin magnitudes. The joint two-dimensional posterior probability distribution of the dimensionless spin magnitude and the spin tilt are shown in Figure 13. Regions of high (low) probability are denoted by a darker (lighter) shade.
Figure 13. Two-dimensional posterior probability distributions for the spin magnitude χi and the spin-tilt angle θi for the primary (left) and secondary (right) compact objects at a reference frequency of 20 Hz. A spin-tilt angle of 0° (180°) indicates a spin that is perfectly aligned (antialigned) with the orbital angular momentum . The pixels have equal prior probability, and shading denotes the posterior probability of each pixel of the high-spin prior analysis, after marginalizing over azimuthal angles. The solid (dashed) contours denote the 90% credible region for the high-spin (low-spin) prior analyses. The probability distributions are marginalized over the azimuthal spin angles.
Download figure:
Standard image High-resolution imageWhen the spins are misaligned with the orbital angular momentum, relativistic spin–orbit and spin–spin couplings drive the evolution of the orbital plane and the spins themselves (Apostolatos et al. 1994; Kidder 1995). The leading-order effect can be captured by an effective precession spin parameter, 0 ≤ χp ≤ 1, which approximately measures the degree of in-plane spin and can be used to parameterize the rate of precession of the orbital plane (Schmidt et al. 2015)
We find the constraints on χp to be uninformative, and we are not able to make any significant statements on precession. The uninformative nature of these results is corroborated by the Bayes factor between the precessing and nonprecessing phenomenological waveform models, .
E.2. Source Location and Distance
As GW230529 was a single-detector observation, the sky localization is poor and covers a sky area of ≈24,100 deg2 at the 90% credible level. The luminosity distance is inferred to be , corresponding to a redshift of computed using the Planck 2015 cosmological parameters (Ade et al. 2016). The luminosity distance has a degeneracy with the inclination angle of the binary's total angular momentum with respect to the line of sight, θJN (e.g., Cutler & Flanagan 1994; Nissanke et al. 2010; Aasi et al. 2013; Vitale & Chen 2018). The posterior distribution on θJN is broadly unconstrained, showing no preference for a total angular momentum vector that is pointed toward or away from the line of sight.
E.3. Tidal Deformability
The tidal deformability of a neutron star is defined by
where k2 is the gravitational Love number of the object and R is its radius (Damour et al. 1992; Mora & Will 2004; Flanagan & Hinderer 2008; Damour & Nagar 2009). However, it is often convenient to introduce a dimensionless tidal deformability
where C = m/R is the compactness of the neutron star in geometrized units.
Using the IMRPhenomPv2_NRTidalv2 model, which allows for tidal deformability on both objects but does not model tidal disruption, we infer Λ1 ≤ 1462 at 90% credibility with the χ1 < 0.99, χ2 < 0.05 spin prior; the posterior peaks at Λ1 = 0. We do not constrain Λ2 relative to the prior. Constraints on the tidal deformability Λ of both components of GW230529 are shown in Figure 14. We also show the tidal deformability predicted by the Λ(m) relation of the specific neutron star EOS model BSK24 (Goriely et al. 2013; Pearson et al. 2018) weighted by the m2 posterior distribution. BSK24 is chosen as an illustrative EOS that is thermodynamically consistent and falls within the range of support of constraints from nuclear physics and astrophysics, including GW170817 (Perot et al. 2019). Even under the assumption that the EOS is known (BSK24), Λ2 is not well constrained, due to the relatively wide m2 posterior. The Λ2 prediction of the set of GP-EOS constraints (red; Section 6) is consistent with both the inferred value of Λ2 (green solid) and the BSK24 value of Λ2 (navy). In addition to the tidal deformability posteriors, we use two distinct likelihood interpolation schemes (Landry & Essick 2019; Ray et al. 2023a) to directly constrain the neutron star EOS using both spectral (Lindblom 2010; Lindblom & Indik 2012, 2014) and GP (Landry & Essick 2019; Essick et al. 2020b) representations. Both methods incorporate consistency with the observed heavy pulsars PSR J0740+6620 (Fonseca et al. 2021) and PSR J0348+0432 (Antoniadis et al. 2013) as priors on the EOS and enforce thermodynamic stability and causality. Unlike other analyses where the GP representations are conditioned on both GW and pulsar data, here we only include the pulsar constraints in the prior, in order to distinguish the EOS information provided by GW230529 alone. These direct EOS inferences are also uninformative, returning their respective priors.
Figure 14. Constraints on dimensionless tidal deformability Λ for the primary (green dashed) and secondary (green solid) components of GW230529. We also show tidal deformabilities predicted for the secondary object under the BSK24 neutron star EOS (navy; Goriely et al. 2013; Pearson et al. 2018; Perot et al. 2019) and constraints obtained from BNS and pulsar observations using the GP model (GP-EOS; red).
Download figure:
Standard image High-resolution imageAppendix F: Testing General Relativity
We perform several analyses to verify whether GW230529 is consistent with general relativity (Sänger et al. 2024). Specifically, we perform parameterized tests searching for deviations in the PN coefficients that determine the phase evolution of the GW signal (Blanchet & Sathyaprakash 1994, 1995; Arun et al. 2006b, 2006a; Yunes & Pretorius 2009; Mishra et al. 2010; Li et al. 2012a, 2012b). We search for parametric deviations to the GW inspiral phasing applied to the frequency-domain waveform models IMRPhenomXP_NRTidalv2 (Colleoni et al. 2023) and IMRPhenomXPHM (García-Quirós et al. 2020; Pratten et al. 2020, 2021) using the TIGER framework (Agathos et al. 2014; Meidam et al. 2018) and SEOBNRv4_ROM_NRTidalv2_NSBH (Matas et al. 2020) and SEOBNRv4HM_ROM (Bohé et al. 2017; Cotesta et al. 2018, 2020) using the FTI framework (Mehta et al. 2023a). These waveform models each capture different physical effects (precession, higher-order multipole moments, and neutron star tidal deformability) to determine whether their absence in the model leads to inferred inconsistencies with general relativity.
For all waveform models and PN orders whose analyses have been completed, we find that GW230529 is consistent with general relativity within the inferred uncertainties on the deviation parameters. We do not yet include results about 0PN deviations, which are being examined separately owing to technical issues related to the degeneracy between the 0PN deviation parameter and the chirp mass (Payne et al. 2023). The constraints obtained at −1PN are an order of magnitude tighter than previously reported bounds for NSBH and BBH (Abbott et al. 2021b). The previously reported bounds using GW170817 remain the tightest constraints on −1PN deviations obtained with GWs (Abbott et al. 2019c).
Appendix G: Compact Binary Merger Rate Methods
A key component of the event-based rate estimation described in Section 5 is the sensitive time–volume of our GW searches. The sensitivities to GW200105-, GW200115-, and GW230529-like events are computed from O1 through the first 2 weeks of O4a according to the corresponding mass and spin posteriors from the IMRPhenomXPHM high-spin analyses of these signals. The posterior samples chosen for this analysis are consistent with previous event-based rate estimates presented for NSBH mergers (Abbott et al. 2021a). Simulated signals, whose binary parameters are drawn from the mass and spin posteriors for each signal, are distributed uniformly in comoving volume and following the other extrinsic parameter distributions given in Appendix D, and they are added to simulated detector noise characterized by representative PSDs for each observing run. The detectability of these injections is then calculated semianalytically to determine the sensitive time–volume by calculating network responses to the simulated signals and applying a threshold on the network optimal S/N of >10. This choice of threshold was previously tuned to the results of matched-filter searches (Abbott et al. 2021c) and is comparable to threshold statistics calculated for semianalytic sensitivity estimates (Essick 2023).
For the population-based approach, we estimate the merger rate of three astrophysical populations (BBH, BNS, and NSBH) by aggregating triggers found by GstLAL from O1 through the first 2 weeks of O4, while accounting for the possibility that some of these triggers are of terrestrial origin. Our population analyses do not, however, include data from the engineering run preceding the start of O4, during which another NSBH candidate was identified (LIGO Scientific, Virgo, & KAGRA Collaboration 2023b). As outlined and implemented in previous LVK results (Abbott et al. 2021a, 2023b), we construct the joint likelihood on the Poisson parameters of the astrophysical populations (ΛBBH, ΛBNS, ΛNSBH) and of the terrestrial triggers (Λbackground; Farr et al. 2015; Kapadia et al. 2020). We then extrapolate the sensitive time–volume to each astrophysical population (〈VT〉BBH, 〈VT〉BNS, 〈VT〉NSBH) obtained at the end of O3 (Abbott et al. 2023b) to account for the additional 2 weeks of observation time during O4 during which GW230529 was detected. Using a uniform prior on the merger rates , we infer their posterior distributions. The three astrophysical populations are defined by dividing up the space of compact binary component masses and spins into disjoint regions. We consider all components with masses between 1 and 3 M⊙ to be neutron stars and assume a maximum dimensionless spin of 0.05 for such components; all components that are more massive are assumed to be black holes with a maximum dimensionless spin of 0.99. The distribution of component masses within each region is assumed to be a power law in primary mass with index −2.35 and uniform in secondary mass.
Even though we infer the posterior distributions of the merger rates of all three populations, we only present given that GW230529 contributes negligibly to and .
Appendix H: Additional Details of Population Analyses
We use three different models to analyze the population of compact-object binaries with and without GW230529.
H.1. Binned Gaussian Process
The Binned Gaussian Process method models the merger rate density per log component masses as a piecewise constant function. By inferring the merger rate density in each bin, we reconstruct the shape of the CBC mass spectrum up to the resolution limit imposed by our choice of binning. A GP prior with an exponential quadratic kernel is imposed on the logarithmic rate densities to smooth out the inferred shapes over sparse regions of the parameter space. To assess the impact of GW230529 on the shape of the CBC mass spectrum, we use the same bin locations and priors on the GP hyperparamters as the GWTC-3 analysis (Abbott et al. 2023b). The means (μ), correlation length (l), and covariance amplitude (σ) of the GP are drawn from normal, lognormal, and half-normal priors, respectively. Further details of these priors are summarized in Table 5.
Table 5. Summary of Binned Gaussian Process Model Parameters (Mandel et al. 2017; Mohite 2022; Ray et al. 2023b)
| Parameter | Description | Prior |
|---|---|---|
| μ | Mean in each bin | N(0, 10) |
| σ | Amplitude of the covariance kernel | HN(0, 10) |
| of the covariance kernel | N(−0.085, 0.93) |
Note. Arguments in the priors specify the mean and standard deviation in a normal (N) or half-normal (HN) distribution.
Download table as: ASCIITypeset image
Unlike more recent implementations of the Binned Gaussian Process model that also fit for the redshift distribution (Ray et al. 2023b), we assume a redshift distribution such that the overall merger rate of compact binaries is uniform in comoving volume and source-frame time. This facilitates a direct comparison with the findings of the GWTC-3 analysis (Abbott et al. 2023b) and avoids any systematic biases in Binned Gaussian Process–based redshift distribution models originating from the inclusion of low-mass events, which remains an active area of study. As in previous analyses (Abbott et al. 2023b; Ray et al. 2023b), we fix the spin distributions for each component to be isotropic in direction and uniform in spin magnitude.
H.2. Power law + Dip + Break
The Power law + Dip + Break model is designed to search for a separation in masses between neutron stars and black holes by employing a broken power law with a dip at the location of the power-law break. The dip is modeled by a notch filter with depth A, which is fit along with other model parameters in order to determine the existence and depth of a potential mass gap (Farah et al. 2022). A value A = 0 corresponds to no gap, whereas A = 1 corresponds to zero merger rate over the interval of the gap, i.e., a maximally deep gap. Power law + Dip + Break also employs a low-pass filter at high black hole masses to allow for a tapering of the mass spectrum, which has the effect of adding a smooth second break to the power law.
The component mass distributions in this model are both fit by the same broken power law with exponents α1 between and and α2 between and . The model additionally includes a power-law pairing function in mass ratio (Fishbach & Holz 2020), assumed to be the same for all component masses. The parameters for this mass model are summarized in Table 6. Like the Binned Gaussian Process model, the Power law + Dip + Break model additionally assumes that CBCs are uniformly distributed in comoving volume and source-frame time and that component spins are isotropically oriented and uniformly distributed in magnitude. These assumptions are made for simplicity and consistency with GWTC-3 analyses (Abbott et al. 2023b). Components with m < 2.5 M⊙ are limited to spin magnitudes <0.4, while components with m > 2.5 M⊙ can have spin magnitudes up to 0.99.
Table 6. Summary of Power law + Dip + Break Model Parameters (Farah et al. 2022)
| Parameter | Description | Prior |
|---|---|---|
| α1 | Spectral index for the power law of the mass distribution below | U(−8, 2) |
| α2 | Spectral index for the power law of the mass distribution above | U(−3, 2) |
| A | Lower mass-gap depth | U(0, 1) |
| (M⊙) | Location of the lower end of the mass gap | U(1.4, 3) |
| (M⊙) | Location of the upper end of the mass gap | U(3.4, 9) |
| ηlow | Parameter controlling how the rate tapers at the low end of the mass gap | 50 |
| ηhigh | Parameter controlling how the rate tapers at the high end of the mass gap | 50 |
| η | Parameter controlling tapering the power law at high black hole mass | U(−4, 12) |
| β | Spectral index for the power law-in-mass ratio pairing function | U(−2, 7) |
| mmin (M⊙) | Minimum mass of the mass distribution | U(1 , 1.4) |
| mmax (M⊙) | Maximum mass of the mass distribution | U(35, 100) |
| Maximum allowed component spin for objects with mass <2.5 M⊙ | 0.4 | |
| Maximum allowed component spin for objects with mass ≥2.5 M⊙ | 0.99 |
Note. Arguments in the Prior column specify the lower and upper bounds of a uniform (U) distribution. The first several entries describe the mass distribution parameters, and the last two entries describe the spin distribution parameters.
Download table as: ASCIITypeset image
H.3. NSBH-pop
The NSBH-pop model (Biscoveanu et al. 2022) is designed to capture the mass and spin distributions of the NSBH population. The black hole mass distribution is fit by a truncated power law, and the conditional mass ratio distribution p(q∣m1) is fit by a truncated Gaussian between and . The black hole spin magnitude distribution is modeled as a Beta distribution (including singular distributions that peak at the edges of the χ1 ∈ [0, 0.99] space), while the neutron star spin magnitude distribution is assumed to be uniform over χ2 ∈ [0, 0.7] to encapsulate the effect of the neutron star breakup spin at the mass-shedding limit (Most et al. 2020a; Shao et al. 2020). We assume that the redshift distribution is uniform in comoving volume and source-frame time and that spins are isotropically oriented. The parameters for this model are summarized in Table 7.
Table 7. Summary of NSBH-pop Model Parameters (Biscoveanu et al. 2022)
| Parameter | Description | Prior |
|---|---|---|
| α | Black hole mass power-law index | U (−4, 12) |
| (M⊙) | Minimum black hole mass | U(2, 10) |
| (M⊙) | Maximum black hole mass | U(8, 20) |
| (M⊙) | Maximum neutron star mass | U(1.97, 2.7) |
| μ | Mass ratio mean | U(0.1, 0.6) |
| σ | Mass ratio standard deviation | U(0.1, 1) |
| αχ | Beta distribution shape parameter (α) for black hole spin distribution | U(0.1, 10) |
| βχ | Beta distribution shape parameter (β) for black hole spin distribution | U(0.1, 10) |
Note. Arguments in the Prior column specify the lower and upper bounds of a uniform (U) distribution.
Download table as: ASCIITypeset image
H.4. Population Reweighting
We use the population-level mass and spin distributions inferred using all three population models to reweight (e.g., Payne et al. 2019) the high-spin combined parameter estimation samples from the default priors described in Appendix D to three different astrophysically motivated priors given by the one-left-out posterior predictive distributions (Galaudage et al. 2020; Callister 2021; Essick & Fishbach 2021) inferred for each model. A comparison of the posteriors under the default and astrophysically motivated priors is shown in Figure 7.
Because of the q–χeff degeneracy, the spin distribution assumptions also impact the inferred component masses, and hence classification, for this source. However, the fact that the Binned Gaussian Process and NSBH-pop model priors recover similar mass posteriors indicates that the different spin distribution assumptions have a subdominant effect compared to the different mass distribution assumptions.
H.5. Selection Effects
For our population analyses, we account for search selection bias to measure the underlying astrophysical mass and spin distributions of the NSBH and CBC populations (Loredo 2004; Farr 2019; Mandel et al. 2019; Vitale et al. 2020). To estimate the sensitivity of the searches over the data recorded by the detector network across the binary parameter space, we use a large suite of injections and impose a significance threshold on the FAR. In contrast to the two methods for calculating merger rates described in Appendix G, we use the same set of injections used for GWTC-3 population studies recovered using matched-filter searches (Abbott et al. 2023b). By using the combined injection sets from the first three observing runs (Abbott et al. 2021c), we are effectively assuming that GW230529 occurred at the end of O3. Without accounting for the extra time–volume provided by the first 2 weeks of O4, this introduces a bias in our inferred estimates of the merger rate in the mass gap (Figure 5), rate of gamma-ray bursts with NSBH progenitors, and total contribution of NSBH mergers to the production of heavy elements (Section 7). However, this effect is negligible given that GW230529 occurred in the first 2 weeks of O4a and the detector sensitivity during this period was not drastically different from that during O3 (see Section 2), meaning that the additional time–volume unaccounted for in the injection sets is small.
To generate samples from the hierarchical likelihood for our population analyses, we use the Dynesty nested sampler (Speagle 2020) as implemented in GWPopulation (Talbot et al. 2019) for the Power law + Dip + Break and NSBH-pop analyses and the Hamiltonian Monte Carlo sampler implemented in the PyMC package (Salvatier et al. 2016) for the Binned Gaussian Process analysis.
Appendix I: Equation of State
To assess the effect of the EOS constraint used for source classification (Section 6) and calculation of fEM-bright (Section 7), we repeat these analyses using GP-EOS constraints additionally conditioned on NICER observations (Legred et al. 2021, 2022) of J0030+0451 (Miller et al. 2019; Riley et al. 2019; Vinciguerra et al. 2024) and J0740+6620 (Miller et al. 2021; Riley et al. 2021; Salmi et al. 2022, 2023).
We find that including the NICER observations only changes the source classification results by at most a few percent. This is because the constraint on the maximum neutron star mass used in the source classification analysis is primarily driven by the observation of massive pulsars. Changing the information included within the EOS constraint leads to a smaller effect than the choice of astrophysical mass and spin distribution, as discussed in Section 6.
The inference on the fraction of NSBH systems that may have EM counterparts changes more significantly when NICER observations are included in the EOS constraint. This is because the GW and pulsar-only results provide much more support for compact neutron stars, which inhibits the probability that the neutron star is disrupted and suppresses the probability of an EM counterpart. The constraints including the NICER observations favor stiffer EOSs, enhancing the inferred EM-bright fraction, , and pulling the posterior peak away from zero.
