PyStoch: Efficient Pipeline for SGWB Mapmaking

• PyStoch is an open-source Python framework for gravitational-wave searches, using sidereal-day folding and HEALPix sky maps to enable fast, directional cross-correlation analyses.
• It leverages matrix multiplications and data compression with seed matrices to reduce computational time from years to minutes while lowering storage demands.
• The unified pipeline supports both pixel-based and spherical-harmonic estimators, allowing multicomponent analyses in anisotropic SGWB and directed continuous-wave searches.

Searching arXiv for PyStoch-related papers to ground the article in the cited literature. PyStoch is an open-source, Python-based pipeline for gravitational-wave radiometry and anisotropic stochastic gravitational-wave background (SGWB) mapmaking that was introduced to exploit sidereal-day data folding, matrix-based computation, and HEALPix-based sky representations in order to make directional cross-correlation analyses computationally practical at scale (Ain et al., 2018). In subsequent work, its role broadened from fast pixel-basis radiometer mapping to a unified framework that also supports spherical-harmonic (SpH) analyses, all-sky all-frequency directional searches, joint multicomponent anisotropy inference, and directed continuous-wave (CW) searches from known sky locations (Suresh et al., 2020).

1. Historical development and problem setting

PyStoch was developed in the context of SGWB searches, where the observable is a weak excess correlation between spatially separated detectors with approximately uncorrelated instrumental noise. Conventional anisotropic SGWB analyses had already used cross-correlation techniques, but they were computationally expensive and, in practice, implemented with distinct numerical pipelines for different sky bases. The 2018 presentation of PyStoch emphasized three design decisions: folding data into one sidereal day, replacing explicit loops with matrix multiplications, and incorporating HEALPix pixelization tools; together, these made the radiometer analysis “a few thousand times” faster and allowed all-sky maps to be produced “in just a few minutes on an ordinary laptop” (Ain et al., 2018).

The pipeline’s scope then expanded along three directions. First, it became part of the anisotropic SGWB workflow used with Advanced LIGO and Advanced Virgo folded data for broadband radiometer, narrow-band radiometer, and spherical-harmonic decomposition searches over the first three observing runs (Collaboration et al., 2021). Second, the algebra of pixel-basis and SpH mapmaking was unified within PyStoch, eliminating the need for manifestly redundant pipelines and making per-frequency SpH searches computationally feasible (Suresh et al., 2020). Third, the same folded cross-correlation machinery was repurposed for directed CW searches, where PyStoch serves as a low-cost first-pass filter for nearly monochromatic signals from known sky positions, with more sensitive CW pipelines reserved for follow-up (Salvadore et al., 20 Jul 2025).

A common misconception is that PyStoch is only a broadband stochastic-background mapmaker. The published literature instead describes it as a general directional cross-correlation engine that produces maps at every frequency bin, supports both pixel and SpH representations, and has been used in both stochastic-background and directed CW searches (Collaboration et al., 2021).

2. Statistical estimator and sky representations

