Jun 15, 2026
A Cross-Correlation Pipeline for Multimessenger Tomography of Inhomogeneous Dark Matter via Gravitational Wave / Electromagnetic Arrival Time Differentials
- Carlos Cruz1,
- Carlos Cruz1
- 1Independent
- Carlos Cruz: Please refer to the attached PDF for precise mathematical formulation that may be compromised due to the interpreter used by the publishing host.
Protocol Citation: Carlos Cruz, Carlos Cruz 2026. A Cross-Correlation Pipeline for Multimessenger Tomography of Inhomogeneous Dark Matter via Gravitational Wave / Electromagnetic Arrival Time Differentials. protocols.io https://dx.doi.org/10.17504/protocols.io.bp2l6on45lqe/v1
License: This is an open access protocol distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Protocol status: In development
Please download the attached PDF to avoid confusion caused by online interpreters
Created: June 15, 2026
Last Modified: June 15, 2026
Protocol Integer ID: 319148
Keywords: multimessenger tomography of inhomogeneous dark matter, specific dark matter scenario, regarding specific dark matter scenario, gravitational wave catalog, dark matter model constraint, generation gravitational wave catalog, dark matter, inhomogeneous distribution of dark matter, cosmological lines of sight, derivation of dark matter model constraint, inhomogeneous dark matter, cosmological line, gravitational wave, em multimessenger events with measured arrival time differential, galaxy density field, multimessenger tomography, sky localization, measured arrival time differential, electromagnetic arrival time differential, electromagnetic arrival time differentials this protocol, em multimessenger event, coincident gw, arrival time differential, reproducible computational pipeline, lensing survey
Disclaimer
This protocol intentionally describes the measurement methodology only. The following are explicitly outside its scope and remain the independent intellectual contribution of the user:
- The physical model relating \(\hat{C}\) to dark matter column density and forward scattering amplitude
- Derivation of constraints on specific dark matter models
- Scientific interpretation and conclusions
- Any results obtained by executing the pipeline
Publishing this protocol establishes a timestamped, citable (DOI-bearing) record of the analysis methodology while leaving the scientific results unpublished and available for separate publication. License: CC BY 4.0, which permits reuse with attribution and preserves the authors' priority over the method.
Abstract
This protocol describes a reproducible computational pipeline for measuring the cross-correlation between gravitational wave / electromagnetic (GW/EM) arrival time differentials and large-scale structure, for the purpose of probing the inhomogeneous distribution of dark matter along cosmological lines of sight. The pipeline takes as input a catalog of coincident GW/EM multimessenger events with measured arrival time differentials and sky localizations, together with large-scale structure maps (weak lensing convergence or galaxy density fields), and produces as output a cross-correlation statistic and its associated significance. The protocol specifies data preparation, estimator construction, systematic controls, and significance evaluation. It is detector- and survey-agnostic, applicable to current and next-generation gravitational wave catalogs and Stage III/IV lensing surveys. This protocol describes the analysis methodology only; scientific interpretation of results, derivation of dark matter model constraints, and conclusions regarding specific dark matter scenarios are outside the scope of this protocol.
Attachments
Protocols_GW.pdf
183KB
Guidelines
This protocol assumes the user has access to:
- A catalog of coincident GW/EM events with arrival time differentials and sky/redshift localization
- Large-scale structure maps overlapping the GW event sky footprint
- A computing environment with standard scientific Python (NumPy, SciPy, Astropy, HEALPix/healpy)
This protocol describes the measurement procedure. It does not prescribe the interpretation of the resulting cross-correlation statistic, nor the mapping from that statistic to dark matter model parameters. Those steps constitute independent scientific analysis and are intentionally left to the user. Background rationale for the method—including the physical motivation for treating gravitational waves as vacuum reference signals and photons as medium-sensitive probes—is summarized in the Background section but not derived here.
Materials
INPUT 1 — Multimessenger event catalog. A table of N coincident GW/EM events, each with: event identifier; measured arrival time differential Δt_obs (EM arrival minus GW arrival), with uncertainty; sky localization (right ascension, declination), with uncertainty region; source redshift or luminosity distance, with uncertainty; EM counterpart type and observed frequency/energy band.
INPUT 2 — Large-scale structure map(s). One or more of: weak lensing convergence (κ) map covering the GW event footprint; galaxy number density field; reconstructed projected matter density map. Each with associated redshift bin / lensing kernel information.
INPUT 3 — Astrophysical source-delay model. A prescription for the expectation value and variance of the source timing uncertainty Δt_source for the relevant class of EM counterparts (e.g., short GRB jet breakout timescale distribution).
SOFTWARE. Scientific Python environment: NumPy, SciPy, Astropy, healpy (or equivalent HEALPix tools), and a statistical resampling library for bootstrap/jackknife.
Safety warnings
CRITICAL: Localization precision directly determines the achievable cross-correlation signal-to-noise. Events with localization regions large compared to the correlation length of the large-scale structure map contribute primarily noise and should be down-weighted or cut.
Procedure
Assemble the multimessenger event catalog (INPUT 1). Verify that each event has a measured arrival time differential, a sky localization, and a redshift/distance estimate. Exclude events lacking any of these.
For each event, define the sky localization as a probability map over the celestial sphere (e.g., HEALPix format). Record the redshift posterior for each event.
Apply quality cuts. Recommended cuts: minimum localization precision (e.g., area of 90% credible region below a threshold), minimum signal-to-noise in both GW and EM channels, and unambiguous association between GW and EM counterpart. Document all cuts and the resulting effective sample size \(N_{\text{eff}}\).
Compute the expected astrophysical source timing delay (\(\Delta t_{\text{source}}\)) from the source-delay model (INPUT 3), appropriate to the EM counterpart class and any available source properties.
Define the residual timing estimator for each event \(i\):
\[
\hat{T}_i = \Delta t_{\text{obs},i} - (\Delta t_{\text{source},i} - \Delta t_{\text{geo},i})
\]
where \(\Delta t_{\text{geo},i}\) is the residual differential lensing delay (Stage 2b). Under the null hypothesis (no dark matter–induced photon delay), \(\hat{T}_i\) contains only source-delay residual scatter, lensing-model residual, and instrumental noise.
Propagate the uncertainty in \(\Delta t_{\text{source}}\) into the per-event weight \(w_i\) used in subsequent steps. Events with poorly constrained source delays receive lower weight.
NOTE: The absolute calibration of \(\Delta t_{\text{source}}\) does not need to be perfect. Any constant offset or large-scale-structure-uncorrelated error in the source-delay model contributes to the null and does not bias the cross-correlation estimator, which is sensitive only to the structure-correlated component of \(\hat{T}\).
Compute the residual differential lensing delay arising from the spatial offset between the GW emission point (merger location) and the EM emission region (jet breakout front or photosphere). To first order in the angular separation:
\[
\Delta t_{\text{geo},i} \approx \frac{(1 + z_{d,i})}{c} \frac{D_d D_s}{D_{ds}} [\nabla \psi_i \cdot \Delta \theta_i(t)],
\]
where \(D_d, D_s, D_{ds}\) are the angular diameter distances to the lens, source, and between lens and source; \(\nabla \psi_i\) is the deflection angle from the foreground potential, reconstructed from the weak-lensing convergence and shear \((\kappa, \gamma)\) along the line of sight; and \(\Delta \theta_i(t)\) is the angular offset between the GW and EM emission regions.
Propagate the lens-model covariance and the uncertainty in \(\Delta \theta_i\) into the per-event weight.
CRITICAL: Unlike the source-delay term, \(\Delta t_{\text{geo}}\) is correlated with the large-scale structure field \(\kappa\), because it is sourced by the same foreground potential. It therefore does not average to zero in the cross-correlation and is not removed by the Stage 6 position/redshift null tests. It must be modeled and subtracted explicitly. The residual factorizes as (lensing response, from \(\kappa/\gamma\), correlated with the target) \(\times\) (\(\Delta \theta(t)\), from merger astrophysics, uncorrelated with the target); error in the latter factor enters \(\hat{C}\) as variance rather than coherent bias to first order, with residual bias appearing only at second order in the lensing deflection.
Project the large-scale structure map (INPUT 2) into the same sky pixelization scheme used for the event localizations.
For each event, extract the large-scale structure value (e.g., convergence \(\kappa\)) at the event sky position, integrated over the event’s localization probability map and weighted by the appropriate lensing/redshift kernel up to the source redshift:
\[
\kappa_i = \int p_i(\theta, \phi) \kappa(\theta, \phi; z_i) \, d\Omega,
\]
where \(p_i\) is the event localization probability map.
Subtract the mean of the large-scale structure field over the survey footprint so that the cross-correlation estimator is built on the fluctuation field \(\delta \kappa = \kappa - \langle \kappa \rangle\).
Construct the weighted cross-correlation estimator:
\[
\hat{C} = \frac{\sum_i w_i \hat{T}_i \delta \kappa_i}{\sum_i w_i},
\]
where \(w_i\) are the per-event weights from Step 2.3.
(Optional, tomographic mode) Bin events by source redshift into shells. Compute \(\hat{C}\) separately within each redshift shell to obtain the redshift dependence of the cross-correlation. Record \(\hat{C}(z)\) for each shell.
Record the point estimate of \(\hat{C}\) (and \(\hat{C}(z)\) if tomographic mode is used).
Significance Evaluation
Estimate the variance of \(\hat{C}\) under the null hypothesis using a resampling procedure. Recommended: randomly shuffle the assignment of \(\hat{T}_i\) values to event sky positions a large number of times (≥ 10^4), recomputing \(\hat{C}\) for each shuffle. This destroys any true correlation while preserving the marginal distributions of both \(\hat{T}\) and \(\delta k\).
Construct the null distribution of \(\hat{C}\) from the shuffled realizations. Compute the \(p\)-value of the observed \(\hat{C}\) against this null distribution.
As an independent cross-check, perform a jackknife over sky regions: recompute \(\hat{C}\) omitting each sky region in turn, to confirm the result is not driven by a single region or event.
Report the point estimate \(\hat{C}\), its resampling-based uncertainty, the \(p\)-value, and the jackknife stability check.
Systematic Controls
Null test — random sky positions. Repeat Stages 3–5 using random sky positions in place of the true event localizations. The cross-correlation should be consistent with zero. A nonzero result indicates a systematic in the pipeline.
Null test — scrambled redshifts. Repeat using scrambled event redshifts (breaking the true lensing-kernel association). The signal should degrade toward zero.
Frequency-split test. If the EM counterparts span a range of observed frequencies/energies, split the sample by frequency and recompute \(\hat{C}\). A genuine refractive signal carries a characteristic frequency dependence; a frequency-independent result flags a non-dispersive systematic.
Differential lensing cross-check. Verify that the residual differential lensing delay \(\Delta t_{\text{geo}}\) (Stage 2b) has been correctly subtracted. Because \(\Delta t_{\text{geo}}\) is correlated with \(\kappa\), it cannot be caught by the position/redshift null tests above; instead, exploit its frequency dependence. The dark-matter refractive delay scales as \(\omega_{\text{EM}}^{-2}\), whereas \(\Delta t_{\text{geo}}\) is achromatic in \(\omega_{\text{EM}}\) to leading order. Split the sample by EM frequency band and confirm that the surviving cross-correlation carries the \(\omega_{\text{EM}}^{-2}\) dependence of a genuine refractive signal and not the achromatic signature of residual lensing. Where multi-band timing of a single counterpart is available, the inter-band difference annihilates the achromatic lensing residual by construction and provides the strongest control.
Document all systematic test outcomes alongside the primary result.
Troubleshooting
Issue: Null tests (6.1/6.2) show nonzero correlation
Likely cause: Pixelization mismatch or mean-subtraction error
Action: Re-verify Stage 3 map projection and \(\delta \kappa\) construction
Issue: Estimator dominated by single event
Likely cause: Insufficient localization cuts
Action: Tighten Step 1.3 quality cuts; rerun jackknife
Issue: Large resampling variance
Likely cause: Insufficient effective sample size
Action: Report as upper limit; await larger event catalog
Issue: Frequency-split inconsistency
Likely cause: Non-dispersive systematic
Action: Investigate instrumental/source frequency dependence
Protocol references
[1] Abbott BP, et al. GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. Phys Rev Lett. 2017;119:161101.
[2] Lorimer DR, Kramer M. Handbook of Pulsar Astronomy. Cambridge University Press; 2004.
[3] Takahashi R, Nakamura T. Arrival Time Differences Between Gravitational Waves and Electromagnetic Signals Due to Gravitational Lensing. Astrophys J. 2017;835:103.
[4] Hu W. Weak Lensing Tomography. Astrophys J Lett. 1999;522:L21.