Distinguishing the Demographics of Compact Binaries with Merger Entropy Index

Based on the Hawking Area Theorem (Hawking, 1971), the entropy of a black hole correlates to the surface area of its event horizon. By expressing the Bekenstein-Hawking Formula (Bekenstein, 1973) in terms of mass (m) and spin ($\chi$), we can compute the entropy of an individual black hole as:

According to the second law of thermodynamics, mergers between BBHs will increase the entropy of the universe (Bekenstein, 1973). Based on this theorem, we introduced a new variable called the Merger Entropy Index ($\mathcal{I}_{\mathrm{BBH}}$) in an earlier study (Hu et al., 2021), which measures the efficiency of entropy transfer during the BBH mergers:

where $\Delta S_{BBH}=S_{\mathrm{f}}-(S_{1}+S_{2})$ measures the increase in entropy during the merger ($S_{\mathrm{f}}$ is the entropy of the post-merger black hole, and $S_{1},~{}S_{2}$ are entropy of the pre-merger black holes). General relativity bounds this parameter between $[0,~{}1]$.

For all the GW events reported by LVK, we utilize the parameter estimation results made available in the GW Open Science Center (Abbott et al., 2021c, 2023c). For BBH events, we use posteriors from SEOBNRv4PHM waveform model (Buonanno & Damour, 1999, 2000; Damour et al., 2008; Ossokine et al., 2020; Gadre et al., 2022), while for NSBH we use the IMRPhenomXPHM posteriors (Khan et al., 2016, 2019; Pratten et al., 2020, 2021a; D’Ambrosio et al., 2024). The files also include the mass and spin associated with the corresponding remnant black hole, which has been computed using phenomenological fits to NR results (Santamaría et al., 2010; Hofmann et al., 2016; Healy & Lousto, 2017; Jiménez-Forteza et al., 2017; Estellés et al., 2022). We input the posteriors of primary, secondary, and remnant black holes in Eqs. 1-2 and compute the corresponding $\mathcal{I}_{\mathrm{BBH}}$ posteriors for each LVK event, as highlighted in Fig. 1.

To test the $\mathcal{I}_{\mathrm{BBH}}$ formalism for upper and lower mass gap events, we employ six different BBH population models, inspired by astrophysics and detection parameter space. These models are treated as uncorrelated and independent observation channels.

Uniform Prior - we adopt the injection file for all LVK observing runs to estimate the searching sensitivity (Abbott et al., 2023a, b). From the injection, we use the parameters of detectable signals without any $p_{\mathrm{astro}}$ filtering to generate the uniform prior. This model assumes a uniform distribution for spin magnitude $(\chi_{1,2})$ and polar angle $(\theta_{1,2})$ of the spin vectors across the entire range, i.e. $\chi_{1,2}\in[0,~{}1]$ and $\cos{\theta_{1,2}}\in[-1,~{}1]$ (Bhagwat et al., 2021). The distribution of mass ratio $m_{2}/m_{1}\equiv q$ is nearly flat across [0,1], with slightly higher probability on the equal-mass region.

Isolated BBHs - in this model, BBHs are expected to have little spin misalignment, so we strictly constrain $\theta_{1,2}\in\{0,~{}\pi\}$ (Steinle & Kesden, 2021). For mass ratio and spin distributions, we assume uniform distributions across the range $q\in[0.5,~{}1]$, $\chi_{1,2}\in[0,~{}0.5]$, as shown in Zevin et al. (2021).

Dynamical BBHs - this model includes dynamical capture and hierarchical merger scenario (1G+1G, 1G+2G, 2G+2G) involving second-generational black holes (2G), which are the remnant black holes of 1G mergers (Yang et al., 2019; Kimball et al., 2020; Fragione & Rasio, 2023). We adopt mass ratio and spin distributions derived in Kimball et al. (2021) and assume a uniform distribution for polar angle of the spin vectors. Both 1G+1G and 2G+2G populations favor the equal-mass regime, while q distribution of 1G+2G shows a sharp Gaussian distribution peaking at q=0.5. Spin distribution of 1G+1G resembles an exponentially decaying curve with MODEL TRUNCGAUSS, with the highest probability at $\chi=0$. The spin distribution of 2G+2G is a Gaussian distribution with mean of 0.7. $\chi_{1},\chi_{2}$ have different distributions in 1G+2G population, overlapping with the spin distribution of 1G+1G and 2G+2G respectively.