PyStoch implements the standard radiometer maximum-likelihood mapmaker. Let $s_1(t;f)$ and $s_2(t;f)$ denote Fourier components of the strain in two detectors for baseline $I$. In the pixel basis, the “dirty” map $X_p$ and Fisher matrix $\Gamma_{pp'}$ are

$$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$

$$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$

Here $\tau$ is the total observation time, $H(f)$ is the assumed SGWB frequency spectrum, $P_{I_1,2}(t;f)$ are the one-sided noise PSDs, and $s_2(t;f)$0 is the overlap reduction function (ORF) for baseline $s_2(t;f)$1 and sky pixel $s_2(t;f)$2 (Suresh et al., 2020). The maximum-likelihood estimate in pixel space is

$s_2(t;f)$3

In the SpH basis, the anisotropy is expanded as

$s_2(t;f)$4

with dirty coefficients and Fisher matrix

$s_2(t;f)$5

The corresponding maximum-likelihood estimator is

$s_2(t;f)$6

This formulation supports both localized and diffuse searches: pixel maps are used to probe pointlike hot spots, while SpH moments constrain diffuse anisotropies (Suresh et al., 2020).

The same formalism also underlies isotropic searches. In the unified PyStoch treatment, the isotropic limit is recovered by setting $s_2(t;f)$7, which yields the standard isotropic SGWB estimator (Suresh et al., 2020).

3. Folded-data architecture and computational organization

The key algorithmic acceleration in PyStoch is sidereal folding. Because antenna patterns, time-delay phase factors, and related baseline-dependent quantities repeat every sidereal day, months of data can be co-added into a single “folded day.” In the all-sky, all-frequency (ASAF) implementation, ordinary analyses step through segments with $s_2(t;f)$8, whereas folding reduces the time axis to $s_2(t;f)$9 sidereal-time slices; subsequent ORF evaluation, noise weighting, and mapmaking are then performed only on the folded data (Collaboration et al., 2021).

PyStoch expresses the mapmaking operations as linear algebra. In the original folded-data implementation, the baseline-dependent sky response is compressed into two “seed” matrices: an antenna-pattern product matrix and a time-delay seed matrix. For each frequency, the full ORF is generated on the fly by multiplying the antenna-pattern seed by the phase factor $I$0, so the full time-frequency-pixel ORF table never needs to be stored (Ain et al., 2018). This compression is central to the claim that PyStoch can stream over frequency while keeping memory requirements modest.

HEALPix is integrated throughout the pipeline. PyStoch stores sky-dependent quantities as HEALPix maps, uses standard angular-to-pixel conversions, and exploits HEALPix spherical-harmonic transforms. In the unified anisotropy pipeline, the antenna-pattern product $I$1 and the time delay $I$2 are stored as frequency-independent “seed” maps for each folded time segment, while conversion to the SpH basis uses HEALPix’s fast map2alm routines; because the ORFs are complex, a thin wrapper around map2alm was added (Suresh et al., 2020). In the ASAF search, PyStoch further factors the ORF into a “time-dependent modulation” tensor and a “frequency-dependent phase” tensor, avoiding a naïve nested loop over $I$3 (Collaboration et al., 2021).

The sky resolution depends on the analysis. The ASAF search used HEALPix $I$4, corresponding to $I$5 equal-area pixels, while the O1–O3 broadband radiometer analysis used $I$6 (Collaboration et al., 2021, Collaboration et al., 2021).

4. Unified mapmaking, basis transformations, and multicomponent inference

A central technical result associated with PyStoch is that pixel-basis and SpH estimators are algebraically transformable. The map coefficients satisfy

$I$7

and the Fisher matrices transform analogously,

$I$8

This relation implies that a pixel-basis analysis can be projected into the SpH basis after the fact, or the SpH Fisher matrix can be built directly from $I$9 (Suresh et al., 2020).

Earth’s rotation introduces an especially useful simplification: $X_p$0 As a result, only the $X_p$1 transform needs to be precomputed, with the remaining time dependence obtained by a phase factor. In the reported implementation, this reduces the cost enough that per-frequency ORFs in the SpH basis are generated in $X_p$2, and pixel and SpH maps with full Fisher matrices for one observing run can be built in $X_p$3 on a modern multi-core CPU (Suresh et al., 2020).

The same output objects support multicomponent anisotropy inference. For multiple assumed spectral shapes $X_p$4, PyStoch constructs dirty maps $X_p$5 and Fisher blocks $X_p$6. The formal joint clean-map solution is

$X_p$7

with covariance $X_p$8, where $X_p$9 (Suresh et al., 2021). At present detector sensitivities, the off-diagonal pixel-pixel terms are reported to be very small compared to the diagonal, and the standard analyses therefore drop them and invert only the much smaller $\Gamma_{pp'}$0 matrices at each pixel. The published argument is that this changes neither upper limits nor point-spread recovery at any level significant for O1–O3, while making joint-index anisotropic mapmaking tractable (Suresh et al., 2021).

This unified architecture also underlies later downstream products. The literature explicitly notes clean maps via SVD or norm regularization, decorrelated $\Gamma_{pp'}$1 spectra, polarization studies, and lensing as analyses that can be built on the common PyStoch output (Suresh et al., 2020).

5. Role in anisotropic SGWB searches

PyStoch has been used in the major directional SGWB searches with Advanced LIGO and Advanced Virgo data. In the O1–O3 anisotropic search, the pipeline was run on data folded over one sidereal day and used for broadband radiometer, narrow-band radiometer, and spherical-harmonic decomposition analyses. No evidence for gravitational-wave signals was found. The reported $\Gamma_{pp'}$2 confidence-level upper limit sky maps on the gravitational-wave energy flux from broadband point sources ranged from $\Gamma_{pp'}$3, while upper limits on the normalized gravitational-wave energy density spectrum from extended sources ranged from $\Gamma_{pp'}$4, depending on direction and spectral index; these limits improved on previous limits by factors of $\Gamma_{pp'}$5 (Collaboration et al., 2021).

In that analysis, the radiometer output was translated into physically interpretable quantities through

$\Gamma_{pp'}$6

with the first expression used for point-source broadband radiometer maps and the second for extended-source SHD maps (Collaboration et al., 2021).

PyStoch was also the enabling pipeline for the first all-sky all-frequency anisotropic SGWB search. That ASAF search analyzed every $\Gamma_{pp'}$7 frequency bin from $\Gamma_{pp'}$8 to $\Gamma_{pp'}$9, except bins likely to contain instrumental artefacts and hence notched, over $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$0 HEALPix sky pixels. No statistically significant evidence for narrowband gravitational-wave signals was found, and $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$1 confidence upper limits on the gravitational-wave strain were set for each pixel-frequency pair, in the range $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$2 (Collaboration et al., 2021). The same work stressed that a broadband analysis can miss narrowband signals because the signal-to-noise ratio of a narrowband signal can be significantly reduced when combined with detector output from other frequencies; this is the main methodological motivation for per-frequency mapmaking (Collaboration et al., 2021).

Candidate identification in ASAF used a noise-background estimation based on unphysical time shifts and a per-frequency maximum-SNR scan over contiguous $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$3 sub-bands. Hardware injections in O2 were recovered as outliers by this procedure, providing an internal validation of the follow-up logic (Collaboration et al., 2021).

6. Directed continuous-wave searches

PyStoch has also been repurposed for directed CW searches from known sky locations, where it functions as a radiometer-style coherent cross-correlation statistic over folded data. In the O3 study targeting four supernova remnants—Vela Jr., G347.3-0.5, Cassiopeia A, and the neutron star associated with the 1987A supernova remnant—the pipeline scanned the $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$4–$$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$5 band with default settings $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$6 and $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$7. No interesting candidates were identified in the real O3 analysis, and $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$8 confidence-level upper limits on the CW strain amplitude $$X_p = \frac{4}{\tau} \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),$$9 were set. The most stringent limit was $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$0 at $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$1 for Cassiopeia A; the best limits for the other targets were $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$2 at $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$3 for G347.3-0.5, $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$4 at $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$5 for Vela Jr., and $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$6 at $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$7 for the 1987A target (Salvadore et al., 20 Jul 2025).

The O3 CW implementation used a “bincomb” strategy to cope with spin-down broadening: adjacent $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$8 bins were combined in a sliding window of width $$\Gamma_{pp'} = 4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.$$9, with candidates flagged when the empirical $\tau$0-value was $\tau$1, roughly $\tau$2. Follow-up pipelines were 5-vector Resampling and Band-Sampled Data Frequency-Hough (Salvadore et al., 20 Jul 2025).

In the O4a supernova-remnant search, PyStoch was one of five directed pipelines and was described as a cross-correlation-based method that uses $\tau$3 and $\tau$4 folded Stochastic Intermediate Data. The O4a workflow included target-specific sky demodulation, adaptive bin grouping over window sizes from 1 bin up to an empirical maximum of 48 bins ($\tau$5), and empirical $\tau$6-value thresholding with $\tau$7, equivalently $\tau$8. In that study, none of the PyStoch candidates survived coherence follow-up. The tightest O4 Pan-PyStoch constraint was $\tau$9 near $H(f)$0 for Vela Jr., corresponding to a sensitivity depth $H(f)$1 (Collaboration et al., 26 Mar 2026).

These CW applications clarify PyStoch’s niche. The O3 supernova-remnant study reported that, relative to standard CW directed searches, PyStoch is approximately $H(f)$2–$H(f)$3 times less sensitive but approximately $H(f)$4–$H(f)$5 times faster, making it suitable as a first-pass filter for targets with poorly constrained spin parameters (Salvadore et al., 20 Jul 2025).

7. Validation, scaling behavior, and limitations

Validation results reported for PyStoch are unusually explicit. In the unified pixel/SpH mapmaking study, SpH moments obtained by direct cross-spectral-density projection and by pixel-to-SpH transformation matched to $H(f)$6 rms for $H(f)$7, while Fisher matrices agreed to $H(f)$8 between the conventional and PyStoch implementations. On that basis, the study concluded that a single skymap is sufficient to describe anisotropies in a stochastic background (Suresh et al., 2020).

Performance claims are analysis-dependent but consistently large. In the original folded-data mapmaking paper, a representative O1-like benchmark was reduced from $H(f)$9 CPU-years for the conventional unfolded pipeline to $P_{I_1,2}(t;f)$0 CPU-days with folding alone and to $P_{I_1,2}(t;f)$1 CPU-minutes with folding plus PyStoch; intermediate storage also shrank from $P_{I_1,2}(t;f)$2 to a $P_{I_1,2}(t;f)$3 folded dataset (Ain et al., 2018). In the ASAF search, PyStoch plus folding was reported to run $P_{I_1,2}(t;f)$4–$P_{I_1,2}(t;f)$5 faster than the original Matlab implementation without folding; each baseline-plus-run required $P_{I_1,2}(t;f)$6 on a single modern 2-core laptop, and the full O1+O2+O3 HLV combination took $P_{I_1,2}(t;f)$7–$P_{I_1,2}(t;f)$8 wall time on one core (Collaboration et al., 2021). In the O3 CW supernova-remnant study, the PyStoch stage ran in $P_{I_1,2}(t;f)$9 on a 4-core CPU (Salvadore et al., 20 Jul 2025).

The literature also records concrete limitations. In the 2025 CW implementation, PyStoch had no built-in demodulation for Doppler or spin-down, relied on the bincomb approximation, used fixed $s_2(t;f)$00, and folded data precluded segment-by-segment vetoes of non-Gaussian artifacts (Salvadore et al., 20 Jul 2025). In anisotropic multicomponent inference, the approximation of neglecting off-diagonal pixel-pixel covariance is justified only because present detector sensitivities render those terms very small relative to the diagonal (Suresh et al., 2021). These are not changes to the basic cross-correlation estimator, but they define the practical envelope within which PyStoch has been deployed.

Taken together, the published record places PyStoch at the intersection of radiometry, folded-data compression, and fast spherical analysis on the two-sphere. Its defining feature is not a new detection statistic, but an implementation strategy that makes frequency-resolved directional inference routine enough to support broadband SGWB searches, all-sky all-frequency scans, multicomponent anisotropy studies, and computationally economical CW candidate generation within a single family of workflows (Ain et al., 2018).