© The Authors 2026

1. Introduction

With the delivery of results and observations from stage IV surveys, such as the data release 1 of the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016; Adame et al. 2025) and the first observations by Euclid (Euclid Collaboration: Mellier et al. 2025), we have entered the era of precision large-scale structure (LSS) cosmology. The key objective of new-generation surveys, such as DESI and Euclid, but also SPHEREx (Crill et al. 2020), as well as the Vera C. Rubin Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2009), is to stress test the standard ΛCDM cosmological model, at late times, measuring the dark energy equation of state, and at early times, detecting primordial features that would serve as the smoking gun of inflation. In particular, inflationary models predict small deviations from Gaussianity in the primordial density distribution, called primordial non-Gaussianities (PNG), whose footprint we ought to be able to observe in the late-time distribution of galaxies.

Different inflation models produce different types of PNG, as described in Baumann (2011). In the present work, we are interested in the PNG of local type, which we parametrised with f_NL. This parameter controls the amplitude of the non-linear relation between a mean-zero Gaussian field, φ, and the primordial gravitational potential, Φ_P = φ + f_NL(φ² − ⟨φ²⟩). We know that local PNG should be very close to zero in the case of single-field inflation (Maldacena 2003; Creminelli & Zaldarriaga 2004; Cabass et al. 2017), whereas we see f_NL ∼ 𝒪(1) for multi-field inflation models (Senatore & Zaldarriaga 2012; Alvarez et al. 2014). Therefore, a detection or non-detection of local PNG with σ_fNL ∼ 1 would rule out single-field inflation or constraint multi-field models, respectively. For this reason, intensive efforts have taken place in the cosmological community to measure and improve our constraining power over this observable using different probes.

Currently, the most stringent constraints on f_NL come from the measurements of the three-point function of the cosmic microwave background (CMB) by Planck, which are consistent with zero at σ_fNL ∼ 5 (Planck Collaboration IX 2020). We expect the new generation of CMB surveys to tighten these bounds by a factor of 2, but without reaching σ_fNL ∼ 1 (Abazajian et al. 2016). On the other hand, local PNG affect the late-time LSS both at the two- and three-point function level. In Fourier space, this corresponds to the power spectrum and the bispectrum (Dalal et al. 2008; Slosar et al. 2008; Desjacques & Seljak 2010). The most stringent bounds from LSS measurements come from the power spectrum analysis of the first data release of DESI, expressed as σ_fNL ∼ 9 (Chaussidon et al. 2025). This corresponds to a 40% improvement over the previous power spectrum-based measurement that used the quasar sample of data release 16 of the extended Baryon Oscillation Spectroscopic Survey (Cagliari et al. 2024) and an improvement of over 50% with respect to previous measurements based on the power spectrum and bispectrum (e.g. D’Amico et al. 2025; Cabass et al. 2022). This offers a first glimpse of the constraining power of stage IV surveys. These results will improve even more when the complete surveys become available.

While we are in the process of gathering new data, we can attempt to improve our constraints on f_NL (or other cosmological parameters) by devising new estimators that retain more information than the standard approach. For example, the authors of Castorina et al. (2019) computed the optimal quadratic estimator to measure the cosmological f_NL signal from the LSS power spectrum multipoles. In addition, many works have combined the power spectrum and the bispectrum (Cagliari et al. 2025), utilising effective field theory (EFT; see Cabass et al. 2023, for a review) to extend the scale range we can probe with these statistics to smaller scales. Alternatively, we could use field-level algorithms that bypass the data compression into n-point statistics. These methods are based on either forward modelling (Andrews et al. 2023, 2026) or machine learning algorithms applied to the late-time galaxy distribution (Kvasiuk et al. 2025), while also reconstructing the initial conditions (Chen et al. 2025). In addition, they can also be applied to the CMB (Nagarajappa & Ma 2024). Several detailed simulation studies have explored the choice (or combinations) of summary statistics to extract PNG information (Coulton et al. 2023b; Jung et al. 2023a; Coulton et al. 2024; Jung et al. 2024) and possibly guide the community effort in this area of study.

