Shear 4PCF: Fourth-Order Weak Lensing Insights

• Shear 4PCFs are fourth-order statistics that quantify non-Gaussian features in the weak-lensing shear field beyond conventional two-point measures.
• They decompose the shear into natural components, separating connected from disconnected parts and ensuring robust E-mode extraction.
• Efficient multipole estimators and tree-based algorithms make practical 4PCF computation feasible for surveys like DES Y3 and improve covariance modeling.

Shear four-point correlation functions (4PCFs) are fourth-order correlation functions of the weak-lensing shear field and encode non-Gaussian information of the projected matter distribution that is not captured by second-order statistics. In cosmic shear, the field is a complex spin-2 quantity, so the 4PCF is not a single scalar object but a structured set of components whose rotational properties, parity behavior, and connected versus disconnected decomposition are central to both theory and measurement. Recent work has established a theoretical framework linking the shear 4PCF to fourth-order aperture statistics and to the convergence trispectrum, has developed efficient estimators with quadratic scaling, and has carried out an application to DES Y3 data, including a measurement of the connected part of the fourth-order aperture mass (Porth et al., 9 Sep 2025).

1. Formal definition and component structure

The cosmic shear field, being a complex spin-2 field, admits a rich set of higher-order correlation functions. The 4PCF captures information about the connected four-point statistics of the field that is not accessible with lower-order moments. Given four positions $$\mathbf{X}_0,\mathbf{X}_1,\mathbf{X}_2,\mathbf{X}_3$$ on the sky and projection angles $\zeta_\mu$, a natural component can be written, for instance, as

$$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$

where the arguments $\vartheta_i$ are the distances from $\mathbf{X}_0$ to the other points and $\phi_{ij}$ are the enclosed angles between legs of the quadrilateral (Porth et al., 9 Sep 2025).

Owing to the polar nature of shear, there are 16 real components of the 4PCF, but these are reorganized into 8 complex-valued natural components, denoted $\Gamma_\mu$, each transforming in a simple way under rotation. This reorganization is the analogue, at fourth order, of the use of spin-adapted quantities in lower-order shear statistics, and it is the natural starting point for any estimator or compression scheme. The same formalism also makes explicit that the observable depends on the full quadrilateral configuration, not only on separations but also on angular variables.

A useful distinction runs through the entire subject: the total four-point statistic contains both disconnected terms, which are products of lower-order correlations, and connected terms, which are irreducible fourth-order information. For a Gaussian random field, all connected four-point correlations vanish, leaving only the disconnected part. In the aperture-mass language this becomes

$$\langle M_\mathrm{ap}^4\rangle = 3\left[\langle M_\mathrm{ap}^2\rangle\right]^2,$$

so any deviation from this relation diagnoses non-Gaussian structure (Silvestre-Rosello et al., 9 Sep 2025).

2. Relation to aperture statistics and the trispectrum

A central development in the shear-4PCF literature is the explicit relation between the 4PCF, fourth-order aperture statistics, and convergence polyspectra. The $n$-th order aperture statistic $\langle M_\mathrm{ap}^n\rangle$ is defined as an ensemble average of products of aperture masses at radii $\zeta_\mu$0, and for $\zeta_\mu$1 it can be written as a five-dimensional integral over the shear 4PCF and analytically derived filter functions. For equal-scale apertures, one representative form is

$\zeta_\mu$2

where $\zeta_\mu$3 are the natural components of the 4PCF and $\zeta_\mu$4 are the derived fourth-order filters (Silvestre-Rosello et al., 9 Sep 2025).

The aperture-mass representation is important because $\zeta_\mu$5 is a filtered projection of the local shear field designed to separate E- and B-modes robustly against masking and systematics. In this setting, the 4PCF becomes the input for fourth-order aperture mass statistics such as $\zeta_\mu$6, which serve as compressed, E-mode-only summaries of non-Gaussian information. The filters decay rapidly for distances much larger than the aperture scale and encode the phase structure required for E/B-mode separation.

The same fourth-order information has a Fourier-space expression through the convergence trispectrum. For equal-scale apertures,

$\zeta_\mu$7

with $\zeta_\mu$8 the Fourier transform of the real-space filter and $\zeta_\mu$9 the convergence trispectrum (Silvestre-Rosello et al., 9 Sep 2025). This makes explicit that the shear 4PCF is the configuration-space counterpart of fourth-order polyspectra.

