Anisotropic Stochastic Mapping

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

An anisotropic stochastic search uses Earth’s rotation and a detector network as a synthetic aperture. Rather than estimate one isotropic amplitude, it reconstructs the sky components to which the baselines are sensitive. The result is an inverse problem with a nontrivial null space: a displayed map is always a response-filtered and regularized estimate, not a direct photograph of gravitational-wave intensity.

Helpful background. Cross-Correlation for Stochastic Backgrounds derives the baseline cross power and isotropic overlap function. Here the overlap becomes a time-, frequency-, and direction-dependent response.

From directional response to a linear map

Assume a separable tensor-background model

\mathcal P(f,\boldsymbol\Omega) = H(f)\mathcal P(\boldsymbol\Omega), \qquad \mathcal P(\boldsymbol\Omega) = \sum_\alpha \mathcal P_\alpha e_\alpha(\boldsymbol\Omega),

where $e_\alpha$ are real sky pixels or real spherical harmonics, and the coefficients $\mathcal P_\alpha$ are real. For a short segment centered at sidereal time $t$ , define the normalized complex cross periodogram $C_{btf}=(2/T)\widetilde s_I^*(t,f)\widetilde s_J(t,f)$ for an untapered segment of length $T$ . Its expectation is

\left\langle C_{btf}\right\rangle = H(f)\sum_\alpha \gamma_{btf,\alpha}\mathcal P_\alpha.

The generalized overlap coefficients are

\gamma_{btf,\alpha} = \sum_{A=+,\times} \int d\boldsymbol\Omega\, e_\alpha(\boldsymbol\Omega) F_I^A(t,\boldsymbol\Omega) F_J^A(t,\boldsymbol\Omega) e^{2\pi i f\boldsymbol\Omega\cdot\Delta\boldsymbol x_b/c}.

The normalization of $\mathcal P$ is fixed by this response equation; it must be converted before comparing with a normalized isotropic overlap function. Use $\Delta\boldsymbol x_b=\boldsymbol x_J-\boldsymbol x_I$ and the positive-exponent forward transform of the preceding page. Rotation changes the antenna factors relative to the celestial sky. Separation supplies phase fringes whose angular scale is roughly

\Delta\theta\sim\frac{c}{f\lvert\Delta\boldsymbol x_b\rvert}.

Higher frequency and longer baselines can improve angular resolution, but only where the detector noise and antenna response provide usable weight.

Dirty maps and the Fisher matrix

Let $V_{btf}=\mathbb E[|C_{btf}-\langle C_{btf}\rangle|^2]$ . For one independent Fourier cell with weak signal and independent detector noise, $V_{btf}\simeq P_{I,tf}P_{J,tf}$ ; averaging $m$ independent cells divides this variance by $m$ . In the Gaussian approximation to the averaged cross spectra, the dirty map and Fisher matrix for real sky coefficients are

X_\alpha =2\operatorname{Re}\sum_{btf} \frac{H(f)\gamma_{btf,\alpha}^*C_{btf}} {V_{btf}},

and Fisher matrix

\Gamma_{\alpha\beta} =2\operatorname{Re}\sum_{btf} \frac{H^2(f)\gamma_{btf,\alpha}^* \gamma_{btf,\beta}} {V_{btf}}.

Their expectations satisfy

\langle\boldsymbol X\rangle = \boldsymbol\Gamma\boldsymbol{\mathcal P}, \qquad \operatorname{Cov}(\boldsymbol X) \simeq \boldsymbol\Gamma.

The factor two uses both real quadratures of each positive-frequency complex datum, whose individual noise variances are $V/2$ . Tapers, overlapping segments, correlated baselines, and strong signals need their actual covariance rather than this diagonal weight. These conventions make $\langle X\rangle=\Gamma\mathcal P$ and $\operatorname{Cov}(X)=\Gamma$ consistent, not merely proportional Thrane et al. 2009, §III.