The work of Giri et al. (2023) produced highly informative estimates of f_NL from the non-linear matter field by training a neural network to compute local estimates of σ₈ and modelling auto- and cross-spectra of the σ₈ field and the matter. The methodology relies on a bias model specifically for local non-Gaussianity, which neatly bypasses the need for running training simulations with non-Gaussian initial conditions. However, it is not clear how this approach might be generalised to other forms of primordial non-Gaussianity.

In this work, we propose a method that approximates global field-based analysis by hierarchically combining standard n-point statistics with a local, neural summary that extracts field-level information. The main idea is that we can use the power spectrum and bispectrum optimally to extract large-scale information, along with machine learning-based data compression to extract small-scale information in the non-linear regime. Given a large cosmological volume, we expect that on large scales the majority of the f_NL information will be encoded in the two- and three-point functions; on the other hand, on small scales, a field-level algorithm would be expected to extract more information than the power spectrum and bispectrum alone (Bairagi & Wandelt 2026). Aside from the application to non-Gaussianity, we can go beyond the dark matter field and focus on catalogues of halos. Neither of them are directly observable, but we note that the halo field is a better proxy than the dark matter field in terms of a signal that is achievable with observations.

This hierarchical combination of global statistics with local field-level analysis is theoretically well-motivated for the problem of disentangling non-linear evolution from primordial non-Gaussianity. Recent works have demonstrated that locality fundamentally protects primordial signals from contamination by late-time gravitational nonlinearities (Baumann & Green 2022). Specifically, primordial non-Gaussianity creates genuinely non-local correlations between widely separated points (generated during inflation), while late-time gravitational evolution is constrained by causality to produce only local effects that cannot mimic these long-range correlations.

This locality principle has important implications for data analysis. A Fisher analysis has shown that higher order bias parameters become increasingly orthogonal to primordial signals in the perturbative regime and that is only quadratic nonlinearities that affect the map-level analysis and all higher order non-linearities end up decoupled (Baumann & Green 2022). This suggests that combining two-point and three-point statistics with local field-level information naturally incorporates the constraints from higher order correlations without requiring their explicit computation. The map-level approach thus captures the essential information content on primordial non-Gaussianity more efficiently than computing progressively higher order correlation functions.

In this spirit, our algorithm measures the power spectrum and bispectrum of the whole observational volume, while it runs a field-level analysis on patches. We analyse the patches with a convolutional neural network (CNN; LeCun et al. 1989), named PATCHNET (Bairagi & Wandelt 2026), which summarises the field-level information as a single number and then aggregates over all the patches to get a single-value compression. Then, the final summary statistic is the combination of the two- and three-point summaries with the field-level compression. This combination was shown in Bairagi & Wandelt (2026) to produce state-of-the-art information extraction for cosmological parameters in the non-linear regime.

On a practical level, this approach has the advantages of a machine learning field-level analysis without the shortcoming of being too computationally expensive. The network only acts on small sub-volumes of the full computational volume (dubbed ‘patches’) and it greatly reduces the number of simulations required to train the networks, compared with a full field-based analysis, since dividing the volume into patches greatly increases the training set size.

In this paper, we estimate the f_NL information content of this new estimator applied to the QUIJOTE-PNG dark matter halos (Coulton et al. 2023a) in a cube of 1 Gpc³ h⁻³. We perform a detailed study of the Fisher information of combining the power spectrum with local patches and the combined power spectrum, bispectrum, and patches as a function of redshift, halo mass cut, and scale cut for the standard summary statistics.

The main finding of this work is in line with the theoretical expectations, asserting that patch-based data compression consistently provides additional information compared to standard summary statistics. The most significant improvements occur at the lowest k_max for the n-point summary statistics, with gains ranging between 30% and 45%. Furthermore, when comparing the power spectrum and bispectrum constraints to those derived from patch-combined estimators, the patches consistently enhance the results, demonstrating they capture information beyond the three-point function.

The paper is organised as follows. In Sect. 2, we present the simulated data we employed in the analysis. In Sect. 3, we introduce the Fisher formalism and the n-point statistics we used, along with a description of the architecture of PATCHNET. We discuss our results in Sect. 4 and present our conclusions in Sect. 5.

2. Simulations

