© The Authors 2026

1 Introduction

Direct observations of binary neutron stars (BNSs) come primarily from radio observations of binary pulsars in the Milky Way (Özel & Freire 2016)¹, and thus probe the local population. The mass distribution of neutron stars (NSs) in these systems is found to be much narrower than that of the overall NS population (Alsing et al. 2018; Antoniadis et al. 2016; Farr & Chatziioannou 2020; Shao et al. 2020), and can be well described by a Gaussian distribution centred at 1.33 M_⊙ with a standard deviation of 0.09 M_⊙ (Özel et al. 2012; Özel & Freire 2016), although there is some evidence of bimodality between recycled and non-recycled pulsars, the former having slightly larger masses (Farrow et al. 2019). This particularity of double radio pulsars could be attributed to the fact that the components of BNSs undergo much less accretion than neutron stars in other types of binaries, or to different masses at birth (Tauris et al. 2017).

Gravitational waves (GWs) provide a new means to directly observe these systems. The LIGO/Virgo/KAGRA (LVK) collaboration has so far confidently reported the detection of two BNSs: GW170817 (Abbott et al. 2017) and GW190425 (Abbott et al. 2020), at z ~ 0.009 and z ~ 0.03, respectively. Next-generation ground-based observatories such as the Einstein Telescope (ET) (Branchesi et al. 2023; Abac et al. 2026) and Cosmic Explorer (CE) (Evans et al. 2021, 2023), expected to come online by the end of the next decade, should detect thousands to tens of thousands of BNSs out to redshift z ~ 3. The redshift distribution of BNS mergers encodes the convolution of the cosmic star formation rate (SFR) with the distribution of time delays between binary formation and coalescence. As such, it provides a direct probe of the physical processes that govern BNS formation. The orbital separation at the formation of the BNS, and therefore the time delay to coalescence, is determined by various stellar evolution processes, including mass transfer episodes, common envelope phases, and supernova explosions with associated natal kicks (e.g. Chruslinska et al. 2018; Giacobbo & Mapelli 2018; Vigna-Gómez et al. 2018; Neijssel et al. 2019; Santoliquido et al. 2021). Observations across a broad red-shift range with next-generation detectors will therefore enable a reconstruction of the merger rate density and place strong constraints on the underlying time-delay distribution, offering a powerful means to discriminate between formation channels. These observations will also enable tests for the presence of multiple BNS sub-populations, as current data may already suggest. While the component masses of GW170817 are consistent with those of the Galactic population, GW190425 appears significantly heavier. This discrepancy is particularly evident in terms of the total mass: ${3.4}_{-0.01}^{+0.03}{M}_{\odot }$ for GW190425, compared to 2.66 ± 0.12 M_⊙ for Galactic BNSs, which potentially points to the existence of a distinct, more massive sub-population.

This intriguing observation naturally raises the question of how GW190425 formed and why such heavy BNSs are not observed in the Milky Way. Dynamical formation in stellar clusters could, in principle, produce heavier BNSs if a massive neutron star exchanges its stellar companion with another neutron star; however, this channel is expected to contribute negligibly to the overall BNS merger rate (Ye et al. 2020; Belczynski et al. 2018). Focusing on isolated binary evolution, some studies suggest that GW190425 could have formed through pathways similar to those of Galactic BNSs (Qin et al. 2024; Nair & Stevenson 2025), and that the spins and associated magnetic fields of heavy BNSs may render them effectively radio-invisible (Safarzadeh et al. 2020; Chu et al. 2025). Alternatively, Romero-Shaw et al. (2020) proposed that massive binaries like GW190425 may originate from a distinct evolutionary channel characterised by much shorter delay times between the formation of the progenitor stars and BNS merger, potentially explaining their scarcity in the Galaxy. In this scenario, the system undergoes unstable ‘case BB’ mass transfer from a helium star onto the first-formed neutron star (Belczynski et al. 2002; Dewi & Pols 2003; Tauris et al. 2017; Vigna-Gómez et al. 2018; Iorio et al. 2023), triggering a second common-envelope phase between the neutron star and the remaining CO core. The resulting tighter orbital separation allows the binary to survive the second supernova kick and merge within a few million years. Following this idea, Galaudage et al. (2021) jointly fitted the radio and GW BNS populations and found mild evidence that GW190425 could belong to a fast-merging population comprising 8–79% of all BNSs, with associated delay times between 5 and 401 Myr (assumed to be the same for all heavy BNSs). We stress that other studies also find that such a discrepancy between the Milky Way and extragalactic populations can naturally arise in population synthesis, without the need for a distinct time-delay distribution (Kruckow 2020). As the exact pathway to produce GW190425 is still open to debate, in this work we explore the consequence of the scenario proposed by Romero-Shaw et al. (2020) for future observations. If this hypothesis is correct, we could expect to observe two sub-populations of BNSs with GWs: a light population and a heavy one, potentially exhibiting distinct redshift evolutions, which ET and CE could be able to probe.

