# Constraining Primordial Non-Gaussianity from Large Scale Structure with the Wavelet Scattering Transform

Matteo Peron, Gabriel Jung, Michele Liguori, Massimo Pietroni

► To cite this version:

Matteo Peron, Gabriel Jung, Michele Liguori, Massimo Pietroni. Constraining Primordial Non-Gaussianity from Large Scale Structure with the Wavelet Scattering Transform. Journal of Cosmology and Astroparticle Physics, 2024, 07, pp.021. 10.1088/1475-7516/2024/07/021. hal-04541235

HAL Id: hal-04541235

https://hal.science/hal-04541235v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Constraining primordial non-Gaussianity from large scale structure with the wavelet scattering transform

To cite this article: Matteo Peron et al JCAP07(2024)021

View the [article online](#) for updates and enhancements.

You may also like

- *Primordial Non-Gaussianity Systematics from Redshift Mismatch with SPHEREX*
  Chandra Shekhar Saraf and David Parkinson

- *A STUDY OF HIGH-ORDER NON-GAUSSIANITY WITH APPLICATIONS TO MASSIVE CLUSTERS AND LARGE VOIDS*
  Sirichai Chongchitnan and Joseph Silk

- *Quijote-PNG: The Information Content of the Halo Mass Function*
  Gabriel Jung, Andrea Ravenni, Marco Baldi et al.
RECEIVED: March 28, 2024

REVISED: May 15, 2024

ACCEPTED: June 10, 2024

PUBLISHED: July 8, 2024

# Constraining primordial non-Gaussianity from large scale structure with the wavelet scattering transform

Matteo Peron<sup>ID,a</sup> Gabriel Jung<sup>ID,b</sup> Michele Liguori<sup>ID,c</sup> and Massimo Pietroni<sup>ID,a</sup>

E-mail: matteo.peron@unipr.it, gabriel.jung@universite-paris-saclay.fr, michele.liguori@unipd.it, massimo.pietroni@unipr.it

**ABSTRACT:** We investigate the Wavelet Scattering Transform (WST) as a tool for the study of Primordial non-Gaussianity (PNG) in Large Scale Structure (LSS), and compare its performance with that achievable via a joint analysis with power spectrum and bispectrum (P+B). We consider the three main primordial bispectrum shapes — local, equilateral and orthogonal — and produce Fisher forecast for the corresponding $f_{\text{NL}}$ amplitude parameters, jointly with standard cosmological parameters. We analyze simulations from the publicly available QUIJOTE and QUIJOTE-PNG N-body suites, studying both the dark matter and halo fields. We find that the WST outperforms the power spectrum alone on all parameters, both on the $f_{\text{NL}}$'s and on cosmological ones. In particular, on $f_{\text{NL}}^{\text{local}}$ for halos, the improvement is about 27%. When B is combined with P, halo constraints from WST are weaker for $f_{\text{NL}}^{\text{local}}$ (at ~ 15% level), but stronger for $f_{\text{NL}}^{\text{equil}}$ (~ 25%) and $f_{\text{NL}}^{\text{ortho}}$ (~ 28%). Our results show that WST, both alone and in combination with P+B, can improve the extraction of information on PNG from LSS data over the one attainable by a standard P+B analysis. Moreover, we identify a class of WST in which the origin of the extra information on PNG can be cleanly isolated.

**KEYWORDS:** cosmological parameters from LSS, galaxy clustering, cosmological simulations

ARXIV EPRINT: 2403.17657

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup>Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma and INFN Gruppo Collegato di Parma, Parco Area delle Scienze 7/A, I-43124, Parma, Italy</footnote>

<footnote><sup>b</sup>Institut d'Astrophysique Spatiale, Université Paris-Saclay, CNRS, Bâtiment 121, F-91405, Orsay, France</footnote>

<footnote><sup>c</sup>Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova and INFN, Sezione di Padova, Via Marzolo 8, I-35131, Padova, Italy</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Summary statistics</td><td>3</td></tr><tr><td></td><td>2.1 Power spectrum and bispectrum</td><td>3</td></tr><tr><td></td><td>2.2 Wavelet scattering transform</td><td>4</td></tr><tr><td>3</td><td>Simulation-based Fisher matrices</td><td>7</td></tr><tr><td>4</td><td>Simulations</td><td>8</td></tr><tr><td>5</td><td>Results</td><td>9</td></tr><tr><td></td><td>5.1 First order scattering coefficients and P(k)</td><td>9</td></tr><tr><td></td><td>5.2 Matter at z = 1</td><td>10</td></tr><tr><td></td><td>5.3 Halos at z = 0</td><td>10</td></tr><tr><td>6</td><td>Discussion</td><td>12</td></tr><tr><td>7</td><td>Conclusions</td><td>18</td></tr><tr><td>A</td><td>Comparison with marked statistics</td><td>20</td></tr><tr><td>B</td><td>Mother wavelets in Kymatio</td><td>20</td></tr><tr><td>C</td><td>Tests of convergence, saturation and gaussianity of the WST</td><td>20</td></tr></table>

## 1 Introduction

An important goal of observational cosmology in the coming years will be the development of methodologies to extract maximum information from non-linear scales in galaxy surveys. Clearly, this will involve going beyond power spectrum analysis, since the power spectrum alone cannot capture phase information and is therefore not sensitive to structures such as, for example, filaments in the cosmic web. A natural choice in this respect is to consider the bispectrum (three-point function of the Fourier modes of the density field), as this is the lowest order non-Gaussian (NG) correlator. Indeed, many joint power spectrum and bispectrum (P+B) Large Scale Structure (LSS) forecasts and analyses are already present in the cosmological literature [1–5] and this kind of approach to the analysis of future galaxy survey data presents some advantages, notably the availability of several well established and validated P+B analysis pipelines and the theoretical interpretability of bispectrum measurements — at mildly non-linear scales — in cosmological perturbation theory. On the other hand, the bispectrum clearly does not exhaust the available information at fully non-linear scales and a straightforward extension to higher order correlators becomes increasingly more cumbersome and difficult to implement, especially when considering the
need for accurate evaluation of covariance matrices in the non-linear regime of interest and for inclusion of realistic observational effects (e.g., window functions). This motivates the search for different NG summary statistics. The ideal properties that we would require from them would be ease of implementation, efficiency at extracting all the available NG information using a limited number of modes and clear physical interpretability. With such goals in mind, several NG summaries have already been considered in LSS analysis. Among those, let us mention for example marked spectra [6–10], power spectra from cosmic web environments [11, 12], skew spectra [13–16], Minkowski functionals [17, 18], k-nearest neighbour cumulative distributions functions [19, 20], probability distribution function (PDF) of late-time cosmic density fluctuations [21, 22], halo mass function probes [23–28], and the void abundance [29, 30].

Here, we focus on another statistic that has received significant recent attention in observational cosmology, namely the Wavelet Scattering Transform (WST) [31, 32]. The WST has the structure of a three-layer convolutional network, in which the operations performed at each stage are, in sequence, non-linear activation via the modulus function, convolution by a wavelet filter and low-pass spatial averaging of the data. The WST has recently been applied to simulations in [33–37], showing that it could efficiently capture NG information and produce significant improvements over power-spectrum-only cosmological constraints, especially on neutrino mass. The main goal of this work is to investigate the WST as a tool for the study of primordial non-Gaussianity (PNG) and to compare in detail its performance with that achievable via a P+B analysis, along with providing a general physical interpretation for this comparison. Non-Gaussian features in the primordial cosmological density field are a distinctive prediction of inflationary models beyond standard single-field slow-roll. Primordial non-Gaussian signatures are encoded in the primordial bispectrum of cosmological perturbations: different deviations from the single-field slow-roll scenario produce specific functional dependencies of the bispectrum on triangular configurations of Fourier modes. This model dependence makes PNG a very powerful tool to constrain inflation. At the same time, the expected PNG signal is very small and strongly degenerate, at non-linear scales, with the bispectrum produced at late times by gravitational evolution of cosmic structures and galaxy bias. Beyond-bispectrum statistics, such as the WST, can help breaking this degeneracy, which is one of the main motivations of the present analysis. Let us note here that, instead of relying on summary statistics, an alternative, more ambitious approach consists in performing forward modeling field level inference. Field level analysis is more and more employed in cosmological data analysis and carries the promise of extracting all the available information in the data (see, e.g., [38–51]). The more traditional approach, based on preliminary compression of the data in a set of summaries, albeit lossy, remains nonetheless a valid alternative, often leading to easier physical interpretation of the results and comparison with analytical models; it is also relevant as it gives insight on where (i.e., in which domain, at which scales, and so on) the bulk of information about the relevant signals lies and on how it can be efficiently recovered.

In our work, we consider the three main primordial bispectrum shapes — local, equi-lateral and orthogonal — and produce Fisher forecast for the corresponding $f_{\text{NL}}$ amplitude parameters, jointly with standard cosmological parameters ($\sigma_8, \Omega_m, n_s, h$).
The main aim of this kind of study is to check whether future LSS observations, in combination with suitable estimators, will be able to improve over current primordial NG constraints from Planck CMB data [52], namely $f_{\text{NL}}^{\text{local}} = -0.9 \pm 5.1$; $f_{\text{NL}}^{\text{equil}} = -26 \pm 47$; and $f_{\text{NL}}^{\text{ortho}} = -38 \pm 24$ (68 % CL, statistical). These are the best current cosmological bounds on $f_{\text{NL}}$. For galaxy clustering, the P+B analysis of BOSS data performed in [1] yields $f_{\text{NL}}^{\text{local}} = 7 \pm 31$; $f_{\text{NL}}^{\text{equil}} = 245 \pm 293$; and $f_{\text{NL}}^{\text{ortho}} = -60 \pm 72$ (68 % CL, statistical). While LSS constraints are at present not competitive with CMB, large improvements are expected in the near future from surveys such as, e.g., *Euclid*.

To build the Fisher matrix up to a scale $k_{\text{max}} = 0.5 h \text{ Mpc}^{-1}$, we compute covariance matrices and summary responses to changes in parameters by means of Monte Carlo averaging and numerical differentiation, using the publicly available QUIJOTE and QUIJOTE-PNG N-body simulation suites. This a set of simulations with 512³ dark matter particles in boxes of size $1 h^{-1}\text{Gpc}$.

The paper is structured as follows. In section 2 we introduce the summary statistics considered in our analysis: power spectrum, bispectrum and WST; in section 3 we describe our methodology to compute Fisher matrices; in section 4 we describe in more detail the input QUIJOTE simulation dataset; in section 5 we describe our analysis and present the main results; in section 6 we discuss the interpretation of our results and draw a comparison between the P+B and WST performance; in section 7 we draw our main conclusions. Finally, in appendix A we briefly discuss the comparison between WST and marked statistics, whereas in appendix B we give the explicit form of the mother wavelet for the transform we consider and in appendix C we check numerical convergence and Gaussianity of the WST coefficients used in the analysis.

# 2 Summary statistics

## 2.1 Power spectrum and bispectrum

To extract the cosmological information of the density contrast field $\delta(\vec{x})$, we consider its two and three-point correlators in Fourier space, the power spectrum and bispectrum.

From a Fourier-transformed density field $\delta(\vec{k})$ defined on a three-dimensional grid, the power spectrum $P(k)$ is easily obtained using a binned estimator (see e.g., [53]):

$$ \hat{P}(k_i) = \frac{1}{VN_i} \sum_{k \in \Delta_i} \delta(\vec{k}) \delta^*(\vec{k}), \qquad (2.1) $$

where $V$ is the surveyed volume, $\Delta_i$ defines the binning of the relevant $k$-range, $N_i$ is the count of modes $k$ in the $i$-th bin and $k_i$ the average modulus.

Rather than computing the binned bispectrum $B(k_1, k_2, k_3)$, we use a modal approach based on its expansion on a basis of separable functions (i.e. symmetrized product of one-dimensional functions in $k$) [54–57]. Then, evaluating the information content of the bispectrum only comes to estimating modal modes of the following form:

$$ \hat{\beta}_n = \frac{1}{V} \int d^3 x M_p(\vec{x}) M_q(\vec{x}) M_r(\vec{x}), \quad \text{with } M_p(\mathbf{x}) \equiv \int \frac{d^3 k}{(2\pi)^3} \frac{q_p(k) \delta(\vec{k})}{\sqrt{k P(k)}} e^{i\vec{k} \cdot \vec{x}}, \qquad (2.2) $$
where the $q_p(k)$ are the chosen one-dimensional functions for the bispectrum expansion. As was shown in [58–60], only a small number of such functions, constructed from simple polynomials (e.g. Legendre polynomials) or from the tree-level matter bispectrum, leading to typically only 100 modal modes $\beta_n$, is necessary to study the matter and halo bispectrum from linear to non-linear scales ($k_{max} = 0.5 h \text{ Mpc}^{-1}$).

## 2.2 Wavelet scattering transform

The observables in the WST are the scattering coefficients, which are built from non-linear operations on $\delta(\vec{x})$, and are organized in different orders, or layers. In this paper, we adopt the implementation of the WST encoded in the KYMATIO<fnref n="1" /> package [32]. The zeroth order coefficient is just the spatial average of the modulus of $\delta(\vec{x})$ raised to some power $q$,

$$S_0^q = \langle |\delta(\vec{x})|^q \rangle. \qquad (2.3)$$

To obtain the higher order coefficients a mother wavelet function is defined, which, for KYMATIO, is

$$\psi_{0l}^m(\vec{x}, \sigma) = \frac{C_l}{(2\pi\sigma^2)^{3/2}} \exp\left[-\frac{x^2}{2\sigma^2}\right] R_l^m\left(\frac{\vec{x}}{\sigma}\right), \qquad (2.4)$$

where $x = |\vec{x}|$ and the $R_l^m(\vec{r})$ are the regular solid harmonics,

$$R_l^m(\vec{r}) = \sqrt{\frac{4\pi}{2l+1}} r^l Y_l^m(\hat{r}), \qquad (2.5)$$

with, $\hat{r} = \vec{r}/r$, and $Y_l^m(\hat{r})$ are the Laplacian Spherical Harmonics. The coefficients $C_l$ are given in appendix B. From the mother wavelet, a family of rescaled ones, $\psi_{jl}^m$, with $j = 1, 2, \dots$, is generated by successive doublings of the length scale $\sigma$,

$$\psi_{jl}^m(\vec{x}, \sigma) = \psi_{0l}^m(\vec{x}, 2^j \sigma) = 2^{-3j} \psi_{0l}^m(2^{-j} \vec{x}, \sigma). \qquad (2.6)$$

Following [33], in this analysis based on simulation data we will set

$$\sigma = 0.8 \frac{L}{N^{1/3}}, \qquad (2.7)$$

where $L$ is the side of the cubic box and $N$ the number of grid points. So, the radial dependence is a gaussian of width 0.8 times the grid step for $j = 0, 1, 2, \dots$ times the grid step for $j = 1, 2, \dots$, and so on, multiplied by the power law $r^l$. The parameter $l$, together with $m$, controls also the angular dependence, via the spherical harmonics. The 0.8 coefficient in eq. (2.7) was empirically determined in ref. [33] as the optimal choice in terms of the final precision on the parameters determination.

In Fourier space, the wavelets read,

$$\tilde{\psi}_{jl}^m(\vec{p}, \sigma) = \int d^3 r e^{i\vec{p} \cdot \vec{r}} \psi_{jl}^m(\vec{x}, \sigma) = i^l C_l \exp\left[-\frac{(2^j \sigma p)^2}{2}\right] R_l^m(2^j \sigma \vec{p}), \qquad (2.8)$$

with $p = |\vec{p}|$.

<footnote n="1">¹https://www.kymat.io.</footnote>
First order coefficients are then defined as

$$S_1^q(j, l) = \left\langle \left( \sum_{m=-l}^{l} |\delta(\vec{x}) * \psi_{jl}^m(\vec{x}, \sigma)|^2 \right)^{\frac{q}{2}} \right\rangle, \quad (2.9)$$

where the symbol '*' indicates convolution in configuration space. Due to the sum over $m$, the coefficients depend, for fixed $q$, on $j$ and $l$ only. In the following, we will also show results for the special case $q = 2$, where the interpretation of these coefficients simplifies considerably. Indeed, in this case, the first order coefficients are directly related to the non-linear power spectrum $P(p)$,

$$ \begin{aligned} S_1^2(j, l) &= \int \frac{d^3 p}{(2\pi)^3} \frac{d^3 p'}{(2\pi)^3} e^{-i\vec{x} \cdot (\vec{p} + \vec{p}')} \langle \tilde{\delta}(\vec{p}) \tilde{\delta}(\vec{p}') \rangle \psi_{jl}^m(\vec{p}) \psi_{jl}^{\tilde{m}}(\vec{p}') , \\ &= C_l^2 \int \frac{d^3 p}{(2\pi)^3} W_l^2[2^j \sigma p] P(p) \\ &= C_l^2 \langle \delta^2(\vec{x}; j, l) \rangle, \end{aligned} \quad (2.10) $$

