Mock Catalogs of Strongly Lensed Gravitational Waves via A Halo Model Approach with Ground-based Detectors

Abstract: As plans for the construction of third-generation gravitational wave (GW) detectors advance, research into strongly lensed GWs has become increasingly critical. It is anticipated that hundreds of multi-image lensed GWs will be detected annually. We present a comprehensive suite of lensed GW mock catalog derived from a composite lens mass model incorporating dark matter halos, galaxies, and subhalos. We analyze three source populations with four detector network configurations considering the earth rotation. Our simulations encompass not only conventional doublets and quadruplets but also subhalo-lensed events, highly magnified systems, and complete three or five image systems with a detectable central image, a feature distinct from optical lensing. For the joint ET+CE network, we forecast an annual detection rate of approximately 400 doublets and 36 quadruplets. Notably, this population includes roughly 107 events lensed by subhalos and 20 complete systems with detectable central images. Furthermore, we analyze high-magnification events ($μ> 3$), predicting approximately 360 such cases. Under a more relaxed selection criterion that requires only at least one lensed signal to exceed the detection threshold, we estimate a total of approximately 617 lensed events. We also investigate the impact of variations in lens mass models and stellar evolution models on event rates, as well as the distributions of SNR pairs and time delays. These results establish a more physically grounded statistical prior for the future identification and authentication of lensed GW signals. The Gravitational Waves-Lensing Mock Catalog (GW-LMC) have been made publicly available.

• The paper develops the public GW-LMC mock catalog using composite galaxy, group, cluster, and subhalo lenses with GWTC-constrained BBH, BNS, and NSBH populations.
• For ET+CE, it forecasts about 395 doublets, 36 quadruplets, 107 subhalo-lensed events, 20 detectable central-image systems, and 617 triggered events annually.
• The catalog provides actionable time-delay and SNR-ratio priors, including threshold-dependent flux inversion, while highlighting binary-evolution and subhalo-model uncertainties exceeding order-of-magnitude rate differences.

Overview and motivation