Using the latest constraints from the LVK on the local BNS merger rate (The LIGO Scientific Collaboration 2025), we constructed mock catalogues of BNS observations on which we performed hierarchical Bayesian analyses to reconstruct the astrophysical distributions. First, considering the case of a single population, we illustrate the accuracy with which the time-delay distribution can be reconstructed. We then turn to the possibility of multiple sub-populations. We find that approximately 500 events are sufficient to identify the presence of two sub-populations, while a few thousand detections allow us to distinguish their redshift distributions; this highlights the potential of next-generation facilities to constrain BNS formation channels.

2 Binary neutron star population

Our population model builds upon the work of Farrow et al. (2019) and Galaudage et al. (2021). We emphasise that it is not directly informed by binary stellar evolution calculations, but is designed to capture the phenomenology relevant to the case studied in this work. In Galaudage et al. (2021), BNS systems are described as pairs composed of a recycled and a slowly spinning NS, each drawn from either a heavy (H) or a light (L) sub-population. Here, we further assumed that BNSs form preferentially as light and light or heavy and heavy, as appears to be the case for GW170817 and GW190425. Under this assumption, the individual NS masses, m₁ and m₂, associated with the heavy/light sub-populations are sampled from Gaussian distributions with means and standard deviations of (μ_1,H/L,σ_1,H/L) and (μ_2,H/L,σ_2,H/L), respectively. The numerical values of these parameters, adopted from Galaudage et al. (2021), are listed in Table 1.

In this formulation, no explicit ordering between m₁ and m₂ is imposed. To prevent labelling ambiguities that may arise in the GW context, we reparametrised the model in terms of the total mass m_t = m₁ + m₂ and the mass ratio q = min(m₁, m₂)/max(m₁, m₂) ≤ 1. The total mass is assumed to follow a normal distribution with a mean and standard deviation of μ_t,H/L = μ_1,H/L + μ_2,H/L and ${\sigma }_{t,\text{H/L}}=\sqrt{{\sigma }_{1,\text{H/L}}^2+{\sigma }_{2,\text{H/L}}^2}$, respectively. The mass ratio is modelled as a truncated lognormal distribution between 0 and 1, characterised by parameters μ_q,H/L. and σ_q,H/L. The parameters of the distribution were chosen to qualitatively reproduce the population of Galaudage et al. (2021), as we show that it is indeed the case in Appendix A. If we denote the fraction of light BNS with λ and use Λ = (λ, Λ_t,L, Λ_q,L, Λ_z,L, Λ_t,H, Λ_q,H, Λ_z,H) to denote the full set of population parameters, i.e. the hyperparameters, our population model reads $$\begin{array}{lll}p\left(m_t,q,z|\Lambda \right)\hfill & =\hfill & \lambda \times p\left(m_t|{\Lambda }_{t,\text{L}}\right)p\left(q|{\Lambda }_{q,\text{L}}\right)p\left(z|{\Lambda }_{z,\text{L}}\right)\hfill \\ \hfill & \hfill & +\left(1-\lambda \right)\times p\left(m_t|{\Lambda }_{t,\text{H}}\right)p\left(q|{\Lambda }_{q,\text{H}}\right)p\left(z|{\Lambda }_{z,\text{H}}\right).\hfill \end{array}$$(1)

The redshift distribution is given by $$p\left(z|{\Lambda }_{z,\text{H/L}}\right)\propto \frac{1}{1+z}\frac{\text{d}{V}_c}{\text{d}z}{\Re }_{\text{H/L}}\left(z\right).$$(2)

Distributions were normalised between z_min = 0 and z_max = 3, as we fond that almost no BNS event is detected past redshift 3 with ET. In the above equation, $\frac{\text{d}{V}_c}{\text{d}z}$ is the differential comoving volume, and ℜ_H/L(z) is the volumetric rate of events. The latter depends on the cosmic SFR, the BNS formation efficiency, η_H/L, and the time-delay distribution, p_H/L(Δt): $${\Re }_{\text{H/L}}\left(z\right)={\displaystyle \int {\eta }_{\text{H/L}}}\left(t\left(z\right)+\Delta t\right)\text{SFR}\left(t\left(z\right)+\Delta t\right)p_{\text{H/L}}\left(\Delta t\right)\text{d}\Delta t.$$(3)

We considered two distinct scenarios for the time-delay distribution. First, we considered the case of a single BNS population, where we assumed a unique time-delay distribution that follows a power law with index α, bounded between Δt_min = 10 Myr and Δt_max = 13 Gyr. Although our mass model assumed a bimodal distribution to capture the high-mass end, our conclusions on the reconstruction of the time-delay distribution are largely insensitive to the precise form of the mass model. For instance, adopting a skewed distribution with an extended high-mass tail would yield comparable results, provided that the abundance of heavy systems is similar. For simplicity, we therefore adopted a single model for the mass distribution, namely the one described above.

In the second scenario, following Romero-Shaw et al. (2020) and Galaudage et al. (2021), we assumed that the heavy subpopulation merges on systematically shorter timescales than the light one. In this case, heavy binaries were assigned a fixed, short delay time of Δt_Η = 30 Myr, modelled as p_Η(Δt) = δ(Δt – Δt_Η), while light binaries were assumed to follow a power-law delay-time distribution with index α_L, bounded between Δt_L,min = 10 Myr and Δt_L,max = 13 Gyr. In both cases we fond that values of the lower cutoff of the power-law time-delay distribution, Δt_L,min and/or the fixed delay of the heavy population, Δt_Η, below ~ 100 Myr have a negligible impact on the resulting red-shift distributions. Hence, the results presented below are largely insensitive to the precise choice of these parameters.