The QUIJOTE simulations (Villaescusa-Navarro et al. 2020) are a set of N-body simulations run with GADGET-III, which was originally developed for the Aquarius project (Springel et al. 2008) and is an improved version of the GADGET-II code (Springel et al. 2005). They are specifically designed to test machine learning applications on cosmological data. They start from initial conditions at z = 127 generated with the 2LPTIC code¹ and simulate five redshift snapshots (z = 0, 0.5, 1, 2, and 3) with a co-moving volume of 1000³ Mpc³ h⁻³ and a fiducial resolution of 512³ particles². The main data products of the QUIJOTE simulations are the particle snapshots over which two halo finders were run. For this work, we used the halo catalogues produced with the friend-of-friends (FoF; Davis et al. 1985) halo finder.

The original QUIJOTE dataset spans over a range of ΛCDM and νwCDM cosmologies. Recently, a new set of simulations has become available, the QUIJOTE-PNG (Coulton et al. 2023a) simulations, which cover four different types of PNG: LSS orthogonal, CMB orthogonal, equilateral, and local PNG. The initial conditions of the QUIJOTE-PNG simulation were generated with the 2LPTPNG code (Coulton et al. 2023a)³. For this work, we only used the datasets related to local PNG. For the training of PATCHNET, we used the 1000 f_NL local Latin hypercube simulations. In this set, f_NL uniformly spans the interval ( − 300, 300). Then, we used 5000 realisations of the fiducial cosmology for the Fisher information study, where f_NL = 0. Finally, we also used the two sets of 500 simulations with f_NL displaced from the fiducial value. These simulations have f_NL = −50 or f_NL = 50 and are used to estimate the derivative of the summary statistics as a function of f_NL. For all these simulations, except for f_NL, all the other cosmological parameters are fixed to the fiducial values, Ω_m = 0.3175, Ω_b = 0.049, h = 0.6711, n_s = 0.9624, σ₈ = 0.834, w = −1, and M_ν = 0 eV.

Considering the redshift range probed by current surveys (e.g. the Euclid Wide Survey; Euclid Collaboration: Scaramella et al. 2022; Euclid Collaboration: Mellier et al. 2025), we are mainly interested in the results at z = 1, which should better represent the clustering observed by this experiment. In addition, we studied the redshift dependence of the algorithm performance at the other four redshifts available in the Quijote dataset, z = 0,  0.5,  2, and 3. We tested two scenarios: the first is the case where the same mass cut is applied at all redshifts (i.e. we selected halos with M_h > 3.12 × 10¹³ M_⊙ h⁻¹) and the second case is where the number density is constant at different redshifts. In this latter case, at every redshift we chose a halo mass cut, M_halo > M_min, that would fix the number density of the halos in the fiducial cosmology to the number density at z = 3 (n(z)∼5 × 10⁻⁶ h³ Mpc⁻³), which is the least populated snapshot. We present the threshold masses, M_min, in Table 1.

3. Methods

3.1. Fisher Information

It is possible to evaluate the amount of information and the constraining power of an estimator through the Fisher information analysis. The Fisher information matrix is defined as

$$\begin{array}{c}\hfill F_{\mathit{ij}}={\left(\frac{\partial \mathbf{x}}{\partial {\theta }_i}\right)}^T{\mathbf{C}}^{-1}\left(\frac{\partial \mathbf{x}}{\partial {\theta }_j}\right),\end{array}$$(1)

where x is a summary statistic of the optimal unbiased estimator of the parameters θ_i and θ_j, and C⁻¹ is the inverse of the covariance matrix of the summary statistics, x. The expected variance of parameter i is σ_i² = (F⁻¹)_ii. Therefore, to estimate the constraining power of a given summary statistic, we first need to compute the covariance matrix of the estimator and its derivatives with respect to the parameters of interest, which in this study is only f_NL. We estimated both C⁻¹ and $\frac{\partial \mathbf{x}}{\partial f_{\mathrm{NL}}}$ numerically. We used 5000 QUIJOTE simulations in the fiducial cosmology to estimate the covariance matrix, ${\widehat{\mathbf{C}}}^{-1}$, which we also corrected with the Hartlap factor (Hartlap et al. 2007),