where we have used the property of spherical harmonics,

$$\frac{4\pi}{2l+1} \sum_{m=-l}^{l} Y_l^m(\hat{p}) Y_l^m(-\hat{p}) = 1, \quad (2.11)$$

and defined the Fourier space filter,

$$W_l[z] = z^l \exp(-z^2/2), \quad (2.12)$$

showing that, for fixed variance of the Gaussian filter (set by $j$), higher values of $l$ probe larger wavenumbers. The last line of eq. (2.10) is the space averaged variance of the inverse Fourier transform of the filtered density field $\tilde{\delta}(\vec{p}; j, l) \equiv W_l[2^j \sigma p] \tilde{\delta}(\vec{p})$.

To define second order coefficients, we first define a new non-linear field,

$$U(j_1, l)(\vec{x}) = \left( \sum_{m=-l}^{l} |\delta(\vec{x}) * \psi_{j_1 l}^m(\vec{x}, \sigma)|^2 \right)^{\frac{1}{2}}, \quad (2.13)$$

and then convolve it again with the wavelet at scale $j_2$, using the same operation as (2.9),

$$S_2^q(j_2, j_1, l) = \left\langle \left( \sum_{m=-l}^{l} |U(j_1, l)(\vec{x}) * \psi_{j_2 l}^m(\vec{x}, \sigma)|^2 \right)^{\frac{q}{2}} \right\rangle, \quad (2.14)$$