We adopt the cosmic SFR from Madau & Dickinson (2014). The formation efficiency of each sub-population was assumed to be constant with redshift and thus acted as a normalisation factor. The overall normalisation was constrained by the total BNS merger rate inferred by the LVK Collaboration (The LIGO Scientific Collaboration 2025), while the relative contributions of each sub-population were set by the choice of the mixing fraction λ. We emphasise that the mixing fraction entering the probability density function, λ, differs from the one entering the volumetric rate, λ_ℜ. If we define the latter such that ℜ(z) = λ_ℜℜ_L(z) + (1 − λ_ℜ)ℜ_H(z), the two are related by $${\lambda }_{\Re }=\frac{\lambda {\displaystyle {\int }_{z_{\text{min}}}^{z_{\text{max}}}\frac{{\Re }_{\text{H}}\left(z\right)}{1+z}}\frac{\text{d}{V}_c}{\text{d}z}\text{d}z}{{\displaystyle {\int }_{z_{\text{min}}}^{z_{\text{max}}}\frac{\lambda {\Re }_{\text{L}}\left(z\right)+\left(1-\lambda \right){\Re }_{\text{H}}\left(z\right)}{1+z}\frac{\text{d}{V}_c}{\text{d}z}\text{d}z}}.$$(4)

Here, we assumed an equal contribution from both populations, i.e. λ = 0.5. We summarise the parameters of our population in Table 1. For completeness, we also report there the corresponding values of λ_ℜ. Our values are compatible with the estimates of Galaudage et al. (2021) that 8–79% of BNS at birth are rapidly merging (in terms of λ_ℜ).

Since the time-delay distribution is poorly constrained, we explored different assumptions for the delay-time distribution of the slowly merging (light) population. Specifically, we fixed the minimum and maximum values of the time-delay range and vary the power-law index, both in the single population and distinct sub-population cases. Population synthesis models typically find α = −1 (e.g. Chruslinska et al. 2018), although at specific progenitor metallicities the distribution can be much shallower (de Sá et al. 2024; Pellouin et al. 2025), while observational-based estimates that use gamma-ray bursts or Galactic pulsars find much steeper values of α = −1.5 to −2 for some sub-populations (Zevin et al. 2022; Maoz & Nakar 2025). In this work we considered α = −0.5, −1, and −1.5. The corresponding volumetric merger rates are shown in Fig. 1.

In each case, the total rate is normalised to the mean of the 90% credible interval on the local BNS merger-rate range reported by the LVK (The LIGO Scientific Collaboration 2025), between 7.6 and 250 Gpc⁻³ yr⁻¹, i.e. 129 Gpc⁻³ yr⁻¹. Using the value of 1.16 × 10⁻² Mpc⁻³ for the number density of Milky Way-equivalent galaxy (MWEG) (Kopparapu et al. 2008), this rate corresponds to a merger rate per MWEG of 0.65–22 MWEG⁻¹ Myr⁻¹. The grey band at z = 0 shows the credible interval reported by the LVK. The short-delay population is rescaled by a different normalisation factor for each model, since for a fixed λ the corresponding volumetric-rate mixing fraction λ_ℜ depends on α_L. Larger values of assign greater weight to long delays, thereby enhancing the contrast with the short-delay population.

3 Mock catalogues and data analysis

We considered a triangular configuration for ET with a 10 km arm length and adopt the ET-10 (10 km-xylophone) sensitivity curve (Danilishin & Zhang 2023) implemented in GWBENCH (Borhanian 2021). Our mock detection and parameter estimation pipeline was based on the approach outlined in Fishbach et al. (2020) and Farah et al. (2023). Differences arising in our implementation are discussed in Appendix B.

Figure 2 shows the cumulative number of detections as a function of signal-to-noise ratio (S/N) threshold for the three values of α in the single population case. The shaded regions represent the uncertainty associated with the local BNS merger rate reported by the LVK (The LIGO Scientific Collaboration 2025), corresponding to the lower and upper bounds of the 90% credible interval. If we assume a S/N threshold of eight in the single population scenario, we estimate that ET will detect between approximately 680–22310, 1710–56280, and 3090–101630 BNS mergers per year for α = −0.5, −1, and −1.5, respectively. In the case of two sub-populations, we estimate these numbers to be 1180–38 920, 2370–77 950, and 3320–109150 BNS mergers per year for α_L, = −0.5, −1, and −1.5, respectively. In all cases, ~93% of detected BNS are within redshift 2. These estimates are compatible with Branchesi et al. (2023) after accounting for the fact that they assumed a local merger rate of 250 Gpc⁻³ yr⁻¹, which lies at the higher end of our range.

Based on these projections, we generated mock catalogues containing 100, 500, 1000, and 5000 detected events, with ten realisations for each case. We then performed hierarchical Bayesian inference to recover the population hyperparameters. These catalogue sizes are consistent with the expected ET detection rate, towards its lower end, while remaining computationally tractable. Current computational resources, in terms of both processing time and memory requirements, prevent us from extending the analysis to the tens of thousands of detections that ET is expected to observe. For this reason, we restricted our study to catalogues of up to 5000 events. The hierarchical Bayesian framework is described in detail in Appendix C; here, we summarise the population modelling choices relevant to this work.