GWTC-3 PowerLaw + Peak (PP) - we adopt the PowerLaw + Peak model (Talbot & Thrane, 2018), which demonstrates the highest Bayes factor in characterizing the Gravitational-Wave Transient Catalog (GWTC) dataset from the Third Observing Run. The most notable feature of this model is the small peak at $\sim 35~{}M_{\odot}$ in the primary mass distribution due to the PISN mass gap (Farmer et al., 2019, 2020; Edelman et al., 2022). The detailed parameter distributions are described in Abbott et al. (2023b).

For each population model, we randomly sample $i=10^{4}$ binaries with inspiral parameters $\Lambda^{\mathrm{ins}}_{i}:\{q,~{}\chi_{1,2},~{}\theta_{1,2}\}$. The bounds on $\Lambda^{\mathrm{ins}}_{i}$ are dictated by the population model. To obtain the corresponding parameters of the remnant black hole, $\Lambda^{\mathrm{rem}}_{i}:\{m_{\mathrm{f}},~{}\chi_{\mathrm{f}}\}$, we utilize the phenomenological fit NRSur7dq4Remnant (Blackman et al., 2017; Boschini et al., 2023) from the SurfinBH package (Varma et al., 2018, 2019; Zertuche et al., 2024). For each $\Lambda_{i}$, we compute $\mathcal{I}_{\mathrm{BBH},i}$ using Eq. 1 - 2. The high sampling rate provides sufficient accuracy to compare the posteriors of $\mathcal{I}_{\mathrm{BBH}}$ between different population models and to test them with those from LVK events. The probability distribution derived from the random sampling for each population model is shown in Fig. 2.

The catalog of GW events observed by the LVK detectors is affected by the selection effects (Singer et al., 2014; Abbott et al., 2016; Chen et al., 2017). This introduces a bias when constructing the $\mathcal{I}_{\mathrm{BBH}}$ posteriors. To correct for the selection effect, we utilize the formalism described in Doctor et al. (2021) by introducing a weight factor for each $\Lambda_{i}$:

We obtain $p(\Lambda_{i}|z_{\mathrm{det}},z_{\mathrm{ref}})$ from the injection file of combined three LVK observing runs. To compute $p(\Lambda_{i}|z_{\mathrm{det}})$, we adopt the PYCBC BBH detection pipeline and compute $p_{\mathrm{astro}}$ for each inspiral parameter (Tiwari, 2018; Farr, 2019; Abbott et al., 2023b). We only take detectable $\Lambda_{i}$ with $p_{\mathrm{astro}}\geq 0.5$. We then apply $\mathrm{w}_{i}$ to the post-merger parameters $\Lambda^{\mathrm{rem}}_{i}$ derived from numerical relativity in Sec. 2.2. All population models used have already accounted for selection effect. Since the selection effect is relevant when comparing population instead of individual event (Thrane & Talbot, 2019), we only correct when grouping LVK events together.

After correcting for the selection effect, we perform Kolmogorov–Smirnov (KS) tests (Plowman et al., 2011) among population models as well as their comparison with individual LVK events. The KS tests compute $\mathcal{D}$ statistics, which is a measure of the inherent difference between the two distributions.

In this study, we quote favored population model for individual and group of LVK events based on $(1-\mathcal{D})$. These values have been quoted in Table 2 and Fig. 3.