where, again following [33], we make the choice of using the same value for $l$ both for the internal ($j_1$) and for the external ($j_2$) convolution.

The $q = 2$ case now gives

$$ \begin{aligned} S_2^2(j_2, j_1, l) &= C_l^2 \int \frac{d^3 p}{(2\pi)^3} W_l[2^{j_2} \sigma p]^2 P_U(p; j_1, l), \\ &= C_l^2 \langle U^2(\vec{x}; j_2, j_1, l) \rangle, \end{aligned} \quad (2.15) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Radial profile of the wavelets at $l = 0$ (solid lines) and $l = 4$ (dashed lines) at $m = 0$ and for different values of $j$, where $j = 0$ is equivalent to the mother wavelet. As discussed in the text, we cutoff the field beyond $k_{max} = 0.5 h \, \text{Mpc}^{-1}$, this is highlighted by the shaded region.</figcaption></figure>

where the power spectrum in the integral now refers to the field $U(j_1, l)(\vec{x})$,

$$ \langle \tilde{U}(j_1, l)(\vec{p}) \tilde{U}(j_1, l)(\vec{p}') \rangle \equiv (2\pi)^3 \delta_D(\vec{p} + \vec{p}') P_U(p; j_1, l), \quad (2.16) $$

and the filtered field $U(\vec{x}; j_2, j_1, l)$ is the inverse Fourier transform of $W_l[2^{j_2} \sigma p] \tilde{U}(j_1, l)(\vec{p})$.

The $q = 2$ case turns out to be particularly useful in order to compare the performance of the WST with that of the PS, as all the extra information beyond the non-linear PS is encoded in scattering coefficients beyond first order. The second order coefficients are just the variances of the filtered fields $U(\vec{x} : j_2, j_1, l)$, which therefore contain information on the bispectrum and higher order correlators.

Fixing $\sigma$ as in eq. (2.7), we must decide how many coefficients to compute at first and second order. For the first order, we take $j = 0, 1, \dots, J$, with $J = 4$, corresponding to Gaussian smoothing filters in real space of widths from $\sigma$ to $2^4 \sigma$. In figure 1, we show the Fourier space filters use to analyze the QUIJOTE suite of simulations presented in section 4. As for the angular dependence of the wavelets, we take $l = 0, 1, \dots, L$, with $L = 4$. When considering second order scattering coefficients, we take the external convolution, controlled by $j_2$, at scales larger than the internal one, therefore $j_1 < j_2 = 1 \dots, J$, with $J = 4$ as for the first order coefficients. This gives a total of 76 coefficients, one of zeroth order, 25 of first order and 50 of second order. As we do not consider scattering coefficients at higher orders, these will form our data vector.

Since the most computationally heavy operation for each WST coefficient is a number of convolutions corresponding its order, numerical evaluation of the full set of 76 coefficients for one realization of $\delta(\vec{x})$ requires $O(1)$ seconds once the operation is cast in Fourier space. This level of performance has been observed on a high-end GPU through the PYTORCH<fnref n="2" /> back-end of KYMATIO.

<footnote n="2">²https://pytorch.org.</footnote>
### 3 Simulation-based Fisher matrices

To benchmark the amount of information on PNG amplitudes, and more generally on cosmological parameters, extracted from the non-linear density field using the different summary statistics presented in section 2, a standard method is to use the Fisher matrix formalism.

Under the reasonable assumption of normally distributed summary statistics $\mathbf{s}$, depending on a set of parameters $\theta$ through their mean $\bar{\mathbf{s}}$, the Fisher matrix is given by:

$$F_{ij} = \left( \frac{\partial \bar{\mathbf{s}}}{\partial \theta_i} \right)^{\mathrm{T}} \mathbf{C}^{-1} \left( \frac{\partial \bar{\mathbf{s}}}{\partial \theta_j} \right), \qquad (3.1)$$

where $\mathbf{C}$ is the covariance matrix. For each parameter $\theta_i$, the quantity $\sqrt{(F^{-1})_{ii}}$ is the Cramer-Rao bound, the expected error bar from an optimal estimator.

To evaluate eq. (3.1), we choose the simulation-based approach (see section 4 for a description of the simulations used in this work). Unbiased estimates of inverse covariance matrices are obtained using the Hartlap/Anderson correction factor [61]:

$$\mathbf{C}^{-1} = \frac{n_r - n_s - 2}{n_r - 1} \hat{\mathbf{C}}^{-1}, \quad \text{with } \hat{\mathbf{C}} = \frac{1}{n_r - 1} (\mathbf{s} - \bar{\mathbf{s}})(\mathbf{s} - \bar{\mathbf{s}})^{\mathrm{T}}, \qquad (3.2)$$

where $n_r$ is the number of realizations used and $n_s$ the number of summary statistics. Derivatives can be computed at a chosen fiducial point (denoted by the superscript *) by central finite difference:

$$\frac{\partial \bar{\mathbf{s}}}{\partial \theta_i} \Big|_{\theta^*} = \frac{\bar{\mathbf{s}}(\theta_i^* + \delta\theta_i) - \bar{\mathbf{s}}(\theta_i^* - \delta\theta_i)}{2\delta\theta_i}. \qquad (3.3)$$

However, a direct application of this method requires a very large number of simulations to obtain numerically converged results, or can otherwise lead to spuriously super-optimistic Fisher error bars (for a discussion of this issue in the context of PNG, see [62, 63]). This issue can be mitigated using the method of [64], consisting in computing the Fisher matrix from score-compressed statistics instead of the raw summary statistics. This score-compression is optimal, in the sense it keeps the full information about parameters, under the same conditions as the ones used to write eq. (3.1) for the Fisher matrix [65, 66]. From $n_s$ summary statistics, $n_p$ score-compressed statistics $\tilde{\mathbf{s}}$, where $n_p$ is the number of parameters, are obtained using:

$$\tilde{s}_i = \frac{\partial \bar{s}}{\partial \theta_i} \Big|_{\theta^*} \mathbf{C}^{-1} (\mathbf{s} - \bar{\mathbf{s}}^*), \qquad (3.4)$$

which are the same ingredients as for the Fisher matrix computation, and thus can be evaluated from the same sets of simulations with the advantage of a much faster numerical convergence. If numerical convergence is not reached, the Fisher bounds computed from the score-compressed statistics tend to be sub-optimal. As discussed in [64], this issue can be significantly mitigated by computing the geometric mean of the Fisher error bars from the compressed and raw summary statistics, which we do for the results reported in this work.
<table><caption>**Table 1.** The parameters of the QUIJOTE and QUIJOTE-PNG N-body simulations and halo catalogues used in this work.</caption><thead><tr><th></th><th>N<sub>sims</sub></th><th>σ<sub>8</sub></th><th>Ω<sub>m</sub></th><th>Ω<sub>b</sub></th><th>n<sub>s</sub></th><th>h</th><th>f<sub>NL</sub><sup>local</sup></th><th>f<sub>NL</sub><sup>equil</sup></th><th>f<sub>NL</sub><sup>ortho</sup></th><th>M<sub>min</sub>(M<sub>☉</sub>/h)</th></tr></thead><tbody><tr><td>Fiducial</td><td>15000</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>σ<sub>8</sub><sup>+</sup></td><td>500</td><td>0.849</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>σ<sub>8</sub><sup>-</sup></td><td>500</td><td>0.819</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>Ω<sub>m</sub><sup>+</sup></td><td>500</td><td>0.834</td><td>0.3275</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>Ω<sub>m</sub><sup>-</sup></td><td>500</td><td>0.834</td><td>0.3075</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>n<sub>s</sub><sup>+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9824</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>n<sub>s</sub><sup>-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9424</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>h<sup>+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6911</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>h<sup>-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6511</td><td>0</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>local,+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>+100</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>local,-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>-100</td><td>0</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>equil,+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>+100</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>equil,-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>-100</td><td>0</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>ortho,+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>+100</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>f<sub>NL</sub><sup>ortho,-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>-100</td><td>3.2 × 10<sup>13</sup></td></tr><tr><td>M<sub>min</sub><sup>+</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.3 × 10<sup>13</sup></td></tr><tr><td>M<sub>min</sub><sup>-</sup></td><td>500</td><td>0.834</td><td>0.3175</td><td>0.049</td><td>0.9624</td><td>0.6711</td><td>0</td><td>0</td><td>0</td><td>3.1 × 10<sup>13</sup></td></tr></tbody></table>

## 4 Simulations

Evaluating with precision the statistical properties of the different observables introduced in section 2 up to non-linear scales, as well as their dependence on cosmological parameters, requires a large number of simulations. In this work, we use the QUIJOTE<fnref n="3" /> and QUIJOTE-PNG<fnref n="4" /> suites of N-body simulations [67, 68] which were designed exactly for this purpose.