Our population prior is given by a mixture model as in Eq. (1). The component distributions in total mass and mass ratio are modelled as Gaussian and truncated log-normal distributions, respectively. We adopted flat priors on the mixture fraction, as well as on the means and standard deviations of these distributions. For the redshift dependence, we explore two complementary approaches.

We assumed that the redshift evolution of each subpopulation follows Eqs. (2) and (3). In the inference, we modelled the delay-time distribution of both sub-populations as a power law, i.e. we do not assume any of them to have a fixed delay time. The delta-function limit is recovered when Δt_max Δt_min. The priors adopted are:

• Δt_min: log-flat between 1 Myr and 1 Gyr;

• Δt_max– Δt_min: log-flat between 10⁻³ Myr and 13.5 Gyr²;

• α: flat between −3 and 0.

When investigating the ability of ET to identify distinct BNS sub-populations, after performing the hierarchical Bayesian inference, we computed the probability that the mass distribution is bimodal, p_bimodal, as follows. For each hyperposterior sample, we drww 2000 (m_t, q) samples and apply Hartigan’s dip test (Hartigan & Hartigan 1985) using the implementation of Urlus (2025). This test quantifies whether the empirical probability distribution function of the samples exhibits a ‘dip,’ which would indicate the presence of two distinct modes. The test returns a p value that can be interpreted as the probability that the distribution is unimodal. Repeating this for all hyperposterior samples yields a distribution of p values for each hierarchical Bayesian analysis. We then computed p_bimodai as the fraction of p values below 0.05, which provides an estimate of the probability that the mass distribution is bimodal with more than 95% confidence. We verified that the number of (m_t, q) samples drawn per hyperparameter sample was sufficient, in the sense that increasing it does not affect the results. Since Hartigan’s dip test is defined for one-dimensional data, we rotated the (m_t, q) samples for each hyperparameter sample using a rotation matrix, and retained the rotation that yields the lowest p value (i.e. the highest probability of bimodality). This procedure allowed us to identify which linear combination of m_t and q exhibits the strongest bimodality. In practice, we find that the angle that maximises the bimodality was approximately zero, which indicates that the total mass distribution itself is the main source of bimodality.

The probability that the two redshift distributions are different, p_dist, was estimated as follows. We randomly drew pairs of hyperparameters from the posterior samples and computed the Kolmogorov-Smirnov (KS) statistic between the resulting redshift distributions of the two sub-populations. This ‘between-population’ KS distribution captures the typical differences between the two sets of distributions. To construct a background distribution, we also computed two ‘within-population’ KS distributions by comparing different samples from the same subpopulation, which provides a reference for variations expected from measurement uncertainty (detector noise and a finite number of events). We then determined the 5% quantile of the between-population KS distribution and find its corresponding quantiles in the within-population KS distributions. The maximum of these two quantiles is defined as p_dist. It measures whether the difference between the sub-populations exceeds the typical measurement-induced variations in at least one subpopulation, indicating that the distributions can be reliably distinguished.

Fig. 1 Merger rate of the long delays population (dashed lines), short delays population (dotted lines), and total population (full line) for three different hypotheses of the time-delays distribution. The rate is normalised at z = 0 to the mean of the interval reported by the LVK following GWTC-4 (The LIGO Scientific Collaboration 2025). The grey band at z = 0 shows the 90% credible interval reported by the LVK. In the case of a single population, the merger rate follows the long-delay curve, with the appropriate normalisation to match the total event rate.

Fig. 2 Number of events above a given S/N threshold per year in the single population scenario for the three values of α considered in this work.

4 Results

4.1 Single population

We first discuss the results in the case of a common redshift distribution, focusing on how well the parameters of the time-delay distribution can be constrained. Figure 3 shows the uncertainty on the power-law index, α, as a function of the number of events.

We find that, with 1000 observations, α can already be constrained within ~0.5 for an injected value of α = −0.5, while a similar accuracy is achieved with ~5000 events for α = −1.0.

For α = −1.5, however, the power-law index remains essentially unconstrained even after 5000 events. In this case, the redshift distribution closely resembles that obtained for a fixed delay, as shown in the right panel of Fig. 1. This behaviour is reflected in the large uncertainty on the maximum time delay, shown in Fig. 4, which remains uncertain at the level of several Gyr even after 5000 observations.

In contrast, for the other scenarios, the maximum delay time is very well constrained, within a few Myr, as it strongly affects the number of sources at z ~ 0, where the observational constraints are tightest. The minimum time delay is constrained to within a few hundred Myr in all cases, as shown in Fig. 5, consistent with our observation in Sect. 2 that varying Δt_min by a few Myr has little impact on the results.

Fig. 3 Credible intervals on the width of the 90% credible interval for the power-law exponent, as a function of the number of events and for different values of α.

Fig. 4 Credible intervals on the width of the 90% credible interval for the maximum time delay, as a function of the number of events and for different values of α.

