A Quiet Night in the Millihertz Sky? LISA Forecasts from LIGO/Virgo/KAGRA Populations

Introduction

The Laser Interferometer Space Antenna (LISA; Colpi et al. [29]) will observe the millihertz inspiral of binary black holes (BBHs). Some of these systems will merge in the frequency band of the ground-based LIGO–Virgo–KAGRA (LVK) detectors (Aasi et al. [1]; Acernese et al. [17]; Akutsu et al. [19]) within years to decades of their LISA observation, allowing the same binary to be observed by LISA during its inspiral and by ground-based detectors near merger (Sesana [79]; Gerosa et al. [47]).

At millihertz frequencies, measurements of orbital eccentricity can help distinguish how black hole binaries formed, before gravitational radiation further circularizes their orbits (Breivik et al. [26]; Nishizawa et al. [69]; Cutler et al. [32]). For systems observed in both bands, combining LISA’s early inspiral with ground-based observations near merger can improve measurements of binary parameters (Vitale [92]; Toubiana et al. [89]) and strengthen tests of general relativity (Toubiana et al. [90]; Datta et al. [33]). These possibilities motivate forecasts of both the number and masses of binary black holes accessible to LISA.

Forecasts based on the third Gravitational-Wave Transient Catalog (GWTC-3; Abbott et al. [14]) predict only a few individually resolved BBHs in a four-year mission at signal-to-noise ratio (SNR) $\rho> 8$ (Buscicchio et al. [28]). Searches without prior source information may require higher SNR thresholds; Moore et al. [68] forecast fewer than one multiband early warning within a ten-year merger window for most of their modeled populations. A small catalog limits population measurements, but massive BBHs whose inspirals lie in LISA’s sensitive band can be detected over a larger volume than lighter systems (Gerosa et al. [47]; Wagg et al. 2022). They can therefore make up a larger fraction of LISA detections than of the intrinsic population.

The populations inferred from GWTC-3 ([14]), GWTC-4.0 (Abac et al. 2026c), and GWTC-5.0 ([2]) contain several features in the source-frame primary mass distribution, including a peak near $10\,M_{\odot}$ and a change in slope around $35\,M_{\odot}$ ([2]). Other analyses find model-dependent structure at higher masses ([87] (Pierra & Papadopoulos 2026)), including tentative evidence for a subdominant component near $60 - 70\,M_{\odot}$ ([63 b] (Abbott et al. 2021)). Beyond features in the primary mass distribution, detectability will also depend on the secondary mass distribution and on how the black holes are paired ([40, 37]).

GW190521 and GW231123 probe the high-mass black hole population in and around the range where pair instability can suppress black-hole formation through stellar collapse (Farmer et al. 2019; Abbott et al. 2020a,b; Abac et al. 2025). Dedicated searches for intermediate-mass black hole (IMBH) binaries have also constrained the rates of these rare systems ([9, 11]). At intermediate masses, forecasts identify binaries that could be observed by both LISA and next-generation ground-based detectors ([66, 55, 43]).

We first use the released GWTC-3 and GWTC-5.0 population posteriors to update the expected LISA population and compare with earlier forecasts. We then consider alternative fits to GWTC-5.0, with two or three mass peaks, using either a primary-mass distribution with a conditional mass-ratio distribution or a common mass distribution for both black hole masses. A common underlying mass distribution is physically motivated when both black holes are drawn from the same population, as can occur in dynamically assembled binaries ([77, 20]). Analyses of GWTC-3 find that shared mass features in both components are consistent with the data ([37]). The pairing weight allows correlations between the component masses ([40]), as expected when dynamical interactions preferentially pair massive black holes ([88]). Comparing these models tests how the LISA forecasts depend on the assumed mass distributions and pairing prescriptions. To assess the contribution of binaries beyond the fitted mass range, we extrapolate each model from its upper primary-mass limit of $300\,M_{\odot}$ to $1000\,M_{\odot}$, since LISA sensitivity is greater in this range than at lower masses.

For all the LISA forecasts in this work, we adopt a fiducial, flat $\Lambda$CDM cosmology with $H_0 = 67.74\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$ and $\Omega_m = 0.3089$, following Planck Collaboration et al. [73]. § 2 describes the population models and LISA selection. § 3 presents the catalog baselines. § 3.2 and § 3.3 examine the alternative mass models and the less restrictive high-mass extrapolation. § 3.4 considers multiband sources, and § 4 discusses what the forecasts imply for the high-mass population and how they depend on the assumed mass distribution and merger-rate history.

Population Models and LISA

Baseline Populations

We use the released GWTC-3 and GWTC-5.0 population posteriors (LIGO Scientific Collaboration et al. 2024, 2026a) for our reference forecasts. The GWTC-5.0 BrokenPowerLaw+2Peaks model describes the primary mass with a broken power law and two Gaussian peaks, a smooth lower cutoff, and a fixed upper limit of $300\,M_{\odot}$ ([2]). Its inferred parameters include the power-law slopes and break location, the peak locations, widths and weights, and the low-mass cutoff and smoothing scale. The GWTC-3 PowerLaw+Peak model instead has one power-law slope and one Gaussian peak, and infers the upper mass cutoff ([84, 14]). In both models, the companion distribution at fixed primary mass follows a power law in $q = m_2/m_1$, with a low-mass taper and normalization over the allowed secondary masses.

Figure 1 shows the inferred mass distributions and merger-rate evolution for both GWTC-3 and GWTC-5.0. We assume that the mass distribution does not evolve with redshift ([2]). For each released catalog population, we use its inferred rate evolution, $\mathcal{R}(z) = \mathcal{R}_0(1+z)^{\kappa}$ ([41, 14, 2]). Each posterior sample contains one complete set of population parameters, including the mass and spin distributions, local merger rate, and redshift exponent. We use these parameters together when calculating a forecast. Appendix A defines the population densities.

Alternative Mass Models

We fit four mass models to a common set of 259 GWTC-5.0 binary black hole events, including GW190521 and GW231123 ([6 e,a]). The primary-mass models combine a broken power law (BPL) with two or three Gaussian peaks (2P or 3P), and are denoted BPL2P and BPL3P. The companion is described by a conditional mass-ratio distribution. The symmetric models, symBPL2P and symBPL3P, use the same underlying distribution for both component masses before pairing, with pairs weighted by $q^\beta$ [40, 37]. Writing the underlying mass distribution as $B_K$, the primary-mass and symmetric models give $p(m_1,q)=B_K(m_1)p(q|m_1)$ and $p(m_1,q)\propto m_1B_K(m_1)B_K(qm_1)q^\beta$, respectively. We express both models in $(m_1,q)$. For the symmetric models, changing variables from $(m_1,m_2)$ introduces the factor $m_1$. Here $K=2$ or $3$, $q=m_2/m_1$, and $m_1\ge m_2$; Appendix § 8.§ 2 defines the normalizations and low-mass tapers.

The two-peak fits describe the structure near 10 and 30 $M_\odot$. For the three-peak fits, we allow the additional peak’s mean to lie between 45 and 90 $M_\odot$ and fit its width and abundance along with the rest of the population. All four fits restrict component masses to 3–300 $M_\odot$, allow the location and width of the smooth low-mass cutoff to vary, and share a Gaussian effective-spin distribution and power-law merger-rate evolution (Appendix § 8.§ 3).

BPL2P uses the same mass-model family as the released LVK Brokenpowerlaw+2Peaks posterior, but we infer its parameters independently. We use the released posterior for the catalog baseline forecasts. The additional peak is motivated by the reported feature near 70 $M_\odot$ [63 b], and we are specifically interested in its addition to explore whether it could significantly contribute to the LISA detectable population. We note that an additional high mass peak has also been used for LVK cosmological analyses (Abac et al. 2026b). These fits compare mass models with two or three peaks and different prescriptions for pairing the component masses.

The observed number of mergers and the LVK detection efficiency determine the local merger rate for each fitted population. Populations that are harder to detect require a higher rate to explain the same catalog. We propagate both the uncertainty in this rate and its dependence on the mass, spin, and redshift parameters into the LISA forecasts [65, 85].

We fit both the primary-mass and symmetric models to GWTC-5.0 using hierarchical population inference, jointly inferring the mass, spin, and redshift parameters [65, 85]. We evaluate the population likelihood using samples from the event posteriors and account for selection effects using simulated signals recovered by the ground-based searches [39, 34] (LIGO Scientific Collaboration et al. 2026b). To limit numerical error from these finite samples, we require the estimated variance of the log likelihood to satisfy $\hat{\sigma}_{\ln\mathcal{L}}^2<1$ [83] (Abac et al. 2026c). Parameter combinations that fail this requirement are excluded from the fit. This also restricts which populations can be explored. Appendix § 8.§ 3 discusses this limitation.

