Population Synthesis Across Surveys and Detectors

A population model is not tested by placing several observed histograms beside one another. It is tested by passing one latent population through every survey exposure, detector response, association rule, and catalog threshold, while representing a source seen twice only once. The central object is therefore a selection-aware marked point process. This page develops that object, supplies a synthetic recovery benchmark, and states where present black-hole catalogs cannot yet support a common-channel verdict. Instrument-specific gravitational-wave likelihood methods are outside its scope.

Helpful background. Black-Hole Geometry and Mechanics supplies mass and spin parameters; Relativistic Matter and Compact Objects supplies source populations and formation channels; Black-Hole Astrophysics supplies emission and environmental models; Gravitational-Wave Astronomy supplies event posteriors and detector selection; Relativistic Cosmology supplies source-frame rates and distance–redshift transformations.

The population is a marked point process

Let $\theta$ contain source-frame masses and spins, redshift, formation channel, environment, and messenger-emission parameters. Hyperparameters $\Lambda$ describe their normalized population density $p(\theta\mid\Lambda)$. Here $R$ is the expected number of parent events over the modeled observing time and cosmic volume, before selection; it is not a rate per unit volume. A detection mark $s$ records a mutually exclusive outcome such as gravitational-wave only, electromagnetic only, neutrino only, or a named overlap. For unique events $j$, the inhomogeneous Poisson likelihood is

$$\mathcal L(\Lambda,R) = \exp\!\left[-\sum_s N_s(\Lambda,R)\right] \prod_{j\in{\rm unique}} R\!\int p(d_j,s_j\mid\theta)\, p(\theta\mid\Lambda)\,d\theta .$$

The expected count is $N_s=R\int P(s\mid\theta)p(\theta\mid\Lambda)d\theta$, summing only detected marks. Observing time and the cosmic-volume distribution enter the parent intensity; sky coverage, instrument duty cycle, and detection/classification enter $P(s\mid\theta)$. An astrophysical rate model must first be integrated with $dV_c$ and $dt_{\rm obs}/(1+z)$ to construct that intensity. The event factor uses the joint $p(d_j,s_j\mid\theta)$, not an unconditioned data likelihood multiplied again by a detection probability from those same data. If gravitational-wave and electromagnetic selections are conditionally independent given $\theta$, the probability of detection by their union is

$$p_{\rm union}(\theta) =p_{\rm GW}(\theta)+p_{\rm EM}(\theta) -p_{\rm GW}(\theta)p_{\rm EM}(\theta).$$

That last term is not optional bookkeeping: it prevents the shared population from being counted twice. Conditional independence is itself an approximation. Correlated weather, follow-up decisions, sky visibility, or an electromagnetic trigger targeted because of a gravitational-wave alert require a joint injection-and-recovery estimate of $p_{\rm union}$.

This construction reduces to the familiar selection-corrected hierarchical likelihood for a single catalog. Its derivation and the distinction between the detected and astrophysical distributions are developed in Mandel, Farr, and Gair 2019, §§2–4, pp. 1087–1092. A nondetection contributes through the expected count; deleting it changes the inferred rate and often the shape.

Posterior recycling, frames, and release-specific inputs

If event $i$ was sampled under interim prior $\pi_i(\theta)$, its population factor can be estimated from $S_i$ posterior draws:

$$\int \mathcal L_i(d_i\mid\theta)p(\theta\mid\Lambda)\,d\theta \ \propto\ \frac{1}{S_i}\sum_{q=1}^{S_i} \frac{p(\theta_{iq}\mid\Lambda)}{\pi_i(\theta_{iq})}.$$

The proportionality constant cancels only in comparisons where the event evidence is irrelevant. Report effective sample size and the largest importance weight; a few dominant draws indicate that fresh event-level sampling is needed. Convert every catalog to one frame before recycling: $m_{\rm det}=(1+z)m_{\rm src}$, and a comoving source-frame rate carries the detector-time factor $(1+z)^{-1}$. A reference realization may use flat $\Lambda$CDM with $H_0=67.4\ {\rm km\,s^{-1}\,Mpc^{-1}}$ and $\Omega_m=0.315$, but those values and their uncertainty must be varied rather than silently inherited Planck Collaboration 2020, §7.4 and Table 2.