Fig. 5 Credible intervals on the width of the 90% credible interval for the minimum time delay, as a function of the number of events and for different values of α.

Fig. 6 Credible intervals on the confidence that the mass distribution is bimodal as a function of the number of events and for the different values of αL. The horizontal red line corresponds to a probability of 0.95.

4.2 Sub-populations

We start by discussing the recovery of the mass distribution. With 100 events, we can already confidently measure that λ ≠ 0 and λ ≠ 1, which indicates that the population cannot be described by a single component. In Appendix D we show examples of how the total mass and mass ratio distributions are recovered as the number of events increases. Figure 6 displays the 90% confidence interval of p_bimodal (centred on the median) as a function of the number of events. With 500 detections, we can already establish in most cases that the total mass distribution is bimodal, and with 1000 detections this can be determined unambiguously, regardless of α_L.

Next, we turn to the distinction between the two redshift distributions. The upper panel of Fig. 7 shows representative examples of the reconstructed redshift distributions for α_L = − 1 as the number of events increases. The lower panel shows the distribution of KS statistics, both between-population and within-population. The value above the bottom panel indicates p_dist. This figure provides a visual interpretation of our method for estimating the probability that the two redshift distributions differ: as the confidence bands overlap less, the between-population KS distribution is moved towards higher values than the within-population ones, and the computed p_dist increases. We observe that the redshift distribution of the heavy population is typically better constrained than that of the light population, as a consequence of their larger S/N.

Our results for p_dist are shown in Fig. 8. For α_L = −0.5, we can already establish that the two redshift distributions differ with more than 95% confidence using 500 events. For α_L = −1, this level of significance is reached in over 95% of realisations with 5000 events. For α_L = −1.5, the distinction becomes more challenging, and it appears that many more events are needed in order to distinguish the distributions, in agreement with the results found in the single population case.

5 Conclusions

Despite major progress over the last decade, the formation pathways of BNSs remain uncertain. In particular, the number of mass transfer episodes and the criterion for their stability have a strong influence on the properties of the resulting population yet are still not fully understood. The redshift distribution of BNS mergers provides an observational handle on BNS formation processes, since it encodes the distribution of time delays between binary formation and coalescence and therefore the initial BNS orbital separation. Moreover, BNS masses could also be correlated with the formation channel. The detection of GW190425 by LVK highlighted the possibility that a significant fraction of BNSs may form through unstable ‘case BB’ mass transfer, resulting in much shorter delay times than the rest of the population (Romero-Shaw et al. 2020; Galaudage et al. 2021). In this scenario, two BNS sub-populations with distinct redshift evolutions could coexist. In this paper, we explore the ability of future ground-based facilities such as the ET to characterise the time-delay distribution of BNSs and to probe the existence of distinct sub-populations.

Using the latest estimate of the BNS merger rate by the LVK (The LIGO Scientific Collaboration 2025), we generated mock BNS populations consisting of an equal mixture of light and heavy sub-populations, inspired by the parametric models of Farrow et al. (2019) and Galaudage et al. (2021). Assuming that the redshift distribution of each sub-population is given by the convolution of the SFR with a time-delay distribution, we considered two scenarios. In the first, both light and heavy BNSs share a common time-delay distribution described by a power law, with index α_L = −0.5, −1, and −1.5. In the second, heavy BNSs are assigned a fixed, short time delay, while light BNSs follow a power-law distribution with the same set of indices. Independently of this choice, we find that ET should observe at least a few thousand BNS mergers over its lifetime, with estimated detection rates ranging from 1710 to 56 280 events per year in the fiducial α_L = −1 case for a common redshift distribution.

We then considered multiple realisations ofcatalogues ranging from 100 up to 5000 events and performed hierarchical Bayesian analyses to reconstruct the astrophysical distributions. Crucially, when considering the possibility of having two distinct sub-populations, we did not constrain either sub-population to a fixed delay time in the inference, allowing for a more general and flexible model. We find that for α = −0.5 and −1, the power-law index can be constrained within ~0.5 after 1000 and 5000 observations, respectively. The minimum delay time is constrained within a few hundred megayears, while the maximum delay time can be determined to within a few megayears. For α = −1.5, however, the parameters remain poorly constrained, as the resulting redshift distribution closely resembles that obtained for a fixed delay of ~100 Myr.

In the case of two sub-populations, for all α_L values, we find that with 500 events we can identify the bimodality of the BNS total mass distribution in most cases, and with 1000 events we can do so in practically all cases. Moreover, a few thousand events are sufficient to disentangle the redshift distributions of the two sub-populations if α_L = −0.5 or −1. For α_L = −1.5, the difference between the two distributions is very small, and more observations would be required to distinguish them, if at all. We could not explore this further due to limited computational resources. This highlights a crucial point: current population analysis frameworks will face significant challenges when the number of observations reaches several thousand.