Extrapolation to Higher Masses

For each fitted model, we extend the primary-mass range from 300 to 1000 $\mathrm{M}_\odot$ using the same mass dependence and pairing prescription. The merger rate assigned to binaries with $m_1 \le300\,\mathrm{M}_\odot$ is unchanged. Binaries with larger primary masses (while their companions may lie below or above 300 $\mathrm{M}_\odot$) add to that rate. Appendix § 8.§ 4 gives the corresponding rate expressions.

Forecasts using the fitted power-law merger-rate evolution are restricted to $0 \le z \le1.9$. The added massive binaries can be detectable at higher redshifts, so for the extrapolations we use a merger-rate history proportional to the star-formation history of Madau & Fragos [62], scaled to match the inferred local merger rate below 300 $\mathrm{M}_\odot$ at $z=0$ (Appendix § 8.§ 4). We neglect delays between binary formation and merger, which can change the merger-rate evolution with redshift (Fishbach & Kalogera [42]).

We also calculate forecasts restricted to $m_1 \le300\,\mathrm{M}_\odot$ with this same redshift evolution. Comparing these with the extrapolations to $m_1 \le1000\,\mathrm{M}_\odot$ isolates the contribution from binaries with $m_1 > 300\,\mathrm{M}_\odot$. We repeat these calculations across the inferred population parameters and local rates, so that the spread in their predictions gives the quoted credible intervals.

To allow a different abundance of binaries above 300 $\mathrm{M}_\odot$, we also consider a free-slope extrapolation of BPL2P. Its mass-distribution slope above 300 $\mathrm{M}_\odot$ varies independently of the fitted slope. We set $\mathrm{d}R/\mathrm{d}m_1 \propto m_1^{-\alpha}$ up to 1000 $\mathrm{M}_\odot$, with $\alpha\sim U[-2,8]$, and match the rate density at 300 $\mathrm{M}_\odot$ to the fitted model. This range allows both rising and falling high-mass distributions. We keep the same conditional mass-ratio prescription and star-formation-history evolution. We use the O3 IMBH search sensitivity ([11]) to weight the allowed slopes by the probability of no detections with primary mass above 300 $\mathrm{M}_\odot$ within its sensitivity coverage. This update changes the slope distribution conditional on each fitted population and rate, but it does not refit GWTC-5.0 (see Appendix § 12).

LISA Detections

We adopt a 4-yr observing interval and classify binaries with mission-integrated signal-to-noise ratio $\rho> 8$ as resolved. We also report forecasts for $\rho> 15$, motivated by Moore et al. [68], who estimated thresholds for searches without prior source information for binaries with component masses of 5–50 $\mathrm{M}_\odot$ and merger times up to 100 years after mission start. We use the higher threshold as a comparison for all our populations, including the higher-mass extrapolations. We do not calculate a search threshold specific to each mass or merger time.

We integrate the signal over the frequencies traversed during the mission, or until merger if it occurs first. We calculate the SNR using two effective channels. The noise model includes the instrumental contribution of Babak et al. [23] and the foreground from unresolved Galactic double white dwarfs of Robson et al. [76]. We average the response over sky position and polarization and integrate over isotropic inclinations. We assume quasi-circular inspirals and use the same noise spectrum for every population. Appendix § 9 specifies the noise spectrum and the conversion from binary strain to the two-channel SNR. For symBPL3P, which permits the largest numbers of massive binaries, we also compare the total binary black hole background predicted by this model with the assumed noise (Appendix § 9.§ 3). This tests whether neglecting that background could affect our forecasts.

The forecast includes binaries observable during the mission even if they merge long after it ends. For a specified population, we integrate over masses, redshifts, and merger times, weighting each contribution by the fraction of inclinations that exceed the SNR threshold. The mission duration sets the observed frequency interval ([24, 28]). Repeating the integral over the inferred population parameters and rates gives the distribution of expected counts. Appendix § 10 shows the integral and its limits.

The SNR estimation involves comparing the signal with the detector noise at each frequency. At leading order, the inspiral strain amplitude in the frequency domain scales as $|\tilde{h}(f)| \propto\mathcal{M}_{\mathrm{det}}^{5/6}f^{-7/6}/d_L$, where $\mathcal{M}_{\mathrm{det}}=(1+z)\mathcal{M}$ is the detector-frame chirp mass and $\mathcal{M}=(m_1m_2)^{3/5}/(m_1+m_2)^{1/5}$ is the source-frame chirp mass ([31]). Here $d_L$ is the luminosity distance and $f$ is the observed frequency. The accumulated SNR also depends on the portion of the inspiral observed and the detector noise over that interval. For inspirals within LISA’s sensitive band, the larger amplitude of massive binaries can increase their detectable volume. They can therefore make up a larger fraction of LISA detections than they are in the astrophysical merger population ([47, 93]).

Figure 2 compares four reference mass pairs over $0.1 \le z \le1.0$, observed during their final four years before merger. In these examples, larger masses produce stronger signals and shift the observed inspiral toward LISA’s most sensitive frequencies. The purple tracks show the strains and accumulated SNRs of BBHs that remain below our detection threshold over the mission.

For a GW231123-like binary with fixed source-frame masses of 137 + 101 $\mathrm{M}_\odot$, the final four years before merger give $\rho\simeq7.7$ at $z = 0.40$, using the angular averaging in Figure 2. At the lower and upper limits of the reported 90% redshift interval, $z = 0.15$ and $0.67$ ([4]), we obtain $\rho\simeq18$ and $5.2$, respectively. With these assumptions, $\rho= 8$ at $z \simeq0.38$, so that depending on the exact redshift a merger of this kind could be detectable by LISA. This comparison varies redshift alone and does not propagate the event’s correlated uncertainties in masses and orientation.

Results

In what follows, we report our expectations for LISA detections of binary black holes. We quote posterior medians with lower and upper errors spanning the central 90% credible interval (CI). For forecast counts, these intervals describe uncertainty in the expected number of detections from the population parameters and merger rate. An actual catalog has an additional Poisson fluctuation about that expectation, which is not included in the quoted intervals. Numerical integration error is assessed separately in Appendix § 10.

Baseline Population Forecasts

For the released GWTC-5.0 population, the median expected number of resolved binaries is $6.6^{+5.5}_{-2.2}$ at $\rho> 8$ in a 4-yr mission. Using the GWTC-3 population posterior, we obtain $4.3^{+3.6}_{-2.0}$. In the GWTC-5.0 forecast, binaries with $100\,M_\odot< m_1 \leq300\,M_\odot$ contribute $1.4^{+5.0}_{-1.1}$ expected detections, or $21^{+33}_{-15}\%$ of the total. Binaries with $m_1 \leq100\,M_\odot$ contribute $5.0^{+1.8}_{-1.4}$, compared with $4.3^{+3.6}_{-2.0}$ for GWTC-3. The higher GWTC-5.0 forecast is therefore driven mainly by binaries with primary masses between 100 and 300 $M_\odot$, while the forecasts below 100 $M_\odot$ remain similar. Buscicchio et al. [28] report a mean of about five detections for a four-year mission at $\rho> 8$. Appendix § 9.§ 2 specifies the SNR convention used here and its relation to that work.

After averaging over population uncertainty and Poisson count fluctuations, the released GWTC-5.0 population gives a probability of about 0.3% that LISA detects no binaries with $\rho> 8$ in four years. At the higher threshold $\rho> 15$, the expected counts fall to $0.91^{+0.64}_{-0.30}$ for GWTC-5.0 and $0.63^{+0.54}_{-0.29}$ for GWTC-3.

Alternative Mass Model Forecasts

Forecasts to 300 $M_\odot$

Within the fitted mass range and using the inferred power-law merger-rate evolution, BPL2P and BPL3P predict $5.7^{+4.9}_{-1.9}$ and $5.9^{+6.6}_{-1.9}$ detections at $\rho> 8$ in four years, respectively. symBPL2P and symBPL3P predict $6.7^{+12}_{-2.6}$ and $7.4^{+37}_{-3.2}$, with broader credible intervals (See Table 3).