In Fig. 1, we showcase the $\mathcal{I}_{\mathrm{BBH}}$ posteriors for 85 LVK events that broadly qualify as BBH systems. For further investigation, we group these BBH events into four categories based on the median value of the component masses in their cosmologically corrected source frame: (i) 74 stellar BBH events in which $m_{1,2}\in[5,~{}60]~{}M_{\odot}$; (ii) 3 PISN-PISN events with $m_{1},m_{2}>60~{}M_{\odot}$, wherein both the primary and secondary are in upper gap $(60-130~{}M_{\odot})$; (iii) 6 PISN BH-Stellar BH events with $m_{1}>60~{}M_{\odot}$, wherein only the primary object is in upper gap, and (iv) 2 stellar BH-mystery BH events, wherein the secondary object is in the lower gap $(2-5~{}M_{\odot})$. Table 1 highlights the GW events in the special categories (ii)-(iv) and their corresponding $\mathcal{I}_{\mathrm{BBH}}$ values.

For the combined 85 LVK BBH events, we find that $\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.12}_{-0.08}$ with 90% confidence intervals. This is consistent with our expectation for the astrophysical prior of $\mathcal{I}_{\mathrm{BBH}}\leq 0.3$ derived in the earlier work (Hu et al., 2021). For categories describing (i) stellar BBH, we find $\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.11}_{-0.08}$; (ii) PISN-PISN BBH, we find $\mathcal{I}_{\mathrm{BBH}}=0.21^{+0.11}_{-0.08}$; (iii) PISN BH-Stellar BH, we find $\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.16}_{-0.07}$; (iv) stellar BH-mystery BH, we find $\mathcal{I}_{\mathrm{BBH}}=0.06^{+0.02}_{-0.01}$. For each of these cases, the combined posterior of $\mathcal{I}_{\mathrm{BBH}}$ is corrected for the selection effects as described in Sec. 2.3.

Among all LVK events, GW191109 shows the largest $\mathcal{I}_{\mathrm{BBH}}$ = $0.33^{+0.13}_{-0.08}$, followed by GW190521 $\mathcal{I}_{\mathrm{BBH}}$ = $0.26^{+0.09}_{-0.07}$. Both these events belong to the high-mass population (categories ii, iii), with primary masses of $65^{+11}_{-11}$ and $85^{+21}_{-14}~{}M_{\odot}$ respectively. Other most massive LVK BBHs, GW190426 (from GWTC-2.1 (Abbott et al., 2024)) and GW200220, have slightly higher $\mathcal{I}_{\mathrm{BBH}}$ than the combined stellar BBH events.

As shown in our earlier work (Hu et al., 2021), $\mathcal{I}_{\mathrm{BBH}}$ has a strong dependence on the orientations of black hole spins. For example, a binary with two highly spinning black holes can either result in a remnant black hole with a moderate spin or maximal spin depending on the spin evolution during the inspiral. Yet, the first case would result in a significant increase in entropy than the latter. $\mathcal{I}_{\mathrm{BBH}}$ is particularly sensitive to the effective inspiral spin parameter $\chi_{\mathrm{eff}}$ when the binaries have near to moderate mass-ratios. A negative (positive) value of $\chi_{\mathrm{eff}}$ generally results in a higher (lower) entropy configuration and thus a larger (smaller) value for $\mathcal{I}_{\mathrm{BBH}}$. As we demonstrate in Sec. 3.3, this feature of $\mathcal{I}_{\mathrm{BBH}}$ makes for a powerful probe to test astrophysical formation channels for a GW event.

In the case of GW191109, the higher value of $\mathcal{I}_{\mathrm{BBH}}$ can be mainly attributed to significant anti-alignment of black hole spins with respect to the orbital angular momentum ($\chi_{\mathrm{eff}}=-0.29^{+0.42}_{-0.31}$). The spin configuration of GW191109 provides hints to dynamical formation scenario (Romero-Shaw et al., 2022; Antonelli et al., 2023; Zhang et al., 2023). However, GW191109 suffers from data quality issues and glitch (Abbott et al., 2021d; Wang et al., 2022; Udall et al., 2024), which affects the parameter estimation including negative $\chi_{\mathrm{eff}}$ and can be the origin of an inconsistent $\mathcal{I}_{\mathrm{BBH}}$ in comparison to other LVK events.