Our study shows that next-generation ground-based detectors will provide information on the processes leading to BNS coalescence, and could reveal a sub-population of heavy BNSs merging on much shorter timescales than the Galactic population. The most likely scenario behind it would be unstable ‘case BB’ mass transfer, which tightens the binary orbit, allowing the system to remain bound after the second supernova and resulting in shorter merger delays. Detecting such a population would provide insight into BNS formation channels, accessible only through GW observations, since rapidly merging heavy BNSs would be largely absent in the Milky Way, explaining the lack of Galactic binaries as massive as GW190425. We emphasise once again that, although this scenario provides a viable explanation for the discrepancy between GW190425 and the Galactic population, it is not the only one capable of doing so. Alternatively, if heavy BNSs exist in our Galaxy but are radio-invisible (Safarzadeh et al. 2020; Chu et al. 2025), the Laser Interferometer Space Antenna (LISA) should enable them to be found in the Milky Way (Korol & Safarzadeh 2021).

An important caveat of our analysis is that it was limited to two sub-populations. For instance, allowing for a light and heavy coupling would introduce a third sub-population, although the (currently limited) data tend to favour mostly heavy and heavy and light and light coupling. Moreover, several formation pathways may contribute to the BNS population (Tauris et al. 2017; Vigna-Gómez et al. 2018), potentially resulting in a more complex distribution. Neglecting this could severely bias population reconstruction when using parametric or astrophysical models for the population prior (Zevin et al. 2021; Toubiana et al. 2021; Cheng et al. 2023). Future work will explore more comprehensive analyses that combine population synthesis models (Pellouin et al. 2025) with multi-dimensional non-parametric approaches, such as that of Tenorio et al. (2025).

Fig. 7 Upper panel: reconstructions of the redshift distributions for increasing catalogue sizes. Coloured bands show the 90% credible intervals for the light (red) and heavy (golden) populations. Solid lines indicate the median reconstructions, while dashed lines mark the true distributions. For the heavy population, the dashed and full lines superimpose almost perfectly. Lower panel: distribution of within-population and between-population KS statistics. The value on top shows the probability that the two heavy and light population have different redshift distributions.

Fig. 8 Credible intervals on the confidence with which we can determine that the redshift distribution of the light and of the heavy population are different as a function of the number of events and for the different values of αL. The horizontal red line corresponds to a probability of 0.95.

Acknowledgements

We are thankful to S. Borhanian, T. Bruel and M. Quartin for fruitful discussions. A.T. is supported by MUR Young Researchers Grant No. SOE2024-0000125, ERC Starting Grant No. 945155-GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, Italian-French University (UIF/UFI) Grant No. 2025-C3-386, MUR Grant “Progetto Dipartimenti di Eccellenza 2023-2027” (BiCoQ), and the ICSC National Research Centre funded by NextGenerationEU.

References

Appendix A Comparison between population parametrisations

Figures A.1 and A.2 compare our population model with that of Galaudage et al. (2021), with the extra condition that only light/light and heavy/heavy pairings are allowed. The total mass distributions are identical, as both models describe them as the sum of two Gaussians. The corresponding mass-ratio distributions are qualitatively similar. The discontinuity visible in the (m_t, q) correlation for the light population originates from the fact that in Galaudage et al. (2021), the component masses are not ordered; enforcing q ≤ 1 thus introduces this artificial feature. Still, our model has a similar support to the original distribution.

Fig. A.1 Comparison between the parametrisation of Galaudage et al. (2021) for the light population and ours in terms of total mass and mass ratio.

Fig. A.2 Comparison between the parametrisation of Galaudage et al. (2021) for the heavy population of ours in terms of total mass and mass ratio.

Appendix B Mock detection and parameter estimation

In their appendices, Fishbach et al. (2020) and Farah et al. (2023) outline a prescription for generating mock gravitational-wave (GW) catalogues that self-consistently incorporate parameter estimation uncertainty and selection effects, while also reproducing the expected scaling of measurement errors with the inverse of the S/N. Their approach leverages the fact that the S/N, ρ, can be written as $$\rho =\frac{D\left({\theta }^{\prime }\right)}{D_L},$$(B.1)

where θ’ denotes all binary parameters entering the S/N computation except for the luminosity distance D_L (detector frame masses, spins, sky location…). For fixed θ’, this expression allows one to transform straightforwardly between D_L and ρ.

The key idea is to treat the S/N itself as a sampling parameter during mock parameter estimation, rather than the luminosity distance. This makes it possible to enforce that the measured S/N, $\tilde{\rho }$, is normally distributed around the true S/N, with unit variance: $\tilde{\rho }~𝒩\left(\rho ,1\right)$, which reflects the expected behaviour of a matched-filter S/N in the case of a single detector (Finn & Chernoff 1993). Selection effects are then included as follows: for each source, the true S/N is computed and scattered according to the above distribution to obtain $\tilde{\rho }$; the event is retained only if $\tilde{\rho }$ exceeds a chosen S/N threshold, ρ_th.

For retained events, posterior samples of the S/N, $\widehat{\rho }$, are then drawn from a normal distribution centred at $\tilde{\rho }$ with unit variance. To model measurement uncertainty on the remaining parameters θ’, “noisy” values ${\tilde{\theta }}^{\prime }$ are generated by sampling from a normal distribution with variance proportional to ${\tilde{\rho }}^{-1}$, mimicking the typical scaling of parameter errors with inverse S/N. Posterior samples ${\tilde{\theta }}^{\prime }$ are obtained by further drawing from normal distributions centred at ${\tilde{\theta }}^{\prime }$ with variances also proportional to ${\tilde{\rho }}^{-1}$. The $\left\{{\tilde{\theta }}^{\prime },\widehat{\rho }\right\}$ samples are then transformed back into $\left\{{\tilde{\theta }}^{\prime },{\widehat{D}}_L\right\}$ samples via Eq. (B.1). The proportionality constants for the measurement errors on θ^′ were calibrated against mock parameter estimation runs.

