Analytic bispectrum covariance for galaxy survey cosmology

The bispectrum of galaxy distribution remains an underutilized statistic in cosmological data analysis. This can be attributed to mathematical and computational challenges associated with a proper modeling of the signal and its covariance. While recent theoretical advances have made it possible to analytically model the bispectrum signal and covariance using perturbation theory, a fully non-Gaussian analytical treatment of the covariance for a realistic survey with window remains elusive. The full bispectrum covariance of observed galaxies is a non-Gaussian six-point function with each of its legs convolved with a complex window function. Although it has so far proven intractable, the problem is well-posed. In this article, we present the derivation, implementation, and validation of a numerical routine for the analytic calculation of the bispectrum covariance for a spherical window, carried out by a large language model (LLM). The derivation uncovered terms that vanish in the windowless case but that dominate the covariance for triangles with survey scale modes. As a development benchmark for the implementation, we used the sample covariance of the halo bispectrum measured from 6000 fiducial Quijote halo catalogs in real space. Over several iterations of refinement, the LLM converged on a regime-specific numerical scheme for the covariance matrix calculation, mixing exact and approximate approaches, resulting in a highly optimized, tractable and accurate covariance across scales and shapes, using only tree-level models for the input polyspectra. Closing the loop still required human expertise, which illustrates what human-AI collaboration in research may look like in the near term.