GWcosmo Pipeline for Cosmology Inference

• GWcosmo pipeline is a computational framework that integrates gravitational-wave data with galaxy and galaxy-cluster catalogs to infer cosmological parameters such as the Hubble constant H0.
• It employs a Bayesian hierarchical model and advanced completeness corrections to build accurate line-of-sight redshift priors from EM surveys, enhancing cosmological inferences.
• The pipeline is modular and extensible, facilitating robust cross-correlation with cosmological backgrounds and future integration of diverse survey data.

The GWcosmo pipeline is a set of methods and computational tools for inferring cosmological parameters, especially the Hubble constant $H_0$, by combining gravitational-wave (GW) observations with electromagnetic (EM) galaxy or galaxy-cluster catalogs. Originating with works by Gray et al., its scope includes both traditional “dark siren” cosmology (where host redshifts are estimated statistically) and more recent extensions incorporating galaxy cluster information, robust catalogue completeness estimates, and cross-correlation with cosmological backgrounds. The pipeline is architected to handle the selection effects, completeness limits, and likelihood formulations that arise in GW cosmology with incomplete or multi-resolution EM datasets, achieving competitive constraints on $H_0$ through hierarchical Bayesian inference and advanced catalogue handling (Datrier et al., 20 Feb 2025, Beirnaert et al., 20 May 2025, Schulze et al., 2023).

1. Core Bayesian Framework and Likelihood Construction

GWcosmo adopts a Bayesian hierarchical model targeting inference on cosmological parameters ($H_0$, extensions to $\Omega_M$, $w_0$, modified gravity). For a set of $N$ GW events with strain data $\{d_i\}$, the central posterior factorizes as:

$$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$

• $p(H_0)$: user-prior over $H_0$
• $H_0$0: single-event GW likelihood, recast from distances to redshifts via luminosity distance $H_0$1
• $H_0$2: line-of-sight redshift prior from EM catalogue(s) $H_0$3
• $H_0$4: selection-function normalization, accounting for cosmology-dependent GW detectability

Marginalizing over the individual source redshifts leads to

$H_0$5

This formalism is adaptable—handling cluster catalogs, mixed catalogs, or hypothetical future catalogs by adjusting $H_0$6 and the completeness model (Beirnaert et al., 20 May 2025).

2. EM Catalog Integration: Galaxy and Cluster Priors

Initially, GWcosmo utilized full-sky galaxy catalogs (e.g., GLADE, GLADE+) to supply $H_0$7 by associating candidate galaxies (with photometric/spectroscopic $H_0$8) to GW localization regions and marginalizing over possible hosts. With the introduction of robust completeness limits, the pipeline discards galaxies beyond a magnitude threshold $H_0$9, determined per-sky-pixel, for assigning host probabilities (Datrier et al., 20 Feb 2025).

Recent advancements adapt the pipeline to leverage galaxy cluster catalogs such as PSZ2 (Planck) and eRASS (eROSITA). Here, a “ClusterCatalogue” class constructs the LOS-$H_0$0 prior as a sum over clusters in each GW localization pixel, weighting by cluster mass and convolving with $H_0$1-uncertainties (typ. Gaussian, width via physical cluster size). For sky areas or redshift ranges where the catalog is incomplete, the missing prior is supplied from a Press–Schechter mass function, normalized by measured completeness curves. The full prior is:

$H_0$2

with $H_0$3 and $H_0$4 built as in (Beirnaert et al., 20 May 2025).

3. Statistical Treatment of Catalogue Completeness

GWcosmo initially used a “median-magnitude” heuristic to estimate catalogue completeness, often discarding a significant fraction of faint galaxies even when well detected. The robust method introduced in (Datrier et al., 20 Feb 2025) operationalizes the statistical test of Rauzy (2001):

• For each galaxy $H_0$5, the rank-ratio statistic $H_0$6 is constructed, where $H_0$7 (resp. $H_0$8) is the count of galaxies with lower $H_0$9 and $\Omega_M$0 (resp. $\Omega_M$1).
• Under the null hypothesis (catalogue complete to $\Omega_M$2), $\Omega_M$3 are i.i.d. uniform$\Omega_M$4.
• The global $\Omega_M$5 is tracked; the threshold $\Omega_M$6 is where $\Omega_M$7 first drops below $\Omega_M$8.
• No galaxy luminosity function is assumed; only spatial stationarity within band/redshift is needed.
• Implementation subsamples pixels for speed (e.g., $\Omega_M$9 galaxies per B-band HEALPix pixel, averaged over 30 runs), incorporates photometric redshift uncertainties, and produces a per-pixel $w_0$0 map.

The pipeline then integrates this limit into the likelihood, discarding galaxies with $w_0$1 and modifying both discrete sum and missing-mass (out-of-catalog) terms.

4. Pipeline Modularity and Workflow

The “gwcosmo” pipeline is structured for extensibility:

