Statistical Dark Siren Method in Cosmology

• Statistical dark sirens are gravitational-wave standard sirens without electromagnetic counterparts that assign redshifts probabilistically using galaxy catalogs and large-scale structure.
• They utilize hierarchical inference techniques combining event-level host marginalization, catalog-based likelihoods, and population models to constrain the Hubble constant and other cosmological parameters.
• Implementations integrate gravitational-wave skymaps, redshift priors, completeness corrections, and cross-correlation analyses to derive robust constraints on the universe’s expansion.

Statistical dark sirens are gravitational-wave standard sirens without an electromagnetic counterpart. Because the gravitational-wave measurement provides a posterior on luminosity distance and sky position but not a unique source redshift, cosmological inference proceeds statistically: one marginalizes over candidate host galaxies in a catalog, over a reconstructed large-scale-structure field, or over a population model that links detector-frame observables to source redshift. In this sense, the statistical dark siren method is a family of hierarchical inferences connecting gravitational-wave distance information to redshift information, most often to constrain the Hubble constant $H_0$ and related cosmological parameters (Collaboration et al., 2019, Yu et al., 2022, Cheng et al., 13 Mar 2026).

1. Definition, scope, and subclasses

Bright sirens and dark sirens differ only in how the redshift is supplied. Bright sirens have an electromagnetic counterpart that identifies the host galaxy and its redshift directly. Dark sirens do not. The statistical dark siren method therefore replaces event-by-event host identification with a probabilistic redshift assignment derived from galaxy catalogs, large-scale structure, or population information. In the galaxy-catalog formulation, one builds an event-level spatial prior for host galaxies using a catalog, marginalizes over all potential hosts, and accounts for catalog completeness and selection. In the cross-correlation formulation, one compresses gravitational-wave and galaxy information into tomographic angular overdensity maps and their auto- and cross-power spectra. In the spectral or population formulation, one uses the observed distribution of compact-binary events, together with an astrophysical population model and selection effects, to infer cosmology statistically (Cheng et al., 13 Mar 2026).

Within the catalog-based class, “golden dark sirens” are exceptionally well-localized gravitational-wave events without electromagnetic counterparts, typically at $z \lesssim 0.1$, whose three-dimensional localization contains only one plausible $L_\star$ host galaxy. A practical working definition used in forecasts is that events at $z \le 0.1$ with 90% sky area below about $0.04\,\mathrm{deg}^2$ are golden, since they are expected to contain a single $L_\star$ galaxy in their localization volume. They remain “dark” because the host is not confirmed by a counterpart, but they differ from classical statistical dark sirens by having tiny sky areas, narrow distance posteriors, and much stronger single-event constraints (Benetti et al., 16 Feb 2026).

This classification matters because different realizations of the method emphasize different observables. Galaxy-catalog dark sirens use line-of-sight galaxy information directly. Cross-correlation dark sirens emphasize two-point statistics of the galaxy and gravitational-wave fields. Population-based realizations emphasize the shape of the detected-event distribution and the source-frame population model. These are complementary rather than mutually exclusive descriptions (Cheng et al., 13 Mar 2026).

2. Event-level galaxy-catalog formalism

A standard single-event formulation marginalizes over latent host identity and redshift. For a gravitational-wave event with data $d_{\rm GW}$ and a galaxy catalog $G$, the posterior on $H_0$ can be written as

$$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$

where $z \lesssim 0.1$0 encodes any prior weighting of galaxies, $z \lesssim 0.1$1 is the galaxy redshift PDF, $z \lesssim 0.1$2 is the galaxy sky position, and $z \lesssim 0.1$3 is the gravitational-wave likelihood evaluated at the galaxy position and the luminosity distance implied by $z \lesssim 0.1$4 and $z \lesssim 0.1$5 (Benetti et al., 16 Feb 2026).

The same structure is often written in continuous form. With a discrete spectroscopic catalog,

$z \lesssim 0.1$6

and for precise spectroscopic redshifts one typically takes $z \lesssim 0.1$7, giving

$z \lesssim 0.1$8

Here $z \lesssim 0.1$9 is the three-dimensional gravitational-wave skymap, $L_\star$0 incorporates survey weights or host weights, and $L_\star$1 is the selection normalization (Ballard et al., 2023).

