Quijote-PNG: Quasi-maximum likelihood estimation of Primordial Non-Gaussianity in the non-linear halo density field

We study primordial non-Gaussian signatures in the redshift-space halo field on non-linear scales, using a quasi-maximum likelihood estimator based on optimally compressed power spectrum and modal bispectrum statistics. We train and validate the estimator on a suite of halo catalogues constructed from the Quijote-PNG N-body simulations, which we release to accompany this paper. We verify its unbiasedness and near optimality, for the three main types of primordial non-Gaussianity (PNG): local, equilateral, and orthogonal. We compare the modal bispectrum expansion with a $k$-binning approach, showing that the former allows for faster convergence of numerical derivatives in the computation of the score-function, thus leading to better final constraints. We find, in agreement with previous studies, that the local PNG signal in the halo-field is dominated by the scale-dependent bias signature on large scales and saturates at $k \sim 0.2~h\,\mathrm{Mpc}^{-1}$, whereas the small-scale bispectrum is the main source of information for equilateral and orthogonal PNG. Combining power spectrum and bispectrum on non-linear scales plays an important role in breaking degeneracies between cosmological and PNG parameters; such degeneracies remain however strong for equilateral PNG. We forecast that PNG parameters can be constrained with $\Delta f_\mathrm{NL}^\mathrm{local} = 45$, $\Delta f_\mathrm{NL}^\mathrm{equil} = 570$, $\Delta f_\mathrm{NL}^\mathrm{ortho} = 110$, on a cubic volume of $1 \left({ {\rm Gpc}/{ {\rm h}}} \right)^3$, at $z = 1$, considering scales up to $k_\mathrm{max} = 0.5~h\,\mathrm{Mpc}^{-1}$.