On the other hand, GW190403 has the lowest $\mathcal{I}_{\mathrm{BBH}}=0.14^{+0.11}_{-0.08}$ among all the 9 events in the high-mass population. This is consistent with expectations since the GW190403 has highest recorded spin-orbit alignment among all LVK events ($\chi_{\mathrm{eff}}=0.68^{+0.16}_{-0.43}$), thus a stronger preference for smaller $\mathcal{I}_{\mathrm{BBH}}$. The large positive $\chi_{\mathrm{eff}}$ of GW190403 indicates significant alignment between spin and orbit, providing hints to isolated formation scenario (Belczynski et al., 2020; Olejak & Belczynski, 2021; Bavera et al., 2020; Qin et al., 2022; Abbott et al., 2024).

Furthermore, GW190403 is the most asymmetric binary in the high-mass population $(q={0.25}_{-0.11}^{+0.54})$. Highly asymmetric BBH has the least entropy transfer and therefore drives $\mathcal{I}_{\mathrm{BBH}}$ towards its lower bound. As seen in Fig. 1, GW190403 also shows the largest standard deviation in $\mathcal{I}_{\mathrm{BBH}}$ distribution among all the LVK events. This event has a low SNR of ${7.6}^{+0.6}_{-1.1}$ (Abbott et al., 2024), resulting in a much larger uncertainty in the inspiral parameters, and subsequently broader $\mathcal{I}_{\mathrm{BBH}}$ posterior. Yet, the GW190403 does exhibit multiple peaks in its $\mathcal{I}_{\mathrm{BBH}}$ distribution, which is characteristic of moderate mass-ratio binary being impacted by its spin orientations.

The link between asymmetric BBH events and lower $\mathcal{I}_{\mathrm{BBH}}$ is further illustrated by GW190412 (Abbott et al., 2020d), which shows a relatively smaller $\mathcal{I}_{\mathrm{BBH}}\approx 0.14$ in comparison to the combined posterior of stellar BBHs peaking at $\mathcal{I}_{\mathrm{BBH}}\approx 0.20$. The unequal mass ratio is not common in isolated pathways, which points to hierarchical merger origin for these two events (Banerjee, 2022; Ford & McKernan, 2022).

The formalism of $\mathcal{I}_{\mathrm{BBH}}$ is deemed unsuitable for BNS or NSBH events due to the uncertainty in the entropy calculation for their pre- and post-merger objects. However, the formalism can be applied to GW events with secondary objects in the lower mass gap (category-iv), such as GW190814 and GW200210 (Abbott et al., 2020a, 2021a, 2023a), by treating them as BBH. As shown in Fig. 1, both these events in the lower mass gap display sharp spikes around $\mathcal{I}_{\mathrm{BBH}}\sim 0.06$. Since the mass-ratio of these binaries is $q\sim 0.1$, their spin orientations $(\chi_{\mathrm{eff}}\sim 0)$ have no substantial impact on $\mathcal{I}_{\mathrm{BBH}}$ unless the black holes are highly anti-aligned (Hu et al., 2021). Since the secondary mass for GW190814 is measured and constrained more precisely than GW200210, their $\mathcal{I}_{\mathrm{BBH}}$ posterior is narrower.

In Fig. 2, we highlight the distribution of $\mathcal{I}_{\mathrm{BBH}}$ for all six population models described in Sec. 2.2. Note that across all the models, $\mathcal{I}_{\mathrm{BBH}}\gtrsim 0.1$. This lower bound is dictated by the fact that all six models predict near-equal to moderate mass-ratio BBH systems. On the other bound, none of the models predict $\mathcal{I}_{\mathrm{BBH}}\gtrsim 0.6$. This upper bound results from our assumption that binaries in all of the models are quasi-circular, which is a simplification since the dynamical channel could have some residual eccentricity event in the LVK frequency range (Romero-Shaw et al., 2020, 2022). As discussed previously in Hu et al. (2021), highly eccentric binaries with spins are more efficient than circular binaries for entropy transfer, thus resulting in a larger $\mathcal{I}_{\mathrm{BBH}}$.