Several implementations compute a galaxy-association probability from an overlap integral. In one form, defining $L_\star$2 and $L_\star$3 for a specific galaxy,

$L_\star$4

and normalizing these over candidate galaxies gives the association probability. The final single-event $L_\star$5 posterior is then a weighted sum over galaxy-specific posteriors, with weights given by association probabilities (Benetti et al., 16 Feb 2026).

The cosmological mapping is supplied by the luminosity-distance relation. In flat $L_\star$6CDM,

$L_\star$7

with

$L_\star$8

At the low redshifts relevant for many nearby dark sirens, the linear approximation $L_\star$9 is often adequate (Benetti et al., 16 Feb 2026).

3. Hierarchical and population-level formulations

A broader statistical dark siren formulation treats the detected catalog as an inhomogeneous Poisson process. If $z \le 0.1$0 denotes event observables in detector space and $z \le 0.1$1 collects cosmological and population parameters, the log-likelihood is

$z \le 0.1$2

where $z \le 0.1$3 is the intensity in observed space. In this framework, cosmology is inferred by adjusting $z \le 0.1$4 so that the observed-space distribution induced by cosmology matches the source-space population model after selection effects are applied (Yu et al., 2022).

For compact-binary populations, one convenient detected-event density is

$z \le 0.1$5

with $z \le 0.1$6 an overall normalization, $z \le 0.1$7 the detection fraction, and $z \le 0.1$8 the astrophysical population model. The Fisher information is then

$z \le 0.1$9

which yields Cramér–Rao forecasts for both cosmological and population parameters (Yu et al., 2022).

In joint dark-siren and population inference, the event-level likelihood is explicitly nested inside a hierarchical model. One implementation writes

$0.04\,\mathrm{deg}^2$0

with normalization by a selection function

$0.04\,\mathrm{deg}^2$1

This formulation was used in a GWTC-4 analysis that jointly inferred cosmology and the compact-binary mass spectrum, including a heavy–black-hole feature at $0.04\,\mathrm{deg}^2$2 (Pierra et al., 6 Jan 2026).

A plausible implication is that the phrase “statistical dark siren method” now covers several inference layers: event-level host marginalization, catalog-level point-process modeling, and joint population–cosmology inference. The common structure is the replacement of a unique host redshift by a modeled redshift distribution conditioned on selection, population, and cosmology (Yu et al., 2022).

4. Redshifts, completeness, weighting, and large-scale-structure priors

The precision and robustness of statistical dark sirens depend strongly on the quality of the galaxy-side prior. Analyses differ in whether they use spectroscopic or photometric redshifts, whether they impose volume-limited cuts, how they model incompleteness, and whether they treat host probability as uniform or proportional to astrophysical tracers such as stellar mass or star-formation rate (Alfradique et al., 24 Mar 2025).

For redshift uncertainties, one line of work compared spectroscopic-like and photometric-like catalogs in four-detector-era simulations. It found a precision gain from spectroscopic-like redshifts compared to photometric-like redshifts, with the greatest improvements for smaller localization areas. The same study found that redshift outliers in realistic photometric catalogs do not introduce bias into the measurement of $0.04\,\mathrm{deg}^2$3, and that at a completeness of 50% the benefit of spectroscopic redshift precision is outweighed by the degradation from incompleteness. In all three scenarios considered there, the inferred $0.04\,\mathrm{deg}^2$4 remained unbiased (Cross-Parkin et al., 25 Feb 2025).

For host weights, the discrete event likelihood is commonly written with weights $0.04\,\mathrm{deg}^2$5 multiplying each galaxy contribution. One study of flux-limited catalogs showed that an unbiased estimate of $0.04\,\mathrm{deg}^2$6 can be obtained when the corrected weighting scheme is applied to a complete or volume-limited catalog, and proposed

$0.04\,\mathrm{deg}^2$7

where $0.04\,\mathrm{deg}^2$8 is a completeness correction. The same work found that equal host probability per galaxy degrades the attainable precision relative to tracer-based weighting (Alfradique et al., 24 Mar 2025).