$$\begin{array}{c}\hfill {\mathbf{C}}^{-1}=\frac{n_{\mathrm{s}}-n_{\mathrm{x}}-2}{n_{\mathrm{s}}-1}\widehat{}\}{C}^{-1},\end{array}$$(2)

to account for the non-infinite number of data we use to estimate the covariance matrix. In Eq. (2), n_s is the number of simulations used to compute the covariance and n_x is the length of the summary statistic. Finally, we estimate the derivatives numerically from the finite difference,

$$\begin{array}{c}\hfill \frac{\partial \mathbf{x}}{\partial f_{\mathrm{NL}}}=\frac{\overline{\mathbf{x}}(f_{\mathrm{NL}}^{\text{fid}}+\delta f_{\mathrm{NL}})-\overline{\mathbf{x}}(f_{\mathrm{NL}}^{\text{fid}}-\delta f_{\mathrm{NL}})}{2\delta f_{\mathrm{NL}}},\end{array}$$(3)

where, $\overline{\mathbf{x}}(f_{\mathrm{NL}}^{\text{fid}}\pm \delta f_{\mathrm{NL}})$ are the mean summary statistics of the simulations with f_NL = ±50. Then, as $f_{\mathrm{NL}}^{\text{fid}}=0$, δf_NL = 50⁴.

3.2. Summary statistics: Power spectrum and Bispectrum

In this work, we aim to combine large- and small-scale information to improve the constraining power over local f_NL measurements. A first solution is to use effective field theory to increase the minimum scale that we can probe with standard summary statistics. The alternative approach we explored is to use field-level data compression based on machine learning to extract small-scale information, while the large-scale information still comes from standard summary statistics.

The summary statistics we utilised are the power spectrum and the bispectrum. Given the halo over-density, δ_h, the power spectrum corresponds to the two-point correlation function in Fourier space and we define it as

$$\begin{array}{c}\hfill ⟨{\delta }_{\mathrm{h}}(\mathbf{k}){\delta }_{\mathrm{h}}^{\ast }({\mathbf{k}}^{\prime })⟩={(2\pi )}^3{\delta }^{\mathrm{D}}(\mathbf{k}-\mathbf{k}\prime )P(k),\end{array}$$(4)

where δ^D is the 3D Dirac’s delta function. Similarly, we can define the halo bispectrum, which is the three-point correlation function in Fourier space, as