For a cutoff analysis, freeze GWOSC GWTC-5.0 data release version 6, published 2026-05-26, DOI 10.7935/bk00-6a89. The associated population analysis uses 267 compact-binary candidates with false-alarm rate below $1\,{\rm yr}^{-1}$, source posteriors, and detectable simulated signals; that sample is not interchangeable with every candidate appearing in the broader catalog LIGO–Virgo–KAGRA 2026, §§2.1–3, Table 2, and Appendix A.2.

Freeze IceCube IceTracks-DR2, covering 2008-04-06 through 2022-05-23, at Dataverse DOI 10.7910/DVN/MMIIZA. Its calibration and event selection change after 2010-06-01, and its binned response functions and effective areas define the public reproduction boundary IceCube Collaboration 2026, public-release “Introduction” and “Data release”. For an electromagnetic cross-check, Cabrera, Palmese, and Fishbach used 76 GWTC-3 binary-black-hole events with explicit ZTF flare selection and found a 90% upper bound of 3.1% on the detectable bright-flare fraction Cabrera, Palmese, and Fishbach 2026, §§3.1, 3.3, and 4.1. That likelihood cannot be appended to GWTC-5.0 without reconciling event overlap and selection. Cosmic-ray terms similarly require the released directional exposure, rigidity-dependent deflection model, energy calibration, and source horizon; arrival directions alone are not a source catalog.

A synthetic recovery benchmark

The benchmark below tests the most failure-prone operation before any real-catalog fit: recovering a known mixture weight through two complementary selections with shared detections. The latent scalar $x$ is drawn from two Gaussian channels with true weight $w=0.65$. Logistic gravitational-wave and electromagnetic selections act on the same 1,200 sources. A 96-node Gauss–Hermite calculation independently benchmarks the selection integral; repeated injection sets test its Monte Carlo error.

The pinned environment is CPython 3.9.6, NumPy 2.0.2, no other scientific packages, PYTHONHASHSEED=0, and random seed 20260818. This deliberately simpler benchmark conditions on the detected count, assumes $x$ is measured exactly, and uses a uniform prior on $0.05\leq w\leq0.95$. The known selection factor for each observed $x$ is independent of $w$ and cancels from the event term, while the population normalization $\alpha(w)^{-N_{\rm obs}}$ remains. It tests mixture recovery, not the full Poisson rate inference. After saving the following source as population_benchmark.py, the program writes population_benchmark.csv and prints the unique/shared counts, quadrature estimate, and central credible interval.

import csv
import numpy as np

SEED, WTRUE, NDRAW = 20260818, 0.65, 1200
LEVELS, REPEATS = (2000, 8000, 32000, 128000), 24

def probabilities(x):
    p_gw = 1.0 / (1.0 + np.exp(-(x + 0.2) / 0.65))
    p_em = 1.0 / (1.0 + np.exp((x - 0.3) / 0.75))
    return p_gw, p_em, p_gw + p_em - p_gw * p_em

def density(x, mean, sigma):
    z = (x - mean) / sigma
    return np.exp(-0.5 * z * z) / (np.sqrt(2.0 * np.pi) * sigma)

rng = np.random.default_rng(SEED)
channel = rng.random(NDRAW) < WTRUE
x = np.where(channel, rng.normal(-1.0, 0.55, NDRAW),
             rng.normal(1.2, 0.70, NDRAW))
p_gw, p_em, _ = probabilities(x)
gw = rng.random(NDRAW) < p_gw
em = rng.random(NDRAW) < p_em
observed, x_obs = gw | em, x[gw | em]

f0 = density(x_obs, -1.0, 0.55)
f1 = density(x_obs, 1.2, 0.70)
grid = np.linspace(0.05, 0.95, 901)
event_term = np.log(grid[:, None] * f0 + (1.0 - grid[:, None]) * f1).sum(axis=1)