Population	$m_{1,\max}/M_{\odot}$	Rate evolution	$\rho> 8$	$\rho> 15$
GWTC-3	Varies	Fitted power law	$4.3^{+3.6}_{-2.0}$	$0.63^{+0.54}_{-0.29}$
GWTC-5.0 (LVK)	300	Fitted power law	$6.6^{+5.5}_{-2.2}$	$0.91^{+0.64}_{-0.30}$
BPL2P	300	Fitted power law	$5.7^{+4.9}_{-1.9}$	$0.79^{+0.57}_{-0.25}$
BPL2P	300	Madau–Fragos	$5.6^{+5.3}_{-1.9}$	$0.78^{+0.59}_{-0.26}$
BPL2P	1000	Madau–Fragos	$5.9^{+9.5}_{-2.2}$	$0.84^{+1.5}_{-0.3}$
BPL3P	300	Fitted power law	$5.9^{+6.6}_{-1.9}$	$0.81^{+0.72}_{-0.26}$
BPL3P	300	Madau–Fragos	$5.8^{+7}_{-2}$	$0.8^{+0.74}_{-0.27}$
BPL3P	1000	Madau–Fragos	$6.2^{+5.3}_{-2.4}$	$0.89^{+5.3}_{-0.33}$
symBPL2P	300	Fitted power law	$6.7^{+12}_{-2.6}$	$0.89^{+1.2}_{-0.32}$
symBPL2P	300	Madau–Fragos	$6.5^{+12}_{-2.6}$	$0.88^{+1.2}_{-0.32}$
symBPL2P	1000	Madau–Fragos	$8^{+65}_{-3.9}$	$1.1^{+16}_{-0.56}$
symBPL3P	300	Fitted power law	$7.4^{+37}_{-3.2}$	$0.97^{+3.6}_{-0.38}$
symBPL3P	300	Madau–Fragos	$7.2^{+35}_{-3.2}$	$0.96^{+3.4}_{-0.38}$
symBPL3P	1000	Madau–Fragos	$14^{+3500}_{-9.4}$	$2.5^{+1500}_{-1.9}$
BPL2P (free slope)	1000	Madau–Fragos	$6.9^{+43}_{-3}$	$1.1^{+16}_{-0.49}$

The inferred local merger rates on the fitted mass range are $19.3^{+5.2}_{-4.2}$, $19.3^{+5.1}_{-4.1}$, $19.8^{+5.4}_{-4.2}$, and $20.0^{+5.3}_{-4.3}\,\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ for BPL2P, BPL3P, symBPL2P, and symBPL3P, respectively. The local rates are similar across the four models. Their LISA forecasts differ because the mass distributions allow different numbers of binaries with two massive components. In the symmetric mass models, the companion distribution depends on the same high-mass structure as the primary, as well as on the pairing weight, and it therefore may produce more LISA detections than the other mass models, which instead describe the companion through a conditional power law in mass ratio.

In BPL3P and symBPL3P, the fitted additional peak has mean $56^{+12}_{-9.1}\,M_\odot$ and $60^{+9.7}_{-11}\,M_\odot$, respectively. Within the fitted mass range, the third Gaussian accounts for $0.45^{+0.81}_{-0.39}$ and $0.52^{+0.65}_{-0.42}\%$ of primary black holes in the intrinsic population, respectively. Its location, width, and abundance are inferred together with the rest of the population. The count intervals overlap within the fitted mass range. The fits also change the continuum, pairing, and low-mass cutoffs, so the forecast differences cannot be attributed to the additional peak alone. The symmetric three-peak fit permits larger expected counts, with the upper end of its 90% credible interval exceeding $\sim40$ detections.

Extrapolation to 1000 $M_\odot$

The GWTC-5.0 population extends to primary masses above $100\,M_\odot$, with GW231123 providing a key observation at its high-mass end ([2, 4]). Beyond component masses of approximately $200\,M_\odot$, the predicted abundance increasingly depends on the assumed high-mass shape of the population model ([2]). For equal-mass binaries and an inclination-averaged response, the four-year $\rho= 8$ horizon rises from approximately $z = 0.1$ at source-frame total mass $M_{\rm tot} = 100\,M_\odot$ to $z = 4$ at $M_{\rm tot} = 1000\,M_\odot$ (Figure 7). For these horizon curves, we choose the gravitational-wave frequency at the start of the mission that gives the largest four-year SNR. The dependence of the SNR on time to merger is described by Gerosa et al. [47].

The extrapolated BPL2P and BPL3P models give $5.9^{+9.5}_{-2.2}$ and $6.2^{+24}_{-2.4}$ expected detections, respectively. Their median counts remain close to the forecasts restricted to 300 $M_\odot$, with the same rate evolution, although the uncertainty increases. Binaries with $m_1 > 300\,M_\odot$ contribute $0.22^{+4.5}_{-0.21}$ and $0.38^{+18}_{-0.38}$ expected detections, respectively.

The symBPL2P and symBPL3P extrapolations give a larger change, with $8^{+37}_{-3.9}$ and $14^{+3500}_{-9.4}$ expected detections, respectively. The high-mass slopes of the broken power law in the primary-mass distribution are $\alpha_2 = 5.5^{+1.5}_{-1.6}$ and $5.0^{+2.1}_{-2.3}$ for BPL2P and BPL3P. For the two- and three-peak symmetric models, the slopes of the common component-mass distribution before pairing are $2.7^{+1.2}_{-1.1}$ and $1.9^{+1.7}_{-2.1}$, respectively. These models favor shallower high-mass distributions; the three-peak model also allows a rising distribution within its 90% credible interval. Extrapolating those shapes can add many binaries with two massive components, to which LISA is particularly sensitive. This is especially important at the upper end of the symBPL3P forecast interval (Figure 3). Figure 4 shows how primary mass is related to mass ratio, redshift, and merger time among the detectable binaries.

At high primary mass, the 90% contours extend to lower mass ratios for the symmetric models, particularly symBPL3P. Higher-mass detections extend to larger redshifts in all four models, consistent with the mass dependence of LISA’s reach (Figure 7). We discuss the merger-time dependence in § 3.4.

The reported forecasts use the same instrumental noise and unresolved Galactic double-white-dwarf foreground. Adding the total binary black hole background predicted by symBPL3P to the noise reduces the expected counts for populations near the upper end of its forecast interval (Appendix § 9.§ 3).

Free High-Mass Slope

The free-slope calculation allows the abundance of binaries above 300 $M_\odot$ to differ from an extrapolation of the fitted power law. The two distributions meet at the same rate per unit primary mass at 300 $M_\odot$.

For BPL2P, the free-slope extrapolation gives $6.9^{+43}_{-?}$ expected detections at $\rho> 8$, compared with $5.9^{+9.5}_{-2.2}$ when the fitted slope is retained. The median changes little, while the upper end of the 90% credible interval increases. At $\rho> 15$, the free-slope forecast falls to $1.1^{+16}_{-0.49}$. The upper-left panel of Figure 3 compares the fitted-slope and free-slope mass distributions.

Before applying the O3 IMBH constraint, the added population predicts $(5.9^{+72}_{-5.5}) \times10^{-5}$ O3 detections within the mass and mass-ratio region covered by the nonspinning injections used here. Weighting the slopes by the probability of this non-detection (see Appendix § 12) therefore leaves the slope distribution and LISA forecast essentially unchanged. The count interval reflects the range of high-mass slopes we allow, together with uncertainty in the fitted population and merger rate.

The O3 sensitivity study includes binaries with total source-frame masses up to 800 $M_\odot$, but its sampled configurations have primary masses no larger than 400 $M_\odot$ (Abbott et al. [11]). Here we use its nonspinning injections as a reference calculation. For the free-slope model, $70^{+28}_{-44}\%$ of the added LISA detections with $m_1 > 300\,M_\odot$ and $\rho> 8$ lie outside the mass and mass-ratio region covered by these injections. This fraction accounts for the predicted primary masses and mass ratios, and it does not include spin-dependent sensitivity.

Figure 5 summarizes the expected counts for the released GWTC-3 and GWTC-5.0 populations, the four alternative mass models, and the free-slope extrapolation. For each fitted mass model, it compares the extrapolations with the corresponding 300 $M_\odot$ calculations using the same merger-rate evolution with redshift.

Multiband Implications

The time to coalescence sets when a LISA binary could be observed near merger by ground-based detectors. We count binaries with $t_{\mathrm{c,det}} \leq T_{\mathrm{obs}} + \Delta t$, where $t_{\mathrm{c,det}}$ is the observer-frame time to merger from mission start, $T_{\mathrm{obs}} = 4$ yr is the mission duration, and $\Delta t$ is the time after it ends. We report mergers during the mission and within 10 or 30 years after its end. We first report merger-time fractions, then assess ground-based detectability with a reference Cosmic Explorer (CE) configuration ([75, 36]). Within the 30-year window, the fraction of LISA detections that merge is $55^{+19}_{-9.8}\%$ when assuming the baseline GWTC-5.0 population.

For extrapolations to 1000 $M_\odot$, the fractions merging during the mission or within 30 years after its end are $54^{+26}_{-11}\%$ for BPL2P and $56^{+35}_{-13}\%$ for BPL3P. They increase to $66^{+30}_{-22}\%$ and $80^{+20}_{-36}\%$ for symBPL2P and symBPL3P, respectively. Allowing a free high-mass slope in BPL2P gives a fraction of $61^{+34}_{-16}\%$. Figure 6 compares the merger fractions during the mission and through 10 and 30 years after its end. The broad intervals, especially for the extrapolated models, reflect the uncertain contribution of massive binaries that evolve more rapidly through LISA’s band. Table 4 gives both the expected counts and fractions.