The problem can be simplified by noting that the S/N may be expressed as $$\rho ={\rho }_{\text{opt}}\Theta ,$$(B.2)

where ρ_opt is the optimal S/N for a source located directly overhead the detector with face-on inclination, and Θ encodes the dependence on sky location, polarisation, and inclination. For a single detector, Θ is bounded by 0 ≤ Θ ≤ 1 , while for ET the range becomes $0\le {\text{Θ}}_{\text{ET}}\le \sqrt{3}$. Denoting by F_+,i and F_×,i the antenna pattern functions of the three equivalent ET detectors, as given in Regimbau et al. (2012), the factor Θ_ΕΤ can be written as $${\Theta }_{\text{ET}}=\sqrt{{\displaystyle \sum _{i=1}^3\left(F_{+,i}^2+F_{\times ,i}^2\right).}}$$(B.3)

This expression generalises the single-detector relation derived in Finn & Chernoff (1993), but is strictly valid only for the case of a triangular ET configuration, where all three detectors are co-located. We note that in a network of multiple detectors, the distribution of the measured S/N, $\tilde{\rho }$, deviates slightly from the simple normal form assumed in the single-detector case (Essick 2023; Gerosa & Bellotti 2024). However, for the purposes of this study, we neglect this correction and adopt the single-detector approximation for simplicity.

Our approach differs from that of Fishbach et al. (2020); Farah et al. (2023) in two main respects. First, although we enforce hard prior boundaries when drawing the posterior samples ${\widehat{\theta }}^{\prime }$, we do not impose such constraints when generating the “noisy” parameter values ${\tilde{\theta }}^{\prime }$. This choice reflects the fact that noise can make it appear as the likelihood peaks outside of the prior limits, which manifests as railing of the posterior against the prior boundary. This is also simpler, as it means that the likelihood is truly a Gaussian and not a truncated Gaussian, avoiding the need to adjust the sampling procedure to account for a modified functional form.

Second, Fishbach et al. (2020); Farah et al. (2023) adopt θ′ = (ℳ_c,d, η, Θ) as sampling parameters, where ℳ_c,d is the detector-frame chirp mass and η is the symmetric mass ratio. However, we find that sampling in the symmetric mass ratio leads to large variances in the Monte Carlo estimators used when evaluating the population likelihood (see App. D). This issue arises because our underlying population model favours nearly equal-mass systems, while sampling in η tends to underpopulate the region near q ~ 1, as also discussed in Farah (2022). To mitigate this, we sample directly in the mass ratio q instead of η. We assume that uncertainties in the inferred mass ratio scale as 0.15ρ_th/ρ, which is a conservative estimate based on the measurement uncertainties reported for GW170817 (Abbott et al. 2017) and GW190425 (Abbott et al. 2020), particularly when tight spin priors are used. Although the mapping from detector-frame chirp mass to total mass is less problematic, we likewise choose to sample directly in the detector-frame total mass, adopting a similarly conservative uncertainty scaling of 0.1ρ_th/ρ. For Θ we use the same as Fishbach et al. (2020) and Farah et al. (2023): 0.21ρ_th/ρ.

We compute S/Ns using the ~~GWBENCH~~ package (Borhanian 2021), setting the spins to 0 and using the IMRPhenomXHM waveform (García-Quirós et al. 2020). For the S/N threshold, we take ρ_th = 8.

Appendix C Hierarchical Bayesian analysis

We denote with θ the set of parameters describing a GW event, p(d|θ) the single event likelihood and p(θ|Λ) the population prior on θ, which depends on hyperparameters A that we wish to infer. The population likelihood for observing N_obs events $\left\{d\right\}=\left\{d_1,\dots ,d_{N_{\text{obs}}}\right\}$, marginalised over the event rate is Mandel et al. (2019); Vitale et al. (2022) $$p\left(\left\{d\right\}|\Lambda \right)={\displaystyle \prod _i^{N_{\text{obs}}}{\displaystyle \int \frac{p\left(d_i|\theta \right)p\left(\theta |\Lambda \right)}{p_{\text{det}}\left(\Lambda \right)}\text{d}\theta .}}$$(C.1)

We have introduced the selection function defined via $$p_{\text{det}}\left(\Lambda \right)={\displaystyle \int p_{\text{det}}}\left(\theta \right)p\left(\theta |\Lambda \right)\text{d}\theta ,$$(C.2) $$p_{\text{det}}\left(\theta \right)={\displaystyle {\int }_{d>\text{threshold}}p\left(d|\theta \right)\text{d}d,}$$(C.3)

where the second integral is restricted to realisations d that exceed the detection threshold (defined via the chosen ranking statistic), here the measured matched filter S/N $\widehat{\rho }$.