The population model 2G+2G results in the largest value of $\mathcal{I}_{\mathrm{BBH}}=0.27^{+0.09}_{-0.08}$, while the isolated formation model has the lowest value of $\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.05}_{-0.04}$. We find that value of $\mathcal{I}_{\mathrm{BBH}}$ decreases with the generations of hierarchical merger, with 2G+2G being the highest, followed by 1G+2G ($\mathcal{I}_{\mathrm{BBH}}=0.22^{+0.08}_{-0.07}$) and then 1G+1G $(\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.09}_{-0.01})$. In comparison, the signature GWTC-3 Powerlaw+Peak (PP) model shows $\mathcal{I}_{\mathrm{BBH}}=0.20^{+0.05}_{-0.03}$, which, as expected, is consistent with the combined posterior of all LVK events discussed in Sec. 3.1. When considering the uniform prior for BBH, $\mathcal{I}_{\mathrm{BBH}}$ shows the broadest distribution from $0.24^{+0.16}_{-0.07}$.

We gain further insights into the distinction between population models by conducting KS test of $\mathcal{I}_{\mathrm{BBH}}$ distribution. We summarize the KS test results in Table 2. We find that the dynamical capture 1G+1G and GWTC-3 PP has the highest similarity of 85% between their $\mathcal{I}_{\mathrm{BBH}}$ distributions. Closely followed is the similarity of 79% between isolated formation model and GWTC-3 PP. Both 1G+1G and isolated formation have similar mass ratio distribution that favors the equal mass region (Giacobbo & Mapelli, 2018; Bouffanais et al., 2019; Romero-Shaw et al., 2020). We also note that the hierarchical merger scenario of 2G+2G is the least similar (34%) with the $\mathcal{I}_{\mathrm{BBH}}$ distribution of GWTC-3 PP. These comparisons are well within the expectations (Zevin et al., 2021) and demonstrate that $\mathcal{I}_{\mathrm{BBH}}$ could be an independent probe for the underlying population of the observed BBH events.

Since $\mathcal{I}_{\mathrm{BBH}}$ is mass invariant by definition, it will not recognize features such as the peak in GWTC-3 PP model before the PISN mass-gap or that the 2G+2G scenario is applicable only to events inside the PISN mass-gap. Yet, the index is sensitive to the underlying spin and mass-ratio differences between these models. Therefore, we notice that $\mathcal{I}_{\mathrm{BBH}}$ distribution of 1G+1G (lower mass BBH) has the least similarity with 2G+2G (34%). In the latter model, the spin distribution peaks $\sim$ 0.7 as the binary transports the angular momentum from the earlier first-generation merger to the remnant black hole. Similarly, the $\mathcal{I}_{\mathrm{BBH}}$ distribution of 1G+1G favors the uniform prior only with 55%, which is a reflection of the mass-ratio being near equal for dynamical capture of the first generation mergers (Kimball et al., 2021). The low similarity of $\mathcal{I}_{\mathrm{BBH}}$ for 2G+2G model in comparison to all other models enables us to distinguish within dynamical channels as well as isolated formation channels.

On the other hand, 1G+2G has the most generalizable $\mathcal{I}_{\mathrm{BBH}}$ distribution among all population models in this study. As shown in Table 2, the similarity of 1G+2G model ranges from $81\%$ for uniform prior, $73\%$ for isolated binaries to the lowest of $64\%$ with 2G+2G. While the 1G+2G model is referring to events with primary black hole in the upper gap, the underlying spin and mass-ratio distribution has much larger overlap with other channels that are suitable only for lower or higher masses. For example, as described in Sec. 2.2, isolated has a uniform $q$ distribution within [0.5, 1], 2G+2G population is mostly composed of $q\sim 1$, while 1G+2G prefer $q=0.5$. In the next section, we discuss that 1G+2G model is indeed the most consistent dynamical formation channel with LVK events in categories-ii,iii.