• Inputs: GW parameter-estimation posteriors (luminosity distance, sky position, mass, spins), EM catalogs (galaxy or cluster, with magnitude, $w_0$2, mass).
• Pixelation: Catalogs subdivided into HEALPix pixels ($w_0$3 typical).
• Completeness: Per-pixel $w_0$4 map computed and cached.
• Likelihood Evaluation: For each GW event, the pipeline reads the $w_0$5 map, constructs LOS-$w_0$6 prior, and evaluates the GW likelihood, including completeness corrections. For clusters, additional Press–Schechter terms supplement out-of-catalog regions.
• Inference: The core sampling engine is unchanged, except for reading in updated priors and completeness limits; base posterior sampling over $w_0$7 and ancillary parameters proceeds as standard.

Integration with external frameworks—such as using GW_CLASS and MontePython for anisotropy studies or joint GW+CMB cosmology—is handled via standard data flows and output chaining (Schulze et al., 2023).

5. Quantitative Performance and Improvements

Application of these methods has resulted in quantifiable gains in cosmological parameter inference:

• For GLADE+ (B$w_0$8-band), robust completeness cuts from the Rauzy test resulted in an $w_0$9 narrower $N$0 credible interval for dark sirens and $N$1 narrower interval when combining with GW170817, compared to previous median-cut approaches (Datrier et al., 20 Feb 2025).
• Using cluster catalogs (PSZ2 and eRASS) rather than galaxies, the interval width on $N$2 improved from approximately $N$3 (GLADE+) to $N$4 (PSZ2, $N$5 gain) and $N$6 (eRASS, $N$7 gain), attributable to the deeper redshift reach of modern cluster catalogs (Beirnaert et al., 20 May 2025).
• For GWTC-3 BBH events in K-band (GLADE+), both median and robust methods yielded no improvement due to catalog incompleteness at the relevant distances: posteriors reverted to empty-catalogue limits with $N$8 difference (Datrier et al., 20 Feb 2025).

The table below summarizes relevant $N$9 uncertainties (rounded; see original works for details):

Catalog/Band	Method	$\{d_i\}$0 Posterior Width (km/s/Mpc)	Gain vs. Median
GLADE+, $\{d_i\}$1	Median	$\{d_i\}$2	--
GLADE+, $\{d_i\}$3	Robust	$\{d_i\}$4	$\{d_i\}$5
PSZ2	Cluster	$\{d_i\}$6	$\{d_i\}$7
eRASS	Cluster	$\{d_i\}$8	$\{d_i\}$9

6. Selection Effects, Systematics, and Future Prospects

Selection function normalization $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$0 is included in all likelihoods, reflecting the joint sensitivity to changes in $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$1 and event detection probability. GWcosmo pipelines also adjust for sky coverage of catalogs (PSZ2: $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$2, eRASS: $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$3), redshift uncertainty models (e.g., Gaussian for clusters, scatter from physical size or X-ray scaling), and EM survey incompleteness via analytic or empirical completeness curves.

Planned extensions include:

• Combining both galaxy and cluster catalogues such that well-mapped low-$$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$4 galaxies and high-$$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$5 clusters provide joined $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$6 support.
• Adding a term to the LOS–z prior accounting for mergers outside clusters (voids/filaments), as neither catalog necessarily captures all hosts.
• Generalizing completeness methodology to accommodate both faint and bright magnitude limits (e.g., for DESI, LSST).
• Adapting to multi-parameter cosmology inference ($$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$7, $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$8, modified gravity), achievable by replacing $$p(H_0, \{z_i\} | \{d_i\}, \mathcal{C}) \propto p(H_0) \prod_{i=1}^{N} p(d_i|z_i,H_0) p(z_i|\mathcal{C}) \times \frac{1}{S(H_0)}$$9-only grids with multidimensional cosmology grids.

Robustness checks—varying input mass-function parameters, redshift scatter, completeness curves—have consistently shown negligible impact on $p(H_0)$0 posteriors.

A complete end-to-end mock data challenge, with known injected cosmological parameters, is recommended before application to the next-generation observing runs.

7. Integration with Cosmological GW Background Studies

The pipeline can be embedded in broader frameworks, connecting GW event-based inference with studies of cosmological GW backgrounds (CGWB) and their cross-correlation with CMB anisotropies (Schulze et al., 2023). The mathematical structure allows for inclusion of CGWB angular power spectra $p(H_0)$1, cross-spectra with CMB temperature $p(H_0)$2, and computation of projected constraints with joint CMB + GW datasets. A modified version of the CLASS Boltzmann solver (GW_CLASS) computes these observables, and interfacing with MontePython enables joint posterior sampling over both conventional and GW-specific cosmological parameters.

The GWcosmo pipeline thus constitutes a central computational framework in the emerging field of gravitational-wave cosmology, capable of adapting to new catalogs, survey depths, and parameter inference challenges as the capabilities of both GW detectors and EM surveys expand.