If at First You Don't Succeed, Trispectrum: I. Estimating the Matter Power Spectrum Covariance with Higher-Order Statistics

Abstract: We present a method to estimate non-Gaussian power spectrum covariance matrices by directly measuring the response of the small-scale power spectrum to long-wavelength perturbations via bispectrum and trispectrum estimators. Specifically, we derive estimators for the complete non-Gaussian matter power spectrum covariance, including the super-sample contribution, in terms of the squeezed bispectrum and collapsed trispectrum of the underlying density field. We apply these estimators to the Quijote simulations, and recover unbiased estimates of the small-scale ($k\gtrsim 0.15~h/{\rm Mpc}$) matter power spectrum covariance at the percent level using only 25 simulations - comparable to the precision of the sample covariance estimated using 5,000 simulations. This technique significantly reduces the number of simulations needed to estimate power spectrum covariances and opens the possibility of inferring power spectrum covariances directly from survey data, enabling stringent tests of simulations and, potentially, power spectrum analyses that do not rely on external covariance matrices.

• The paper introduces a trispectrum-based method that estimates the full non-Gaussian covariance using squeezed bispectrum and collapsed trispectrum.
• The methodology achieves percent-level accuracy with only 25 simulations, contrasting the thousands required by traditional sample covariance methods.
• The approach cleanly separates connected non-Gaussian and super-sample covariance components, paving the way for direct data-driven analysis in future surveys.

Estimating the Non-Gaussian Matter Power Spectrum Covariance with Higher-Order Statistics

Introduction

Precise estimation of the matter power spectrum covariance is a fundamental technical requirement in cosmological parameter inference from large-scale structure surveys. While the Gaussian approximation holds on large linear scales, the covariance structure becomes intricate on non-linear scales due to mode coupling and non-Gaussianity, with contributions dominated by the connected trispectrum. Traditional approaches rely on brute-force sample covariance estimation using thousands of $N$-body simulations, incurring prohibitive computational costs. This work develops and demonstrates a method for extracting the full non-Gaussian covariance (including both connected and super-sample components) directly from higher-order statistics—the squeezed bispectrum and collapsed trispectrum—substantially reducing simulation requirements and enabling the prospective empirical inference of covariances from survey data itself ["If at First You Don't Succeed, Trispectrum: I. Estimating the Matter Power Spectrum Covariance with Higher-Order Statistics" (2604.18581)].

Theoretical Background

The small-scale covariance of the matter power spectrum receives significant contributions from the connected four-point function (trispectrum), especially beyond the quasi-linear regime where analytic models based on perturbation theory break down. Two key higher-order statistics are exploited:

• Squeezed bispectrum: $B_m(q, k)$ probes the response of small-scale power to long-wavelength overdensity ($q \ll k$). Its leading contribution is entirely determined by the matter consistency relations, specifically by the response function $\partial P(k) / \partial \delta_L$.
• Collapsed trispectrum: $T_m(q, k, k')$ correlates power spectrum estimates in spatially separated regions modulated by a long-wavelength mode. In the soft limit ($q \rightarrow 0$), the trispectrum decomposes into a $q$-independent “connected non-Gaussian” term and a $q$-dependent super-sample covariance term.

All contributions to the power spectrum covariance can thus be reconstructed as explicit functionals of squeezed and collapsed higher-point measurements, leveraging symmetry-enforced dependences (large-scale structure consistency relations) that persist into the highly non-linear regime.

Figure 1: Summary statistics from the Quijote N-body ensemble: large-scale power spectrum, squeezed bispectrum, and collapsed trispectrum, all highly correlated realization-by-realization and probing the same long-wavelength modes.

Methodology

The estimation pipeline proceeds as follows:

1. Summary statistic computation: Fast Fourier Transform-based estimators are employed to measure the large-scale binned power spectrum, squeezed bispectrum, and collapsed trispectrum from a limited ensemble of $N$-body simulations. Hard modes $k$ span the non-linear regime, while soft modes $B_m(q, k)$0 are chosen strictly within the linear regime to guarantee the validity of consistency relations.
2. Parameter inference: Theoretical modeling, validated by simulation outputs, expresses these statistics as functions of:
   • The leading power spectrum response $B_m(q, k)$1,
   • The $B_m(q, k)$2-independent part of the trispectrum $B_m(q, k)$3. These parameters, once inferred, are mapped into the connected non-Gaussian and super-sample covariance components.
3. Statistical analysis: A full joint likelihood (multivariate $B_m(q, k)$4-distribution to account for finite sample estimation) is constructed for the measured statistics across bins. Sample variance cancellation between the correlated statistics ensures optimal error reduction.

A notable technical point is the estimator construction for the collapsed trispectrum. The disconnected (Gaussian) contributions are subtracted via realization-dependent estimators, which are numerically robust to uncertainty in the fiducial power spectrum even when estimated with a single simulation.

Figure 2: Marginalized posteriors for covariance parameters derived from joint fits to power spectrum, squeezed bispectrum, and collapsed trispectrum, showing unbiased, percent-level constraints with only 25 simulations.

Results

Application to the Quijote N-body simulation suite establishes several key empirical benchmarks:

• Efficiency: Using only 25 $B_m(q, k)$5 simulations, the method achieves percent-level accuracy in the recovered power spectrum covariance for $B_m(q, k)$6, matching the precision of sample covariance estimation with approximately 5,000 simulations.
• Separation of covariance components: The connected non-Gaussian and super-sample contributions are simultaneously and cleanly separated, avoiding the biases inherent in jackknife and other internal estimators which mix these physical sources.
• Sample variance cancellation: The strong correlation between higher-order statistics for the same realization yields dramatic error reduction (factor $B_m(q, k)$715 improvement over direct sample variance estimation for the connected non-Gaussian covariance). The joint fitting of all statistics is essential for optimality.

Figure 3: Comparison of non-Gaussian covariance estimates via joint higher-order statistics (red) and sample covariance (gray) with the "true" result from thousands of simulations; the higher-order approach matches the precision of much larger simulation sets.

Figure 4: Correlation matrices for the matter power spectrum covariance: sample covariance (noisy, right), true result (left), and the collapsed trispectrum method (middle), which reproduces the correlation structure to percent-level accuracy from vastly fewer realizations.

Robustness checks validate the pipeline using Gaussian random fields (where the estimator correctly yields zero trispectrum and matches analytic Gaussian covariance), and estimator bias due to fiducial power spectrum inaccuracies is shown to be negligible for the realization-dependent subtraction.

Practical and Theoretical Implications

• Simulation resource savings: The technique enables full non-Gaussian covariance estimation with two orders of magnitude fewer simulations, relaxing computational constraints for future high-resolution survey analyses.
• Direct data-driven covariance: For sufficiently large observational volumes, measurement of the power spectrum, bispectrum, and trispectrum from the survey itself allows for covariance estimation without reliance on external simulations or analytic approximations.
• Separation of physical effects: The method explicitly partitions the super-sample and connected non-Gaussian covariance components, in contrast to the mixing intrinsic to partition-based estimates (jackknife, bootstrap), thus eliminating sources of systematic bias.
• Applicability: While the current work is presented for the 3D matter density field in real space, the underlying methodology is structurally general and can extend (modulo additional complications such as shot noise, bias, and redshift-space distortions) to halos, galaxies, and projected fields (e.g., cosmic shear, CMB lensing).

Future Prospects

Immediate extensions would address:

• Generalization to biased tracers (halos, galaxies) and incorporation of redshift space distortions.
• Projection to angular (spherical) statistics relevant for weak lensing and CMB surveys.
• Incorporation of survey geometry and systematics via window function deconvolution.
• Development of analogous approaches for bispectrum and higher-point function covariance estimation.

The formalism is also naturally aligned with approaches exploiting large-scale structure consistency relations for primordial non-Gaussianity estimation at the field level.

Conclusion

This work introduces and validates a collapsed trispectrum-based methodology for accurate, simulation-efficient estimation of the non-Gaussian matter power spectrum covariance. By exploiting the information content of squeezed higher-order statistics and leveraging large-scale structure symmetries, percent-level precision is realized with minimal computational overhead. This framework provides a foundation for data-driven covariance analysis in next-generation surveys and offers a robust pathway for systematic error control in cosmological inference.

Figure 5: Ratio of trispectrum-based and sample-covariance-based estimates to the true connected non-Gaussian covariance, demonstrating the accuracy and superiority of the higher-order method across all relevant scales.