Fig. 3 shows the results of KS convergence tests between $\mathcal{I}_{\mathrm{BBH}}$ posteriors of 6 population models described in Sec. 2.2 and individual LVK events from category ii, iii in Table 1. The consistency between BBH groups and uniform prior serves as a baseline for comparison. For categories describing (i) stellar BBH: $\mathcal{I}_{\mathrm{BBH}}$ shows the highest consistency with 1G+1G (76%) isolated models (79%), while only 48% consistency is observed with the 2G+2G population; (ii) PISN-PISN BBH: $\mathcal{I}_{\mathrm{BBH}}$ favors 2G+2G (83%); (iii) PISN BH-Stellar BH: $\mathcal{I}_{\mathrm{BBH}}$ favors 1G+2G (94%). $\mathcal{I}_{\mathrm{BBH}}$ of GWTC-3 PP model shows an 82% match with category i, but only a 52% match with category ii. The high spin typical of the 2G+2G population contributes to the lower consistency in $\mathcal{I}_{\mathrm{BBH}}$ between the PISN-PISN BBH and the GWTC-3 population. Stellar BBH (category i) also make up the majority of the GWTC-3 population (Abbott et al., 2023b).

GW190521, the second-most massive BBH event reported to date, with both black holes residing in the PISN mass gap, has $\mathcal{I}_{\mathrm{BBH}}=0.26^{+0.09}_{-0.07}$, which shows a 91% match with 2G+2G model. In contrast, $\mathcal{I}_{\mathrm{BBH}}$ of this event has lower consistency to other models, especially for 1G+1G by only 39% similarity. $\mathcal{I}_{\mathrm{BBH}}$ distribution favors the 2G+2G model, suggesting that the high progenitor masses likely result from hierarchical mergers. This conclusion matches with results from many recent literature Romero-Shaw et al. (2020); Abbott et al. (2020b); Graham et al. (2020); Estellés et al. (2022); Holgado et al. (2021); Gamba et al. (2022); Morton et al. (2023). Additionally, $\mathcal{I}_{\mathrm{BBH}}$ of GW190521 has low consistency (41%) with GWTC-3 populations, indicating that a flatter PP model may be required to fit these high-mass events (Kimball et al., 2021).

The most massive BBH detected so far in GW190426 has $\mathcal{I}_{\mathrm{BBH}}=0.21^{+0.11}_{-0.08}$, which is consistent with $\mathcal{I}_{\mathrm{BBH}}$ posteriors of GWTC-3 population by 80%. Unlike GW190521, $\mathcal{I}_{\mathrm{BBH}}$ of this event matches with 2G+2G the least (55%), while shows higher consistency with 1G+2G (91%) and isolated (82%) models. The primary black hole, with a mass of $106^{+45}_{-24}~{}{M}_{\odot}$, can reasonably be second-generational. However, it is unlikely that the secondary black hole would reach $76^{+26}_{-36}~{}{M}_{\odot}$ solely through dynamical capture mechanisms. With our 1G+2G model described in Sec. 2.2, P$(m_{2}=76~{}{M}_{\odot})\sim 0$, suggesting that GW190426 is unlikely to be a 1G+2G binary despite its high consistency in $\mathcal{I}_{\mathrm{BBH}}$. GW190426 shows no orbital eccentricity (Romero-Shaw et al., 2022), which is a signature of dynamical formation channel. This lack of eccentricity, combined with the 82% $\mathcal{I}_{\mathrm{BBH}}$ consistency with the isolated model, suggests that this massive BBH may evolve from isolated formation channels.

GW200220 shows the highest similarity (93%) of $\mathcal{I}_{\mathrm{BBH}}$ with a population model (1G+2G). The secondary black hole $m_{2}=61^{+21}_{-25}~{}{M}_{\odot}$ lies near the boundary of the PISN mass gap ($60~{}{M}_{\odot}$) by our definition. However, the exact bound of the PISN mass gap remains uncertain and may extend to 65 $M_{\odot}$ (Woosley & Heger, 2021). When defining PISN mass gap as $[65,~{}130]~{}M_{\odot}$, the mass distribution of GW200220 matches with 1G+2G binary. The secondary black hole may be first-generational given P($m_{2}\leq$ 65 $M_{\odot}$) $\geq 60\%$. In contrast to GW190426, GW200220 shows explicitly higher consistency in $\mathcal{I}_{\mathrm{BBH}}$ with 1G+2G (93%) than isolated models (75%), so $\mathcal{I}_{\mathrm{BBH}}$ disfavors the isolated formation channel.