Catalog incompleteness has motivated several completion strategies. A robust completeness test has been implemented in gwcosmo to estimate the apparent-magnitude completeness limit of a magnitude–redshift sample without prior knowledge of the luminosity function. Applied to GWTC-1 with GLADE and GLADE+, it improved the dark-siren-only inference of $0.04\,\mathrm{deg}^2$9 by 3.4% and 8.6%, respectively, although the same approach yielded no improvement for GWTC-3 with GLADE+ $L_\star$0-band because the catalog provided little or no coverage in that band for the relevant events (Datrier et al., 20 Feb 2025).

More elaborate approaches replace the usual “uniform out-of-catalog” component with an LSS-informed reconstruction. Variance completion models the missing population through a ratio

$L_\star$1

and uses

$L_\star$2

so that the homogeneous completion is modulated by large-scale structure while approximately preserving normalization. In a separate field-level Bayesian reconstruction, the true galaxy Poisson rate $L_\star$3 is inferred jointly with the absolute-magnitude distribution and the spatial density field, yielding a host prior

$L_\star$4

which is then inserted into the standard dark-siren integral (Dalang et al., 2024, Leyde et al., 16 Jul 2025).

These developments indicate a shift from treating catalog incompleteness as a purely radial nuisance to treating it as a field-level inference problem. That shift is strongest in analyses designed for deep but incomplete surveys, where preserving large-scale structure in the missing population becomes as important as modeling the observed galaxies themselves (Leyde et al., 16 Jul 2025).

5. Biases, systematics, and robustness

A central issue in statistical dark siren cosmology is the competition between shrinking statistical errors and coherent systematic effects. In population-based Fisher analyses, a 1% deviation in the astrophysical model can induce a bias $L_\star$5, and this bias is approximately independent of the number of detections. Because statistical errors shrink as $L_\star$6, the same study concluded that beyond $L_\star$7 binary-black-hole detections, astrophysical-model systematics can dominate the $L_\star$8 error budget unless controlled (Yu et al., 2022).

Host localization errors are another distinct failure mode. In a Bayesian catalog-based study with ET+CE forecasts, the precision of $L_\star$9 was found not to be compromised when the number of well-localized dark sirens is significantly below 300, even in the extreme scenario that all the dark sirens are localized incorrectly. As the number exceeds 300, incorrect spatial localizations begin to produce statistically noticeable effects such as slow convergence of the posterior. In the same framework, simulations indicated that incorrect spatial localizations will dominate a systematic error of $d_{\rm GW}$0 if as much as 10% of a sample of 300 well-localized dark sirens are affected by their environments (Zhu et al., 2023).

Galaxy weighting can bias the inference even when the catalog is complete. A systematic study of weighting prescriptions found two distinct effects: the assumption of an incorrect galaxy redshift distribution, and preferentially weighting incorrect host galaxies during the inference. The magnitudes of these biases depend on the number of galaxies along each line of sight, the measurement uncertainty in the gravitational-wave luminosity distance, and correlations in the parameter space of galaxies. The same work proposed hierarchical inference as a diagnostic of incorrectly weighted prescriptions, because it can simultaneously infer the correct weighting scheme and the cosmological parameters (Hanselman et al., 2024).

Waveform and calibration systematics become especially acute for golden dark sirens in the next-generation detector era. For reference signals injected with NRHybSur3dq8 and recovered with IMRPhenomXPNR or SEOBNRv5PHM, one study of 111 nearby binary black holes found that for IMRPhenomXPNR recovery about 80%, 85%, 95% of true hosts lie within the 50%, 67%, 90% three-dimensional credible regions, respectively, whereas for SEOBNRv5PHM the corresponding fractions drop to $d_{\rm GW}$1, $d_{\rm GW}$2, $d_{\rm GW}$3. The most-probable host matches the true host for $d_{\rm GW}$4 of events with IMRPhenomXPNR, versus $d_{\rm GW}$5 with SEOBNRv5PHM. Yet the same study emphasized that imperfect host identification does not automatically imply biased $d_{\rm GW}$6: when several plausible hosts share the same group or cluster redshift, their $d_{\rm GW}$7 ridges nearly coincide and the combined posterior can remain nearly as informative as that of a true golden event (Benetti et al., 16 Feb 2026).