Population	Count, +30 yr	During mission (%)	+10 yr (%)	+30 yr (%)	+100 yr (%)
GWTC-3	$1.8^{+1.7}_{-0.91}$	$5.8^{+1.5}_{-1.1}$	$25^{+4.2}_{-3.4}$	$42^{+5.3}_{-4.5}$	$67^{+5}_{-5}$
GWTC-5.0 (LVK)	$3.6^{+5.2}_{-1.6}$	$14^{+16}_{-6}$	$39^{+22}_{-11}$	$55^{+19}_{-9.8}$	$75^{+11}_{-6.6}$
BPL2P	$3.2^{+9.1}_{-1.5}$	$15^{+30}_{-8}$	$39^{+33}_{-12}$	$54^{+26}_{-11}$	$74^{+16}_{-7.6}$
BPL3P	$3.4^{+24}_{-1.7}$	$17^{+46}_{-11}$	$40^{+45}_{-14}$	$56^{+35}_{-13}$	$75^{+20}_{-8.3}$
symBPL2P	$5.2^{+65}_{-3.3}$	$28^{+44}_{-20}$	$53^{+40}_{-26}$	$66^{+30}_{-22}$	$81^{+17}_{-14}$
symBPL3P	$11^{+3500}_{-8.8}$	$49^{+38}_{-42}$	$72^{+28}_{-45}$	$80^{+20}_{-36}$	$89^{+11}_{-21}$
BPL2P (free slope)	$4.1^{+43}_{-2.3}$	$23^{+53}_{-15}$	$47^{+45}_{-19}$	$61^{+34}_{-16}$	$78^{+19}_{-10}$

In all four extrapolated models, binaries merging 10–30 years after the LISA mission tend to have lower primary masses than those merging during it (Figure 4).

Binaries merging during the mission and those merging 10–30 years after its end are detected during the same four-year LISA observation. At a fixed observed frequency, binaries with smaller detector-frame chirp masses evolve more slowly (Peters [70]). At fixed primary mass, redshift, and observed frequency, a smaller companion lowers the chirp mass and lengthens the time to merger. The different companion-mass prescriptions in the primary-mass and symmetric models therefore affect the fraction of LISA detections that merge within each follow-up window.

We use these merger time cuts to identify potential multiband sources. We calculate their SNRs in one 40-km CE detector, using a lower frequency cutoff of 5 Hz and requiring $\rho_{\mathrm{CE}} > 8$ (Appendix § 11). Assuming CE is operating at merger, more than 99.9% of the GWTC-3 and GWTC-5.0 systems merging in each time window exceed this threshold, even at the lower end of the central 90% credible interval on that fraction.

Among LISA detections that merge during the mission or within 30 years after its end, the fraction detectable by CE is $55^{+5}_{-21}\%$ for symBPL3P extrapolated to 1000 $M_\odot$ and $99^{+1}_{-17}\%$ for BPL2P with a free high-mass slope. Some binaries in these extrapolations have such large redshifted masses that too little of their signal remains above 5 Hz to exceed the CE threshold. Table 2 gives the fractions for all models. For systems accessible in both bands, LISA’s earlier observation can constrain the merger time and sky position before the merger ([79, 47, 57]).

Population	During mission	Mission + 10 yr	Mission + 30 yr
GWTC-3	$> 99.9\%$	$> 99.9\%$	$> 99.9\%$
GWTC-5.0 (LVK)	$> 99.9\%$	$> 99.9\%$	$> 99.9\%$
BPL2P	$99.3^{+0.7}_{-3.3}\%$	$99.7^{+0.3}_{-2.3}\%$	$99.8^{+0.2}_{-2.1}\%$
BPL3P	$98.7^{+1.3}_{-7.6}\%$	$99.4^{+0.6}_{-6.2}\%$	$99.6^{+0.4}_{-5.9}\%$
symBPL2P	$97.5^{+2.4}_{-8.1}\%$	$98.7^{+1.3}_{-7.1}\%$	$98.9^{+1.1}_{-7.1}\%$
symBPL3P	$92^{+8}_{-21}\%$	$95^{+5}_{-20}\%$	$95^{+5}_{-21}\%$
BPL2P (free slope)	$96^{+4}_{-21}\%$	$98^{+2}_{-19}\%$	$99^{+1}_{-19}\%$

Discussion

Using the released GWTC-3 and GWTC-5.0 population constraints, we forecast $4.3^{+3.6}_{-2.0}$ and $6.6^{+5.5}_{-2.2}$ resolved binaries, respectively, in a four-year LISA mission at $\rho> 8$. The updated population gives a higher median count, while retaining the expectation of a small catalog found in earlier GWTC-3 forecasts ([28]). Binaries with primary masses between 100 and 300 $M_\odot$ account for most of the increase in the posterior mean forecast, while the updated population below 100 $M_\odot$ also contributes (§ 3.1).

LISA’s sensitivity to the high-mass population depends on the companion masses as well as the primary masses. At fixed primary mass, redshift, and observed frequency, a heavier companion increases the chirp mass, giving a larger strain but a faster frequency evolution. Both affect the SNR accumulated during the four-year mission. In symBPL3P, the additional high-mass Gaussian enters the mass distribution used for both black holes. Even at comparable local merger rates, the four f fits can therefore give different LISA counts because they assign different abundances to binaries with two massive components. The continuum, pairing, and low-mass cutoffs are also refitted, so this comparison does not measure the effect of the additional Gaussian alone. The fitted populations are restricted by the numerical precision of the likelihood: combinations that fail the variance requirement are excluded even if they are physically allowed. Its effect on the inferred population parameters and LISA counts has not been quantified (Appendix § 8.§ 3).

Extending the mass range does not by itself imply a large increase in the LISA count. With the merger-rate evolution held the same in the 300 and 1000 $M_\odot$ calculations, the steep BPL2P primary-mass distribution adds $0.22^{+4.5}_{-0.21}$ expected detections above 300 $M_\odot$. The symmetric fits permit shallower high-mass component distributions and hence more binaries with massive companions. For symBPL3P, the expected count changes from $7.2^{+35}_{-3.2}$ to $14^{+3500}_{-9.4}$ under the same rate evolution (Table 3). Its long upper interval comes from extrapolating populations that allow many more massive binaries than those giving counts near the median. These binaries add to the merger rate; the rates assigned to binaries with primary masses below 300 $M_\odot$ are unchanged.

Whether such a high-mass population exists depends on how these black holes form. Pair instability can prevent stellar collapse from producing black holes over a range of masses, while sufficiently massive progenitors can form black holes above that range ([38, 50, 96]). Remnants of earlier black hole mergers provide another route to high masses; their abundance depends in part on whether the remnants remain bound to their environment and merge again ([45, 46]). Such repeated mergers are one possible interpretation of reported high-mass structure ([63, 21, 48]). Our extrapolations do not impose a pair-instability gap or calculate the populations produced by these channels. The 300 $M_\odot$ boundary is the upper limit of the fits, not a physical gap boundary. The possible formation of black holes above the pair-instability gap motivates testing a different high-mass distribution, but does not prescribe its slope. Their contribution to LISA also depends on their merger rate and companion masses. ([46, 21]). The extrapolated counts also assume the fixed Madau–Fragos merger-rate evolution with redshift. Formation delays can change this evolution ([42]), and that uncertainty is not included in our intervals.

The O3 non-detection gives little additional information for the free-slope extrapolation considered here. Matching the rate density to BPL2P at 300 $M_\odot$ gives an expected O3 count of only $(5.9^{+4.7}_{-5.5}) \times10^{-5}$ in the region covered by the nonspinning injections. The probability of zero detections is therefore close to unity, leaving the LISA forecast essentially unchanged. The wider LISA interval reflects the allowed mass distributions above the matching point. This result is specific to the adopted normalization and O3 sensitivity coverage. In particular, we note that the sensitivity used here does not extend to all masses and mass ratios in the LISA extrapolation (Appendix § 12).

The largest symBPL3P forecasts also reach the limits of the signal and noise approximations. At large redshifted masses, the end of the inspiral moves to lower observed frequencies ([54, 56]), whereas our LISA calculation uses a leading-order inspiral without a mass-dependent termination (Appendix § 9.§ 2). The total BBH background can also become comparable to the assumed noise. Adding this entire background to the noise reduces the count by about 13% for a population near the 95th percentile of expected count and 60% for one near the 99th percentile (Appendix § 9.§ 3). The background used in these checks includes signals that could be subtracted, and the selected populations do not give a corrected credible interval. The extreme counts need to be assessed with mass-dependent signal termination and a calculation of detection and source subtraction together; the upper percentile is also only approximate because it is estimated from relatively few samples in the high-count tail (Appendix § 10).