The dirty map is the true sky convolved with the network point-spread function. If $\Gamma$ were well conditioned, the maximum-likelihood clean map and covariance would be

\widehat{\boldsymbol{\mathcal P}} = \boldsymbol\Gamma^{-1}\boldsymbol X, \qquad \boldsymbol\Sigma = \boldsymbol\Gamma^{-1}.

In practice, some sky combinations produce almost no cross power. Small eigenvalues amplify noise catastrophically, so inversion uses a truncated singular-value decomposition, a prior, or another regularizer. If $W$ is the chosen regularized inverse, then $\widehat{\mathcal P}=WX$ , its mean is $W\Gamma\mathcal P$ , and its noise covariance is $W\Gamma W^T$ . A prior-regularized inverse is therefore not automatically the sampling covariance. Every reported map must state which eigenmodes were retained. Setting discarded modes to zero is an estimator choice, not evidence that the true sky has no power there.

Pixel maps, harmonics, and targeted directions

A pixel basis is intuitive but can contain many nearly redundant degrees of freedom. A spherical-harmonic basis makes angular scale explicit through $\ell_{\max}$ , and a radiometer targets pointlike directions with a chosen point-spread function. These analyses test different source morphologies.

Changing $\ell_{\max}$ after viewing a hot spot incurs trials. So does scanning many spectral indices or source directions. Detection statistics must be calibrated over the complete search. Positivity of intensity can be imposed, but it changes the estimator distribution and must be included in upper-limit coverage.

Correlated magnetic or instrumental noise can also project into celestial modes. Sidereal modulation helps, but does not guarantee rejection. Environmental witnesses, alternative baselines, time slides, and maps in detector-fixed coordinates are important null tests.

Computational convergence and reproduction

Record the basis normalization, sky frame, pixelization or $\ell_{\max}$ , sidereal segmentation, window overlap, frequency binning, spectral shape, notches, calibration, baseline set, and singular-value threshold. Save both $\boldsymbol X$ and $\boldsymbol\Gamma$ when feasible; a cleaned map alone cannot reproduce the covariance or the regularization.

A controlled benchmark injects a point source and a low-order spherical-harmonic sky into correlated Gaussian strain. Then:

refine the sky grid and frequency bins until recovered integrated flux stabilizes;
increase $\ell_{\max}$ and identify modes whose uncertainties diverge;
vary the singular-value cutoff and report bias–variance changes;
verify the point-source response against the predicted beam;
perform zero-signal, time-slide, scrambled-sky, and environmental-coupling null tests;
confirm injection coverage using the same regularization and trials as the search.

An apparent fine angular feature that disappears when the basis is refined or the cutoff is moved is not resolved structure.

Observational status through 2026-08-18

The official O4a directional search reports no significant persistent gravitational-wave signal and sets direction-, spectrum-, and morphology-dependent upper limits LVK 2025. Those upper limits are not detections of map pixels. They are statements about the analyzed baselines, frequency bands, sky basis, regularization, calibration, and pipeline recovery.

Exact numerical limits and sky products should be taken from the cited release. Additional observing data, a changed calibration, a new official directional product, or an altered correlated-noise model triggers a refresh.

Exercises

Consider

\boldsymbol\Gamma = \begin{pmatrix} 1&1\\ 1&1 \end{pmatrix}.

Find its measured and null eigenmodes.

Solution

The normalized eigenvectors are

\boldsymbol v_+ = \frac{1}{\sqrt2} \begin{pmatrix}1\\1\end{pmatrix}, \qquad \boldsymbol v_- = \frac{1}{\sqrt2} \begin{pmatrix}1\\-1\end{pmatrix},

with eigenvalues $2$ and $0$ . The data measure only the sum of the two pixels. Their difference is in the exact null space and cannot be reconstructed without prior information.

A baseline has length $D=3000\ {\rm km}$ . Estimate the fringe scale at $f=100\ {\rm Hz}$ .