These are N-body simulations, containing each 512³ dark matter particles in boxes of size 1 Gpc h⁻¹, in which the gravitational evolution of structures up to z = 0 is computed with the TreePM code GADGET-III [69] starting from initial conditions determined at z = 127 using 2LPTIC [70] or 2LPTPNG [68, 71].<fnref n="5" /> The suite contains 15 000 simulations with Planck-like cosmological parameters [72], and subsets of 500 were one of the parameters is slightly shifted above or below its fiducial value. The characteristics of these different sets are recalled in table 1.

In this work, we focus on the parameters $\vec{\theta} = \{\sigma_8, \Omega_m, n_s, h, f_{\text{NL}}^{\text{local}}, f_{\text{NL}}^{\text{equil}}, f_{\text{NL}}^{\text{ortho}}\}$, and study both the matter (at z = 1) and halo (at z = 0) density fields. The definitions of the

<footnote n="3">³https://quijote-simulations.readthedocs.io/en/latest/.</footnote>

<footnote n="4">⁴https://quijote-simulations.readthedocs.io/en/latest/png.html.</footnote>

<footnote n="5">⁵https://github.com/dsjamieson/2LPTPNG.</footnote>
PNG parameters, $f_{\text{NL}}^{\text{local}}$, $f_{\text{NL}}^{\text{equil}}$, and $f_{\text{NL}}^{\text{ortho}}$, characterizing different shapes of the primordial bispectrum, are given in [68]. The halos have been identified using friends-of-friends (FOF) [73]. Note that we only include halos with a mass above $M_{\text{min}} = 3.2 \times 10^{13} M_{\odot}/h$, and include $M_{\text{min}}$ as a nuisance parameter by varying it in the halo analyses. As shown in previous analyses [63, 68], $M_{\text{min}}$ can be thought of as a proxy for the linear bias parameter $b_1$.

We estimate the WST coefficients from density fields computed by interpolating the dark matter particle or halo positions on regular three-dimensional grids of size $N_{\text{grid}} = 256$ using the CIC scheme implemented in the code PYLIANS3.<fnref n="6" /> We apply a sharp, low-pass filter on the input field $\delta(\vec{x})$ that excludes all modes above $k_{\text{max}} = 0.5 h \, \text{Mpc}^{-1}$ (that is, the smallest accessible scale with QUIJOTE simulations), and then compute all 76 scattering coefficients.

For the power spectrum, we consider 79 linearly spaced bins, from the fundamental frequency, in our case $k_{\text{min}} = 2\pi/10^3 h \, \text{Mpc}^{-1}$, up to $k_{\text{max}} = 0.5 h \, \text{Mpc}^{-1}$. The bispectrum modal modes are also studied on the same k-range, using a basis constituted of the tree-level matter bispectrum and polynomial functions up to degree 10 (for a total of 73 bispectrum modes). Both summary statistics were previously computed in [60] and [63] for the matter and halo fields, respectively, where these choices of binning and modal basis were shown to give optimal and numerically stable Fisher error bars.

Note that the matter density field is studied in real space, while the halo density field is transformed to redshift-space applying the velocity correction along the *x* direction.

In order to estimate the Fisher matrix for each statistic as detailed in section 3, we make use of the 15 000 simulations at fiducial cosmology to estimate the covariance matrix, while the sets of 500 simulations are used to estimate the central derivatives.

# 5 Results

What follows are the result of the Fisher matrix analysis presented in section 3 applied to the summary statistics introduced in section 2. For the WST we present results for $q = 2$, which, as we have discussed in section 2.2, provides a point of comparison with the more familiar $P(k)$ at the level of first order coefficients. For $q \neq 2$ one expects that some information from higher order correlators should leak into first order scattering coefficients. Moreover, values of $q < 1$ emphasize underdensities in the field, where most of the cosmological information is still recoverable. Therefore, we also show results for $q = 0.8$, a value chosen to allow a direct comparison with the results of [33]. We also choose to stick to the same values of *J* and *L*, as we find that increasing them does not lead to significant improvements in the bounds.

## 5.1 First order scattering coefficients and $P(k)$

As we have shown in eq. (2.10), first order scattering coefficients with the specific choice $q = 2$, that is, $S_1^2(j, l)$, contain no more information than the (non-linear) power spectrum $P(k)$. To check this explicitly, we run a Fisher matrix analysis on the non-linear dark matter field at redshift $z = 1$. The resulting constraints for the set of cosmological parameters $\vec{\theta} = \{\sigma_8, \Omega_m, n_s, h\}$ are shown in figure 2. The near to perfect overlap between the 2D contours for $P(k)$ and $S_1^2(j, l)$ shows that the information content of the two statistics is

<footnote n="6">⁶https://github.com/franciscovillaescusa/Pylians3.</footnote>
nearly identical. $P(k)$ gives slightly better constraints, which shows that the equivalence between the 25 integrals of $P(k)$ of eq. (2.10) and the 78 $P(k)$ bins is not perfect. This can be probably improved by increasing the maximum $j$ and/or $l$. However, to limit the computational cost, we decide to stick to our choice since, as we will see, including second order coefficients outperforms the power spectrum for all parameters, including PNG ones, both for matter and for halos.

As expected, the contours for the $S_1^{0.8}(j, l)$ coefficients differ from the ones for the $S_1^2(j, l)$ and for the $P(k)$. Moreover, for all the parameters but $n_s$, the performance is actually worse, suggesting that the redistribution of the information in passing from $q=2$ to $q=0.8$ moves it from first order to higher order coefficients.

## 5.2 Matter at $z=1$

We now consider the matter field at $z=1$ and extend the analysis of the previous subsection by enlarging the parameter vector to include the PNG parameters, and by including second order scattering coefficients and the bispectrum. Specifically, in figures 3 and 4 we use the parameter vector $\vec{\theta} = \{\sigma_8, \Omega_m, n_s, h, f_{\text{NL}}^{\text{local}}, f_{\text{NL}}^{\text{equil}}, f_{\text{NL}}^{\text{ortho}}\}$, whereas in the upper table of figure 5, each PNG shape is analyzed jointly with cosmological parameters and independently from the two others PNG shapes. For cosmological parameters, we report the largest bound of the three analyses (differences being small anyway).

In figures 3 and 4 we show the constraints from the combination of all the 76 WST coefficients computed with $q=0.8$ and $q=2$, respectively. We compare them with those obtained by the combination of the power spectrum and the bispectrum ($P+B$), and the combination of all the statistics considered, $P+B+WST$. The corresponding $1\sigma$ Cramer-Rao bounds are listed in figure 5. We do not show results for $P$ alone, as they are outperformed by both versions of WST on all parameters.

As discussed in the previous subsection, for $q=2$ the extra information comes entirely from the second order scattering coefficients. For $q=0.8$, the inclusion of second order coefficients improves the performance, which is now overall comparable with that of $q=2$: it is worse for all cosmological parameters except $n_s$, but slightly better for the PNG parameters. However, when combined with $P+B$, WST at $q=0.8$ is slightly more efficient than $q=2$ at removing degeneracies between parameters.

## 5.3 Halos at $z=0$

We then apply the Fisher matrix analysis on the non-linear dark matter halo field at $z=0$. In order to test the performance of the WST on a biased tracer in a more realistic scenario, the field is transformed in redshift space. The set of cosmological parameters analyzed in figures 6 and 7 is $\vec{\theta} = \{\sigma_8, \Omega_m, n_s, h, M_{\text{min}}, f_{\text{NL}}^{\text{local}}, f_{\text{NL}}^{\text{equil}}, f_{\text{NL}}^{\text{ortho}}\}$. We show the bounds given by all the 76 WST coefficients computed with $q=0.8$ and $q=2$ respectively, compared with those given by the first two power spectrum multipoles ($P_0+P_2$), their combination with the bispectrum ($P_0+P_2+B$), and the combination of all the statistics. The $1\sigma$ Cramer-Rao bounds for all the above statistics are listed in the lower table of figure 5, where each PNG shape is analyzed independently from the other two.

<footnote n="7">⁷We have verified that the hexadecapole does not improve the constraints.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Expected distribution of the parameters around the fiducial QUIJOTE-PNG cosmology, with covariance given by the Fisher matrix computed on the non-linear matter field at $z = 1$. The statistics compared in the figure are 1st order WST coefficients $S_1^{0.8}(j, l)$ (olive), $S_1^2(j, l)$ (grey) and $P(k)$ (aquamarine).</figcaption></figure>

As for matter, the WST outperforms the power spectrum alone for both values of $q$. This is true in particular for the PNG parameters, including for $f_{\text{NL}}^{\text{local}}$ (we do not report the constraints on the other two PNG shapes from $P_0 + P_2$, in figure 5, as they are huge). $P_0 + P_2 + B$ performances are comparable with WST for the cosmological parameters (and for $M_{\text{min}}$), while for the PNG parameters we find a mixed situation: $f_{\text{NL}}^{\text{local}}$ constraints are better (by about 10-15%) for $P_0 + P_2 + B$, $f_{\text{NL}}^{\text{equil}}$ are better (by about 20-30%) for WST, while $f_{\text{NL}}^{\text{ortho}}$ ones are of the same order for WST with $q = 2$ and about 30% better for $q = 0.8$.

