gwcosmo: Bayesian Dark Siren Cosmology

• gwcosmo is a hierarchical Bayesian inference pipeline for gravitational-wave cosmology that jointly samples cosmological parameters and compact-binary population hyperparameters using dark siren events.
• The framework combines gravitational-wave posterior samples with selection effects and redshift priors from galaxy or cluster catalogues to mitigate biases in H0 and Omega_m estimation.
• A GPU-accelerated implementation of gwcosmo achieves dramatic speed-ups, allowing for efficient analyses of large event catalogs in preparation for O5-era gravitational-wave observations.

gwcosmo is a hierarchical Bayesian inference pipeline for gravitational-wave standard-siren cosmology and compact-binary population inference, with particular emphasis on dark sirens, i.e. events without identified electromagnetic counterparts. In this setting, gwcosmo combines gravitational-wave posterior samples, selection effects, and a statistical redshift prior derived from galaxy catalogues, galaxy cluster catalogues, or population-based redshift models, and returns posteriors on cosmological and astrophysical hyperparameters such as $H_0$, $\Omega_m$, merger-rate evolution, and source-mass distributions. The pipeline was developed across the Gray et al. 2019–2023 sequence and has since been extended to joint population–cosmology inference, modified gravitational-wave propagation, deep photometric catalogues, galaxy-cluster tracers, and GPU-accelerated large-catalog analyses (Gray et al., 2023, Papadopoulos et al., 22 May 2026, Beirnaert et al., 20 May 2025).

1. Origins, scope, and position in gravitational-wave cosmology

gwcosmo was originally developed for dark standard sirens with galaxy catalogues and later generalized into a broader framework for joint cosmological and population inference. Earlier gwcosmo applications treated dark sirens with galaxy catalogues and tackled galaxy-catalogue incompleteness via pixelated completeness corrections; Gray et al. (2023) then extended gwcosmo to joint population and cosmology inference with galaxies, and Chen et al. (2023) used gwcosmo to test modified gravitational-wave propagation with dark sirens (Beirnaert et al., 20 May 2025). In the GWTC-3 reanalysis, the enhanced pipeline was described as a “significantly enhanced” version of the Python package GWCOSMO that allows joint estimation of cosmological and compact-binary population parameters, specifically to avoid bias from fixing unknown population distributions during an $H_0$ inference (Gray et al., 2023).

Within the current methodological landscape, gwcosmo is one of several state-of-the-art pipelines for gravitational-wave cosmology. The GPU-acceleration study explicitly places it alongside icarogw and CHIMERA, while distinguishing gwcosmo by its hierarchical Bayesian structure, KDE-based treatment of gravitational-wave posteriors, and ability to incorporate galaxy catalogues, population models, and selection effects (Papadopoulos et al., 22 May 2026). Previous versions were used by the LIGO–Virgo–KAGRA Collaboration for GWTC-4.0 cosmological analyses, and later work adapted the same framework to the Dark Energy Survey Year 6 Gold catalogue, to galaxy cluster catalogues such as PSZ2 and eRASS, and to host-galaxy model comparison in simulated O5-like data (Papadopoulos et al., 22 May 2026, McMahon et al., 4 Feb 2026, Beirnaert et al., 20 May 2025, Li et al., 21 Aug 2025).

The pipeline’s scope is therefore broader than a single $H_0$ estimator. In the papers surveyed here, gwcosmo functions as a unified engine for dark-siren cosmology, spectral-siren inference from mass-spectrum features, joint population hyperparameter estimation, catalog-completeness modeling, and parameterized tests of departures from General Relativity affecting gravitational-wave propagation (Chen et al., 2023, Chen et al., 12 Jun 2026).

2. Hierarchical Bayesian formulation

The central mathematical object in gwcosmo is a population-level hierarchical posterior over cosmological and source-population hyperparameters. In the GPU implementation, letting $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$ denote cosmological and population hyperparameters, the posterior is written schematically as

$$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$

with the selection function

$$\alpha(\Theta)=\int p(D_{\rm GW}\mid \theta,\Theta,I)\,p(\theta\mid \Theta,I)\,d\theta.$$

Under a scale-free Jeffreys prior on the detected rate, the explicit cosmology dependence of $p(N_{\rm det}\mid\Theta)$ can be neglected, leaving the usual $\alpha(\Theta)^{-N_{\rm det}}$ normalization (Papadopoulos et al., 22 May 2026).