Alternatively, if we do not wish to marginalise over the rate, we can write $$p\left(\left\{d\right\}|\Lambda \right)=e^{-N\left(\Lambda \right)p_{\text{det}}\left(\Lambda \right)}{\displaystyle \prod _i^{N_{\text{obs}}}{\displaystyle \int p\left(d_i|\theta \right)}\frac{\text{d}N}{\text{d}\theta }}\left(\Lambda \right)\text{d}\theta .$$(C.4)

where $\frac{\text{d}N}{\text{d}\theta }$ is the differential number of events, ${\displaystyle \int \frac{\text{d}N}{\text{d}\theta }\left(\Lambda \right)\text{d}\theta }$ is the total expected number of events, and we have the relation $\frac{\text{d}N}{\text{d}\theta }-N\left(\Lambda \right)p\left(\theta |\Lambda \right)$.

In practice, the single-event likelihood can be written in terms of the posterior and the parameter estimation prior via Bayes’ theorem. The integral over θ is evaluated using Monte Carlo integration, based on posterior samples obtained as described in App. B.

Similarly, the selection function is estimated via an injection campaign in which events are drawn from an injection distribution π_inj(θ), and only those satisfying the detection criterion are retained, forming the set θ_det. The detection probability is then computed using importance sampling: $$p_{\text{det}}\left(\Lambda \right)={\displaystyle \sum _{{\theta }_{\text{det}}}\frac{p\left({\theta }_{\text{det}}|\Lambda \right)}{{\pi }_{\text{inj}}\left({\theta }_{\text{det}}\right)}}.$$(C.5)

As discussed in Talbot & Golomb (2023), the finite number of samples used to evaluate the integrals in population analyses can introduce significant variance in the Monte Carlo estimators, which may in turn lead to biases. To mitigate this issue in the estimation of the selection function, we choose the injection prior to match the true population, thereby reducing the variance in the importance sampling weights. Additionally, in line with the recommendations of Talbot & Golomb (2023), we impose a threshold on the total variance of the log-likelihood and discard hyperparameters exceeding this limit.

The posterior on A is then obtained through Bayes’ theorem: p(Λ|{d}) ∝ p({d}|Λ)π(Λ), with the priors described in Sec. 2 and in the next one. The sampling is performed with the Eryn sampler³ (Karnesis et al. 2023).

Appendix D Reconstruction of the mass distribution

Figure D.1 shows representative examples of how the measurement of the total mass and the mass ratio distribution improves as the number of events increases. These plots correspond to realisations of the α_L = −1 case; however, the choice of α_L has little impact on this result.

Fig. D.1 Reconstructions of the total mass and mass ratio distributions for increasing catalogue sizes. Coloured bands show the 90% credible intervals for the light (red) and heavy (golden) populations. Solid lines indicate the median reconstructions, while dashed lines mark the true distributions.

All Tables

All Figures

Fig. 1 Merger rate of the long delays population (dashed lines), short delays population (dotted lines), and total population (full line) for three different hypotheses of the time-delays distribution. The rate is normalised at z = 0 to the mean of the interval reported by the LVK following GWTC-4 (The LIGO Scientific Collaboration 2025). The grey band at z = 0 shows the 90% credible interval reported by the LVK. In the case of a single population, the merger rate follows the long-delay curve, with the appropriate normalisation to match the total event rate.

In the text

	Fig. 2 Number of events above a given S/N threshold per year in the single population scenario for the three values of α considered in this work.
In the text

	Fig. 3 Credible intervals on the width of the 90% credible interval for the power-law exponent, as a function of the number of events and for different values of α.
In the text

	Fig. 4 Credible intervals on the width of the 90% credible interval for the maximum time delay, as a function of the number of events and for different values of α.
In the text

	Fig. 5 Credible intervals on the width of the 90% credible interval for the minimum time delay, as a function of the number of events and for different values of α.
In the text

	Fig. 6 Credible intervals on the confidence that the mass distribution is bimodal as a function of the number of events and for the different values of αL. The horizontal red line corresponds to a probability of 0.95.
In the text

Fig. 7 Upper panel: reconstructions of the redshift distributions for increasing catalogue sizes. Coloured bands show the 90% credible intervals for the light (red) and heavy (golden) populations. Solid lines indicate the median reconstructions, while dashed lines mark the true distributions. For the heavy population, the dashed and full lines superimpose almost perfectly. Lower panel: distribution of within-population and between-population KS statistics. The value on top shows the probability that the two heavy and light population have different redshift distributions.

In the text

	Fig. 8 Credible intervals on the confidence with which we can determine that the redshift distribution of the light and of the heavy population are different as a function of the number of events and for the different values of αL. The horizontal red line corresponds to a probability of 0.95.
In the text

	Fig. A.1 Comparison between the parametrisation of Galaudage et al. (2021) for the light population and ours in terms of total mass and mass ratio.
In the text

	Fig. A.2 Comparison between the parametrisation of Galaudage et al. (2021) for the heavy population of ours in terms of total mass and mass ratio.
In the text

	Fig. D.1 Reconstructions of the total mass and mass ratio distributions for increasing catalogue sizes. Coloured bands show the 90% credible intervals for the light (red) and heavy (golden) populations. Solid lines indicate the median reconstructions, while dashed lines mark the true distributions.
In the text