Overall, our results on halos confirm what we found for matter, namely, that $q = 0.8$ and $q = 2$ WST give comparable constraints basically on all parameters. This means that
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Expected distribution of the parameters around the fiducial QUIJOTE-PNG cosmology, with covariance given by the Fisher matrix computed on the non-linear matter field at $z = 1$. The statistics compared in the figure are $P+B$ (olive), WST coefficients $S_n^{0.8}(j, l)$ with $n = 0, 1, 2$ (grey), and the full combination $P+B+WST$ (aquamarine).</figcaption></figure>

taking different exponents does not change significantly the amount of extracted information, once second order coefficients are included.

In appendix A, we compare directly the information content on PNG and cosmological parameters of WST with marked statistics, which is very similar.

# 6 Discussion

In order to understand where the information on PNG parameters is stored in the WST, we plot in figures 8 and 9 the sensitivity of WST coefficients to each parameter. More
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Expected distribution of the parameters around the fiducial QUIJOTE-PNG cosmology, with covariance given by the Fisher matrix computed on the non-linear matter field at $z = 1$. The statistics compared in the figure are $P+B$ (olive), WST coefficients $S_n^2(j, l)$ with $n = 0, 1, 2$ (grey), and the full combination $P+B+WST$ (aquamarine).</figcaption></figure>

specifically, for each parameter $\alpha$, we plot the absolute value of the difference between the WST coefficients computed in the set of simulations with parameter $\alpha^+$ and the one for $\alpha^-$, normalized to the ones for the fiducial cosmology, (see table 1 for the specific values taken for each parameter). Figure 8 is for dark matter at $z = 1$ and real space, and figure 9 for halos at $z = 0$ and redshift space. The upper row of each panels contains PNG parameters, while the lower row all the other ones. Moreover, the left and right columns are obtained by setting the exponent parameter, $q$, to $q = 0.8$ and $q = 2$, respectively.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The Cramer-Rao $1\sigma$ bounds for different combinations of summary statistics (power spectrum, bispectrum and WST, “all” referring to the three of them altogether). Every statistic is measured up to $k_{max} = 0.5\ h\ Mpc^{-1}$ in the QUIJOTE N-body simulations at $z=1$ (top panel), or in the corresponding halo catalogues at $z=0$ (bottom panel), for a volume of $1\ (h^{-1}Gpc)^3$. The colour scale indicates the ratio of each error bar with respect to its power spectrum + bispectrum equivalent, highlighting the significant additional information picked by the WST in most cases. Note that each PNG shape is analyzed jointly with cosmological parameters and independently from the two others PNG shapes. For cosmological parameters, we report the largest bound of the three analyses (differences being small anyway).</figcaption></figure>

The WST coefficients from 0 to 75 are ordered in decreasing length scale of the most external convolution. So, the first WST coefficient on the left is the second order one with $j_2 = 4$, $j_1 = 3$, $l = 0$, then $j_2 = 4$, $j_1 = 3$, $l = 1$, and so on up to $j_2 = 4$, $j_1 = 0$, $l = 4$ (end of the first light brown horizontal segment on the right). Then the first order WST coefficients with $j = 4$, $l = 0, \dots, 4$ (dark brown segment), followed by second order coefficients with $j_1 < j_2 = 3$, and so on.

First of all, we notice no clear qualitative difference between left and right columns of each panel, meaning that the intuition developed in section 2.2 for the $q=2$ case, namely, that first order coefficients are related to the power spectrum, extends, to a large extent, also to $q=0.8$. This explains why the first order coefficients for matter are almost insensitive to PNG parameters (figure 8, upper row, dark brown segments), especially at large scales. Indeed, the $f_{NL}$ dependence of the matter power spectrum is a suppressed one-loop effect,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Expected distribution of the parameters around the fiducial QUIJOTE-PNG cosmology, with covariance given by the Fisher matrix computed on the non-linear halo field at $z = 0$. The statistics compared in the figure are the power spectrum monopole and quadrupole $P_0 + P_2$ (olive), $P_0 + P_2 + B$ (grey), WST coefficients $S_n^{0.8}(j, l)$ with $n = 0, 1, 2$ (aquamarine), and the full combination $P_0 + P_2 + B + \text{WST}$ (violet).</figcaption></figure>

and therefore it manifests itself only at small scales, as seen by the small but nonzero values for the WST coefficients at small $j$'s. The same argument explains why also second order coefficients with $l = 0$ are almost insensitive to PNG parameters, as the sum in $U(j, l)(\vec{x})$ of eq. (2.13) involves just one term and then the second order coefficient is also related to the power spectrum. We conclude that, for matter, the information on PNG parameters comes mainly from second order coefficients with $l \neq 0$. On the other hand, for the other parameters the information is stored also in the power spectrum at linear order, and therefore we expect signal also in first order WST coefficients, as seen in the second row of figure 8.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Expected distribution of the parameters around the fiducial QUIJOTE-PNG cosmology, with covariance given by the Fisher matrix computed on the non-linear halo field at $z = 0$. The statistics compared in the figure are the power spectrum monopole and quadrupole $P_0 + P_2$ (olive), $P_0 + P_2 + B$ (grey), WST coefficients $S_n^2(j, l)$ with $n = 0, 1, 2$ (aquamarine), and the full combination $P_0 + P_2 + B + \text{WST}$ (violet).</figcaption></figure>

Before concluding the analysis of the results for matter, we recall that the information content not only depends on the derivatives of the WST coefficients, but also on their covariances. In the first row of figure 10 we show the correlation coefficients defined as

$$r_{ij} \equiv \frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}}, \qquad (6.1)$$

where $C_{ij}$ is the covariance between the $i$-th and the $j$-th WST coefficient. $r_{ij}$ can take values between -1 (complete anticorrelation) and 1 (complete correlation), with $r_{ij} = 0$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Variation of WST coefficients for matter in real space with respect to different parameters, normalized to the fiducial values, for the matter field at $z=1$. Left: $q=0.8$, Right: $q=2$. Upper row: NG parameters; lower row: all other parameters. The light brown segments on the x-axis identify second order WST coefficients, whereas dark brown segments identify first order ones. There is essentially no NG signal in first order WST coefficients because the power spectrum of the matter field — to which such coefficients are sensitive — does not depend on $f_{NL}$ at lowest order. The residual signal in the dark brown segment can be attributed to higher order power spectrum corrections.</figcaption></figure>

corresponding to no correlation. The correlation patterns of $q=0.8$ and the $q=2$ cases are clearly different. For the $q=2$ case, the correlation between first order coefficients is very small at large scales and increases at smaller ones, reflecting the behaviour of the power spectrum bins. First order coefficients are also correlated with the second order ones for $l=0$, confirming our interpretation that the latter are also related to the power spectrum. The most striking features are, perhaps, the strong correlation among second order coefficients with $l \neq 0$, and the very small correlation between second order coefficients (with the exception of $l=0$ ones) and first order ones, especially at large scales. This separation is absent for $q=0.8$, confirming our conclusion (see section 5.1) that the PS information is confined to 1st order coefficients for $q=2$ and distributed among 1st and higher order ones for $q \neq 2$. Moreover, the strong correlation between second order coefficients for $q=2$ suggests that a smaller set of them may be enough to extract the same amount of information.

When we turn to halos, the differences between the PNG sensitivities of first and second order coefficients, and between the correlation patterns for $q=2$ and $q=0.8$, basically disappear for the local shape. In this case, first order WST coefficients are sensitive to $f_{NL}$ since, as is well known, the power spectrum exhibits scale dependent bias at large scales [74]. This can be understood as a new contribution to the density field

$$ \delta_h = b_1 \delta_m + b_\phi f_{\text{NL}}^{\text{local}} \phi + \epsilon, \qquad (6.2) $$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Variation of WST coefficients for the halo field in redshift space with respect to different parameters, normalized to the fiducial values. Now, first order WST coefficients pick up significant primordial NG signal for the local shape, due to the scale dependent bias signature in the power spectrum on large scales. See the main text for more explanations. Left: $q = 0.8$, Right: $q = 2$.</figcaption></figure>

where $b_1$ and $b_\phi$ are tracer-dependent bias parameters (for matter, $b_1 = 1$ and $b_\phi = 0$) and $\phi$ is the gravitational potential. Finally, a stochastic contribution, $\epsilon$, uncorrelated to $\delta_m$ and $\phi$, describes ‘shot noise’ in the halo field, which, analogously to what happens for $P$ and $P+B$, degrades the cosmological bounds from halos with respect to matter, as can be seen by comparing the upper and lower tables in figure 5. Shot noise is also responsible for the decorrelation between WST coefficients in the $q=2$ case seen in figure 10.

# 7 Conclusions

In this paper we have investigated the WST as a tool for the study of PNG, and compared its performance with that achievable via a P+B analysis. We analyzed the matter and the halo fields, in real and redshift space, measured from the publicly available QUIJOTE and QUIJOTE-PNG N-body simulation suite, producing Fisher forecast for the amplitude parameters, $f_{\text{NL}}^{\text{local}}$, $f_{\text{NL}}^{\text{equil}}$, and $f_{\text{NL}}^{\text{ortho}}$, along with the cosmological parameters.