For dark sirens, gwcosmo separates sky position and redshift explicitly. The sky is pixelized, typically with HEALPix, and the event likelihood becomes a sum over pixels with a line-of-sight redshift prior $p(z\mid \Omega_j,\Theta)$. At the event level, gravitational-wave parameter-estimation samples are reassigned to pixels and converted from luminosity distance to redshift using the trial cosmology. gwcosmo then constructs a kernel density estimate over redshift in each pixel and reweights posterior samples by the ratio between the current population prior and the original parameter-estimation prior,

$\Omega_m$0

thereby implementing hierarchical importance sampling directly on the posterior samples (Papadopoulos et al., 22 May 2026).

A major development in the GWTC-3 analysis was that gwcosmo no longer required the compact-binary mass distribution to be fixed when galaxy information was used. Instead, it jointly sampled $\Omega_m$1, mass-distribution hyperparameters, and merger-rate evolution. This change was introduced precisely because informative priors on the mass distribution and redshift evolution can impact the inferred posterior constraints on the Hubble constant, and fixing them risks bias (Gray et al., 2023).

3. Line-of-sight redshift priors, tracer catalogues, and incompleteness

The distinguishing practical feature of gwcosmo is its treatment of the line-of-sight redshift prior. In the galaxy-catalogue formulation, the prior is split into in-catalogue and out-of-catalogue contributions. The in-catalogue term is built from galaxies in a given sky pixel, each carrying a redshift PDF and an optional host weight. In the DES Y6 Gold implementation, each galaxy contributes a Gaussian photo-$\Omega_m$2 PDF, and the host weight is

$\Omega_m$3

with $\Omega_m$4 in the fiducial luminosity-weighted analysis, so that host probability scales with luminosity (McMahon et al., 4 Feb 2026).

The same framework supports explicit incompleteness corrections. In the DES Y6 analysis, gwcosmo pre-computes a HEALPix line-of-sight prior from the effective number density of catalogued galaxies plus an out-of-catalogue term modeled as uniform in comoving volume and isotropic. That study identified significant features in the galaxy redshift distribution which can introduce biases, and therefore restricted the catalogue to $\Omega_m$5 to maintain consistency with a uniform in comoving volume galaxy distribution (McMahon et al., 4 Feb 2026). In a complementary methodological study, a robust test of apparent-magnitude completeness was implemented in gwcosmo, replacing the earlier median-magnitude threshold. For GWTC-1 this yielded a $\Omega_m$6 improvement on the inference of $\Omega_m$7 using dark sirens only with GLADE $\Omega_m$8-band data and an $\Omega_m$9 improvement with GLADE+, while the same method showed no improvement for GWTC-3 with GLADE+ $H_0$0-band because the catalogue provided little or no coverage in that band for the relevant events (Datrier et al., 20 Feb 2025).

gwcosmo has also been adapted to large-scale-structure-informed incompleteness corrections. The “variance completion” construction introduces a ratio

$H_0$1

which multiplies the standard homogeneous out-of-catalogue prior,

$H_0$2

This makes the missing-galaxy term sensitive to expected large-scale-structure variance while remaining compatible with existing gwcosmo LOS machinery (Dalang et al., 2024).

A further generalization replaces galaxies with galaxy clusters as redshift tracers. In the PSZ2 and eRASS adaptation, the only structural change to gwcosmo is that the in-catalogue term uses cluster-based redshift PDFs and mass weights, while the out-of-catalogue term is tuned using cluster selection functions and a Press–Schechter mass function rather than galaxy luminosity functions. Beirnaert et al. explicitly assume that the gravitational-wave merger rate density is proportional to cluster mass and use clusters as tracers of the matter field rather than as strict host centers (Beirnaert et al., 20 May 2025).

4. Population modeling and theoretical extensions

gwcosmo’s cosmological inference is tightly coupled to population modeling because detector-frame masses satisfy $H_0$3. The pipeline has therefore been used with several mass-distribution families. The GWTC-3 joint inference used the Power Law + Peak model for BBHs together with parametric rate evolution (Gray et al., 2023). Later analyses introduced richer models such as FullPop-4.0 for multi-class compact binaries and MultiPeak for BBHs in O5-like mock studies (McMahon et al., 4 Feb 2026, Papadopoulos et al., 22 May 2026). In the DES Y6 GWTC-4 analysis, gwcosmo jointly inferred cosmological and GW population parameters with the 19-parameter FullPop-4.0 mass model and a Madau–Dickinson-like merger-rate evolution (McMahon et al., 4 Feb 2026).