$$\begin{array}{c}\hfill ⟨{\delta }_{\mathrm{h}}({\mathbf{k}}_1){\delta }_{\mathrm{h}}({\mathbf{k}}_2){\delta }_{\mathrm{h}}({\mathbf{k}}_3)⟩={(2\pi )}^3{\delta }^{\mathrm{D}}({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3)B({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3).\end{array}$$(5)

To measure the power spectra and bispectra of the simulated boxes, we used the public codes Pylians3 (Villaescusa-Navarro 2018) and pySpectrum (Scoccimarro 2015; Hahn et al. 2020), respectively. We measured the power spectra on a 512³ grid from the fundamental wave number of the box, k_f ∼ 0.006 h Mpc ⁻¹, to k_max ∼ 0.8 h Mpc ⁻¹ with a k_f linear binning; on the other hand, we used a 256³ grid for the bispectrum, measured with a 3 k_f linear binning from k_min = 3 k_f to k_max ∼ 0.5 h Mpc ⁻¹.

3.3. Field level compression with PATCHNET

We employed field-level data compression, using a CNN to extract small-scale information following the approach of Bairagi & Wandelt (2026). The network does not analyse the whole simulated box; it only analyses a sub-cube, which we dubbed a ‘patch’. This network, called PATCHNET, offers two main advantages in its approach. First, it reduces the memory requirement for the CNN while maintaining the pixel (or better voxel) resolution, as it processes a smaller volume. Second, the number of training data greatly increases as we can cut out a set of N_sb patches from one simulated box to train the network, effectively increasing the amount of training data by a factor of N_sb. As we will explain below, we trained PATCHNET to compress each patch into one informative statistic (i.e. the value of f_NL) and we then aggregated these statistics by averaging across the N_sb patches of a given simulated box. Finally, we concatenated this aggregated information with the power spectrum or the combination of the power spectrum and bispectrum.

PATCHNET is characterised by the following architecture: first, there are three blocks composed of a convolutional layer and a 3D average pooling layer. Their output is flattened and processed by three dense layers before the final output. The convolutional layers, as well as the average pooling, have a 3 × 3 × 3 kernel with a stride of 1 and 0 padding. The first convolutional layer has eight filters in terms of output and the filters are increased by a factor of 4 in the following layers. After flattening, the latent features are compressed into 512 neurons and subsequently halved in terms of their number by the other dense layers. Both the convolutional and dense layers have a LeakyReLU (Maas et al. 2013) with a negative slope of 0.5 as activation function. We used a modified version of the mean squared error as a loss function,

$$\begin{array}{c}\hfill \text{Loss}(\mathbf{y},\widehat{}\}y)={\displaystyle \sum _i^{n_{\mathrm{b}}}}\frac{1}{n_{\mathrm{b}}}{({\displaystyle \sum _j^{n_{\text{sb}}}}\frac{1}{n_{\text{sb}}}{y}_{j,i}-{\displaystyle \sum _j^{n_{\text{sb}}}}\frac{1}{n_{\text{sb}}}{\widehat{y}}_{j,i})}^2,\end{array}$$(6)

where y represents the outputs of PATCHNET, $\widehat{\mathbf{y}}$ the labels (f_NL), n_b is the batch size (i.e. the number of simulations we loaded), and n_sb is the number of patches we use for each simulation. We used n_b = 32 and n_sb = 64; therefore, a batch actually contains n_b × n_sb different inputs.

Before inputting the halo catalogues to the network, we pre-processed them. First, we divided the 1000 simulations with varying f_NL into the training, validation, and test sets, which we selected based on a standard split of 75%, 15%, and 10% of the data, respectively. Then, we built the patches from the simulated halo catalogue density contrast, which we computed by painting the halos on a 128³ grid with the cloud in cell (CIC) algorithm. We note that during the painting process, we did not use any halo mass information. We divided the 3D number density contrast into 16³ non-overlapping cubes⁵. This resulted in 512 patches for each simulation that have a physical dimension of 125³ Mpc³ h⁻³ and a resolution of 8 Mpc h⁻¹. In Fourier space, these scales correspond to k_min ∼ 0.05 h Mpc ⁻¹ and k_max ∼ 0.4 h Mpc ⁻¹. Finally, before feeding the density contrast patches to the network, we centred and scaled them with the mean and standard deviations computed from the density contrast of the fiducial cosmology realisations. We note that we computed these two quantities for each redshift and mass cut configuration (see details in Sect. 2).

4. Results

In this section, we present and discuss the results of the Fisher information analysis. In particular, we compare the information content of the power spectrum and bispectrum with their combination with the patch information.

Figure 1 shows the results of the Fisher information analysis at redshift z = 1 for two mass cuts as a function of the k_max of the power spectrum and bispectrum. Adding patch information always produces tighter bounds than the standard summary statistics alone. Moreover, we also observe an improvement in the combination of the power spectrum and the patches (solid orange line with triangular markers) with respect to the combined power spectrum and bispectrum (dash-dotted green line with reverse triangular markers), showing that PATCHNET extracts information not captured by the two- and three-point functions. This is also backed up by the fact that even though the improvement due to the patch information diminishes when k_max increases, as the minimum scale probed by the standard summary statistics gets closer to the minimum scale of the patches ($k_{\max }^{\mathrm{patch}}\sim 0.4h{\text{Mpc}}^{-1}$), we still observe an improvement when $k_{\max }^{\mathrm{patch}}<k_{\max }^{P,B}$, which indicates once again that the patches use information beyond the three-point function. Nevertheless, the bispectrum still brings in a small amount of information that the network cannot extract, as the combination of the three statistics (see the sparse dash-dotted purple line with pentagonal markers) is always slightly lower than the combination of the power spectrum and patches alone. This behaviour is not unexpected as the bispectrum has access to the very large-scale squeezed triangle configuration, which contains non-Gaussian information that the patches cannot access due to their limited size.

Fig. 1. Expected uncertainty for the different summary statistics at z = 1 for two halo-mass cuts as a function of the power spectrum and bispectrum kmax. In both cases, the PATCHNET information (patches) improves the constraining power of the power spectrum and the power spectrum combined with bispectrum analyses.

Comparing the left and right panels, we see that the overall information content of the smaller mass cut (M_min = 3.2 × 10¹³ M_⊙ h⁻¹ in the left panel) is higher than the larger mass cut (M_min = 1.05 × 10¹⁴ M_⊙ h⁻¹ in the right panel). This is because when the mass cut is smaller, the halo field has a lower shot noise and contains more information. However, the improvement related to the patch information is greater in the case of the larger mass cut for k_max ≥ 0.2 h Mpc ⁻¹. This suggests that in terms of real objects, this method remains effective for the rarer, more highly biased objects that usually constitute the large volume required for f_NL measurements. The dependence of the performance on the bias would need to be explored with a more detailed tracer-specific study.

We present the constraining power of the different statistics as a function of redshift in Figs. 2 and 3. As shown in Fig. 2, we kept the M_min constant, while Fig. 3 shows how we varied it as a function of redshift to keep the number density of the catalogues equal to the halo number density at z = 3. The feature that is common to both plots is that the addition of the patch information always makes the expected error at k_max = 0.1 h Mpc ⁻¹ comparable to (if not better than) the constraint on the standard summary statistics alone at k_max = 0.5 h Mpc ⁻¹. This is an interesting result, considering that modelling the power spectrum and bispectrum becomes more involved as we enter the non-linear regime, k_max > 0.1 h Mpc ⁻¹, requiring a large number of simulations (or EFT) to model the observed data. The use of patch information can instead leverage the numerical simulation model while requiring a smaller simulation budget than neural estimators acting on the full field. Nevertheless, to provide a definitive statement about the advantage of the PatchNet approach over EFT modelling, we would need to perform additional tests. In particular, we would have to study the robustness of the algorithm against different halo finders and on non-linear and baryonic physics.

Fig. 2. Expected uncertainty for the different summary statistics as a function of the power spectrum and bispectrum kmax at different redshift snapshots. At all redshifts, we applied the same halo-mass cut, Mh > 3.2 × M⊙ h−1. The patch information adds constraining power at every redshift.

Fig. 3. Expected uncertainty for the different summary statistics at fixed number density over redshifts as a function of the power spectrum and bispectrum kmax. At each redshift, we selected a halo-mass cut that has in the fiducial cosmology the same number density of z = 3 (see Table 1). As in Fig. 2, the patch information improves the constraining power of the standard summary statistics.

Finally, the authors of Jung et al. (2023b) have shown that there is a high information content related to f_NL in the halo mass function. As some halo mass information is retained in the density contrast, PATCHNET may be just learning the relation between the halo mass function and f_NL, making the field-level analysis superfluous. To check whether this is indeed the case, we trained PATCHNET keeping the number of objects in the catalogue fixed. The difference between this analysis and the one with constant density is that in the latter, the mass cut was chosen at the fiducial cosmology (with f_NL = 0) and used for all the other simulated cosmologies. Therefore, there was a variation in the number density as a function of f_NL. On the other hand, for this test, we kept the number of objects fixed even when we vary the value of f_NL. In doing so, we removed the halo mass function information. Even in this configuration, PATCHNET learns a relationship between the density contrast field and f_NL. In this case, the Pearson correlation coefficient is r_{w/o HMF} ∼ 0.78, which we compared with r_{w HMF} ∼ 0.99 for the standard configuration. This shows that the HMF helps improve the precision of the network, as r_{w HMF} > r_{w/o HMF}. Nevertheless, the network can still extract PNG information even when we remove the HMF from the field, as r_{w/o HMF} ≫ 0.5.

5. Conclusions

In this work, we propose an alternative approach to EFT up to large k_max or more standard ML-based field-level algorithms to analyse galaxy redshift surveys. The algorithm combines the large-scale information encoded in the n-point summary statistics, such as the power spectrum or bispectrum, with the small-scale information obtained by averaging over the compressed field-level information of small patches of the simulated box, the PATCHNET approach (Bairagi & Wandelt 2026). The main advantage of this approach is that the analysis of small patches is computationally lighter than analysing the whole box; additionally, by cutting the simulated boxes into patches, we can greatly increase the number of training examples for the PATCHNET (Bairagi et al. 2025). We estimated the Fisher information of this combined summary statistics to extract local PNG information.

We tested the algorithm on the FoF halo catalogues of the QUIJOTE-PNG simulations. We compared its performance with the more standard analyses based on the power spectrum or the combination of the power spectrum and bispectrum. We performed this comparison for different values of k_max for the n-point statistics and studied it as a function of redshift and for different halo mass cuts.

Our preliminary study shows that the patch-based data compression always provides additional information relative to the standard summary statistics. We observed the greatest improvements when the power spectrum and the bispectrum were truncated at $k_{\max }=0.1h{\text{Mpc}}^{-1}$, where the small-scale structure in the patches provides the most complementary information. At this scale (depending on the redshift, mass cut, and the statistics combination), the improvement varies between 30% and 45%. At k_max = 0.4 h Mpc ⁻¹, where the patches and the standard summary statistics get to the same minimum scales, the improvement falls to a more moderate range between 15% and 25%, but it is still present. Moreover, when we compare the power spectrum plus bispectrum bounds with the patch-combined bounds (both power spectrum with patches and power spectrum plus bispectrum with patches), the patches once again always bring an improvement, proving that they contain information beyond the three-point function. Nevertheless, we also observed that combining the bispectrum with the PATCHNET compression also provides an improvement, demonstrating that the large-scale bispectrum encodes information that the patches end up missing due to their limited size. In this regard, it is possible that the use of larger volumes (thereby increasing access to the squeezed bispectrum information) could actually reduce the improvement observed when adding the patch information, especially in the case of the local f_NL. On the other hand, local higher order information should help disentangle the non-linear coupling due to non-linear gravity from primordial non-Gaussianity (Baumann & Green 2022). The use of larger volumes will also allow us to test the effect super-sample covariance has on the method, even though we expect this effect to be subdominant at the patch-level, as the patches are already sub-volumes of the bigger simulation box. We would already be observing the effect of super-sample covariance.

A first extension to this work would be to test our method’s performance for other forms of PNG and the effect of marginalising over cosmological parameters. Second, we need to perform robustness tests to check whether the method has the same performance with different and more accurate halo finders. Third, to fully understand the potential of the patch-based approach, it would be of interest to compare its information content with other alternative summary statistics; for instance, a wavelet scattering transform (e.g. Peron et al. 2024) or a marked power spectrum (e.g. Marinucci et al. 2025). In addition, to prepare this approach for applications to real data, we will need to introduce a galaxy model in the simulations and add the lightcone, the survey mask, and systematic effects to the simulations. We plan to study these issues in future work.

Acknowledgments

We thank Francisco Villaescusa-Navarro, William Coulton, and Michele Liguori for their feedback on the first version of this manuscript. MSC thanks ‘Fondazione Angelo della Riccia’ for financial support. The work of MSC is supported by the Agence Nationale de la Recherche (ANR) grant n. ANR-23-CPJ1-0160-01. AB acknowledges support from the Simons Foundation as part of the Simons Collaboration on Learning the Universe. The Flatiron Institute is supported by the Simons Foundation. All neural training and post-analysis have been done on IAP’s Infinity cluster. The authors thank Stéphane Rouberol for his efficient management of this facility.

References

All Tables

All Figures

	Fig. 1. Expected uncertainty for the different summary statistics at z = 1 for two halo-mass cuts as a function of the power spectrum and bispectrum kmax. In both cases, the PATCHNET information (patches) improves the constraining power of the power spectrum and the power spectrum combined with bispectrum analyses.
In the text

	Fig. 2. Expected uncertainty for the different summary statistics as a function of the power spectrum and bispectrum kmax at different redshift snapshots. At all redshifts, we applied the same halo-mass cut, Mh > 3.2 × M⊙ h−1. The patch information adds constraining power at every redshift.
In the text

	Fig. 3. Expected uncertainty for the different summary statistics at fixed number density over redshifts as a function of the power spectrum and bispectrum kmax. At each redshift, we selected a halo-mass cut that has in the fiducial cosmology the same number density of z = 3 (see Table 1). As in Fig. 2, the patch information improves the constraining power of the standard summary statistics.
In the text