In order to make a fair comparison, WST and P+B analyses were implemented on the same fields. In particular, to reduce possible spurious effects picked up from WST at small scales, a sharp cut of wavenumbers above $k_{\text{max}} = 0.5 \, \text{h Mpc}^{-1}$ was applied. Moreover, we compared WST’s with two values of the exponent parameter $q$, namely, $q = 0.8$ and $q = 2$, as the latter allows for a more clear interpretation of the results: first order scattering coefficients are equivalent to the power spectrum, see figure 2, so all information beyond that comes from second order coefficients.
We find that WST's (both for $q = 0.8$ and $q = 2$) outperform the power spectrum alone on both on PNG and cosmological parameter constraints, as well as on the 'bias' parameter $M_{\min}$. Specifically, on $f_{\text{NL}}^{\text{local}}$ for halos, the improvement is about 27%. When B is combined with P, halo constraints from WST are weaker for $f_{\text{NL}}^{\text{local}}$ (at $\sim 15\%$ level), but still stronger for $f_{\text{NL}}^{\text{equil}}$ ($\sim 25\%$) and $f_{\text{NL}}^{\text{ortho}}$ ($\sim 28\%$ for $q = 0.8$).

Our results on the standard cosmological (non-PNG) cosmological parameters are in line with previous analyses [33, 34], which also display a better performance of the WST with respect to the power spectrum alone. Concerning PNG parameters — which are the main focus of this work — we show that the WST can improve the extraction of PNG information from LSS data over that attainable by a standard P+B analysis. The level of improvement for QUIJOTE-PNG simulations is competitive with the best results we obtained in previous works, using other summary statistics, such as the marked power spectrum and bispectrum [75].

In order to apply the methodology presented in this paper to current and future surveys, the current analysis should be developed in two directions. On one hand, the level of realism of data should be increased, in order to match as closely as possible that of a realistic survey, covering volumes at least $O(10)$ times that of the QUIJOTE simulations. Modeling of galaxy bias, survey properties, and so on, should be included in our pipeline. With all the current analysis limitations in mind, it is however also worth pointing out that the authors of [36, 37] have recently performed WST analyses of BOSS data, for standard cosmological parameters, yielding results in line with Fisher forecasts; it is reasonable to assume that these analysis pipelines can be directly extended to PNG parameters, with a comparable performance, in the near future. For the time being, extrapolating our WST non-Gaussian forecasts to expectations for forthcoming next generation data, such as those by *Euclid*, is not always straightforward. If we consider, for example, shapes with considerable contributions from squeezed triangles (local, orthogonal), a simple volume scaling argument does not hold; we can however apply such argument to equilateral type bispectra, which are the main focus of this analysis. For a *Euclid*-like survey, a rough estimate of an expected $f_{\text{NL}}^{\text{equil}}$ bound from WST can be thus obtained by rescaling our QUIJOTE-PNG results by a factor 6–8, leading to possible improvements over current *Planck* measurements.

Further optimizations can be implemented, both for marked statistics and WST, in order to improve the performance of these estimators on PNG parameters. For example, in the WST analysis we can study the impact of different choices of wavelet transform, change the $q$ exponent or drop the angle-average approximation, whereas for marked statistics we can optimize the exponent and filter radius. All these possibilities will be the subject of future studies and a suitable hyperparameter optimization procedure might lead some summary statistics to outperform the others. Regardless of this, we believe that the main message is that both marked statistics and WST provide very powerful estimators to extract PNG information at non-linear scales. Both methodologies should thus be implemented in the PNG analysis of galaxy clustering data and combined for cross-validation and robustness purposes.

On a more fundamental side, advances in the theoretical interpretation of the WST are highly desirable. We have provided a first attempt in this paper, by considering the WST for $q = 2$. This choice enables a cleaner theoretical description, while carrying at the same time basically no information loss with respect to other $q$ values. It will be interesting to
investigate more general families of wavelets, in search for those best tailored for extracting information, especially on PNG parameters, from the LSS. For instance, the angular structure could be exploited to minimize the effect of redshift space distortions, and the underlying symmetries of the problem, such as extended galileian invariance [76–81], could be exploited to reduce the number of independent coefficients to take into account in the analyses.

## Acknowledgments

We thank William Coulton, Francisco Villaescusa-Navarro and Ben Wandelt for useful discussions and feedback on the manuscript. M. Liguori and M. Pietroni acknowledge support by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022 — grant 20228RMX4A. GJ acknowledges support from the ANR LOCALIZATION project, grant ANR-21-CE31-0019 / 490702358 of the French Agence Nationale de la Recherche.

## A Comparison with marked statistics

In [75], the marked power spectrum and marked bispectrum, computed from re-weighted density fields, were identified as promising observables of PNG on non-linear scales. In figure 11, we compare their reported $1\sigma$ Fisher bounds combining standard and marked power spectra and bispectra to our strongest constraints based on WST reported in figure 5. We obtain very similar constraints on PNG by including marked statistics or the WST, indicating that they probe the same information from higher order correlators beyond the bispectrum.

## B Mother wavelets in Kymatio

The $C_l$ coefficients appearing in (2.4) are given by

$$C_l = \begin{cases} \frac{1}{(l+1)!!}, & l \text{ even} \\ \frac{(2\pi)^{1/2}}{2^{(l+3)/2} (\frac{l+1}{2})!}, & l \text{ odd.} \end{cases} \qquad (B.1)$$

## C Tests of convergence, saturation and gaussianity of the WST

In order to rely on the Fisher matrix approach presented in section 3, we need first to verify the assumption of Gaussianity of the statistics used. While this has been extensively done for the power spectrum and bispectrum, the WST is still largely unexplored due to its novelty. Some works attempt to shed light on this [82], but no clear answer is given. Here we test the assumption on WST by simply performing an Anderson-Darling test on the way each WST coefficient distributes across full 15 000 simulations of the QUIJOTE suite at fiducial cosmology. All WST coefficients pass the test, indicating that the statistic is normally distributed. Figure 12 shows the distribution of on coefficients picked at random for dark matter and halos for both values of $q$. Having obtained our results from a simulation-based approach, we must check that our statistics have converged numerically. We test this by comparing the variation in the Cramer-Rao bounds obtained by using a growing number
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The correlation matrix of the WST coefficients. Upper row: matter, lower row: halos. Left: $q = 0.8$, Right: $q = 2$.</figcaption></figure>

$N_{\text{der}} < 500$ of simulations to compute the numerical derivative of eq. (3.3), with the same result given by the full simulation suite (that is, with $N_{\text{der}} = 500$). In figure 13 we report the ratio between the two bounds as a function of $N_{\text{der}}$ for the WST, highlighting the $f_{\text{NL}}$ configurations for clarity. If a line asymptotically converges to 1, then the bound for that parameter is considered to have numerically converged. We see that the bounds appear to converge rather quickly, especially for $f_{\text{NL}}$ in the case of halos. In the case of matter, the situation is reversed, however the absolute difference in the values never exceeds 3%, which implies very high numerical stability. For similar tests, showing the stability of the P + B analysis, see [60, 62, 63, 68]. It is also interesting to see how the Cramer-Rao bounds improve when adding more and more WST coefficients. This is shown in figure 14, where,
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of the 1-σ Fisher bounds on cosmological parameters and PNG amplitudes for different combinations of summary statistics (power spectrum, bispectrum, WST ($q = 0.8$), marked power spectrum and marked bispectrum), measured up to $k_{max} = 0.5 h \text{ Mpc}^{-1}$ in the QUIJOTE halo catalogues at $z = 0$. Numbers and the corresponding light gray lines indicate the scaling of the Fisher bounds, with respect to the largest ones in each case.</figcaption></figure>

as in figures 8, 9, and 10, we order the coefficients starting from $j = 4$, all the way to $j = 0$ (see section 6). The horizontal line at 1 indicates the corresponding bounds from $P + B$. Focusing on PNG parameters, we see that extra information with respect to $P + B$ comes in already at intermediate scales (starting from the 40th coefficient, approximately) and, in the case of halos it shows clear saturation at small scales, supporting the robustness of the bounds with respect to spurious small scale effects. Before concluding this appendix, let us briefly comment on our customary choice of *Planck* best-fit parameters to define our fiducial cosmology in the Fisher analysis. One might wonder here whether the covariance might be parameter dependent and the results could sensibly change for a different parameter choice. We do not expect this, for a number of reasons: first of all, as long as we are interested in reasonable fiducial cosmologies, current *Planck* bounds allow only for small variations with respect to the best-fit model; moreover, not only we tested here that the likelihood is with good approximation Gaussian, but we are also working in a high signal-to-noise regime, for which the parameter dependence in the likelihood can essentially be treated as linear; finally, the PNG parameter dependence is explicitly linear at tree level.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distribution of WST coefficients over the 15 000 simulations at fiducial cosmology. One coefficient at random is shown for matter and halos at each value of $q$. The green line shows a normalized Gaussian distribution with the same mean and variance as the coefficient. Upper row: matter, lower row: halos. Left: $q = 0.8$, Right: $q = 2$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Convergence of the Cramer-Rao $1-\sigma$ bounds on the cosmological parameters as more simulations are used to compute the derivatives of the WST coefficients used in the Fisher matrix. Upper row: matter, lower row: halos. Left: $q = 0.8$, Right: $q = 2$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Change in the Cramer-Rao $1\sigma$ bounds on the cosmological parameters obtained as we add more WST coefficients to the Fisher matrix. The horizontal lines represent the bounds given by $P+B$. Upper row: matter, lower row: halos. Left: $q = 0.8$, Right: $q = 2$.</figcaption></figure>