The framework has also been extended beyond standard $H_0$4CDM inference. In the modified-propagation study, gwcosmo was generalized to sample parameterized deviations from General Relativity through a modified gravitational-wave luminosity distance. One implementation adopts the Belgacem parameterization

$H_0$5

while others use a Horndeski-inspired $H_0$6 parameterization and an extra-dimension model with parameters $H_0$7 (Chen et al., 2023). The GWTC-4 extra-dimension analysis then applied the dark-siren method with GLADE+ and obtained

$H_0$8

for a prior $H_0$9 and $H_0$0, while also showing that the inferred constraint on $H_0$1 depends sensitively on the assumed prior range of $H_0$2 (Chen et al., 12 Jun 2026).

A different extension uses gwcosmo as a Bayesian model-comparison engine for host-galaxy weighting models. In a simulated O5-like LVK scenario with a mock spectroscopic catalogue from MICECATv2, the analysis compared Bayes factors among $H_0$3-band, $H_0$4-band, and uniform host weighting. The Bayes factors showed a minor preference for the true model for a $H_0$5-detection case, and a decisive preference for the true model over the uniform model for a $H_0$6-detection case, with the result strongly driven by a small number of well-localised events (Li et al., 21 Aug 2025). This suggests that gwcosmo is not restricted to cosmological parameter estimation; it can also be used to learn the host–binary connection from the data themselves.

5. Computational implementation and scalability

A major recent development is the GPU-accelerated rewrite of gwcosmo. The upgraded implementation represents posterior samples as a dense three-dimensional tensor of shape $H_0$7, padded where necessary, and evaluates the full gravitational-wave catalogue in parallel at each likelihood call. Tensor operations are implemented in PyTorch, while the most expensive pieces—vectorized KDE construction and pixel-wise interpolation and summation—use custom Numba CUDA kernels (Papadopoulos et al., 22 May 2026).

This reorganization changed the practical scale of the pipeline. The 2026 scalability paper reports a speed-up of 1000 times over the previous version and states that the new implementation facilitates analyses of O5-like numbers of GW events on wall-clock timescales of hours (Papadopoulos et al., 22 May 2026). For a mock catalog of 2000 simulated events, a single full likelihood evaluation dropped from approximately $H_0$8 ms on the legacy CPU code to $H_0$9 s on an H100 GPU, or $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$0 s when posterior samples were capped at 2500 samples per pixel. In the GWTC-4.0 comparison, the CPU gwcosmo analysis on 32 cores required 32 days wall time, whereas the GPU implementation required 96 hours with all samples and 15 hours with the 2500-sample cap; the corresponding energy consumption dropped from $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$1 kWh to $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$2 kWh for the downsampled H100 run (Papadopoulos et al., 22 May 2026).

The same study undertook explicit validation and systematic tests. Posteriors from the GPU code and the legacy CPU code were found to be in excellent agreement for the 141-event GWTC-4.0 dark-siren analysis using FullPop-4.0 and GLADE+ K-band weighting, with marginal Kullback–Leibler divergences below an internally estimated sampling-noise bound (Papadopoulos et al., 22 May 2026). Moderate downsampling to 2500 samples per pixel, single precision, and robust alternatives to Scott’s-rule KDE bandwidth were found to be safe, whereas too-coarse redshift grids ($\Theta=\{\Theta_{\rm cosmo},\Lambda\}$3) or overly strict effective-sample cuts could bias the $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$4 posterior (Papadopoulos et al., 22 May 2026).

6. Representative scientific results

The scientific output of gwcosmo has spanned nearby bright-siren combinations, dark sirens with shallow and deep catalogues, and new tracer classes. In the GWTC-3 reanalysis that first enabled joint population–cosmology inference with galaxies, the combined posterior from BBHs, NSBHs, dark-siren galaxy information, and GW170817 yielded

$\Theta=\{\Theta_{\rm cosmo},\Lambda\}$5

reported as the maximum a posteriori probability and 68% highest density interval (Gray et al., 2023). The point of that analysis was not only the numerical value but also the demonstration that fixing the population model can lead to artificially tight constraints.

In the DES Y6 Gold analysis using 142 compact binary coalescences from GWTC-4.0, gwcosmo jointly inferred cosmological and GW population parameters from a deep photometric catalogue. Restricting the galaxy catalogue to $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$6, the study obtained