On the other hand, GW191109 has the least similarity across all population models. The $\mathcal{I}_{\mathrm{BBH}}$ posterior for this event is only 10% consistent with isolated BBHs and 14% consistent with 1G+1G binaries, but shows a relatively higher 57% match with the 2G+2G population, supporting the hierarchical merger model for the massive primary black hole and spin configuration. Since GWTC-3 PP model is derived with a strong bias towards $\chi_{\mathrm{eff}}\sim 0$ (Abbott et al., 2023b), $\mathcal{I}_{\mathrm{BBH}}$ between GW191109 and GWTC-3 PP only matches by 11%. As suggested in Sec. 3.1 and Table 1, this event has the largest $\mathcal{I}_{\mathrm{BBH}}$ due to significant spin-orbit misalignment. Since this event does not show high similarity (green category) with any population model, we provide an independent confirmation that GW191109 may be affected by a glitch.

For all the six events in category iii, we note their $\mathcal{I}_{\mathrm{BBH}}$ favors 1G+2G model more than 1G+1G model. This is consistent with the simulations from Kimball et al. (2021), which show PISN BH-Stellar BH binaries have mass ratio distributions consistent with 1G+2G population. The $\mathcal{I}_{\mathrm{BBH}}$ distribution of GW190602 is 85% favored towards 1G+2G and about the same as the isolated channel, making it challenging to differentiate between them. However, the remaining events, GW190519, GW190706, and GW190929, show a slight preference for the isolated formation channel. Due to the low SNR and broader $\mathcal{I}_{\mathrm{BBH}}$ posterior discussed in Sec. 3.1, the $\mathcal{I}_{\mathrm{BBH}}$ formalism of GW190403 does not favor any of six population model with medium similarity.

In April 2024, LVK announced GW230529 from O4 (Abac et al., 2024), a potential BNS coalescence or NSBH merger event. The primary compact object has a mass of $3.6^{+0.8}_{-1.3}M_{\odot}$, which sits in the lower mass gap with 99% confidence interval. The primary object is too heavy to fit into our current understanding of neutron star mass distribution (Özel et al., 2012; Kiziltan et al., 2013; Landry & Read, 2021) and the neutron star equation of state (Godzieba et al., 2021). The secondary object ($m_{2}=1.4^{+0.6}_{-1.2}M_{\odot}$) is likely a neutron star.

To understand the demographics of GW230529, we apply the $\mathcal{I}_{\mathrm{BBH}}$ formalism to the available parameter estimation samples from the IMRPhenomXHM model (Pratten et al., 2021b). This model does not incorporate tidal disruptions expected (Zhu et al., 2024a), but instead assumes coalescence between light black holes, thus permitting us to compute the merger entropy. The $\mathcal{I}_{\mathrm{BBH}}$ distribution of GW230529 is shown in Fig. 4. Alongside, we have included $\mathcal{I}_{\mathrm{BBH}}$ distribution from the three known NSBH events GW191219, GW200105, and GW200115. For these NSBH events, we have the samples from IMRPhenomXHM model to compute $\mathcal{I}_{\mathrm{BBH}}$.

When comparing with individual LVK NSBH events, we find that $\mathcal{I}_{\mathrm{BBH}}$ of GW230529 shows the highest similarity with GW200115 (61%), which is a confirmed NSBH event with $q=0.26^{+0.35}_{-0.10}$. In contrast, GW230529 favors the other two asymmetric NSBH events, GW191219 and GW200105, only with 1.5% and 17% respectively.

Furthermore, we derive the $\mathcal{I}_{\mathrm{BBH}}$ distribution for an NSBH population model described in Biscoveanu et al. (2022). This model features the parameter distributions for binaries with secondary compact object $\lesssim 5~{}M_{\odot}$. Here, the mass ratio distribution peaks at $q\sim 0.1$. The spin distribution of neutron stars is assumed to be uniform between 0 and 0.7. The $\mathcal{I}_{\mathrm{BBH}}$ posterior of NSBH population model peaks $\sim 0.06$. We find that the $\mathcal{I}_{\mathrm{BBH}}$ of GW230529 favors the GWTC-3 PP model (60%) over the NSBH population model (33%).