nodes, weights = np.polynomial.hermite.hermgauss(96)
q0 = probabilities(-1.0 + np.sqrt(2.0) * 0.55 * nodes)[2]
q1 = probabilities(1.2 + np.sqrt(2.0) * 0.70 * nodes)[2]
alpha0 = np.sum(weights * q0) / np.sqrt(np.pi)
alpha1 = np.sum(weights * q1) / np.sqrt(np.pi)
alpha_ref = WTRUE * alpha0 + (1.0 - WTRUE) * alpha1

logpost_ref = event_term - len(x_obs) * np.log(grid * alpha0 + (1.0 - grid) * alpha1)
post_ref = np.exp(logpost_ref - logpost_ref.max())
post_ref /= post_ref.sum()
cdf_ref = np.cumsum(post_ref)
interval_ref = grid[np.searchsorted(cdf_ref, (0.05, 0.95))]

rows = []
for n_inj in LEVELS:
    errors, estimates = [], []
    for repeat in range(REPEATS):
        trial = np.random.default_rng(SEED + 10000 * (repeat + 1) + n_inj)
        a0 = probabilities(trial.normal(-1.0, 0.55, n_inj))[2].mean()
        a1 = probabilities(trial.normal(1.2, 0.70, n_inj))[2].mean()
        alpha = grid * a0 + (1.0 - grid) * a1
        estimates.append(grid[np.argmax(event_term - len(x_obs) * np.log(alpha))])
        error = WTRUE * a0 + (1.0 - WTRUE) * a1 - alpha_ref
        errors.append(error * error)
    rows.append((n_inj, np.median(estimates), np.sqrt(np.mean(errors))))

with open("population_benchmark.csv", "w", newline="") as stream:
    writer = csv.writer(stream)
    writer.writerow(("N_inj_per_channel", "w_hat", "selection_RMSE"))
    writer.writerows(rows)
print(int(gw.sum()), int(em.sum()), int((gw & em).sum()), int(observed.sum()))
print("quadrature estimate and central 90% interval:",
      grid[np.argmax(logpost_ref)], interval_ref)
print(rows)

Run it in a clean directory:

python3 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.0.2
PYTHONHASHSEED=0 python3 population_benchmark.py

The program emits 570 gravitational-wave detections, 754 electromagnetic detections, 246 shared detections, and 1,078 unique detections. Using the quadrature selection rather than injection noise gives $\hat w=0.672$ and a central 90% grid interval $0.646\leq w\leq0.696$, which contains the injected truth. The convergence record is:

Injection resolution $N$ per channel	Mixture result $\hat w$	Selection error RMSE
2,000	0.672	6.56e-4
8,000	0.672	2.78e-4
32,000	0.672	1.72e-4
128,000	0.672	7.51e-5

The stable mixture estimate says only that selection Monte Carlo noise is subdominant here; its displacement from 0.65 is finite-catalog fluctuation. Repeat the benchmark with correlated detection draws, shifted channel means, and a third messenger. Recovery should fail visibly if the analysis keeps the product selection while the simulator introduces correlated selection.

Results and limitations through 2026-08-18

Results after 2026-08-18 are out of scope. The matrix separates released inputs from channel interpretations.

Evidence or result	Population claim it can test	Scope and limitations	Independent check
GWTC-5.0 version 6 and its 267-event population sample	Mass, spin, redshift, rate, and subpopulation structure of merging compact binaries	Threshold, waveform family, injections, posterior recycling, and assumed evolution remain model inputs	Recover the reported rate with the released selection injections
76 GWTC-3 events plus selected ZTF flares	Fraction of detected mergers producing bright AGN flares	Applies to the adopted flare cuts and coverage; no confirmed one-to-one population identity follows	Refit after removing shared candidates and varying delay and magnitude cuts
IceTracks-DR2, 2008–2022	Steady or transient neutrino-source population under its public response	A submitted data-release paper and binned response do not select a unique black-hole source class	Reproduce the release benchmark before stacking a black-hole catalog
Synthetic two-channel program above	Selection and overlap implementation recovers known $w$	One scalar, Gaussian channels, known measurements, and conditionally independent selections are deliberately simplified	Compare injections to Gauss–Hermite quadrature and introduce correlated selection

