4 610 universes — that's exactly how many cosmological models astrophysicists needed to show that the familiar Gaussian mask no longer works for next-generation data. Future sky surveys, such as the Euclid mission and the Vera C. Rubin Observatory, promise maps of dark matter via gravitational lensing — the same dark matter whose existence was discovered by Vera Rubin through the anomalous rotation of galaxies. These maps are built from tiny distortions in the shapes of distant galaxies — weak lensing, which senses the entire large-scale structure of the Universe. From them, we can measure the matter density Ωm and the clumpiness parameter S8, and thereby probe the nature of dark energy and the expansion of the universe discovered by Edwin Hubble and theoretically justified by Georges Lemaître. However, extracting these parameters requires statistical inference. The classical method — explicit likelihood inference (ELI) — puts a 'Gaussian mask' on the data: it assumes all uncertainties follow a normal distribution, and the connection to cosmology is described by a simple emulator trained on numerical simulation. The alternative — likelihood-free inference (LFI) — removes the mask: using Machine Learning, specifically Convolutional Neural Networks, it learns the shape of the distribution directly from examples, without requiring Gaussian dogma.

To settle the debate, scientists built an ideal test bed. Inside the GLASS simulation pipeline, they generated 4 610 independent universes with different combinations of Ωm and σ8 within the ΛCDM model, where gravity is described by general relativity, and the Friedmann equation \[ H^2(a) = H_0^2 \left( \Omega_m a^{-3} + \Omega_\Lambda \right) \] sets the expansion rate. Each universe is a full celestial sphere, divided into pixels with a resolution of about 7 arcminutes. Realistic galaxy shape noise was added to the maps, as in the Euclid final release (30 galaxies per square arcminute), and a patch of 100 square degrees was cut out — a typical size for early surveys. For each universe, two types of statistics were measured: classical shear correlation functions, which, according to gravitational lensing theory, are related to the three-dimensional dark matter power spectrum as \[ C_\ell^{\gamma\gamma} = \int_0^{\chi_H} d\chi \frac{W^2(\chi)}{\chi^2} P_\delta(\ell/\chi) \], and a convolutional neural network estimator trained directly on the maps. The data were then compressed to two numbers: either by the linear MOPED method, optimal for Gaussian fields, or by a nonlinear neural network. Then two inference procedures were run. ELI trained an emulator using a Gaussian process, estimated the covariance matrix, and applied Gaussian likelihood. LFI trained a neural density estimator on most of the universes to learn the shape of the likelihood without assumptions. Both methods underwent Bayesian inference — a fair duel on a single test bed.

The first comparison revealed a drama. For ELI with linear compression, the Gaussian mask cracked: constraints on Ωm were twice as wide as those from LFI, and the distribution shifted away from the truth. The reason was a poorly emulated second compression component, which carried information about matter density and had mild non-Gaussianity. When that component was removed or the same model was fed to LFI, the difference disappeared — both methods converged. With nonlinear compression by a neural network, the Gaussian mask shattered: the distribution of compressed data was no longer normal, and ELI produced errors 2,2 times larger for S8 and 1,3 times larger for Ωm, while LFI retained accuracy. The conflict became even more pronounced when comparing correlation functions with CNN. In an ideal Gaussian world, two-point statistics should contain all information, but ELI, due to its mask, showed that CNN measures Ωm 2 times better and S8 1,6 times worse. LFI, freed from the mask, confirmed theoretical equality: CNN and compressed correlations give consistent results with a difference of only 30%. Calibration tests, usually applied only to LFI, became the main diagnostic tool: they clearly indicated where ELI loses calibration. Traditional tests for Gaussianity (Kolmogorov-Smirnov test, χ²) turned out to be too weak and did not detect dangerous deviations.

The moral is clear: the Gaussian mask can no longer be the standard in the era of high-precision surveys. The work proved that methods without explicit likelihood are not just a trendy toy but a necessary tool for honest data analysis. At the same time, it rehabilitates ELI, showing that with careful emulation and absence of non-Gaussianity, it performs no worse. The main practical takeaway: any future cosmological analysis must be accompanied by rigorous calibration tests; otherwise, we risk mistaking method artifacts for new physical effects. The choice of data compression method is critical for ELI but almost irrelevant for LFI — this gives freedom in designing optimal statistics. The next steps are clear: full N-body simulations with non-Gaussian evolution, tomographic slices, accounting for systematics. The stakes — the ability to distinguish dark energy from modified gravity. And if the mask has fallen on these 4 610 universes, then without it we will see the true structure of the cosmos in the future.