Frequentist Regularised GW Mapping Analysis

• Frequentist regularised gravitational-wave mapping is a method that reconstructs the sky's angular power by decomposing the GW signal and stabilising inversion with regularisation techniques.
• It employs a maximum-likelihood framework along with Tikhonov regularisation and truncated SVD to handle ill-conditioned inversions and improve point-source recovery.
• The approach unifies data from various detector networks, enabling researchers to quantify anisotropy and optimize angular resolution under instrumental and geometric constraints.

Frequentist regularised gravitational-wave mapping analysis constitutes a foundational methodology for reconstructing the angular power distribution of gravitational-wave (GW) signals on the sky using data from detector networks such as pulsar timing arrays (PTAs) or ground-based interferometers. Based on a maximum-likelihood formalism, it utilizes explicit sky decompositions, detector response modeling, and sophisticated regularisation techniques to stabilize inversion in the presence of ill-posedness. This framework is crucial for producing GW sky maps, quantifying anisotropy, and resolving point sources, as well as for understanding the effective angular resolution achievable given instrumental and geometric limitations (Cornish et al., 2014, Grunthal et al., 20 Jan 2026, 0708.2728).

1. Mathematical Foundations: Sky Decomposition and Linear Data Model

The analysis begins with a decomposition of the GW field $h_A(\hat\Omega)$ (for polarizations $A\in\{+,\times\}$) on the sphere into an orthonormal basis $\{\phi_i(\hat\Omega)\}$, satisfying

$$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$

Both pixel-space and spherical-harmonic expansions can be employed, with $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$. All coefficients are collected into a $2N$-vector $$s = (s_{+,1},...,s_{+,N}; s_{\times,1},...,s_{\times,N})^T$$ (Cornish et al., 2014).

For PTAs, each pulsar $p$ yields timing-residual data $\delta t_p$ modeled as

$\delta t_p = \sum_{A,i} F_{p,A,i} s_{A,i} + n_p,$

where $A\in\{+,\times\}$0 is the detector response matrix projected onto the chosen basis, incorporating geometric quantities (antenna patterns, pulsar terms, light travel delays). The vector form is $A\in\{+,\times\}$1, with noise $A\in\{+,\times\}$2 having Gaussian covariance $A\in\{+,\times\}$3 (Cornish et al., 2014, 0708.2728).

Alternatively, one can construct cross-correlated "dirty maps" from pairwise detector data (e.g., for LIGO/Virgo or PTA pairs), integrating over frequency and time with optimal filters, producing a linear model $A\in\{+,\times\}$4 noise, where $A\in\{+,\times\}$5 is the response matrix linking spherical-harmonic powers $A\in\{+,\times\}$6 to observations (Grunthal et al., 20 Jan 2026, 0708.2728).

2. Likelihood Formulation and Maximum-Likelihood Estimation

The frequentist analysis is rooted in a composite Gaussian likelihood: $A\in\{+,\times\}$7 which, analogously, holds for the cross-correlation vector in radiometer or PTA map-making formalisms. Maximization yields normal equations: $A\in\{+,\times\}$8

$A\in\{+,\times\}$9

with Fisher or normal matrices $\{\phi_i(\hat\Omega)\}$0 ($\{\phi_i(\hat\Omega)\}$1 or $\{\phi_i(\hat\Omega)\}$2) and "dirty map" vectors $\{\phi_i(\hat\Omega)\}$3 or $\{\phi_i(\hat\Omega)\}$4 (Cornish et al., 2014, Grunthal et al., 20 Jan 2026).

The formal unregularised maximum-likelihood (ML) map solution is: $\{\phi_i(\hat\Omega)\}$5 However, in practice, $\{\phi_i(\hat\Omega)\}$6 is typically ill-conditioned due to limited detector number, non-uniform sky sensitivity, and incomplete coverage, necessitating further regularisation (Cornish et al., 2014, Grunthal et al., 20 Jan 2026, 0708.2728).

3. Regularisation: Tikhonov and Truncated SVD

Direct inversion of $\{\phi_i(\hat\Omega)\}$7 amplifies noise in weakly constrained sky modes. Two primary regularisation strategies are employed:

• Tikhonov (Ridge) Regularisation: Addition of a penalty term $\{\phi_i(\hat\Omega)\}$8 (typically $\{\phi_i(\hat\Omega)\}$9) yields $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$0, with the regularised solution $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$1. The parameter $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$2 is chosen via L-curve, cross-validation on noise, or empirical Bayes (Cornish et al., 2014, 0708.2728).
• Truncated Singular Value Decomposition (SVD): $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$3 is decomposed as $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$4 with eigenvalues $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$5. Only the $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$6 largest singular values are retained, forming a pseudo-inverse $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$7 with $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$8 for $$\int_{S^2} d\Omega\, \phi_i(\hat\Omega)\, \phi_j(\hat\Omega) = \delta_{ij}.$$9, else $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$0 (Grunthal et al., 20 Jan 2026). This procedure discards noise-dominated sky modes, providing a stable "clean map" $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$1 and corresponding uncertainties (Cornish et al., 2014, Grunthal et al., 20 Jan 2026, 0708.2728).

Table: Key Regularisation Approaches

Regularisation Type	Modification	Typical Use Case
Tikhonov	$$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$2	General ill-posed linear inversion
Truncated SVD	Keep $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$3 largest $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$4	Sky mapping with limited baselines

These regularisation approaches are robust, tractable, and adaptable to detector geometry and noise properties.