The released catalogs and null constraints are durable. A universal latent population, particular formation-channel weights, and claimed cross-messenger tension remain unresolved because emission mappings, selection completeness, overlap, and cosmic evolution differ. No current evidence establishes that one stationary source population explains electromagnetic, gravitational-wave, neutrino, and cosmic-ray catalogs simultaneously.

Competing positions, parameter regimes, and failure modes

Credible competing positions are: one latent population with messenger-dependent emission; several physically distinct formation channels; apparent tension created by inconsistent selection or cosmology; double-counted catalog overlap; and redshift or environmental evolution invalidating a stationary model. The important parameter regime axes are low versus high redshift, isolated versus dense-environment formation, prompt versus delayed emission, weak versus strong beaming, and sparse versus selection-complete follow-up.

Failures have recognizable signatures. Duplicate source identifiers make the inferred rate high and uncertainty too small. Treating shared events as independent multiplies the same astrophysical draw twice. Omitting nondetections favors bright emitters. Mixing detector- and source-frame masses deforms redshift evolution. Inconsistent $H_0$ values shift volumes and rates. Posterior recycling fails when effective sample size collapses. Too few injections make the hyperposterior track the random seed. A mismatch between predicted and observed overlap counts is especially diagnostic: it can expose an incorrect joint selection even when every marginal histogram looks acceptable.

Reproducing the result and testing its robustness

A real analysis must publish the exact catalog IDs and dates, file checksums, event posteriors and interim priors, GW injections, neutrino effective areas, electromagnetic cadence and masks, cosmic-ray exposure, association windows, cosmology, source/detector frame transformations, overlap table, random seeds, environment lock, and one command rebuilding rates, hyperposteriors, evidences, and posterior-predictive distributions. Report statistical intervals separately from calibration, population-model, emission, cosmology, and selection systematics. Each catalog must pass its own rate and distribution benchmark before the joint likelihood is evaluated.

What could change the conclusion. The interpretation could change when any participating catalog or selection release changes, a previously unknown overlap is identified, source evolution or cosmology is revised materially, importance weights become unstable, or a catalog-specific or joint posterior-predictive check fails despite adequate injection convergence.

Formulation-choice diagnostic

The benchmark reports 570 gravitational-wave detections and 754 electromagnetic detections, including 246 sources present in both. How many unique sources enter the event product, and when is multiplication of the two catalog likelihoods valid?

Solution

Inclusion–exclusion gives

$$N_{\rm unique}=570+754-246=1078.$$

The 246 shared sources each require one joint event block, with both data products conditioned on the same latent $\theta$. Multiplying two catalog likelihoods as if they contained disjoint sources would create 1,324 independent draws and double count every overlap. A product is valid only after partitioning the catalogs into mutually exclusive marks, or when independence and zero overlap have been established. The selection normalization must use the same partition: under conditional independence it contains $p_{\rm GW}+p_{\rm EM}-p_{\rm GW}p_{\rm EM}$; under triggered or otherwise correlated observing, it must come from joint injections.

References

• Cabrera, T., Palmese, A., and Fishbach, M. (2026), “Multimessenger Constraints on LIGO/Virgo/KAGRA Gravitational-Wave Binary Black Holes Merging in Active Galactic Nucleus Disks,” The Astrophysical Journal 1000, 234, doi:10.3847/1538-4357/ae4d37.
• IceCube Collaboration et al. (2026), “IceCube Second Track Data Release IceTracks-DR2: Data from 2008–2022 for Neutrino Source Searches,” arXiv:2605.19040v1, arXiv; public data release, doi:10.7910/DVN/MMIIZA.
• LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (2026), “GWTC-5.0: Population Properties of Merging Compact Binaries,” arXiv:2605.27226, arXiv; GWTC-5.0 data release version 6, doi:10.7935/bk00-6a89.
• Mandel, I., Farr, W. M., and Gair, J. R. (2019), “Extracting Distribution Parameters from Multiple Uncertain Observations with Selection Biases,” Monthly Notices of the Royal Astronomical Society 486, 1086–1093, doi:10.1093/mnras/stz896.
• Planck Collaboration (2020), “Planck 2018 Results. VI. Cosmological Parameters,” Astronomy & Astrophysics 641, A6, corrigendum 652, C4, doi:10.1051/0004-6361/201833910.