INTRODUCTION
The coming generation of spectroscopic and photometric galaxy surveys -e.g., Euclid, DESI, Spherex, Rubin Observatory, Roman (Laureijs et al. 2011;DESI Collaboration et al. 2016;Doré et al. 2014;LSST Science Collaboration et al. 2009) -will allow us to study galaxy clustering with an unprecedented level of accuracy and precision, shedding further light on many open questions in cosmology. Among the many exciting possibilities, an interesting prospect, which we mainly focus on in this work, will be that of improving our under-standing of Early Universe physics, via high precision tests of Primordial non-Gaussianity (PNG).
Cosmic Microwave Background (CMB) measurements (Akrami et al. 2020), in agreement with theoretical expectations, have constrained the primordial cosmological perturbation field to be at most weakly non-Gaussian. This implies, for a large majority of Early Universe scenarios, that most of the PNG information is contained in the primordial bispectrum. For this reason, the bispectrum of dark matter tracers (e.g., galaxies) in Large Scale Structure (LSS) can be a powerful probe of PNG. Crucially, the 3D galaxy bispectrum gives us also access, in principle, to a larger number of modes with respect to the 2D (angular) bispectrum of CMB anisotropies. Therefore, LSS bispectrum analyses can potentially lead to significant improvements in PNG constraints over current, CMB-based, results. Achieving such improvements will require however to include non-linear scales in the analysis, carrying strong non-Gaussian (NG) signatures which are not primordial, but arise from late-time, non-linear evolution of cosmic structures. Disentangling the NG late time component from the subdominant primordial one is therefore a crucial challenge in this kind of studies. It can be addressed either by analytical modeling of the bispectrum -via a suitable perturbative approach at mildly non-linear scales (Cabass et al. 2022a,b;D'Amico et al. 2022) -or by relying on fully numerical approaches, which evaluate the bispectrum (and/or other summary statistics, see, e.g., Biagetti et al. 2021;Friedrich et al. 2020; Valogiannis & Dvorkin 2022) using large mock datasets; field-level inference on large scales, not relying on specific statistical summaries, has also been recently considered, see Andrews et al. (2023). In this work, which is the fourth in a series of papers, following Jung et al. (2022); Coulton et al. (2023a,b), we base our study on the Quijote-png simulation suite, recently presented in Coulton et al. (2023a). Our main goal is to quantify the accuracy with which f NL can be constrained using both the power spectrum and bispectrum of the dark matter halo field, up to strongly non-linear scales (k max = 0.5 h Mpc −1 ), using simulations with different kinds of PNG. More precisely, we consider the three main PNG bispectrum shapes, namely, the local, equilateral, and orthogonal shapes, which are predicted in a large variety of inflationary scenarios.
This work extends our initial analysis presented in Jung et al. (2022), where we worked at the level of the matter field. As in the previous analysis, we derive forecasts for PNG and standard cosmological parameters, by combining power spectrum and bispectrum measurements at non-linear scales; our main focus is then on building an unbiased and nearly optimal quasi-maximum likelihood estimator, based on applying a MOPED-like compression algorithm to a modal decomposition of the data bispectrum. In our companion paper (Coulton et al. 2023b) we independently perform a similar analysis at a different redshift (z = 0 in Coulton et al. 2023b, vs. z = 1 in this work), but we employed a binned decomposition in k-space instead of the modal approach adopted here; in that work we studied in detail the PNG information content in the halo field while focusing on important numerical convergence issues. Therefore, the two analyses are comple-mentary to study the robustness of our approach and to cover a full range of crucial issues, from numerical stability, to the precise quantification of the information gain obtained from different observables at different scales and the demonstration of nearly optimal and unbiased bispectrum data compression for parameter estimation. Taken together, we think these works represent an important development in the effort to build a data analysis pipeline to be applied to observations. We release the halo catalogues of the Quijote-png 1 suite used in these works.
Since in the current work we consider tracers of the underlying density field, an additional signature of PNG arises -in comparison to the previous matter field analysis -in the form of a scale-dependency in the tracer bias. Such feature has a power law behaviour, with degree determined by the squeezed limit of the PNG bispectrum shape under study: it is most prominent for local NG, with a ∼ 1/k 2 behaviour and absent in the equilateral case. Scale-dependent bias has been the object of significant study in the literature, see, e.g., Dalal et al. (2008) Giri et al. (2023), and was used to extract local PNG constraints from BOSS data Slosar et al. (2008); Ross et al. (2013); Leistedt et al. (2014); Mueller et al. (2021); Cabass et al. (2022a);D'Amico et al. (2022). Recently, it has been however pointed out that accurate modeling of scale-dependent bias from PNG also depends on details of galaxy formation, making its use as a tool to measure the PNG amplitude f NL significantly more challenging than previously thought Barreira (2020Barreira ( , 2022a. The effect of scale-dependent bias is automatically incorporated in our analysis, where we are mostly concerned with assessing its relative constraining power on different NG shapes, as compared to the bispectrum, and verifying agreement with both our analysis in Coulton et al. (2023b) and theoretical expectations (see, e.g. de Putter 2018; Karagiannis et al. 2018).
The paper is structured as follows. In section 2 we briefly review the NG models considered in the analysis. In section 3 we illustrate our methodology for data compression and parameter estimation. In section 4 we discuss our Fisher matrix analysis, showing expected parameter constraints on different scales, and describe the application of our quasi-maximum likelihood, joint power spectrum and bispectrum estimator to simulated data. In section 5 we summarize our main results and draw our final conclusions. Finally, in appendix A we provide more details about the implementation of shotnoise modes in the bispectrum estimator. In appendix B we show a comparison between the modal and binned approaches to bispectrum estimation, and in appendix C we discuss the results of a preliminary study aimed at the application of the CARPool technique to the evaluation of covariances and numerical derivatives.

BISPECTRUM SHAPES
Violating any condition of the standard inflationary model induces a deviation from the perfect Gaussian initial conditions, which leads to non-zero high-order correlators. The largest of them, in most inflationary models, is the bispectrum, i.e. the three-point correlation function of Fourier modes, defined as: (1) The primordial bispectrum is generally written as where f NL is the dimensionless amplitude parameter corresponding to a given primordial bispectrum shape F (k 1 , k 2 , k 3 ), which encompasses the dependence of the bispectrum on triplets of Fourier space modes.
In this work, we focus on building estimators to measure f NL of three of the most common primordial shapes, namely the local, equilateral and orthogonal 2 bispectra (see Coulton et al. 2023a, and references therein for the complete description of these templates).
For dark matter tracers (e.g. halos), the presence of PNG has a significant impact, due to the introduced coupling between large and small scale modes, with the most known example being that of the local type. In this case, the halo overdensity on large scales will no longer depend only on the matter overdensity, but also on the primordial gravitational potential (see Desjacques et al. 2018, for a review). This results in a scale-dependent term that introduces an important PNG signature on the large scales of a correlator. This is of particular importance to the power spectrum of the observed dark matter tracers, since it enhances the PNG signal within the two-point correlation function, which otherwise would have been very limited, e.g. in the case of a dark matter field analysis (see e.g. Coulton et al. 2023a;Jung et al. 2022).
The effect of the scale dependent term on the power spectrum has been extensively studied in the literature, especially for the local PNG type. However, recent developments have made the measurement of f NL , by such a term in the power spectrum, challenging, due to the perfect degeneracy between f NL and the scale dependent bias coefficient b φ (Barreira 2020(Barreira , 2022b. For the halo bispectrum, a significant amount of the PNG signal is located within the primordial part (eq. 2), while the scale dependent terms, studied at a theoretical level e.g. in Karagiannis et al. (2018), that could carry a notable amount of information on local PNG, suffer from the same limitations as the power spectrum (Barreira 2022a).
The effect of the scale-dependent bias will be taken into account within the framework of the forwardmodeling. In a simulation-based approach we assume tight priors on the scale dependent bias coefficient b φ , in order to focus on the f NL constraints (see also Coulton et al. 2023b, for an extensive discussion).

METHOD
In this section we review the main aspects of our methodology for data compression and quasi-maximum likelihood estimation of cosmological and PNG parameters, starting from the evaluation of power spectrum and modal bispectrum summary statistics.

Quasi maximum-likelihood estimator
Starting from a given data vector d (a given set of summary statistics, like the power spectrum and/or the bispectrum) that depends on some parameters of interest denoted θ (e.g. f NL ), one can write the following quasi maximum-likelihood estimator for the value of the parameters (see Alsing & Wandelt 2018, for details): where the subscript * denotes that the quantities are evaluated at some chosen fiducial point, and µ and C are, respectively, the mean and the covariance of d.
The two key ingredients of this estimator, which are the Fisher information F and the compressed score statistic t, will be detailed below. Note also that in this expression we assume a Gaussian likelihood and a dependence on parameters through the mean only, a reasonable assumption as verified in Jung et al. (2022). The Fisher matrix, a standard method to evaluate the information content of some observables, is given by: This requires knowledge of the derivatives ∇ θ µ and the covariance C, which can be both evaluated from a large set of simulations, as we do in this work (see section 4). However, reaching numerical convergence for the joint analysis of multiple parameters may be very challenging and therefore prone to wrong results. In Jung et al. (2022) we checked, by studying the matter field, that if the covariance matrix is not converged, it would typically induce suboptimal error bars, and that non-converged derivatives could bias the estimated parameters. In Coulton et al. (2023b) we showed that noisy derivatives could lead to overconfident error bars when working with the halo field.
To tackle this problem, the alternative compressed Fisher method described in Coulton et al. (2023b) (see also Coulton & Wandelt 2023, for details) can provide conservative bounds. It consists of two steps, the first of which is to compress the data to the score function using (see Alsing & Wandelt 2018) that is equivalent to the MOPED compression scheme of Heavens et al. (2000). This operation reduces the data vector d of size n down to only p numbers, where p is the number of parameters of interest, while keeping all relevant information about these parameters. Then, to compute the compressed Fisher matrix one has only to apply the standard expression of eq. (4) to the compressed data. An important subtlety of this scheme is that it requires to use two separate sets of simulations for the two different steps. The first part is used for the compression step, to build a new summary statistics, which will be suboptimal if derivatives are noisy. The second part is then compressed and used to estimate the Fisher matrix from the compressed statistics, which is suboptimal if the compression step is suboptimal, but is also a lot less noisy due to the much lower dimensionality of the compressed statistics. 3

Summary statistics
In this work, we use the same observables as in Jung et al. (2022), based on the power spectrum and bispectrum statistics as they contain significant and complementary information about both the ΛCDM cosmological parameters and the PNG amplitudes f NL .
The standard power spectrum estimator of a field δ(k) defined on a grid of fundamental mode k f is given bŷ where V is the survey volume, and a binning of the krange has been introduced with each bin ∆ i having a width k f and containing N i independent vectors of k.
As was initially shown for CMB NG analysis in Fergusson et al. (2010Fergusson et al. ( , 2012a, and extended later to LSS in Fergusson et al. (2012b); Regan et al. (2012); Schmittfull et al. (2013) (see also Lazanu et al. 2016Lazanu et al. , 2017Hung et al. 2019a,b;Byun et al. 2021;Byun & Krause 2022), the bispectrum information can be efficiently extracted from data by measuring the following modal coefficients: where for a well-chosen basis of one-dimensional functions q p (k) and mode triplets n ↔ (p, q, r). We refer the reader to Jung et al. (2022) for the details of the exact setup we use for the analyses presented in section 4, as they are almost identical (the only change is the addition of two special modes, introduced in Byun et al. (2021) and recalled in appendix A, describing the shot-noise component of the bispectrum expected from halos).

Specifications
For our analysis we use the publicly available Quijote 4 and Quijote-png 5 suites of N-body simulations Coulton et al. 2023a). Each simulation represents a periodic cubic box of length 1 h −1 Gpc, which contains 512 3 particles, run with the Gadget-III code (Springel 2005). Initial conditions are generated at z i = 127 with the codes 2LPTIC (Crocce et al. 2006) in the Gaussian case, and 2LPTPNG 6 in the non-Gaussian case (Scoccimarro et al. 2012;Coulton et al. 2023a); linear matter power spectra and transfer functions are obtained from CAMB (Lewis et al. 2000). Finally, dark matter halos are identified using the Friends-of-Friends algorithm (Davis et al. 1985) with a value of the linking length equal to b = 0.2; we select those with a mass larger than M min = 3.2 × 10 13 M /h, corresponding to a number densityn ∼ 5.1 × 10 −5 h 3 Mpc −3 at z = 1 (see Hahn et al. 2020, for a power spectrum and bispectrum analysis of these halo catalogues focused on cosmological parameters). We construct the halo density field in redshift-space at z = 1 by depositing the halo positions, displaced radially by the velocity, on a grid of size N grid = 256, using a fourth-order interpolation scheme implemented in the Pylians3 code 7 (Villaescusa-Navarro 2018). We then measure the power spectrum and modal bispectrum monopoles using the estimators (6) and (7) including modes up to k max = 0.5 h Mpc −1 .
For the numerical computation of the covariance matrix, we have 15, 000 simulations at fiducial cosmology, whereas smaller sets of 500 realizations, with varying input parameters, are used to evaluate the derivatives in eq. 5. In particular the analysis is focused on the PNG amplitudes of the three shapes considered here, f local NL , f equil NL , f orth NL ; four cosmological parameters Ω m , n s , σ 8 and h; and one parameter related to the halo bias, M min . The variation of the minimum halo mass generates distinct catalogs with M fid min ± ∆M min (see table 1), which consequently leads to a variation of the halo number density. This is roughly equivalent to a variation of the linear bias contribution (see e.g. Desjacques et al. 2018, for details), while it propagates, to a minor extent, to higher order terms. Although, this bias model is quite simplistic, it is still useful within the framework of a first-order analysis presented in this work. A thorough investigation on the impact of the bias parameters, within a simulation-based approach, requires the population of halos with a HOD and the variation of the HOD parameters, which is left for future work. More details about the specifications of these simulations, concerning all the parameters considered in our analyses, can be found in table 1.

Fisher constraints
We aim to evaluate the information content on ΛCDM parameters and PNG amplitudes contained in the power spectrum and bispectrum of the halo field at redshift z = 1 using a Fisher matrix formalism. This analysis 7 https://github.com/franciscovillaescusa/Pylians3 complements the work of Coulton et al. (2023a), as we focus on a different redshift and make use of a different bispectrum estimator. We explore the dependence on the number of simulations used, the chosen k max or the role of the different summary statistics. We show the results in figures 1, 2, 3 and table 2.
As highlighted in Coulton et al. (2023a), a main difficulty of this simulation-based approach is to accurately compute numerical derivatives of both the power spectrum and bispectrum with respect to the different parameters considered. This is illustrated in figure 1, where we show that using smaller subsets of the 500 available pairs of simulations per parameter leads to spurious smaller 1-σ uncertainties when analyzing jointly Instead of the computationally intensive possibility of producing many more simulations, we use here an alternative method, described briefly in section 3.1 to compute conservative constraints from a lower number of simulations. As expected, the resulting 1-σ error bars decrease when we use more simulations to calculate numerical derivatives, and using the full set they are only between 5% and 25% larger than the unconverged standard Fisher constraints. The results are stable when using 250 pairs of simulations or more to compute each derivative.
In the simpler situation where we consider only the three following parameters {σ 8 , Ω m , f local NL } in the analysis, the two methods give very similar results when numerical convergence is reached (using at least 100 pairs of simulations per derivative). This is why in the rest of this work we always use the conservative approach, knowing it is equivalent to the standard Fisher approach in the cases where numerical accuracy can be reached with the available simulations, and otherwise only leads to a reasonable overestimation of order 10% as verified. Note that we manage to keep this overestimation small here due to the use of the modal bispectrum, rather than a standard "binned" approach, because it compresses the original data more efficiently leading to more stable numerical derivatives (we need less than 50 modes to extract the full information of the bispectrum up to k max = 0.5 h Mpc −1 ). This is shown explicitly in appendix B.
In figure 2 we study the dependence of the constraints on k max , considering values from 0.1 to 0.5 h Mpc −1 . The largest improvement (for both ΛCDM cosmological and PNG parameters) is obtained between k max = 0.1 and 0.2 h Mpc −1 , at which point error bars on f local NL become saturated (as well as for h). However, for the equilateral and orthogonal shapes considering smaller scales yields better constraints (a few percent for each additional increase of 0.1 h Mpc −1 ). For other parameters, the gain can even be larger, justifying the need to probe these nonlinear scales. Note also that all these improvements are computed using the conservative error bars obtained from the compressed summary statistics. Including smaller scales in the analysis typically leads to less converged numerical derivatives, and the less converged these derivatives are the more suboptimal the conservative approach becomes. This means that we may be underestimating slightly the constraining power of the small scales (in any case, this effect should not be large, as for k max = 0.5 h Mpc −1 this overestimation is ∼ 10% as can be seen in figure1).
In figure 3 we show the information content of the halo power spectrum, the halo bispectrum, and their combination. The bispectrum is a much more efficient probe of the equilateral and orthogonal shapes than the power spectrum, while for other parameters they yield constraints of the same order separately. Their combination always helps to reduce degeneracies, although to a lesser extent than for the matter field studied previously in Coulton et al. (2023a); Jung et al. (2022).
In table 2, we present the 1-σ conservative constraints on ΛCDM parameters and PNG amplitudes using jointly the power spectrum and bispectrum and including small scales up to k max = 0.5 h Mpc −1 . Unlike the matter field case discussed in Coulton et al. (2023a); Jung et al. (2022), including PNG shapes in the analysis increases slightly error bars on ΛCDM cosmological parameters. The different PNG shapes are also less correlated, as analyzing them jointly increases only slightly their own error bars.

Parameter estimation
As was shown in Jung et al. (2022) for the matter field, the simple quasi maximum-likelihood estimator (see eq. 3) built from the Fisher matrix at some chosen fiducial cosmology is very efficient to measure ΛCDM cosmological parameters and PNG amplitudes using the power spectrum and bispectrum information. Here we extend this conclusion to the halo field.
The key ingredient of the estimator is the Fisher matrix, which in this work is fully evaluated from a very large set of simulations. As discussed in the previous section, we use a two-step conservative approach for its computation leading to slightly suboptimal results, because numerical convergence is difficult to reach with the standard method. We verify that this leads nonetheless to unbiased and near-to-optimal measurements of parameters, by estimating jointly σ 8 , Ω m , n s , h, f local NL , f equil NL and f ortho NL in the Quijote simulations using both the power spectrum and the bispectrum.
In figure 4, we study the effect of varying k max on the 1-σ error bars of the quasi-maximum likelihood estima-  Table 2. Joint 1-σ error bars on cosmological parameters and PNG from the power spectrum and the modal bispectrum of the halo field at z = 1, at kmax = 0.5 h Mpc −1 .
In the first part, we report the Fisher constraints described in section 4.2 and in the second part the corresponding error bars of the quasi-maximum likelihood estimator used in section 4.3. We analyzed 15000 Quijote halo catalogues of 1 (Gpc/h) 3 volume at fiducial cosmology, and sets of 500 simulations with one adjusted parameter. tor. To compute these error bars, we use a set of 1000 simulations at fiducial cosmology and analyze it with the estimator calibrated using all other simulations (using 14000 simulations instead of 15000 to calculate the covariance has been verified to have no impact on the results). We repeat the procedure for different sets of 1000 simulations, and compute the standard deviation of the results. As expected, this highlights a very similar behaviour as for the Fisher constraints discussed in the previous section. Concerning PNG, there is no im-provement for f local NL above k max = 0.2 h Mpc −1 , while for the other two shapes there is no clear saturation yet (although the gain between k max = 0.4 and 0.5 h Mpc −1 is only a few percent). For every other parameter considered (except h), the decreasing of error bars is significant up to k max = 0.5 h Mpc −1 . In table 2, we report the corresponding error bars at k max = 0.5 h Mpc −1 , considering cosmological parameters only or jointly with the PNG shapes. For all parameters, the error bars of the quasi-maximum likelihood estimator are close to, or slightly larger than the Fisher constraints reported in the same table (less than 10% difference).
In figure 5, we compare the estimated parameters to their input values for different cases, focusing here on changes of PNG amplitudes. We first study the mildly nonlinear regime (k max = 0.2 h Mpc −1 ) and then include also nonlinear scales (k max = 0.5 h Mpc −1 ). The measured parameters match their expected values for both ranges of scales when studying datasets at fiducial cos-mology or having PNG of the equilateral or orthogonal types (f equil NL = +100 or f ortho NL = +100). There are however large statistical deviations on several parameters for the simulations with local NG (in particular f local NL , several datasets giving a value more than 5-σ away from the expected one). This difference of behaviour between this specific set and the others can be explained, by the fact that f local NL = 100 is more than 2-σ away from the fiducial value of f NL = 0 (based   Figure 5. Relative difference of measured cosmological parameters and PNG amplitudes using the quasi maximum-likelihood estimator (eq. 3) with their expected value. We use the power spectrum and the bispectrum of the halo field jointly, for kmax = 0.2 h Mpc −1 in the top row and kmax = 0.5 h Mpc −1 below. Each column corresponds to a given parameter (cosmological or PNG). Each panel corresponds to a different input cosmology of the data samples (i.e. one with Gaussian initial conditions and the three types of PNG). For each input cosmology, we analyze five independent datasets of 100 realizations, each being indicated by its own colour and marker. The dark and light grey bands represent, respectively, the 2 and 1-σ intervals around the expected deviation (0).  on error bars given in table 2) while f ortho NL = 100 and f equil NL = 100 are respectively smaller and a few times smaller than a 1-σ deviation from f NL = 0. The NG simulations of the three shapes correspond to different regimes where a parameter is more or less displaced from the model we use to calibrate the estimator. This is confirmed in figure 6, where we check that simulations with f local NL = 50 (thus roughly a 1-σ deviation) give this time the expected results.
These tests confirm the unbiasedness of the quasi maximum-likelihood estimator, with the caveat that the estimator must be calibrated relatively close to the actual parameter values. This, of course, is due to the fact that the entire method is based on a linear approxi- mation of the likelihood around the fiducial parameters. For the same reason, however, it is clear that the issue can be immediately addressed -at the computational cost of producing new sets of simulations -by implementing a standard recursive procedure, in which the estimated parameters at the previous step generate the new fiducial model for the following step, until convergence. Note that this scenario is not bound to occur in practice, since current cosmological parameter constraints from, e.g., CMB datasets such as Planck produce already quite narrow priors.
While it was shown in section 4.2 that using a lower number of simulations to compute derivatives leads to more suboptimal Fisher matrices, it is also important to verify the effect of changes in the number of simulations used to compute the covariance matrix. We explore this in figure 7, where we show the increase of error bars due to using fewer simulations. Above 1000 simulations, error bars of the quasi-maximum likelihood estimator are stable (variations at the percent level) and close to the Fisher estimates (10% difference at most).

CONCLUSIONS
In this paper, we have developed a joint power spectrum and bispectrum quasi-maximum likelihood estimator of cosmological and PNG parameters and applied it to the study of the halo field in the Quijote-png simulation suite. The data analysis pipeline applies the optimal data compression methodology developed in Alsing & Wandelt (2018); Heavens et al. (2000) to a set of power spectrum and modal bispectrum summary statistics, efficiently extracted from the input mock realizations. In this way, we extended our previous analysis (Jung et al. 2022), which considered the matter field in the same dataset.
The main arising technical complication was related to the convergence of numerical derivatives that are used to compute the Fisher information and to perform the final compression step. This turns out to be much slower now, with respect to the previous matter field analysis, now leading to potential problems such as spurious "superoptimal" error bars in the final estimator. Interestingly, though, we have also found that our modal decomposition of the bispectrum makes derivative convergence much faster with respect to the binning approach we implemented in Coulton et al. (2023b). Although still not sufficient for a brute force computation with the available realizations, such faster convergence suggests that more investigation should be done in the future to find the optimal bispectrum decomposition scheme, for the best numerical stability. In the meantime, to circumvent the issue, we have implemented the method first described in Coulton et al. (2023b), which is based on computing the Fisher matrix of MOPED-compressed statistics, extracted from an independent simulation set. This approach leads to stable, robust results, at the price of slight suboptimality in the final estimator. Despite such small suboptimality, we have verified that the forecasted errors significantly improve after including nonlinear scales up to k max = 0.5 h Mpc −1 (see figure 2 as a summary of our main results), in agreement with our findings in Coulton et al. (2023b). Given the significant contribution provided by small scale, shot-noise dominated, bispectrum triangles, further improvements could be in principle achieved in a future galaxy density analysis, by selecting higher-density tracers. In contrast to other parameters, we have observed a saturation of the f local NL error at a scale k ∼ 0.2 h Mpc −1 ; this is again consistent with our previous findings and with other forecasts, such as those in Karagiannis et al. (2018), where it was shown that the f local NL signal is dominated by the scale-dependent bias signature, on large scales, both in the power spectrum and in squeezed bispectrum configurations.
After investigating the power spectrum and bispectrum information content on non-linear scales, the final step of our analysis consisted in testing our quasimaximum likelihood estimator on the simulated dataset. We have verified that we can recover unbiased results, deep into the non-linear regime, up to k max = 0.5 h Mpc −1 (see figures 5 and 6). Unbiasedness is of course verified only provided the starting fiducial parameter values in the estimator are close enough to the real ones. We studied this in more detail by varying the input value of f local NL in the analyzed simulations and verifying that biased results are obtained when the true f local NL in the data is ∼ 2σ away from the fiducial choice in the estimator. In a realistic observational scenario, this issue can of course always be addressed by implementing a recursive estimation procedure, which however becomes more and more expensive, by requiring new mock realizations and re-calibration of the estimator weights at each step. This suggests to investigate the possibility to reduce the overall computational cost of simulations. We have started a preliminary analysis in this direction, using the CARPool method , which is further discussed in Appendix C. Another possibility is to use machine-learning-augmented simulations, see Kaushal et al. (2022); Jamieson et al. (2022); Piras et al. (2022) for examples. Making use of these different techniques will play a key role in enabling simulationbased inference with the upcoming generation of galaxy surveys, which will have a much higher tracer density.
The recovered error bars are, as expected, slightly larger than the optimal Fisher bound. This is a direct consequence of the fact that, to secure unbiasedness and robustness of the results, we have calibrated the estimator weights using the stable, yet conservative approximation of the Fisher matrix described above. Also in this case though, the slight suboptimality does not prevent us from obtaining large improvements in precision for the final parameter estimates, when we include nonlinear scales in the analysis (see figure 4). By extending our previous analysis to the halo field in redshift space, we have made a significant step forward toward the final development of an efficient, joint power spectrum and bispectrum estimation pipeline, able to extract cosmological and PNG parameters at strongly non-linear scales from actual observations. In a follow-up work we will further extend the current analysis, by looking at the galaxy density field, simulated via a suitable Halo Occupation Distribution (HOD), following Hahn & Villaescusa-Navarro (2021). Marginalization over HOD parameters will also allow us to significantly improve the accuracy of our bias model, which is currently defined by a single parameter which describes the leading order contribution and only to a minor extent captures higher order effects.
Our conclusions are in full agreement with those in our companion work, Coulton et al. (2023b), where we performed an independent analysis at a different redshift (z = 0 in Coulton et al. 2023b, vs. z = 1 in this paper) and used a standard binning scheme for the bispectrum, rather than the modal approach developed here. Besides increasing the robustness of our conclusions via cross-validation of independent data analysis pipelines, the two works complement each other in several ways and together cover a significant range of crucial issues: Coulton et al. (2023b) focused on addressing numerical stability issues, on assessing the information content of our observables at different scales and on evaluating in detail all possible contributions to the error budget (such as, e.g., shot noise and super-sample covariance effects), whereas the present study, while cross-checking the previous Fisher matrix results, is more centered on optimal data compression and on the development and testing of related statistical estimators. LV acknowledges ERC (BePreSySe, grant agree-ment 725327), PGC2018-098866-B-I00 MCIN/AEI/10.13039/501100011033 y FEDER "Una manera de hacer Europa", and the "Center of The left column is obtained with the modal bispectrum estimator used throughout this paper, while the two others use a standard "binned" approach for different widths of bin (3k f in the middle panels and 2k f in the right panels). Otherwise, this figure is similar to 1, for kmax = 0.2 h Mpc −1 .

A. SHOT NOISE MODAL MODES
The shot-noise contribution to the matter bispectrum at tree-level is given by wheren is the halo number density and P L (k) is the linear matter power spectrum. As introduced in Byun et al. (2021), this can be fully described in the modal way by using the two triplets (0, 0, 1) and (0, 0, 0) combining the following one-dimensional basis functions B. COMPARISON WITH THE STANDARD "BINNED" BISPECTRUM ESTIMATOR A key ingredient to compute the Fisher matrix (eq. 4), which is used both for constraint forecasts and to build estimators, is to have accurate derivatives of the summary statistics with respect to the different parameters considered. As discussed in section 4.2, even the large sets of 500 paired simulations for each parameter of the Quijote and Quijote-png collections are not sufficient to reach the necessary numerical convergence for the power spectrum and bispectrum derivatives. This typically leads to an underestimation of 1-σ error bars.
On the other hand, a two-step computation, consisting first on compressing optimally the data, and then computing the Fisher matrix from this compressed data (using different datasets for the two steps), yields slightly overestimated error bars. Combining the power spectrum and the bispectrum information of the halo field at z = 1, we have verified that even when we include nonlinear scales up to k max = 0.5 h Mpc −1 , the difference between the lower and upper bounds of constraints is at most of order 20% on the different parameters, a very reasonable difference. In this appendix, we show that the modal estimator, which by construction compresses the bispectrum information in the data, is a necessary ingredient for the efficiency of the method.
In figure 8, we compare the convergence of standard and conservative Fisher 1-σ uncertainties obtained with the modal bispectrum (as in the rest of this paper) and a standard "binned" bispectrum estimator. 8 A main result is that

P(k)
Power spectrum . The CARPool method applied to the power spectrum (left column) and modal bispectrum (right column) of the halo field, at z = 1. On the top row, the black dashed lines correspond to the averages from the 15000 Quijote simulations at fiducial cosmology with f local NL = 0 (note that in the bispectrum case, all modal coefficients are normalized by dividing by the modes from these 15000 simulations at fiducial cosmology). The blue lines correspond to the average from 500 simulations with f local NL = +100. The red dotted lines have been computed using the CARPool method (see eq. C5), using 10 simulations at f local NL = +100 as the high-fidelity simulations and the 15000 simulations at fiducial cosmology as surrogates. The blue areas and red vertical lines show the respective error bars from the two cases (they correspond to standard errors for sets of 10 simulations, and have been multiplied by a factor 10 for visibility in the power spectrum case). In the bottom row, we show the difference between the CARPool estimates and averages from 500 simulations, normalized by the standard deviation. The blue areas correspond to the standard error for 10 simulations, and error bars. Error bars on the CARPool estimates, shown in red, are computed by applying the CARPool method to many different sets of 10 simulations at f local NL = +100. the constraints obtained with the modal estimator are the most stringent, with a difference of order 10% with the standard estimator with bins of width 3k f , and even more with smaller bins of width 2k f . Indeed the estimator using the smallest bins gives here the largest error bars, despite the fact that in principle it should keep more information, due to the greater difficulty of computing sufficiently accurate numerical derivatives. This lack of convergence is also very clear when we compare the lower and upper bounds on error bars for all three methods. Using the full sets of simulations, the lower bounds are 10-20% smaller than the upper limits for the modal estimator, 30% for bins of width 3k f , and as much as two times smaller for bins of width 2k f . The modal estimator gives more stringent constraints, which are proven to be closer to the actual Fisher uncertainties, and should converge totally with a smaller number of simulations, as shown in the first row in the simple situation, where the modal estimator has fully converged and the other two have not.

C. APPLICATION OF CARPOOL
As we have verified in this paper, the quasi-maximum likelihood estimator is a powerful method to infer cosmological parameters and PNG amplitudes from halo catalogues using information beyond the mildly non-linear regime, which however, as other simulation-based methods, can require a large number of costly forward simulations. Therefore, a key component of future applications will be to include the variance reduction CARPool technique, developed in ; , 2022, into the full analysis pipeline. The basic idea behind CARPool is to use a relatively small number of high fidelity simulations combined with a large number of less accurate simulations, or surrogates, to measure some chosen summary statistics with much smaller error bars. In , these surrogates were computed using much faster, but less precise, N-body solvers like COLA (Tassev et al. 2013). This could for example be applied to the case of numerical derivatives, for which reaching numerical convergence typically requires thousands of costly simulations. We leave this application for future work, and instead focus here on the use of CARPool to speed-up the iteration process of quasi-maximum likelihood estimation.
To obtain unbiased estimates of cosmological parameters or f NL 's, it is important that the fiducial cosmology where we evaluate the covariance and numerical derivatives is not too far from the actual parameter values. For example, in section 4, we have seen that with a fiducial cosmology at f local NL = 0, the estimator is unbiased in measuring f local NL in simulations with an input of f local NL = 50, but not for f local NL = 100. Note that even if the measured bias were large when averaging from hundreds of simulations, it would still be smaller than the 1-σ error bar, making it a good first estimate. Then, working by iteration and choosing a new fiducial cosmology at these roughly measured parameters should yield unbiased results. To avoid producing a completely new large set of simulations at the new fiducial cosmology, we can consider the original simulations as the surrogates of the CARPool method (the idea of combining simulations at different cosmology was also explored in Ding et al. 2022).
The main ingredients of the CARPool method are: • A set of N paired high-fidelity simulations and surrogates, sharing the same random seeds to produce their initial conditions, from which we measure some chosen summary statistic denoted y or c (simulation or surrogate respectively) and the corresponding sample covariance given bŷ • A separate set of M surrogates, to compute the mean of c with the standard expression Then, the key quantity to compute is which by construction has the same ensemble average as y (i.e.x =ȳ). The variance of x is minimized when the control matrixβ is given byβ = Σ yc Σ −1 cc .
In , it was shown that a very efficient choice, using only the diagonal elements of Σ yc and Σ cc , is the following diagonal control matrix: β diag = diag cov(y 1 , c 1 ) σ(c 1 ) 2 , cov(y 2 , c 2 ) σ(c 2 ) 2 , ..., cov(y n , c n ) σ(c n ) 2 , where n is the size of the vectors y and c.
In figure 9, we show the results obtained with the CARPool technique applied to the Quijote-png set of halo catalogues. We use the large set of 15000 Quijote simulations at fiducial cosmology and with no PNG in their initial conditions as the surrogates, and a small set of 10 non-Gaussian simulations (f local NL = +100) with the same ΛCDM cosmological parameters as the high-fidelity simulations, the goal being to predict the power spectrum and bispectrum more accurately outside of fiducial cosmology. We compare the CARPool results to the 500 simulations with f local NL = +100 at our disposal and verify that they are indeed unbiased, as expected. We repeat the procedure to many 10 simulation subsets among the 500 to check that the result is not spurious, and to derive error bars on the CARPool averages. For all power spectrum and bispectrum modes, the error bars are significantly smaller than the standard errors on the average from 10 simulations alone. The effect is the strongest on linear scales (small k for the power spectrum, and the first few bispectrum modes which describes the tree-level matter bispectrum), but is also present in the non-linear regime.
In figure 10, we follow a similar procedure to study the derivatives of the power spectrum and bispectrum with respect to f local NL . We use the derivatives evaluated at f local NL = 0 by finite difference applied to the 500 f local NL = ±100 simulations as surrogates, to compute the derivatives at f local NL = 50 using only a few simulations with f local NL = 0 or 100. For the power spectrum the improvement is small, or even negligible in some cases, outside of the largest scales. For the bispectrum, there is a significant improvement of the first few modes describing the tree-level matter bispectrum. All error bars are reduced by the CARPool method, although the improvement is very small for some modal coefficients. One issue here is the small number of surrogates compared to the previous application (only 500 instead of 15000), adding to the fact that we know that the surrogate derivatives are not even fully converged numerically (as discussed thoroughly in section 4).
These examples illustrate briefly the possibilities of the CARPool technique. We leave its full implementation in the pipeline for future works, where we will also include the powerful "CARPool Bayes" introduced in Chartier & Wandelt (2022) for the fast and accurate estimation of the covariance matrix and its inverse.