References

[1] G. D'Amico, M. Lewandowski, L. Senatore and P. Zhang, *Limits on primordial non-Gaussianities from BOSS galaxy-clustering data*, arXiv:2201.11518 [INSPIRE].

[2] G. D'Amico et al., *The BOSS bispectrum analysis at one loop from the effective field theory of large-scale structure*, JCAP 05 (2024) 059 [arXiv:2206.08327] [INSPIRE].

[3] O.H.E. Philcox et al., *Cosmology with the redshift-space galaxy bispectrum monopole at one-loop order*, Phys. Rev. D 106 (2022) 043530 [arXiv:2206.02800] [INSPIRE].

[4] M.M. Ivanov et al., *Cosmology with the galaxy bispectrum multipoles: optimal estimation and application to BOSS data*, Phys. Rev. D 107 (2023) 083515 [arXiv:2302.04414] [INSPIRE].
[5] SIMBIG collaboration, *Cosmological constraints from the nonlinear galaxy bispectrum*, Phys. Rev. D 109 (2024) 083534 [arXiv:2310.15243] [INSPIRE].

[6] R.K. Sheth, *The halo-model description of marked statistics*, Mon. Not. Roy. Astron. Soc. 364 (2005) 796 [astro-ph/0511772] [INSPIRE].

[7] R.K. Sheth, A.J. Connolly and R. Skibba, *Marked correlations in galaxy formation models*, astro-ph/0511773 [INSPIRE].

[8] O.H.E. Philcox, E. Massara and D.N. Spergel, *What does the marked power spectrum measure? Insights from perturbation theory*, Phys. Rev. D 102 (2020) 043516 [arXiv:2006.10055] [INSPIRE].

[9] E. Massara et al., *Using the marked power spectrum to detect the signature of neutrinos in large-scale structure*, Phys. Rev. Lett. 126 (2021) 011301 [arXiv:2001.11024] [INSPIRE].

[10] E. Massara et al., *Cosmological information in the marked power spectrum of the galaxy field*, Astrophys. J. 951 (2023) 70 [arXiv:2206.01709] [INSPIRE].

[11] T. Bonnaire, N. Aghanim, J. Kuruvilla and A. Decelle, *Cosmology with cosmic web environments — I. Real-space power spectra*, Astron. Astrophys. 661 (2022) A146 [arXiv:2112.03926] [INSPIRE].

[12] T. Bonnaire, J. Kuruvilla, N. Aghanim and A. Decelle, *Cosmology with cosmic web environments — II. Redshift-space auto and cross-power spectra*, Astron. Astrophys. 674 (2023) A150 [arXiv:2212.06838] [INSPIRE].

[13] D. Munshi and A. Heavens, *A new approach to probing primordial non-Gaussianity*, Mon. Not. Roy. Astron. Soc. 401 (2010) 2406 [arXiv:0904.4478] [INSPIRE].

[14] A. Moradinezhad Dizgah, H. Lee, M. Schmittfull and C. Dvorkin, *Capturing non-Gaussianity of the large-scale structure with weighted skew-spectra*, JCAP 04 (2020) 011 [arXiv:1911.05763] [INSPIRE].

[15] M. Schmittfull and A. Moradinezhad Dizgah, *Galaxy skew-spectra in redshift-space*, JCAP 03 (2021) 020 [arXiv:2010.14267] [INSPIRE].

[16] J. Hou et al., *Cosmological constraints from the redshift-space galaxy skew spectra*, Phys. Rev. D 109 (2024) 103528 [arXiv:2401.15074] [INSPIRE].

[17] G. Pratten and D. Munshi, *Non-Gaussianity in large scale structure and Minkowski functionals*, Mon. Not. Roy. Astron. Soc. 423 (2012) 3209 [arXiv:1108.1985] [INSPIRE].

[18] W. Liu, A. Jiang and W. Fang, *Probing massive neutrinos with the Minkowski functionals of large-scale structure*, JCAP 07 (2022) 045 [arXiv:2204.02945] [INSPIRE].

[19] A. Banerjee and T. Abel, *Nearest neighbour distributions: new statistical measures for cosmological clustering*, Mon. Not. Roy. Astron. Soc. 500 (2020) 5479 [arXiv:2007.13342] [INSPIRE].

[20] W.R. Coulton, T. Abel and A. Banerjee, *Small-scale signatures of primordial non-Gaussianity in k-nearest neighbour cumulative distribution functions*, arXiv:2309.15151 [INSPIRE].

[21] O. Friedrich et al., *Primordial non-Gaussianity without tails — how to measure $f_{NL}$ with the bulk of the density PDF*, Mon. Not. Roy. Astron. Soc. 498 (2020) 464 [arXiv:1912.06621] [INSPIRE].

[22] C. Uhlemann et al., *Fisher for complements: extracting cosmology and neutrino mass from the counts-in-cells PDF*, Mon. Not. Roy. Astron. Soc. 495 (2020) 4006 [arXiv:1911.11158] [INSPIRE].
[23] S. Matarrese, F. Lucchin and S.A. Bonometto, *A path integral approach to large scale matter distribution originated by non-Gaussian fluctuations*, *Astrophys. J. Lett.* **310** (1986) L21 [INSPIRE].

[24] S. Matarrese, L. Verde and R. Jimenez, *The abundance of high-redshift objects as a probe of non-Gaussian initial conditions*, *Astrophys. J.* **541** (2000) 10 [astro-ph/0001366] [INSPIRE].

[25] J. Robinson, E. Gawiser and J. Silk, *Constraining primordial non-Gaussianity with the abundance of high redshift clusters*, *Astrophys. J.* **532** (2000) 1 [astro-ph/9906156] [INSPIRE].

[26] M. LoVerde and K.M. Smith, *The non-Gaussian halo mass function with $f_{NL}$, $g_{NL}$ and $\tau_{NL}$*, *JCAP* **08** (2011) 003 [arXiv:1102.1439] [INSPIRE].

[27] A.E. Bayer et al., *Detecting neutrino mass by combining matter clustering, halos, and voids*, *Astrophys. J.* **919** (2021) 24 [arXiv:2102.05049] [INSPIRE].

[28] G. Jung et al., *Quijote-PNG: the information content of the halo mass function*, *Astrophys. J.* **957** (2023) 50 [arXiv:2305.10597] [INSPIRE].

[29] M. Kamionkowski, L. Verde and R. Jimenez, *The void abundance with non-Gaussian primordial perturbations*, *JCAP* **01** (2009) 010 [arXiv:0809.0506] [INSPIRE].

[30] G. D'Amico, M. Musso, J. Norena and A. Paranjape, *Excursion sets and non-Gaussian void statistics*, *Phys. Rev. D* **83** (2011) 023521 [arXiv:1011.1229] [INSPIRE].

[31] S. Mallat, *Group invariant scattering*, *Commun. Pure Appl. Math.* **65** (2012) 1331.

[32] M. Andreux et al., *Kymatio: scattering transforms in python*, arXiv:1812.11214.

[33] G. Valogiannis and C. Dvorkin, *Towards an optimal estimation of cosmological parameters with the wavelet scattering transform*, *Phys. Rev. D* **105** (2022) 103534 [arXiv:2108.07821] [INSPIRE].

[34] M. Eickenberg et al., *Wavelet moments for cosmological parameter estimation*, arXiv:2204.07646 [INSPIRE].

[35] G. Valogiannis and C. Dvorkin, *Going beyond the galaxy power spectrum: an analysis of BOSS data with wavelet scattering transforms*, *Phys. Rev. D* **106** (2022) 103509 [arXiv:2204.13717] [INSPIRE].

[36] SIMBIG collaboration, *Galaxy clustering analysis with SimBIG and the wavelet scattering transform*, *Phys. Rev. D* **109** (2024) 083535 [arXiv:2310.15250] [INSPIRE].

[37] G. Valogiannis, S. Yuan and C. Dvorkin, *Precise cosmological constraints from BOSS galaxy clustering with a simulation-based emulator of the wavelet scattering transform*, *Phys. Rev. D* **109** (2024) 103503 [arXiv:2310.16116] [INSPIRE].

[38] A. Andrews, J. Jasche, G. Lavaux and F. Schmidt, *Bayesian field-level inference of primordial non-Gaussianity using next-generation galaxy surveys*, *Mon. Not. Roy. Astron. Soc.* **520** (2023) 5746 [arXiv:2203.08838] [INSPIRE].

[39] J. Stadler, F. Schmidt and M. Reinecke, *Cosmology inference at the field level from biased tracers in redshift-space*, *JCAP* **10** (2023) 069 [arXiv:2303.09876] [INSPIRE].

[40] J. Alsing, T. Charnock, S. Feeney and B. Wandelt, *Fast likelihood-free cosmology with neural density estimators and active learning*, *Mon. Not. Roy. Astron. Soc.* **488** (2019) 4440 [arXiv:1903.00007] [INSPIRE].