The same work derived accuracy requirements for waveform and calibration errors. If the unfaithfulness is $d_{\rm GW}$8 and the network SNR is $d_{\rm GW}$9, the indistinguishability criterion

$G$0

implies that unfaithfulness must scale as $G$1 to keep biases subdominant at high SNR. For quasi-circular BBHs with $G$2, requiring recovery within the nominal 90% credible region gives

$G$3

which yields $G$4; at $G$5, the threshold is $G$6. Calibration errors must satisfy an analogous criterion, and for representative golden sirens in Cosmic Explorer the median tolerance envelope reaches $G$7 in parts of the CE40 band, far below current Advanced LIGO calibration uncertainties of order $G$8 (Benetti et al., 16 Feb 2026).

6. Implementations, measurements, and future regimes

The method has progressed from single-event demonstrations to multi-event analyses with survey-specific catalogs.

Study	Data	Reported $G$9
GW170814 with DES (Collaboration et al., 2019)	first dark standard siren measurement	$H_0$0
GW190412 with DESI (Ballard et al., 2023)	single BBH dark siren with DESI galaxies	$H_0$1
10 O3 dark sirens with DELVE (Alfradique et al., 2023)	10 well-localized dark sirens	$H_0$2
15 O1–O4a dark sirens (Bom et al., 2024)	catalogue method only	$H_0$3

These measurements remain prior-sensitive and catalog-dependent at current precision, but they establish the observational workflow: ingest a three-dimensional gravitational-wave skymap, combine it with galaxy redshift information and completeness modeling, evaluate a per-event likelihood on an $H_0$4 grid or within a hierarchical sampler, and then combine events multiplicatively (Collaboration et al., 2019, Bom et al., 2024).

Beyond direct host marginalization, joint population–cosmology analyses have begun to use dark sirens together with spectral information. Using 142 GWTC-4 compact-binary events with false alarm rate smaller than $H_0$5 and the GLADE+ catalog, one study reported

$H_0$6

from the dark siren method, and attributed an approximately 36.2% reduction in uncertainty to the inclusion of a heavy–black–hole mass feature at $H_0$7 (Pierra et al., 6 Jan 2026).

Future regimes split naturally into golden-event analyses and field-level cross-correlation analyses. For golden dark sirens, single-event $H_0$8 widths scale roughly as $H_0$9, and for optimally oriented golden sirens the error propagation in one next-generation study gave $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$0 at $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$1–$$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$2 for chirp mass $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$3 when peculiar velocities are included. The same study emphasized that forecasts of order $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$4 with A+ golden sirens and order $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$5 with a CE40+CE20+ET network in two years will only be realized if waveform and calibration systematics are reduced with roughly $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$6 scaling (Benetti et al., 16 Feb 2026).

For cross-correlation dark sirens, a unified harmonic framework treats galaxies and gravitational-wave sources as biased Poisson tracers of the matter field and performs inference from tomographic auto- and cross-spectra rather than explicit host assignment. In a forecast with 2 Einstein Telescopes and 1 Cosmic Explorer, the gravitational-wave–galaxy cross-correlation part alone was found capable of jointly measuring $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$7 and $$P(H_0 \mid d_{\rm GW}, G) \propto p(H_0)\sum_{g\in G} p_g \int dz\, p(z\mid g)\, p\!\left(d_{\rm GW}\mid \Omega_g, d_L[z;H_0,\Omega]\right),$$8 to 1% and 5% precision with two years of data. That result was explicitly contrasted with current-detector prospects, for which the same study remained pessimistic because the method implicitly requires large-number statistics (Cheng et al., 13 Mar 2026).

Across these formulations, the statistical dark siren method is therefore best understood not as a single algorithm but as a hierarchy of related inferences. At one end are event-by-event catalog sums over candidate hosts; at the other are point-process, spectral, and cross-correlation analyses that trade explicit host assignment for population or field statistics. What unifies them is the same cosmological task: reconstructing redshift information probabilistically from incomplete and structured data, then combining it with gravitational-wave luminosity distances to infer the expansion history.