$\Theta=\{\Theta_{\rm cosmo},\Lambda\}$7

from dark sirens alone and

$\Theta=\{\Theta_{\rm cosmo},\Lambda\}$8

when combined with the bright siren GW170817 (McMahon et al., 4 Feb 2026). That work also showed that extending the catalogue to higher redshift without controlling redshift-distribution features could shift the $\Theta=\{\Theta_{\rm cosmo},\Lambda\}$9 posterior systematically downward.

The cluster-catalogue adaptation established that gwcosmo’s line-of-sight prior abstraction can be retargeted to more distant tracers. Using GWTC-3 binary black holes together with the PSZ2 and eRASS cluster catalogues, Beirnaert et al. obtained

$$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$0

respectively, corresponding to precision improvements of $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$1 and $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$2 over the traditional galaxy-catalogue result (Beirnaert et al., 20 May 2025). An important qualitative result of that study was that many individual events acquired visible peaks in their $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$3 posteriors attributable to in-catalogue clusters.

On the purely methodological side, the GPU paper used a 2000-event O5-like spectral siren mock catalog with the MultiPeak mass model and recovered

$$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$4

demonstrating approximately $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$5 precision under its assumptions and, more importantly, the practical feasibility of population-level analyses with $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$6 events (Papadopoulos et al., 22 May 2026).

7. Limitations, systematics, and future directions

The gwcosmo literature emphasizes that the main remaining limitations are now as much astrophysical and catalog-related as numerical. Galaxy catalogues are incomplete and anisotropic; host weighting by luminosity is convenient but uncertain; photometric-redshift systematics can introduce artificial features in the redshift distribution; and in cluster-based analyses the assumption that cluster bias and GW source bias are equal is explicitly identified as a simplification (McMahon et al., 4 Feb 2026, Beirnaert et al., 20 May 2025). In the DES Y6 analysis, the authors advocated explicit diagnostics of $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$7 and conservative cuts where deviations exceed $$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$8 (McMahon et al., 4 Feb 2026). In the cluster analysis, they suggested that a realistic future gwcosmo implementation should include a third LOS-$$p(\Theta \mid \{x_{\rm GW}\}, \{D_{\rm GW}\}, I) \propto p(\Theta \mid I)\, p(N_{\rm det}\mid \Theta,I)\, \alpha(\Theta)^{-N_{\rm det}} \prod_{i=1}^{N_{\rm det}} \int p(x_{{\rm GW},i}\mid \theta,\Theta,I)\, p(\theta\mid \Theta,I)\, d\theta,$$9 term for mergers hosted in galaxies not associated with clusters (Beirnaert et al., 20 May 2025).

Population uncertainty remains structurally important. The 2023 GWTC-3 study made the case that dark-siren analyses necessarily depend on the mass distribution of compact objects and the evolution of their merger rate with redshift, and that informative priors on these quantities will impact the inferred posterior constraints on $$\alpha(\Theta)=\int p(D_{\rm GW}\mid \theta,\Theta,I)\,p(\theta\mid \Theta,I)\,d\theta.$$0 (Gray et al., 2023). The host-weighting study sharpened that point by showing that host weighting and rate evolution are entangled in the full-sky LOS prior, implying that future cosmological analyses may need to marginalize over classes of host models rather than fixing a single prescription (Li et al., 21 Aug 2025).

The computational side is more favorable. Internal approximations such as moderate posterior downsampling, single precision, and alternative robust KDE bandwidth estimators have now been validated for current catalog sizes, and the vectorized GPU architecture appears capable of handling O5-like and eventually much larger data sets (Papadopoulos et al., 22 May 2026). The papers surveyed here repeatedly point toward the same next steps: joint galaxy-plus-cluster LOS priors, richer cosmological models including $$\alpha(\Theta)=\int p(D_{\rm GW}\mid \theta,\Theta,I)\,p(\theta\mid \Theta,I)\,d\theta.$$1, $$\alpha(\Theta)=\int p(D_{\rm GW}\mid \theta,\Theta,I)\,p(\theta\mid \Theta,I)\,d\theta.$$2, and modified gravity parameters, deeper spectroscopic and photometric surveys, and third-generation detectors whose event counts and localization accuracy will turn gwcosmo from a precision-limited dark-siren framework into a high-statistics cosmological pipeline (Beirnaert et al., 20 May 2025, Papadopoulos et al., 22 May 2026).