This paper (Li et al., (2603.09289) presents a comprehensive mock catalog of strongly lensed gravitational wave (GW) events, constructed using a composite halo-model lens population rather than the simplified galaxy-scale SIS/SIE models that have dominated previous event-rate forecasts. The motivation is twofold: third-generation (3G) detectors such as the Einstein Telescope (ET) and Cosmic Explorer (CE) are expected to detect hundreds of lensed multi-image systems annually, and the identification of such systems requires statistically grounded priors on time delays, SNR ratios, and image multiplicities. By extending the SL-Hammocks framework of Abe et al. [(2603.09289) references therein], the work unifies lensing predictions across galaxy, group, and cluster scales, and explicitly includes dark matter subhalos—components largely absent from prior GW lensing simulations such as Yang et al. (2021), which considered only galaxy-scale doublets.

The paper's central quantitative results are striking: for a combined ET+CE network, the authors forecast approximately 400 doublets and 36 quadruplets of BBH events per year, roughly 107 subhalo-lensed events annually, about 20 systems per year in which the typically demagnified central image is detectable, and roughly 360 high-magnification ($\mu > 3$) events. Under the looser criterion that at least one image exceeds SNR 8, the total triggered population is approximately 617 events per year. The resulting catalog, GW-LMC, is publicly released on GitHub.

Source population modeling

The source model covers BBH, BNS, and NSBH populations. For BBHs, the authors adopt the LVK's fiducial "Broken Power Law + 2 Peaks" mass model from the GWTC-4.0 population analysis, which captures the observed overdensities at approximately $10\,M_\odot$ and $35\,M_\odot$, in place of the sharp-cutoff power laws used in earlier studies. Mass ratios follow a power law modulated by a low-mass smoothing function; spin magnitudes are drawn from truncated normals, with spins assumed perfectly aligned with the orbital angular momentum (precession neglected). BNS systems use the Peak mass model with a low-spin assumption ($\chi \le 0.05$); NSBH systems follow the GWTC-3 hierarchical framework of Biscoveanu et al. (2022), with the BH mass capped at $20\,M_\odot$ to exclude systems like GW190814. Redshift distributions derive from the Dominik et al. (2013) population synthesis, which the authors cross-validate against the GWTC-4.0 power-law evolution model and find consistent within current uncertainties—an assumption that carries forward the caveat that observational constraints currently extend only to $z \lesssim 1.5$, while the simulated catalog reaches $z \sim 21.5$.

A methodological strength is the use of the full posterior predictive distribution: population hyperparameters are drawn dynamically from the LVK posterior chains rather than fixed at best-fit values, propagating population-inference uncertainties into the synthetic catalog. Sampling uses rejection sampling for the primary mass and a hybrid inverse-transform/rejection scheme for the mass ratio. The synthetic population comprises roughly $4.0 \times 10^5$ BBH, $5.7 \times 10^4$ BNS, and $3.1 \times 10^4$ NSBH events.

Composite lens model and the equivalent luminosity function

The deflector model follows the halo-model approach of Abe et al., comprising five components: the host NFW halo (with mass-dependent ellipticity and concentration from the Diemer relations), a central galaxy with an elliptical Hernquist profile and a stellar mass–halo mass relation calibrated on JWST size–mass measurements, a subhalo population based on the Oguri analytic subhalo mass function (including tidal stripping and dynamical friction), satellite galaxies, and redshift-dependent external shear. A notable approximation is that subhalo lensing is computed with non-truncated elliptical NFW profiles using accretion-epoch masses rather than stripped masses; the authors justify this on the grounds that strong lensing near the Einstein radius is dominated by the central mass distribution, but this remains an approximation whose impact on subhalo-lensed rates is not quantified.

To interface with SL-Hammocks, which was designed for optical surveys, the authors construct an "equivalent luminosity function": since SNR scales as $10\,M_\odot$0 while optical flux scales as $10\,M_\odot$1, they map the joint redshift–SNR distribution into an equivalent absolute-magnitude space via Pogson's equation, with a deliberately conservative intrinsic SNR threshold of 0.1 to retain events that could be strongly magnified. This is presented strictly as an operational bridge; the mapping itself introduces no additional physics but enables reuse of the validated lensing pipeline.

Event rates and detector configurations

Simulations cover four networks—A+ (LIGO Hanford/Livingston 4 km plus AdV+ 3 km plus KAGRA at design sensitivity), ET (triangular, 10 km, ET-D PSD), CE (L-shaped, 40 km, CE1 ASD), and ET+CE—over a ten-year baseline with Earth's rotation properly treated, so that each lensed image experiences a distinct antenna-pattern response at its arrival time. The headline rates for BBH sources are:

Network	Doubles	Quads	Subhalo	Central image ($10\,M_\odot$2)	High-mag $10\,M_\odot$3	Total triggered
ET+CE	395.4	36.1	107.1	20.0	359.3	617.2
ET	263.1	30.8	74.0	10.5	327.5	556.8
CE	260.8	33.6	75.2	11.7	328.7	546.8
A+	1.7	1.1	0.9	0.0	15.9	17.3

BBHs constitute roughly 95% of detected lensed events in the ET+CE configuration, with lensed BNS and NSBH rates of approximately 4.8 and 13.4 per year, respectively. Doublets account for 91.5% of multiple-image events. Three results deserve emphasis. First, the subhalo-lensed population (~110 events per year) is a direct consequence of the composite lens model and exploits the point-like nature of GW sources, which are sensitive to small-scale perturbations that finite-size optical sources wash out; accumulating this sample would constrain the low-mass halo mass function and test the CDM paradigm. Second, the detection of complete three- or five-image systems including the central image (~20 per year for $10\,M_\odot$4) is unique to the GW band, where the central image is not obscured by lens light; such events probe the innermost mass distribution of lens galaxies and can discriminate cuspy (NFW-like) from cored profiles. Third, the A+ network yields only marginal lensed detections, confirming that this science is essentially a 3G endeavor.

Model dependencies

The rate forecasts depend strongly on binary evolution assumptions. Across four StarTrack-based evolution models (Standard, Optimistic CE, Delayed SN, High BH Kicks) at both high and low metallicity, the Optimistic CE model yields roughly three times the Standard-model lensed rate, while the High BH Kicks model produces only about 1/20 of it. The authors conclude that the total intrinsic merger rate is the dominant driver of the lensed population, and that metallicity effects are subdominant to the choice of evolutionary pathway. This is a significant systematic: the headline ET+CE rates should be read as conditional on the low-metallicity Standard model, with an uncertainty band spanning more than an order of magnitude across plausible evolution scenarios.

Regarding mass models, the comparison against the power-law assumptions of Yang et al. (2021) shows that the Broken Power Law + 2 Peaks model is robust: event-rate fluctuations remain within about 10 events per year as the mass range is widened, because approximately 99% of the BBH population lies within $10\,M_\odot$5–$10\,M_\odot$6. Simple power-law models, by contrast, are highly sensitive to the mass cutoff and tend to overestimate rates by extending the high-mass tail. Under matched power-law configurations, the authors' rates come out 18–20% below the earlier predictions, which they attribute to the composite lens model and updated treatment.

Statistical distributions for lens identification

A major deliverable is a set of priors for lens identification algorithms. For the 6,172-event triggered catalog over ten years, the mean source and lens redshifts are $10\,M_\odot$7 and $10\,M_\odot$8, with density peaks at $10\,M_\odot$9 and $35\,M_\odot$0. Magnification bias shifts the high-magnification subsample to higher redshift (median $35\,M_\odot$1 versus 5.20), with tails extending to $35\,M_\odot$2—supporting the use of lensed events as probes of early-universe binary formation. The angular-separation distribution peaks at $35\,M_\odot$3 arcseconds, is dominated by galaxy-scale lenses (83.5%, with ~15% group and 1.5% cluster scale), and host halos contribute 75% versus ~25% from subhalos.

The SNR-ratio and time-delay statistics contain several results with direct observational consequences. For doublets, the late-arriving image is fainter in 88.3% of cases (median $35\,M_\odot$4), and time delays peak at 75.9 days with a median of 65.5 days but a mean of 157 days—implying that complete recovery of a lensed system requires at least six months of continuous monitoring. For quadruplets, the situation inverts: because highly magnified pairs form near caustics, the second and third arrivals are louder than the first in 77.8% and 63.0% of cases respectively (median ratios 1.33 and 1.16), while the fourth image is demagnified in 82.5% of cases and falls below threshold in about 14.4% of systems. The authors draw the operational conclusion that quadruplet identification requires both prospective monitoring and retrospective ("pre-trigger") searches of archival data for the fainter first image. Globally, the median time delay across all multiple-image pairs is approximately 47 days, which the authors propose as a prior peak for pair-identification Bayes factors combining time delay and SNR ratio.

The analysis of selection effects yields a notable and potentially counterintuitive result: near the detection threshold ($35\,M_\odot$5), a "flux ratio anomaly" appears, with late-arriving images louder than early ones in 58.1% of detected doublets and 78.7% of quads—opposite to the intrinsic population trend. This arises because systems whose first image barely exceeds threshold are detectable only if the second image is intrinsically louder. The implication is that identification criteria must be adaptive across SNR regimes; applying intrinsic-distribution priors uniformly would misclassify threshold events. As SNR increases, the ratio distribution converges to its intrinsic form, plateauing at a peak of $35\,M_\odot$6 for doublets. The joint SNR distributions show moderate positive correlation ($35\,M_\odot$7 for doublets), decorrelated by Earth-rotation modulation of the antenna response and the diversity of magnification ratios; approximately 12.2% of doublets exhibit "flux inversion" ($35\,M_\odot$8).

Limitations and open questions

Several limitations are stated or implicit. The redshift distribution relies on Dominik et al. (2013) population synthesis validated only against observations at $35\,M_\odot$9, while the catalog extends to $\chi \le 0.05$0; the high-redshift tail is therefore extrapolative. The subhalo lensing treatment uses non-truncated NFW profiles with accretion-epoch masses, an approximation whose effect on the ~110 per year subhalo rate is not quantified. The spin model assumes perfect alignment and neglects precession, and the NSBH model is inherited from GWTC-3 because GWTC-4.0 provided no new NSBH constraints. The rate forecasts are conditional on the low-metallicity Standard evolution model, with the evolution-model uncertainty spanning more than an order of magnitude. Wave-optics diffraction effects are acknowledged as relevant for subhalo-lensed events but are not modeled in the geometric-optics pipeline. Finally, the paper restricts attention to detectable systems; the companion analyses of identification efficiency and space-borne detector predictions remain forthcoming, leaving open how many of the forecast 617 annual ET+CE triggers would actually be confirmed as lensed by existing Bayesian identification methods.

Conclusion

This work provides the most physically complete mock catalog of strongly lensed GW events to date, unifying galaxy-, group-, and cluster-scale lensing with explicit subhalo contributions and realistic GWTC-4.0-constrained source populations. Its quantitative forecasts—roughly 400 doublets and 36 quadruplets per year for ET+CE, with ~107 subhalo-lensed events, ~20 detectable central images, and ~360 high-magnification events annually—establish concrete benchmarks for 3G-era lensing searches. The derived priors on time delays, SNR ratios, and their SNR-dependent selection biases, particularly the threshold-region flux ratio anomaly and the caustic-dominated brightness ordering of quadruplets, are directly actionable for identification algorithms and observational strategy. The dominant residual uncertainties lie in the binary evolution model and the simplified subhalo lensing treatment, and the catalog's public release enables the community to test identification pipelines against a realistic, complete population of lensed signals.