A recurrent practical implication is that aperture statistics are not identical to the full 4PCF. They are a compression of it. This suggests a trade-off: compression can improve interpretability and E/B-mode control, but it can also reduce access to configuration-level information.

3. Multipole decomposition and efficient estimation

Direct estimation of the shear 4PCF by summing over all quadruplets scales as $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$0 and is intractable for large modern surveys. The recent estimator for fourth-order shear statistics builds on the multipole decomposition of the shear 4PCF, expanding the angular dependence in a basis of complex exponentials. This recasts the computation as a sum over quadratic rather than quartic operations, with scaling

$$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$1

rather than brute-force $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$2 (Porth et al., 9 Sep 2025).

The estimator is formulated in terms of pre-computed catalogue sums with angular weightings and per-object weights. Symmetries of the 4PCF under permutation of legs are exploited to reduce redundant computation, and multiple-counting corrections are applied to ensure that only unique quadruplets are counted. In practice, the multipole components are computed over binned radial and angular configuration space, reflecting the finite sampling imposed by real data.

Further acceleration is obtained from a tree-based implementation. Trees are built by dividing the survey into hierarchical grids, enabling spatially coarse representation of the data for large-scale separations and the use of approximated tree contributions whenever fine detail is not needed. The “DoubleTree”/“BaseTree” strategy is implemented so that both memory and CPU demands are minimized without sacrificing accuracy. Low-memory and high-performance implementations are provided, exploiting threading and spatial partitioning, and the code is publicly available as part of the orpheus package (Porth et al., 9 Sep 2025).

The computational point is substantive rather than merely technical. The estimator’s quadratic scaling is what makes direct application to datasets with $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$3, such as DES Y3 and beyond, practically feasible. A plausible implication is that the observational bottleneck for fourth-order shear statistics shifts from raw tuple enumeration toward validation, covariance modeling, and data compression.

4. Validation, numerical precision, and covariance modeling

Validation has proceeded along two complementary routes. First, the estimator has been tested on mock ellipticity catalogues obtained from Gaussian random fields. In this regime, the measured 4PCF can be compared to the theoretically expected disconnected part derived from the two-point functions using Wick’s theorem. The multipole expansion converges rapidly, with $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$4 sufficing, and the agreement is at the percent level (Porth et al., 9 Sep 2025).

Second, the estimator has been applied to realistic $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$5-body simulations, including the SLICS suite. There, validation uses two estimation pathways: integration of the 4PCF in multipole space and a direct estimator using aperture-mass moments for unmasked fields. Good agreement is found at the 2–10% level over a wide range of scales, both for the connected and disconnected contributions to the fourth-order aperture mass. The estimator is also shown to produce negligible B-mode or parity-violating artifacts (Porth et al., 9 Sep 2025).

A separate line of work addresses the numerical integration pipeline that converts 4PCFs into fourth-order aperture statistics. The calculation is a five-dimensional integral over quadrilateral configurations, performed via Riemann sums over logarithmically spaced radial bins and linearly spaced angular bins. By analyzing and mitigating numerical effects within the integration pipeline, two-percent-level precision on the fourth-order aperture statistics for a Gaussian random field is achieved, remaining well below the noise budget of Stage IV surveys (Silvestre-Rosello et al., 9 Sep 2025).

Covariance estimation is a major issue because fourth-order data vectors are high-dimensional and their sampling distributions are not generically Gaussian. In the Gaussian Random Field approximation, analytic covariance matrices for galaxy $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$6-point correlation functions have been derived in the isotropic basis of Cahn & Slepian, and the formalism is described as particularly powerful for shear and large-scale-structure analyses because only the isotropic, physically meaningful degrees of freedom appear (Hou et al., 2021). Beyond the GRF approximation, two analytical templates for 4PCF covariance include 1-loop corrections from second-order and third-order density fields, second-order and third-order galaxy biasing, and decompositions into a finite number of structural classes; these papers state that the isotropic basis formalism and structural decomposition are directly useful for shear 4PCF calculations, although adjustments are required for the spin nature of shear (Leonard et al., 5 Sep 2025, Leonard et al., 5 Sep 2025).

The practical significance of this covariance work is underscored by the measurement side: the sampling distribution of the fourth-order aperture mass is found to be significantly skewed, and internally estimated uncertainties can underestimate true uncertainties when large apertures are cosmic-variance dominated (Porth et al., 9 Sep 2025).