[41] D.H. Ribli et al., *Weak lensing cosmology with convolutional neural networks on noisy data*, *Mon. Not. Roy. Astron. Soc.* **490** (2019) 1843 [arXiv:1902.03663] [INSPIRE].
[42] M. Ntampaka, D.J. Eisenstein, S. Yuan and L.H. Garrison, *A hybrid deep learning approach to cosmological constraints from galaxy redshift surveys*, arXiv:1909.10527 [DOI:10.3847/1538-4357/ab5f5e] [INSPIRE].

[43] N. Jeffrey and B.D. Wandelt, *Solving high-dimensional parameter inference: marginal posterior densities & moment networks*, in the proceedings of the 34th conference on neural information processing systems, (2020) [arXiv:2011.05991] [INSPIRE].

[44] P. Villanueva-Domingo and F. Villaescusa-Navarro, *Learning cosmology and clustering with cosmic graphs*, Astrophys. J. 937 (2022) 115 [arXiv:2204.13713] [INSPIRE].

[45] T.L. Makinen et al., *The cosmic graph: optimal information extraction from large-scale structure using catalogues*, Open J. Astrophys. 5 (2022) 18 [arXiv:2207.05202] [INSPIRE].

[46] H. Shao et al., *Robust field-level inference of cosmological parameters with dark matter halos*, Astrophys. J. 944 (2023) 27 [arXiv:2209.06843] [INSPIRE].

[47] N.S.M. de Santi et al., *Robust field-level likelihood-free inference with galaxies*, Astrophys. J. 952 (2023) 69 [arXiv:2302.14101] [INSPIRE].

[48] A. Roncoli et al., *Domain adaptive graph neural networks for constraining cosmological parameters across multiple data sets*, in the proceedings of the 37th conference on neural information processing systems, (2023) [arXiv:2311.01588] [INSPIRE].

[49] C.H. Hahn et al., *A forward modeling approach to analyzing galaxy clustering with SIMBIG*, Proc. Nat. Acad. Sci. 120 (2023) e2218810120 [arXiv:2211.00723] [INSPIRE].

[50] SIMBIG collaboration, *Field-level simulation-based inference of galaxy clustering with convolutional neural networks*, Phys. Rev. D 109 (2024) 083536 [arXiv:2310.15256] [INSPIRE].

[51] N.-M. Nguyen et al., *How much information can be extracted from galaxy clustering at the field level?*, arXiv:2403.03220 [INSPIRE].

[52] PLANCK collaboration, *Planck 2018 results. IX. Constraints on primordial non-Gaussianity*, Astron. Astrophys. 641 (2020) A9 [arXiv:1905.05697] [INSPIRE].

[53] H.A. Feldman, N. Kaiser and J.A. Peacock, *Power spectrum analysis of three-dimensional redshift surveys*, Astrophys. J. 426 (1994) 23 [astro-ph/9304022] [INSPIRE].

[54] J.R. Fergusson, M. Liguori and E.P.S. Shellard, *General CMB and primordial bispectrum estimation I: mode expansion, map-making and measures of $f_{\text{NL}}$*, Phys. Rev. D 82 (2010) 023502 [arXiv:0912.5516] [INSPIRE].

[55] J.R. Fergusson, M. Liguori and E.P.S. Shellard, *The CMB bispectrum*, JCAP 12 (2012) 032 [arXiv:1006.1642] [INSPIRE].

[56] J.R. Fergusson, D.M. Regan and E.P.S. Shellard, *Rapid separable analysis of higher order correlators in large scale structure*, Phys. Rev. D 86 (2012) 063511 [arXiv:1008.1730] [INSPIRE].

[57] M.M. Schmittfull, D.M. Regan and E.P.S. Shellard, *Fast estimation of gravitational and primordial bispectra in large scale structures*, Phys. Rev. D 88 (2013) 063512 [arXiv:1207.5678] [INSPIRE].

[58] J. Hung, J.R. Fergusson and E.P.S. Shellard, *Advancing the matter bispectrum estimation of large-scale structure: a comparison of dark matter codes*, arXiv:1902.01830 [INSPIRE].

[59] J. Byun, A. Oddo, C. Porciani and E. Sefusatti, *Towards cosmological constraints from the compressed modal bispectrum: a robust comparison of real-space bispectrum estimators*, JCAP 03 (2021) 105 [arXiv:2010.09579] [INSPIRE].
[60] G. Jung et al., *Quijote-PNG: quasi-maximum likelihood estimation of primordial non-Gaussianity in the nonlinear dark matter density field*, Astrophys. J. 940 (2022) 71 [arXiv:2206.01624] [INSPIRE].

[61] J. Hartlap, P. Simon and P. Schneider, *Why your model parameter confidences might be too optimistic: unbiased estimation of the inverse covariance matrix*, Astron. Astrophys. 464 (2007) 399 [astro-ph/0608064] [INSPIRE].

[62] W.R. Coulton et al., *Quijote-PNG: the information content of the halo power spectrum and bispectrum*, Astrophys. J. 943 (2023) 178 [arXiv:2206.15450] [INSPIRE].

[63] G. Jung et al., *Quijote-PNG: quasi-maximum likelihood estimation of primordial non-Gaussianity in the nonlinear halo density field*, Astrophys. J. 948 (2023) 135 [arXiv:2211.07565] [INSPIRE].

[64] W.R. Coulton and B.D. Wandelt, *How to estimate Fisher information matrices from simulations*, arXiv:2305.08994 [INSPIRE].

[65] A. Heavens, R. Jimenez and O. Lahav, *Massive lossless data compression and multiple parameter estimation from galaxy spectra*, Mon. Not. Roy. Astron. Soc. 317 (2000) 965 [astro-ph/9911102] [INSPIRE].

[66] J. Alsing and B. Wandelt, *Generalized massive optimal data compression*, Mon. Not. Roy. Astron. Soc. 476 (2018) L60 [arXiv:1712.00012] [INSPIRE].

[67] F. Villaescusa-Navarro et al., *The Quijote simulations*, Astrophys. J. Suppl. 250 (2020) 2 [arXiv:1909.05273] [INSPIRE].

[68] W.R. Coulton et al., *Quijote-PNG: simulations of primordial non-Gaussianity and the information content of the matter field power spectrum and bispectrum*, Astrophys. J. 943 (2023) 64 [arXiv:2206.01619] [INSPIRE].

[69] V. Springel, *The cosmological simulation code GADGET-2*, Mon. Not. Roy. Astron. Soc. 364 (2005) 1105 [astro-ph/0505010] [INSPIRE].

[70] M. Crocce, S. Pueblas and R. Scoccimarro, *Transients from initial conditions in cosmological simulations*, Mon. Not. Roy. Astron. Soc. 373 (2006) 369 [astro-ph/0606505] [INSPIRE].

[71] R. Scoccimarro, L. Hui, M. Manera and K.C. Chan, *Large-scale bias and efficient generation of initial conditions for non-local primordial non-Gaussianity*, Phys. Rev. D 85 (2012) 083002 [arXiv:1108.5512] [INSPIRE].

[72] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[73] M. Davis, G. Efstathiou, C.S. Frenk and S.D.M. White, *The evolution of large scale structure in a universe dominated by cold dark matter*, Astrophys. J. 292 (1985) 371 [INSPIRE].

[74] N. Dalal, O. Dore, D. Huterer and A. Shirokov, *The imprints of primordial non-Gaussianities on large-scale structure: scale dependent bias and abundance of virialized objects*, Phys. Rev. D 77 (2008) 123514 [arXiv:0710.4560] [INSPIRE].

[75] G. Jung et al., *Quijote-PNG: optimizing the summary statistics to measure primordial non-Gaussianity*, arXiv:2403.00490 [INSPIRE].

[76] M. Peloso and M. Pietroni, *Galilean invariance and the consistency relation for the nonlinear squeezed bispectrum of large scale structure*, JCAP 05 (2013) 031 [arXiv:1302.0223] [INSPIRE].

[77] A. Kehagias and A. Riotto, *Symmetries and consistency relations in the large scale structure of the universe*, Nucl. Phys. B 873 (2013) 514 [arXiv:1302.0130] [INSPIRE].
[78] P. Creminelli, J. Noreña, M. Simonović and F. Vernizzi, *Single-field consistency relations of large scale structure*, JCAP **12** (2013) 025 [arXiv:1309.3557] [INSPIRE].

[79] M. Peloso and M. Pietroni, *Ward identities and consistency relations for the large scale structure with multiple species*, JCAP **04** (2014) 011 [arXiv:1310.7915] [INSPIRE].

[80] M. Peloso and M. Pietroni, *Galilean invariant resummation schemes of cosmological perturbations*, JCAP **01** (2017) 056 [arXiv:1609.06624] [INSPIRE].

[81] G. D'Amico, M. Marinucci, M. Pietroni and F. Vernizzi, *The large scale structure bootstrap: perturbation theory and bias expansion from symmetries*, JCAP **10** (2021) 069 [arXiv:2109.09573] [INSPIRE].

[82] C.F. Park, E. Allys, F. Villaescusa-Navarro and D.P. Finkbeiner, *Quantification of high-dimensional non-Gaussianities and its implication to Fisher analysis in cosmology*, Astrophys. J. **946** (2023) 107 [arXiv:2204.05435] [INSPIRE].