4. Sky-Map Statistics: Isotropy, Anisotropy, and Point-Source Recovery

Isotropic Limit and Hellings–Downs Cross-Correlation

For a statistically isotropic GW background, the covariance structure simplifies: $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$5 leading to the classic Hellings–Downs overlap reduction function (for PTAs, $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$6), and just the monopole sky mode survives regularisation. Marginalisation over sky coefficients analytically recovers the standard cross-correlation statistic (Cornish et al., 2014).

Anisotropic Backgrounds and Point Sources

No assumption of isotropy is necessary. Anisotropy is modeled by retaining higher spherical harmonic modes or pixel basis elements. For point-source recovery, localized delta-function templates $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$7 are included in the basis, with their amplitudes estimated via the same regularised normal equation machinery (Cornish et al., 2014). Broader prior covariances can be assigned for more flexibility in capturing extended anisotropy (Cornish et al., 2014, Grunthal et al., 20 Jan 2026).

5. Quantifying and Optimizing Angular Resolution

A defining innovation is the characterization and optimization of angular resolution using the point spread function (PSF) and distortion matrix formalism (Grunthal et al., 20 Jan 2026). The PSF directly quantifies the sky-patch area to which a true point source's recovered flux is spread due to the measurement process and regularisation. The distortion matrix is constructed as: $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$8 where $$h_A(\hat\Omega) = \sum_{i=1}^N s_{A,i}\, \phi_i(\hat\Omega)$$9 relates the pixel and spherical-harmonic bases. The effective PSF area at each sky location $2N$0 is given by summing pixels above half-max: $2N$1 yielding a PSF map $2N$2 across the sky. In regions with dense detector coverage, $2N$3 approaches the geometric limit set by the nearest-neighbor spacing, expressed as

$2N$4

Simulations show that local resolution closely tracks the underlying detector geometry, and optimal $2N$5 is chosen such that the scale $2N$6, with $2N$7 (Grunthal et al., 20 Jan 2026).

A significant improvement over earlier methods is the implementation of an adaptive, variable local-resolution scheme, allowing finer resolution where the detector network is more densely populated, as confirmed by MeerKAT PTA simulations (Grunthal et al., 20 Jan 2026).

6. Numerical Implementation and Computational Aspects

The map-making pipeline follows a sequence:

1. Choosing a sky basis (pixels, spherical harmonics).
2. Computing frequency-dependent detector response matrices.
3. Measuring data cross-spectra, assembling noise or covariance matrices.
4. Whitening data and response matrices ($2N$8 such that $2N$9).
5. SVD or eigendecomposition of the compactified response.
6. Regularisation via Tikhonov or truncation.
7. Solving for sky coefficients in the low-dimensional subspace.
8. Inverting transforms to reconstruct sky maps ($$s = (s_{+,1},...,s_{+,N}; s_{\times,1},...,s_{\times,N})^T$$0).
9. Quantifying angular resolution with PSF maps. 10. Examining higher-order or localized coefficients for anisotropy or individual source detection.

Table: Computational Steps in Frequentist Regularised GW Mapping

Step	Major Operation	Reference
1	Basis choice, response modeling	(Cornish et al., 2014, Grunthal et al., 20 Jan 2026)
2–3	Data/covariance assembly	(Cornish et al., 2014, Grunthal et al., 20 Jan 2026)
4	Whitening, low-rank projection	(Cornish et al., 2014)
5–6	SVD/truncation, regularisation	(Grunthal et al., 20 Jan 2026)
7–8	Sky map solution and inverse transform	(Cornish et al., 2014)
9	Local PSF/area calculation	(Grunthal et al., 20 Jan 2026)

The dominant computational costs arise in beam and response calculation, typically $$s = (s_{+,1},...,s_{+,N}; s_{\times,1},...,s_{\times,N})^T$$1, which are handled by exploiting matrix symmetry, parallelisation, and fast transform techniques (e.g., FFT interpolation) (0708.2728). The rank of the relevant matrices never exceeds $$s = (s_{+,1},...,s_{+,N}; s_{\times,1},...,s_{\times,N})^T$$2, allowing efficient low-rank representation and inversion (Cornish et al., 2014).

7. Extensions, Applications, and Performance in Practice

The frequentist regularised framework is applicable across detector types, including GW radiometers for ground-based interferometers (0708.2728) and PTAs (Cornish et al., 2014, Grunthal et al., 20 Jan 2026). It extends naturally to:

• Multiple-baseline, multi-detector networks by stacking response and noise matrices (0708.2728).
• Polarised backgrounds by including amplitude coefficients for each polarisation (0708.2728).
• Empirical-Bayes embedding for hyperparameter selection (e.g., regularisation strength) (Cornish et al., 2014).
• Adaptive-resolution mapping to match local detector density, yielding significant gains (e.g., the clean-map S/N for a simulated continuous GW increased by ~factor 2 and hotspot area shrank by ~40% in MeerKAT simulations) (Grunthal et al., 20 Jan 2026).

In the isotropic limit, the method reduces to the standard Hellings–Downs cross-correlation statistic. For anisotropic or point-like signals, the framework provides a rigorous basis for statistical significance quantification and limits spurious anisotropy artifacts. Interpretation of map coefficients, PSF structure, and detection statistics are geometrically transparent within this formalism.

This approach provides a unified and robust methodology for extracting GW sky maps, probing anisotropy, and quantifying the true resolving power of current and near-future GW detector networks (Cornish et al., 2014, Grunthal et al., 20 Jan 2026, 0708.2728).