5. Application to DES Y3 and empirical findings

The most developed observational application so far is to DES Y3. The DES Y3 shape catalogue, containing approximately 100 million galaxies over $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$7, is divided into 100 overlapping patches; 4PCF multipoles are computed with $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$8 up to large angular separations, and tree-based approximations are used where justified to improve efficiency (Porth et al., 9 Sep 2025).

Covariance matrices for the aperture-mass statistics are estimated using 864 DES-like mock footprints from the Takahashi et al. 2017 ray-tracing simulation suite, matched to the real data properties and geometry, including realistic shape noise, weights, and responses. In this framework, the effect of non-Gaussian sampling distributions is quantified, and the fourth-order aperture-mass distribution becomes strongly skewed for large apertures (Porth et al., 9 Sep 2025).

The observational result is a robust, high-significance detection of the connected part of the fourth-order aperture mass in DES Y3. The E-mode is consistently detected, while B-mode, parity-violating, and cross-modes are consistent with zero, implying low systematics. At the same time, the observed E-mode amplitude is significantly below that predicted by the Takahashi et al. simulations, in qualitative agreement with previous findings that DES Y3 displays reduced higher-order lensing signals compared to typical $$\Gamma_0(\vartheta_1,\vartheta_2,\vartheta_3,\phi_{12},\phi_{13}) = \left\langle \gamma(\mathbf{X}_0;\zeta_0)\, \gamma(\mathbf{X}_1;\zeta_1)\, \gamma(\mathbf{X}_2;\zeta_2)\, \gamma(\mathbf{X}_3;\zeta_3) \right\rangle,$$9CDM expectations (Porth et al., 9 Sep 2025).

This empirical success should not be conflated with large immediate gains in parameter constraints. A Fisher forecast for a DES-Y3-like setup reports a minimal improvement in the constraining power of aperture statistics when fourth-order statistics are added to a joint analysis of $\vartheta_i$0 and $\vartheta_i$1, when the analysis is restricted to non-tomographic equal-scale aperture statistics (Silvestre-Rosello et al., 9 Sep 2025). The same study attributes this to correlations and degeneracies between second-, third-, and fourth-order moments, the limited extra sensitivity of equal-scale non-tomographic fourth-order statistics, and the lack of tomographic information.

6. Interpretation, parity, and related higher-order programs

Several interpretive points are now well established. First, fourth-order shear statistics contain part of the non-Gaussian information of the projected matter field and can provide additional constraints on cosmological parameters when combined with second-order statistics, but the amount of additional constraining power depends strongly on the chosen compression, scale dependence, and tomography (Silvestre-Rosello et al., 9 Sep 2025). Second, the connected and disconnected pieces must be treated separately: only the connected 4PCF isolates genuinely new fourth-order information.

Third, parity and mode structure are diagnostically important. In the DES Y3 application, B-mode, parity-violating, and cross-modes are consistent with zero, whereas in broader 4PCF literature parity-odd modes are emphasized as clean tests for parity violation or systematics. In three-dimensional large-scale structure, the parity-odd 4PCF is presented as a null test because standard $\vartheta_i$2CDM contributes no parity-odd signal, and efficient estimators based on isotropic basis functions make such tests practical (Philcox, 2022). For shear, this suggests that parity-odd fourth-order channels are both a systematics diagnostic and a possible target for extensions beyond standard assumptions.

Fourth, the fourth-order program now extends beyond pure cosmic shear. Fourth-order galaxy-galaxy-lensing has been formulated through the galaxy-galaxy-galaxy-shear 4PCF, its relation to the trispectrum, and aperture statistics such as $\vartheta_i$3; an FFT-based direct estimator on pixelized data is reported to detect the connected part of this aperture statistic in a realistic Stage IV mock setup with a signal-to-noise ratio of roughly nine on small aperture scales (Oel et al., 20 Apr 2026). This indicates that fourth-order weak-lensing observables form a broader family of shear-based probes rather than a single statistic.

A common misconception is that any fourth-order detection automatically implies a correspondingly large cosmological information gain. The current results do not support that simplification. What they do show is that fourth-order shear statistics are now measurable, that their sampling properties are non-Gaussian, that they require dedicated covariance modeling, and that their most promising uses may include both cosmological inference and internal consistency tests. A plausible implication is that Stage IV surveys will make the distinction between these roles increasingly important.