For the released GWTC-5.0 population, $39^{+22}_{-11}\%$ of LISA detections merge during the four-year mission or within ten years after its end. More than 99.9% of those mergers exceed our reference CE threshold, even at the lower end of the central 90% credible interval on that fraction, assuming CE is operating at merger. For this population, the number of joint observations is therefore limited mainly by which binaries LISA detects and when they merge, rather than by CE’s sensitivity. The follow-up interval also selects different masses among LISA detections: in all four extrapolated models, binaries merging 10–30 years after mission end tend to have lower primary masses than those merging during the mission (Figure 4). At a fixed observed frequency, smaller detector-frame chirp masses take longer to merge.

to merge. Waiting longer for mergers therefore adds lower-mass systems to the potential multiband sample.

Even a small joint sample would follow the same binaries from their early inspiral to merger. LISA can provide information about the orbit before it has further circularized, while ground-based observations measure the late inspiral and remnant ([26, 69, 27]). Combining the two bands can improve binary-parameter measurements and tests of general relativity ([92, 33]). Eccentricity, masses, and spins can help distinguish formation channels, although they need not uniquely identify a binary containing a previous merger remnant ([46]). Our circular-inspiral forecasts estimate the available sources, not the precision of those measurements. Further ground-based detections of systems such as GW190521 and GW231123 will constrain the abundance and companion masses of massive binaries, improving the forecasts for both LISA and joint observations.

Conclusion

We have used ground-based population constraints to forecast the number and masses of binary black holes detectable by LISA. For a 4-year mission at $\rho> 8$, we obtain $6.6^{+5.2}_{-2.2}$ expected detections using the baseline GWTC-5.0 constraints, compared with $4.3^{+3.6}_{-2.0}$ using GWTC-3. The quoted errors enclose central 90% intervals from population and merger-rate uncertainty.

We find that the assumed mass model can have a significant impact on the forecast. For our GWTC-5.0 fits, extrapolating the primary-mass models from 300 to 1000 $M_\odot$ changes their median expected counts little, while widening their count intervals. The symmetric mass models allow a larger contribution from binaries with two massive components. The broadest forecasts arise when a weakly constrained high-mass distribution is extrapolated for both components. For the symmetric broken power law plus three peaks mass model, symBPL3P, the expected count changes from $7.2^{+35}_{-3.2}$ to $14^{+3500}_{-9.4}$, using the same Madau–Fragos merger-rate evolution in both cases and preserving the merger rate below 300 $M_\odot$. The broad upper interval reflects the uncertain high-mass abundance, but it is worth mentioning that our inspiral and noise approximations also limit the reliability of its extreme counts.

For the broken power law plus two peaks model, BPL2P, allowing the mass distribution above 300 $M_\odot$ to differ from its fitted form raises the upper end of the expected-count interval. When this extrapolation is matched to the fitted rate per unit primary mass at 300 $M_\odot$, it predicts negligible detections within the O3 search’s sensitivity coverage. The O3 IMBH non-detection therefore leaves this forecast essentially unchanged.

For the released GWTC-5.0 population, $39^{+22}_{-11}\%$ of LISA detections merge during the mission or within ten years after its end. Essentially all of these mergers would be detectable by our reference CE configuration if it is operating at merger. In all four extrapolated models we consider, binaries merging 10–30 years after mission end tend to have lower primary masses than those merging during it. Further ground-based observations will constrain how often massive black holes merge and the masses of their companions, improving these forecasts. Observations in both frequency bands would follow the same binary from its early inspiral to merger, combining LISA’s information about the orbit with ground-based measurements of the late inspiral and remnant.

Data Availability

The lisarates analysis code and data needed to reproduce the forecasts and figures will be made publicly available with this work.

Acknowledgments

We thank the LISA Consortium and the LIGO–Virgo–KAGRA Collaboration for public data products. This work made use of publicly released GWTC-3 and GWTC-5.0 population posterior samples and the O3 IMBH search sensitivity data.

This material is based upon work supported by the National Aeronautics and Space Administration under Grant No. 22-LPS22-0025. Computations used the Vera cluster at the Pittsburgh Supercomputing Center, dedicated to the McWilliams Center for Cosmology and Astrophysics at Carnegie Mellon University. This work used Jetstream2 at Indiana University through allocation PHY240308 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (AC-CESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

This research has made use of data or software obtained from the Gravitational Wave Open Science Center (https://gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA [14] (Abac et al. 2026f). This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation, as well as 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 GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.

Claude Fable 5.1 (Anthropic) and Codex GPT-6 Astra (OpenAI) assisted with code development, numerical checks, figure preparation, and manuscript editing. The authors are responsible for the analysis and manuscript.

Facilities: LIGO, Virgo, KAGRA

Software: lisarates (this work), dynesty [81], JAX [25], LEGWORK [93], LALSuite [58, 95], Astropy [22], NumPy [49], SciPy [91], Matplotlib [53], Seaborn [94]

Population Definitions and Merger Rates

We use source-frame component masses $m_1 \ge m_2$ and mass ratio $q = m_2/m_1$. The population parameters, collectively denoted by $\Lambda$, specify the intrinsic mass and spin distributions and the redshift evolution of the merger rate. The density $p(m_1, q \mid\Lambda)$ describes all binaries in the model, whether or not they are detectable. For compactness, the equations we write show $p(m_1, q)$ in place of $p(m_1, q \mid\Lambda)$, but the mass distributions still depend on the fitted population parameters. The local merger rate $R_0$ is the number of mergers per unit comoving volume and source-frame time at $z = 0$, in $\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$. Thus $R_0p(m_1,q \mid\Lambda)\,dm_1\,dq$ gives the local rate in a small primary-mass and mass-ratio interval.

Baseline Populations

For GWTC-3, the primary-mass distribution is a mixture of a truncated power law and a Gaussian peak, multiplied by a low-mass taper (Talbot & Thrane [84]; Abbott et al. [14]). The power-law slope, peak location and width, relative peak weight, mass cutoffs, and taper width are inferred parameters.

For a component mass $m$, let $x = m - m_{\min}$ be the distance above the lower mass cutoff. The taper rises from zero to one over a mass interval of width $\delta$ (Talbot & Thrane [84]):

For $\delta= 0$, we use a sharp lower cutoff. The same taper is applied to the secondary mass $m_2 = qm_1$. At fixed primary mass, the normalized mass-ratio density is (Talbot & Thrane [84])

Its support is $m_{\min}/m_1 \le q \le1$; it vanishes outside this interval. The exponent $\beta_q$ controls the relative abundance of unequal- and equal-mass binaries. The normalization $Z_q(m_1)$ depends on the primary mass because the allowed secondary-mass range changes with $m_1$.

The released GWTC-5.0 Brokenpowerlaw+2Peaks posterior uses a broken power law and two Gaussian peaks, with separate low-mass tapers for the primary and secondary masses (Abac et al. [2]). We write the normalized broken power law as (Abbott et al. 2021)

Here $m_b$ is the break mass, and $\alpha_1$ and $\alpha_2$ set the slopes below and above it. The two power laws meet continuously at $m_b$. The lower edge is $m_{\mathrm{low}}$, the fitted upper limit is $M_{\mathrm{fit}} = 300\,M_\odot$, and $Z_b$ is the integral of the piecewise shape over this interval.

We evaluate these distributions using the joint posterior samples released with each catalog, as described in § 2.1.

Alternative Mass Models

The four new fits use the same broken power law, with either $K = 2$ or $K = 3$ Gaussian peaks. Before tapering, each Gaussian is normalized on the fitted mass interval:

This density is zero outside the integration limits. The parameters $\mu_j$ and $\sigma_j$ are the location and width of the underlying Gaussian, before truncation and tapering. The normalized mass distribution is

The coefficients $\lambda_j$ weight the component densities before the low-mass taper is applied: $\lambda_0$ weights the broken power law, $\lambda_j$ with $j \ge1$ weights peak $j$, and all weights are nonnegative with $\sum_{j=0}^{K} \lambda_j = 1$. The normalization $Z_K$ is the integral of the numerator over $m_{\mathrm{low}} \le m \le M_{\mathrm{fit}}$, so that $\int B_K(m)\,\mathrm{d}m = 1$. The taper changes the relative integrals of the components; the $\lambda_j$ are therefore not generally the fractions remaining in each component after tapering.

For BPL2P and BPL3P, $B_K$ is the primary-mass density. Multiplying it by the companion’s conditional mass-ratio density gives

where $p(q\mid m_1)$ follows Eq. (2), replacing $m_{\min}$ by the secondary-mass cutoff $m_{\mathrm{low},2}$ and $\delta$ by its taper width $\delta_{m,2}$. The secondary cutoff satisfies $3\,M_\odot\le m_{\mathrm{low},2} \le m_{\mathrm{low}}$.

For symBPL2P and symBPL3P, the same $B_K$ describes both component masses before pairing. We weight each pair by $q^\beta$. The resulting probability density in primary mass and mass ratio is

on the ordered domain $3\,M_\odot\le m_2 \le m_1 \le300\,M_\odot$. Here $Z_{\mathrm{comp}}$ normalizes the distribution over ordered binary masses, so that $\int p_{\mathrm{comp}}(m_1,q)\,\mathrm{d}m_1\,\mathrm{d}q = 1$ on the fitted mass range:

When extrapolating beyond $300\,M_\odot$, we keep this normalization fixed so that the rate below this mass is unchanged (Appendix § 8.§ 4).

The factor $m_1$ in Eq. (7) follows from $\mathrm{d}m_2 = m_1\,\mathrm{d}q$ at fixed primary mass. Positive $\beta$ gives greater weight to pairs with similar masses, while negative $\beta$ favors unequal masses. In the symmetric models, $B_K$ describes both component masses before pairing. The primary-mass distribution also depends on the abundance of allowed companions and their pairing weights, as we obtain it by integrating $p_{\mathrm{comp}}(m_1,q)$ over $q$ at a fixed $m_1$. It is therefore not simply $B_K(m_1)$. At fixed primary mass, the companion distribution is proportional to $B_K(qm_1)q^\beta$, so the peaks also affect which companions are favored.

These distinctions matter when comparing the peak fractions between models. We quote the fraction of primary black holes belonging to the third Gaussian in the intrinsic population, rather than its mixture weight. We calculate this fraction by integrating the third Gaussian’s contribution to the normalized primary-mass distribution over the fitted mass range, including the low-mass taper in both models and pairing in the symmetric model. For the symmetric model, we use the third-Gaussian term of $B_K(m_1)$ in Eq. (7), but sum over all companion components through $B_K(qm_1)$ and keep the full population normalization $Z_{\mathrm{comp}}$. The companion need not belong to the third Gaussian.

Population Inference and the Local Merger Rate

All four fits use the same 259 GWTC-5.0 events and the cumulative O1–O4b injection set ((LIGO Scientific Collaboration et al. 2026b)), which contains simulated signals used to estimate the probability of detection by the ground-based searches. We infer the mass, spin, and redshift parameters with a hierarchical likelihood that integrates over the local merger rate using the prior $p(\mathcal{R}_0) \propto1/\mathcal{R}_0$ and accounts for selection ([65, 85, 39]). The effective spin is $\chi_{\mathrm{eff}} = (\chi_1 + q\chi_2)/(1+q)$ ([74, 18]), where $\chi_1$ and $\chi_2$ are the dimensionless component spins projected along the orbital angular momentum. Its density is a Gaussian with location $\mu_\chi$ and width $\sigma_\chi$, normalized on $[-1,1]$ ([78, 67]). Within each population, this spin distribution is independent of mass and redshift, and the mass distribution does not evolve with redshift. The merger-rate density follows $\mathcal{R}(z) = \mathcal{R}_0(1+z)^\kappa$ over the fitted range $0 \le z \le1.9$; $\kappa$ controls its growth or decline. The fits use nested sampling ([80]) with dynesty ([81]). Table 1 lists the mass, spin, and redshift priors.

Parameter	Prior
Power-law slopes $\alpha_1,\alpha_2$	$U[-4,12]$
Break mass $m_b/M_\odot$	$U[20,50]$
Primary or common lower edge $m_{\mathrm{low}}/M_\odot$	$U[3,10]$
Taper widths $\delta_m/M_\odot$	$U[0,10]$
Secondary edge $m_{\mathrm{low},2}/M_\odot$	$U[3,m_{\mathrm{low}}]$
First peak mean $\mu_1/M_\odot$	$U[5,20]$
Second peak mean $\mu_2/M_\odot$	$U[25,60]$
Third peak mean $\mu_3/M_\odot$	$U[45,90]$
Peak widths $\sigma_j/M_\odot$	$U[0,10]$
Mixture weights $(\lambda_0,\ldots,\lambda_K)$	$\operatorname{Dirichlet}(1,\ldots,1)$
Mass-ratio or pairing exponent $\beta_q$ or $\beta$	$U[-2,7]$
Spin mean $\mu_\chi$	$U[-1,1]$
Log spin width $\ln\sigma_\chi$	$U[-5,0]$
Rate exponent $\kappa$	$U[-10,10]$

**Note—**The primary-mass models have separate lower mass cutoffs and smoothing widths for the primary and secondary. The symmetric models use the same mass distribution, cutoff, and smoothing width for both black holes. The Dirichlet prior is uniform over nonnegative mixture weights that sum to one. The third-peak parameters apply only to the three-peak models. All fits have a fixed upper mass limit of $300M_\odot$.

The hierarchical likelihood is computed by reweighting single event-posterior samples and simulated signals. Let $I_i$ denote the population-weighted likelihood integral for event $i$ and $\xi$ the population-averaged detection efficiency. After marginalizing over the merger rate, the log likelihood contains $\sum_i \ln I_i-N_{\mathrm{obs}}\ln\xi$, where $N_{\mathrm{obs}}=259$ is the number of events. Assuming independent Monte Carlo estimates, we estimate the numerical variance as [83]

During inference, we set the likelihood to zero wherever this estimate is at least one. This avoids using poorly determined likelihood values, but also restricts the populations the fit can explore (Abac et al. 2026c).

Many posterior samples have estimated log-likelihood variance close to the imposed limit, $\widehat{\sigma^2_{\ln\mathcal{L}}}=1$. The forecast intervals are conditional on this requirement; we have not quantified how it changes the inferred populations or LISA counts. The variance requirement controls Monte Carlo error in the standard hierarchical likelihood, arising from both the selection-function injections and the event-level parameter-estimation samples. Following Talbot & Golomb [83] and Essick & Farr [35], these sample sets must provide sufficiently precise likelihood estimates to avoid biasing population constraints. The imposed limit of 1 on the estimated log-likelihood variance is a numerical reliability criterion, not a physical population constraint, and it does not by itself establish that the posterior is unbiased [51]. Testing populations excluded by this requirement requires more accurate likelihood estimates, rather than simply relaxing the threshold.

We reconstruct a local merger-rate distribution for each draw of the population parameters. A population that is harder to detect requires a larger intrinsic rate to account for the same observed catalog. This dependence is described by the sensitive spacetime volume [86, 39]

Here $T_{\mathrm{GW}}$ is the ground-based observing duration and $dV_{\mathrm{c}}/dz$ is the full-sky comoving volume per unit redshift Hogg [52]. The detection probability $\overline{P}_{\mathrm{det}}$ is averaged over the population’s masses and spins, source orientation, and the observing periods. The factor $(1+z)^\kappa$ describes the rate evolution, while $(1+z)^{-1}$ converts source-frame time to observer time. We estimate $\langle VT\rangle_\Lambda$, in $\mathrm{Gpc}^3\mathrm{yr}$, by weighting the recovered injections by the population density relative to the injection density, normalizing by the total number generated, and including the observing duration [86, 34].

For fixed $\Lambda$, the expected number of ground-based detections is $\mathcal{R}_0\langle VT\rangle_\Lambda$. The Poisson count likelihood and the prior $p(\mathcal{R}_0)\propto1/\mathcal{R}_0$ give the conditional rate posterior Mandel et al. [65]

Here $\{d_i\}$ denotes the event data, $N_{\mathrm{obs}}=259$ is the number of events, and $\Gamma$ is the Gamma function. This is a Gamma density on $\mathcal{R}_0>0$, with shape $N_{\mathrm{obs}}$, rate parameter $\langle VT\rangle_\Lambda$, and conditional mean $N_{\mathrm{obs}}/\langle VT\rangle_\Lambda$. Both the expected-count intervals and the shaded mass-distribution bands in the Figures include the conditional rate uncertainty in Eq. (11). We integrate over this distribution for the count intervals. For the bands in Figure 3, we pair each population sample with one rate drawn from its own conditional distribution and calculate the spread of the resulting merger-rate and detected-count densities.

Extrapolated Merger Rates and Redshift Evolution

Let $\tilde{p}(m_1,q\mid\Lambda)$ denote the mass density obtained by extending the upper mass limits of the fitted functions to $1000\,M_\odot$, retaining the original normalization of each mixture component and of $B_K$. Below $300\,M_\odot$, $\tilde{p}$ equals the fitted density exactly. Above this boundary, we retain the fitted slopes, Gaussian parameters, and pairing prescription. For BPL2P and BPL3P, the conditional mass-ratio density is normalized over the allowed companion masses for each primary mass, including primaries above $300\,M_\odot$. For symBPL2P and symBPL3P, we use the value of $Z_{\mathrm{comp}}$ calculated with both masses below $300\,M_\odot$. Recalculating it over the enlarged mass range would reduce the rate of the original binaries, contrary to the extrapolation defined here. The extrapolated local rate is

Unlike the fitted density, $\tilde{p}$ is not normalized to one over the enlarged mass range. The total local rate becomes

Thus $\mathcal{R}_0$ still denotes the rate on the original mass range, and $\Delta\mathcal{R}$ is the rate of binaries with a primary above $300\,M_\odot$. The rate for every pair with $m_1\leq300\,M_\odot$ is unchanged, as is the rate per unit primary mass below this boundary. In all four models, an added binary may have one or both components above $300\,M_\odot$. The abundance of secondaries below this boundary can therefore increase when they pair with newly added primaries.

For all extrapolated populations and their matched $300\,M_\odot$ comparisons, we obtain the rate at redshift $z$ by multiplying the local mass-dependent rate by $g(z)$, where

Here $\psi(z)$ is a dimensionless function with the star-formation-history shape of Madau & Fragos [62], parameterized by its peak redshift as in Babak et al. [24]. Only the ratio $\psi(z)/\psi(0)$ enters the merger rate so no absolute star-formation-rate normalization is required. The exponents $d$ and $r$ give the low- and high-redshift slopes in $1+z$, and $z_p$ is the redshift of the maximum. These parameters are fixed, not inferred from the LVK catalog. The 300 and 1000 $M_\odot$ forecasts use the same $g(z)$, so their difference is due to the contribution of the added masses. We note that this assumed merger-rate history does not include a distribution of delays between formation and merger, and the forecast intervals do not include uncertainty in that history.

LISA Noise and Signal-to-Noise Ratio

Instrumental Noise and Galactic Foreground

We use the equal-arm effective strain sensitivity in Eq. (56) of Babak et al. [23], with arm length $L = 2.5 \times10^{9}\ \mathrm{m}$, acceleration-noise level $3 \times10^{-15}\ \mathrm{m}\ \mathrm{s}^{-2}\ \mathrm{Hz}^{-1/2}$, and optical-metrology noise level $15 \times10^{-12}\ \mathrm{m}\ \mathrm{Hz}^{-1/2}$. The effective strain-noise spectrum $S_{n,\mathrm{inst}}$ includes the sky- and polarization-averaged response for one channel, without the improvement from combining two channels. We add the four-year foreground from unresolved Galactic double white dwarfs of Robson et al. [76], using the single-channel amplitude $A = 1.80 \times10^{-44}\ \mathrm{Hz}^{-1}$ of Cornish & Robson [30]. The Galactic foreground is expressed in the same single-channel convention. Writing its spectrum as $S_c$, we use $S_n = S_{n,\mathrm{inst}} + S_c$. The factor of two in squared SNR for two effective channels enters in Eq. (19), $\rho^2 = 2 \times4 \Xi^2 \rho_0^2$. No additional channel factor is applied to $S_n$ or $\rho$.

Signal-to-Noise Ratio

For source-frame chirp mass $\mathcal{M}$ and source-frame time to coalescence $\tau$, the leading-order circular inspiral gives the source-frame frequency (Peters [70]; Cutler & Flanagan [31])

A source at redshift $z$ evolves for $T_{\mathrm{obs}}/(1 + z)$ in its own frame during the mission, or until merger if that occurs first. We convert to observed frequency using $f_{\mathrm{det}} = f/(1 + z)$ and integrate over the frequencies traversed during the mission, restricted to the adopted LISA band, $3 \times10^{-5}$–$0.5\ \mathrm{Hz}$. We use the leading-order inspiral waveform up to the upper frequency of this band without a mass-dependent termination. In the integral below, $f$ denotes the observed frequency. We define

Here the limits are the observed frequencies at the start and end of the observation, restricted to the adopted band. The dimensionless SNR factor is

At inclination $\iota$,

These expressions average over sky position and polarization and include both effective channels [23]. Isotropic inclinations give $\langle\Xi^2\rangle= 4/5$ and $\langle\rho^2\rangle_\iota= (32/5)\rho_0^2$, while the face-on value is $\rho^2 = 16\rho_0^2$. No additional channel factor is applied.

For matched source parameters, single-channel noise curve, and frequency limits, our face-on $\rho^2$ is a factor of two below the approximate face-on $\rho^2$ in Eq. (2.2) of Buscicchio et al. [28]. They use that expression for preliminary catalog selection and recompute SNRs from simulated instrument data for their reported detection forecasts. As an independent check of the SNR normalization, we compare with LEGWORK, a code for gravitational waves from stellar-origin binaries [93]. For circular binaries whose frequency changes little during the mission, the SNRs agree to better than 1% when the source parameters and noise curves are matched.

For the four reference binaries in Figure 2, we use effective characteristic strain to display the same SNR calculation. Writing Eq. (18) as $\rho_0 = \mathcal{A}\sqrt{I}$, we define

The ratio of these curves gives the SNR accumulated during the observation:

The reference curves and bands include both channels and the root-mean-square response over isotropic inclinations. For each fixed source-frame mass pair, the bold curve uses $z = 0.40$ and the band spans $0.1 \le z \le1.0$. These bands compare redshifts; they are not uncertainty intervals. We recalculate the luminosity distance, detector-frame chirp mass, and starting frequency at each redshift, keeping the observation in the final four observer years before merger.

For the purple sample, we choose one BPL2P population near the posterior medians of its mass-distribution parameters. We extrapolate it to $1000\,M_\odot$ with the Madau–Fragos merger-rate history. We sample masses, redshifts, merger times, and isotropic inclinations, weighting sources by their expected abundance during the mission. Each track covers the frequencies traversed in four observer years, or until the signal cutoff if reached earlier; merger need not occur during the mission. We consider sources with total four-year SNR $0.1 \le\rho\le8$ and a track intersecting the displayed frequency range. The lower cut selects the louder part of the undetected population for illustration. The shading shows the smoothed, population-weighted distribution of individual tracks in frequency and strain in the upper panel, and in frequency and accumulated SNR in the lower panel. Each panel is scaled to its own maximum for display. The lower-panel SNR is accumulated from the start of each source’s observation. No BBH contribution is added to the noise in Figure 2, but the summed BBH background is treated separately in Appendix § 9.§ 3.

Figure 7 (left) instead chooses the starting frequency that maximizes the four-year SNR for each mass and redshift. Its $\rho= 8$ contour is therefore a reach estimate with the averaged response, whereas the four reference binaries in Figure 2 are observed during their final four years before merger. The population forecasts integrate over merger times and inclinations. The horizon search extends to $z = 14$; the plot stops at $z = 10$. This numerical search limit is separate from the redshift integration used for the population forecasts.

Binary Black Hole Background

To check whether the binary black hole (BBH) background could affect our largest detection forecasts, we calculate the total circular-inspiral background for symBPL3P extrapolated to $1000\,M_\odot$, with the adopted merger-rate evolution with redshift, using the relation between merger rate and emitted gravitational-wave energy [71 Babak] et al. 2023). This background includes signals from individually detectable BBHs. Subtracting their signals would reduce the remaining background, but we do not estimate the residual background here.

Figure 7 (right) shows the BBH background power spectral density (PSD), $S_{\mathrm{BBH}}$, relative to $S_n$, the sum of instrumental noise and the foreground from unresolved Galactic double white dwarfs. The band includes uncertainty in the fitted mass distribution and merger rate, conditional on the extrapolation and fixed merger-rate evolution with redshift. At 3 mHz, $S_{\mathrm{BBH}}/S_n = 0.08^{+0.55}_{-0.02}$ (median and central 90% interval). The dashed curve is a separate example selected near the 99th percentile of expected detection count; it is not the upper edge of the background interval.

To test the effect on individual BBH detections, we repeat the SNR calculation for populations near the median, 95th, and 99th percentiles of expected count, after adding each population’s total BBH background to the assumed noise. This treats the combined BBH signal as stationary Gaussian noise, including signals from individually detectable BBHs. The expected counts above $\rho= 8$ decrease by about 5%, 13%, and 60%, respectively. These checks evaluate individual populations, not posterior credible intervals or the residual noise after source subtraction. The most extreme forecasts therefore require a calculation of detection and source subtraction together ([24]).

Expected Counts and Population Uncertainty

Let $r(m_1, q, z)$ be the merger rate per unit source-frame time, comoving volume, primary mass, and mass ratio. For each population and rate draw, the expected number observable by LISA is

The factor $P_{\iota}$ is the fraction of isotropic inclinations for which Eq. (19) exceeds the SNR threshold, using the frequencies traversed during the mission. The time to coalescence $t_{c,\mathrm{det}}$ is measured in the observer frame from mission start. We use the merger rate at redshift $z$ to infer the number of binaries that are still inspiraling at that redshift. We assume that this rate changes negligibly over their remaining time to merger, which can extend beyond the fiducial four-year mission ([24, 28]). At fixed masses and redshift, equal intervals of time to coalescence therefore contain the same expected number of BBHs before applying the LISA selection. We assume the redshift dependence of $r(m_1, q, z)$, allowing the merger rate to vary across cosmic time.

We evaluate the mass, redshift, and time integrals numerically. For the catalog forecasts and fitted populations, we integrate the inferred power-law rate evolution over $0 \le z \le1.9$. For the extrapolations and their matched 300 $M_\odot$ comparisons, we use Eq. (15) and include redshifts and merger times throughout the detectable range. The upper time-to-merger bound follows from the low-frequency edge of the LISA band. We choose the redshift bound so that the maximum possible SNR is below threshold beyond it. The mission limits the observed inspiral to at most four years, but we stop earlier if the binary merges or leaves the LISA frequency range used in the calculation.

The final baseline forecasts use all 11,184 released GWTC-3 samples and all 8,200 released GWTC-5.0 samples. The alternative forecasts use 4,000 population samples per model, with conditional rate uncertainty from Eq. (11). We assess numerical stability with finer quadrature, and sensitivity to the finite population sample by comparing subsets and bootstrap resamples of the stored draws. For symPBL3P extrapolated to 1000 $M_\odot$, the upper end of the central 90% interval on the expected count is about 3500 detections. Resampling the inferred populations with replacement gives a central 90% range of approximately 3000–4200 for this endpoint. This measures uncertainty in the endpoint caused by the finite population sample, not an additional interval on the expected count. We therefore treat the quoted upper endpoint as approximate.

For a time $\Delta t$ after mission end, the merger-window count uses $0 \le t_{c,\mathrm{det}} \le T_{\mathrm{obs}} + \Delta t$ in Eq. (22). For each population, the expected fraction merging in that window is $f_{\Delta t} = \mu_{\Delta t}/\mu_{\mathrm{all}}$, where $\mu_{\Delta t}$ and $\mu_{\mathrm{all}}$ are calculated for the same population parameters. The overall rate normalization cancels in this ratio. We report the median and central 90% interval of $f_{\Delta t}$ across the population posterior.

The quoted 90% count intervals describe uncertainty in $\mu$ from the inferred population parameters and merger rate. Neither these intervals nor those on $f_{\Delta t}$ include Poisson fluctuations in a realized LISA catalog. To calculate the probability of an empty LISA binary black hole catalog, we include these fluctuations by averaging $P(N = 0 \mid\mu) = e^{-\mu}$ over the population and rate posterior:

Cosmic Explorer Detectability

We use the published strain sensitivity for a 40-km CE detector (CE-T2000017-v9; [82])¹. The adopted sensitivity curve is evaluated at $15^\circ$ from normal incidence. We use IMRPhenomXHM waveforms, which include inspiral, merger, ringdown, and higher harmonics ([44]). We calculate the CE SNR through:

where $F_+$ and $F_\times$ are the detector antenna factors and $S_{n,\mathrm{CE}}$ is the one-sided strain-noise PSD ([31]). The waveform uses the same luminosity distance as the LISA calculation. We use the long-wavelength antenna response of an L-shaped detector with the published strain curve and do not include the full direction-dependent finite-arm response.

For the four fitted models and the free-slope extrapolation, we set the aligned component spins to $\chi_{1z} = \chi_{2z} = \chi_{\mathrm{eff}}$, with $\chi_{\mathrm{eff}}$ drawn from the fitted distribution. For the released catalog populations, we use the aligned projections of their inferred spin magnitudes and tilts. Our prescriptions neglect precession. We integrate the CE detection probability over isotropic sky positions and polarization angles, assuming the inclination distribution selected by LISA.

Let $\mu_{\mathrm{LISA+CE},\Delta t}$ be the expected count after imposing $\rho_{\mathrm{CE}} > 8$ in addition to the LISA selection and merger-time window in Eq. (22). Table 2 reports $F_{\mathrm{CE}} = \mu_{\mathrm{LISA+CE},\Delta t}/\mu_{\Delta t}$: the fraction detectable by CE among LISA sources that merge in that window.

For the released GWTC-5.0 population, the median CE SNR among LISA sources merging during the mission or within 30 years after its end is $(2.5^{+0.9}_{-1.1}) \times10^3$ (with the spread indicating the population uncertainty on this median). For mergers during the mission or within ten years after its end, the fraction of all LISA detections that also exceed the CE threshold is $39^{+22}_{-11}\%$, unchanged at the quoted precision from the fraction without a CE requirement. The lower ends of the central 90% credible intervals on the CE detection fractions also exceed $99.9\%$ for all four fits restricted to $300\,M_\odot$, with either rate history.

These estimates assume that CE is operating when the binaries merge. The intervals include population uncertainty, conditional on the adopted sensitivity and spin prescription.

Free-Slope Extrapolation and IMBH Limits

For the free-slope extrapolation of BPL2P, let $\mathcal{R}_{\mathrm{ref}}$ denote its fitted rate. The primary-mass rate above $m_j = 300\,M_\odot$ is

We retain the reference model’s conditional mass-ratio prescription and use Eq. (15) for rate evolution. The original population below $m_j$ is unchanged. For each draw of the reference model’s shape and rate, we start with $\alpha\sim U[-2, 8]$ and update it by

Here $\mathrm{O3}$ denotes zero detections with $m_1 > 300\,M_\odot$, and $\mu_{\mathrm{O3}}$ is the expected number of O3 detections from this added population. For a given fitted population and merger rate, we normalize the distribution above over $\alpha$. We average these slope distributions over the posterior from our BPL2P fit to GWTC-5.0, keeping its population and rate sample weights unchanged. This update uses the O3 non-detection to reweight only the extrapolation slope.

GW190521 remains part of the fitted population; the zero count applies only to binaries with $m_1 > 300\,M_\odot$ ([10, 11]). We calculate $\mu_{\mathrm{O3}}$ by integrating the tail rate against the nonspinning sensitive spacetime volume from the released O3 IMBH injections ([11]).² We weight the recovered injections by the redshift dependence in Eq. (15), normalized to its local value, retaining the release’s spacetime-volume normalization and cosmology. The published sensitivity nodes and their upper limits provide a check of this reconstruction. They are limits for specified source populations, not upper bounds on $\mathrm{d}R/\mathrm{d}m_1$ at every mass.

We interpolate the sensitive spacetime volume only within the region enclosed by the sampled mass and mass-ratio combinations in $(\log m_1, \log q)$. Its direct coverage ends at $m_1 = 400\,M_\odot$ and has $q \ge0.1$, with additional restrictions from the node layout. No direct O3 sensitivity is assigned outside this region, because the grid does not supply it. The calculation is conditional on the nonspinning exposure and does not marginalize over the spin dependence of the search. At the rate densities inferred at $300\,M_\odot$, $\mu_{\mathrm{O3}} \ll1$ across the tested slopes. The O3 non-detection therefore has a negligible effect on the LISA forecast.

Forecast Tables

**Note—**The GWTC-3 upper mass cutoff is inferred and differs between posterior draws. The GWTC-5.0 (LVK) row uses the released posterior; BPL2P and BPL3P use our primary-mass fits, and symBPL2P and symBPL3P use our symmetric fits to GWTC-5.0. The free-slope row uses BPL2P with an independent slope above 300 $M_{\odot}$. Rows with fitted power-law evolution integrate over $0 \le z \le1.9$. Madau–Fragos denotes the fixed rate history in Eq. (15). The extrapolations leave the merger rate for every component-mass pair with $m_1 \le300 M_{\odot}$ unchanged. For each assumed population model, entries give the median expected count with errors to its 5th and 95th posterior percentiles. These errors describe uncertainty in the merger rate and population parameters; they do not give the range of integer counts in a realized LISA catalog. The broad upper interval for symBPL3P depends on the approximations in Appendix § 9; its upper endpoint is also approximate because of the finite population sample (Appendix § 10).

**Note—**Entries show median expected counts or merger fractions, with lower and upper errors enclosing their central 90% credible intervals. The GWTC-3 and GWTC-5.0 rows use the released posteriors. Our primary-mass fits (BPL2P, BPL3P) and symmetric fits (symBPL2P, symBPL3P) are extrapolated to 1000 $M_\odot$; the free-slope row uses BPL2P with an independent slope above 300 $M_\odot$. These extrapolations use the Madau–Fragos merger-rate history. No ground-based detection cut is applied.