Perturbative or effective field theory EFTbased fullshape analyses of galaxy clustering data involve nuisance parameters to capture various observational effects such as the galaxydark matter connection galaxy bias. We present an efficient approach to set informative physically motivated priors on these parameters. We extract these priors from simulated galaxy catalogs based on halo occupation distribution HOD models. First we build a joint distribution between EFT galaxy bias and HOD parameters from a set of 10500 HOD mock catalogs. We use the field level EFT technique that allows for cosmic variance cancellation enabling a precision calibration of EFT parameters from computationally inexpensive smallvolume simulations. Second we use neural density estimatorsnormalizing flowsto model the marginal probability density of the EFT parameters which can be used as a prior distribution in full shape analyses. As a first application we use our HODbased priors in a new analysis of galaxy power spectra and bispectra from the BOSS survey in the context of single field primordial nonGaussianity. We find that our priors lead to a reduction of the posterior volume of bias parameters by an order of magnitude. We also find inlineformulammlmath displayinlinemmlmsubsupmmlmifmmlmimmlmrowmmlmiNLmmlmimmlmrowmmlmtextequilmmlmtextmmlmsubsupmmlmommlmommlmn320mmlmnmmlmommlmommlmn300mmlmnmmlmathinlineformula and inlineformulammlmath displayinlinemmlmsubsupmmlmifmmlmimmlmrowmmlmiNLmmlmimmlmrowmmlmtextorthommlmtextmmlmsubsupmmlmommlmommlmn100mmlmnmmlmommlmommlmn130mmlmnmmlmathinlineformula at 68 CL in a combined twotemplate analysis representing a inlineformulammlmath displayinlinemmlmommlmommlmn40mmlmnmmlmommlmommlmathinlineformula improvement in constraints on single field primordial nonGaussianity equivalent to doubling the survey volume.
['Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Fullshape analysis with simulationbased priors Constraints on single field inflation from BOSS
{'Full-shape analysis with simulation-based priors: constraints on single field inflation from BOSS': "Mikhail M. Ivanov, 1, ∗ Carolina Cuesta-Lazaro, 2, 3, 4, † Siddharth  \nMishra-Sharma, 5, 1, 6, ‡ Andrej Obuljen, 7, § and Michael W. Toomey 1, ¶  \n1 Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 2 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, Cambridge, MA 02139, USA 3 Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n4 Center for Astrophysics - Harvard & Smithsonian,  \n60 Garden Street, MS-16, Cambridge, MA 02138, USA  \n5 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions  \n6 Department of Physics, Harvard University, Cambridge, MA 02138, USA  \n7 Department of Astrophysics, University of Zurich,  \nWinterthurerstrasse 190, 8057 Zurich, Switzerland  \nPerturbative, or effective field theory (EFT)-based, full-shape analyses of galaxy clustering data involve 'nuisance parameters' to capture various observational effects such as the galaxy-dark matter connection (galaxy bias). We present an efficient approach to set informative physically motivated priors on these parameters. We extract these priors from simulated galaxy catalogs based on halo occupation distribution (HOD) models. First, we build a joint distribution between EFT galaxy bias and HOD parameters from a set of 10,500 HOD mock catalogs. We use the field level EFT technique that allows for cosmic variance cancellation, enabling a precision calibration of EFT parameters from computationally inexpensive small-volume simulations. Second, we use neural density estimators - normalizing flows - to model the marginal probability density of the EFT parameters, which can be used as a prior distribution in full shape analyses. As a first application, we use our HOD-based priors in a new analysis of galaxy power spectra and bispectra from the BOSS survey in the context of single field primordial non-Gaussianity. We find that our priors lead to a reduction of the posterior volume of bias parameters by an order of magnitude. We also find f equil NL = 320 ± 300 and f ortho NL = 100 ± 130 (at 68% CL) in a combined two-template analysis, representing a ≈ 40% improvement in constraints on single field primordial non-Gaussianity, equivalent to doubling the survey volume.", '1. INTRODUCTION': "The last three decades of precision cosmological observations have yielded significant evidence for new physical phenomena in the form of dark matter, dark energy, and primordial accelerated expansion of the Universe, called cosmic inflation. Future progress in understanding of these phenomena will depend on our ability to extract cosmological information from upcoming large-scale structure surveys, such as DESI [1], Euclid [2], LSST [3], and Roman Space Telescope [4].  \nOn scales larger than ∼ 10 Mpc, large scale structure formation occurs in the mildly non-linear regime. In this regime the data can be described with perturbation theory (its consistent formulation is called effective field theory (EFT) of large-scale structure [5-7]), which provides a systematic framework for modeling galaxy clustering based only on symmetries and scale separation. Currently, this is the only approach to galaxy clustering that offers a sub-percent accuracy for galaxy separations larger than 10 Mpc [8, 9]. Another key advantages of EFT is its efficiency. Namely, EFT model templates can be computed in less than 1 second [10-12]. EFT thus offers an unmatched flexibility in analyzing cosmological models beyond the vanilla ΛCDM model (see e.g. [1214]). All these virtues made EFT a useful tool for the analysis of the galaxy clustering data on large scales, see e.g. [15-18] for a sample of analyses of BOSS data [19].  \nThe main disadvantages of EFT techniques are: (a) the breakdown on small scales (where the statistical power of data can be significant), and (b) the proliferation of free 'nuisance' parameters. These 'nuisance' parameters include, e.g. the classical perturbative bias parameters such as linear, quadratic, and tidal biases b 1 , b 2 , b G 2 etc [20]. These parameters encapsulate effects of smallscale galaxy formation physics and typically are determined from data in actual analyses. However, marginalization over EFT parameters within uninformative priors leads to a significant loss of information in current EFTbased full-shape analyses [21-23]. 1  \nThe loss of information due to free EFT parameters can be avoided by using informative priors. These priors can be extracted from phenomenological or empirical galaxy formation models. In the EFT context, they are UV-complete models available for matching calculations.  \nPhenomenological approaches such as the local Lagrangian approximation, co-evolution of galaxies and dark matter, the peak-background split etc. (see [20, 25] and references therein), predict certain relationships between galaxy bias parameters, which have been often used in actual data analyses, see e.g. [26]. However, accurate determination of these parameters from numerical simulations have shown that commonly used phenomenological analytic relations are not very accurate in practice, see e.g. [25, 27-29]. This suggests that their use in data analysis can bias cosmological parameter recovery. An alternative is to match the EFT parameters from empirical galaxy formation models based on hydrodynamical simulations or the halo occupation distribution (HOD) approach [30-33].  \nIn what follows we will focus on the HOD framework. This approach is motivated by the fact that galaxies fundamentally reside inside dark matter halos. Based on that, the HOD framework naturally assumes that key galaxies' properties are derived from those of the halos. Emulators based on HOD models have been successfully applied to data, see e.g. [34-38], which proves that they can reproduce the observed galaxy clustering with sufficient accuracy even though their robustness on small  \nscales remains to be determined.  \nThere have been notable efforts to determine EFT parameters from hydrodynamical simulations and HOD catalogs [25, 29, 39-41]. In particular, Ref. [25] measured non-linear bias parameters of 3 BOSS-like HOD galaxy samples and 4 halo catalogs using a combination of the power spectrum and bispectrum. Ref. [29] extracted second order galaxy bias parameters from ≈ 30 different samples of galaxies from IllustrisTNG hydrodynamical simulations. The relatively small size of data points in both of these analyses does not allow one to robustly explore the EFT parameter distribution, which is expected to have a complicated correlation structure. In particular, one can estimate that building an accurate distribution for 14 dimensional parameter space of HOD and EFT parameters requires a total of 10 4 samples. 2 Previously, [41] measured Lagrangian bias parameters of a hybrid EFT model for ∼ 8000 galaxy and halo samples. Their results, however, cannot be directly applied to pipelines based on the traditional EFT models.  \nIn this paper, we present a new approach for precise determination of priors for EFT-based full shape analyses. The key object of our study is a precise map between EFT parameters and HOD parameters for BOSSlike galaxy catalogs. We build this map from 10,500 mocks, which represents a significant improvement in size over previous EFT measurements. This map can be used in multiple ways. For instance, we consider joint and conditional distributions of EFT and HOD parameters, p ( θ EFT , θ HOD ) and p ( θ EFT | θ HOD ), respectively, which clearly display the response of galaxy bias parameters to variations of HOD parameters. These distributions give us new insights into the physical meaning of galaxy bias parameters and help us connect galaxy bias models on small and large scales.  \nFinally, we propose to use the marginal density of the EFT parameters from our samples p ( θ EFT ) as a prior distribution in EFT-based full-shape analyses. Similar ideas have been previously presented in [42] in the context of the Halo Zel'dovich perturbaiton theory model. We demonstrate the power of this approach in a new analysis of the BOSS data in the context of inflation-  \nry models with non-local single-field primordial nonGaussianity (PNG) captured by the equilateral and orthogonal templates. This type of non-Gaussianity probes inflaton self-interactions and the propagation speed [4350]. In the past, it has been shown that constraints on single-field PNG from large-scale structure depend significantly on the assumptions about non-linear bias parameters [22, 51-54] (see also [52, 55-60] for recent related work and discussions in the context of multi-field inflation). We find that with our HOD-based priors the constraints improve by ≈ 40%, analogous to a twofold increase of the survey volume. The volume of the posterior distribution of non-linear galaxy bias parameters for each independent chunk of BOSS data shrinks by an order of magnitude.  \nOur paper is structured as follows. We describe the technical aspects of how to produce a map of EFT and HOD parameters in Sec. 2. Sec. 3 presents the map and selected conditional distributions. Our PNG constraints on BOSS with HOD-informed priors are presented in Sec. 4. Sec. 5 compares the optimal values of the HOD parameters implied by our EFT full shape analysis to those based on other techniques involving small scale data. Finally, Sec. 6 draws conclusions and discusses directions of further work. Additional plots are presented in Appendix A. Appendix B summarizes details of our normalizing flow training, while Appendix C studies the residual UV-dependence of our results.", '2. THE METHOD': "The main idea of our method is to create a detailed mapping (joint probability density) between EFT and HOD parameters. From the physical point of view, it will be interesting to consider a conditional version of this mapping, p ( θ EFT | θ HOD ), which will describe the dependence of the EFT parameters on the underlying halo model physics. In the context of EFT-based full shape analyses, one is interested in the marginal distribution of the EFT parameters 'informed' by the HOD models,  \np ( θ EFT ) = ∫ p ( θ EFT | θ HOD ) p ( θ HOD ) dθ HOD , (1)  \nwhere p ( θ HOD ) is the prior distribution of the HOD parameters. Creating p ( θ EFT ) is the main practical goal of our paper.  \nTo build a joint EFT/HOD distribution, we measure EFT parameters from a large set of simulations with varying HOD parameters. This poses several challenges. The first one is cosmic variance. A traditional way to determine the EFT parameters is to extract them by fitting a combination of simulated correlation functions, typically the power spectrum and bispectrum, see e.g. [61]. In this approach, the fits have to be performed on large scales where EFT is applicable. But the large scales are also most affected by cosmic variance. As a consequence, a high precision measurement of the EFT parameters requires computationally expensive simulations with large volume. For example, a measurement of the quadratic bias parameters with absolute errorbars σ ≃ 0 . 02 [23, 28] requires a cumulative simulation volume of 566 h -3 Gpc 3 . We solve this problem by employing a field-level EFT [27, 62-75] that allows for cosmic variance cancellation, and thus enables precision calibration of EFT parameters from computationally cheap smallvolume simulations. The field-level EFT also provides a computationally inexpensive way to take into account information beyond the two-point function, which is important in order to break degeneracies between EFT parameters.  \nIn what follows we describe in detail the creation of HOD catalogs and field level EFT fits.", '2.1. HOD mocks': 'We build a large set of HOD mocks for BOSS-like galaxies based on the AbacusSummit [76] suite of simulations. For the purposes of this work, we use a set of mocks with an underlying Planck 2018 baseline cosmological model [77]. In principle, the EFT parameters depend on cosmology, so a full analysis of the BOSS data with variations of cosmological parameters will require a distribution that samples cosmological parameters. This will be presented in future work.  \nWe produce a suite of 10,500 mock catalogs of galaxies similar to luminous red galaxies of BOSS. Specifically, we generate galaxies whose linear bias is similar to that of BOSS galaxies (around ≈ 2), and number density is less or equal to that of the CMASS sample of BOSS, ¯ n ≈ 3 . 6 · 10 -4 h 3 Mpc -3 . Note that the HOD-to-EFT mapping will have to be calibrated anew for galaxy samples different from that of BOSS, e.g. for emission line  \ngalaxies of DESI.  \nWe use available N -body simulation boxes with periodic boundary conditions from the small boxes AbacusSummit covariance suite. This suite has a total of 1883 independent realizations of dark matter initial conditions. Each box has a site length 500 h -1 Mpc. Our mocks are produced from snapshots at z = 0 . 5, which matches the typical redshift of the CMASS-type galaxies of BOSS. The mocks are fitted with a resolution of 256 pixels and have been corrected for the cloud-in-cell (CIC) window. Thanks to the small volume of our simulations, we can generate galaxy catalogs and perform EFT fits in a very short time, ∼ 20 and ∼ 40 sec., respectively.  \nEach N -body simulation is populated with HOD galaxies [32]. We sample HOD parameters { log M cut , log M 1 , log σ , α , κ , B cen , B sat } from uniform uninformative priors specified in [35]. These are:  \nlog M cut ∈ [12 . 4 , 13 . 3] , log M 1 ∈ [13 . 2 , 14 . 4] , log σ ∈ [ -3 . 0 , 0 . 0] , α ∈ [0 . 7 , 1 . 5] , κ ∈ [0 . 0 , 1 . 5] , B cen ∈ [ -0 . 5 , 0 . 5] , B sat ∈ [ -1 . 0 , 1 . 0] . (2)  \nIn the HOD model, average numbers of central and satellite galaxies are given by  \n⟨ N c ⟩ ( M ) = 1 2 [ 1 + Erf ( log M -log M cut √ 2 σ )] , ⟨ N s ⟩ ( M ) = ⟨ N c ⟩ ( M ) ( M -κM cut M 1 ) α , (3)  \nwhere Erf( x ) is the Gauss error function, M cut is the minimum mass of a halo that can host a galaxy, M 1 is the typical (i.e. most probable) halo mass that hosts one satellite galaxy, 3 κM cut is the minimum mass of a halo that can host a satellite, α is the slope of the satellite probability distribution function.  \nWe focus on PNG analysis in this work. As mentioned in the Introduction, the only relevant parameters in these case are non-linear bias parameters b 2 , b G 2 which can be extracted from real space mocks. Hence, for our purpose, it is sufficient to consider real space clustering only, and hence we ignore the velocity bias parameters that are relevant in redshift space.  \nIn addition, we use B cen , B sat to capture assembly bias for satellites and centrals [78-80]. Specifically, the galaxy  \noccupation enhancement is captured by promoting M cut and M 1 to functions of δ R via  \nlog M cut → log M cut + B cen ( δ R -0 . 5) , log M 1 → log M 1 + B sat ( δ R -0 . 5) , (4)  \nwhere δ R is the (smoothed) dark matter overdensity around halo centers ranked within halo mass bin and normalized to range from 0 to 1, and we use a top-hat filter with radius R = 5 h -1 Mpc for smoothing, following [78]. Note that the main halo itself is excluded when computing δ R , i.e only neighbor halo contributions are taken into account. Positive (negative) values of B sat / cen mean that galaxies of the relevant type preferably form in more (less) dense environments.', '2.2. Field-level EFT': "The standard Eulerian EFT bias model relevant for the one-loop EFT-based full shape analyses is given by [15]  \nδ g = b 1 δ + b 2 2 ( δ 2 -σ 2 )+ b G 2 G 2 + b Γ 3 Γ 3 -b ∇ 2 δ ∇ 2 δ + ε , (5)  \nwhere and ε is the stochastic density component, δ is the nonlinear dark matter field, ∇ ≡ ∂ i ∂ i , the 'non-local' bias operators are given by [81]  \nG 2 = ( ∂ i ∂ j ∇ 2 δ ) 2 -δ 2 , Γ 3 = 4 7 δ G 2 -4 7 ∂ i ∂ j δ ∇ 2 ∂ i ∂ j G 2 ∇ 2 , (6)  \nand σ 2 ≡ ⟨ δ 2 ⟩ , enforcing non-renormalization of the background density. The direct use of the above model in field-level EFT, however, leads to a poor match to simulations because the fully non-linear density field δ above has a strong UV sensitivity [27, 65]. This sensitivity also makes it hard to compare field level results with those based on traditional correlation functions because the latter feature renormalized biased parameters as opposed to 'bare' ones in the field-level calculations. In order to compare field-level and correlation function EFT calculations in a more consistent manner, one would want to evaluate the Eulerian bias model (5) with perturbative matter fields, as is done in the Eulerian EFT loop expansion [7, 82]. The perturbative fields, however, miss large contributions from enhanced Zel'dovich displacements [82, 83], which again, leads to a failure of  \nFIG. 1. A typical HOD mock galaxy distribution from our set (left), field-level EFT fit to it (center), and the residuals (right). The overdensity field has been smoothed with a R = 4 h -1 Mpc 3D Gaussian filter and the depth of each panel is ≈ 60 h -1 Mpc.  \n<!-- image -->  \nnaive field-level EFT models to predict the galaxy density field. A way around is to use shifted operators proposed in [65, 71, 84], which have a well controlled smallscale behavior and retain large displacements at the same time. A shifted operator ˜ O is obtained by shifting the Lagrangian operator O ( q ) by the Zel'dovich displacement ψ 1 :  \n˜ O ( k ) = ∫ d 3 q O ( q ) e -i k · ( q + ψ 1 ( q )) , (7)  \nwhere q are Lagrangian coordinates. Note that in this approach higher order Lagrangian displacements are treated perturbatively, which is appropriate given that they are suppressed just like other operators in perturbation theory [20]. The statistics of shifted operators can be shown to be equivalent to IR-resummed Eulerian EFT [65] (see [85-89] for details of IR resummation).  \nIn the context of real-space clustering that we study here, the field level EFT forward model based on shifted operators takes the form  \nδ g ∣ ∣ EFT = β 1 ˜ δ 1 + β 2 ( ˜ δ 2 1 ) ⊥ + β G 2 ˜ G ⊥ 2 + β 3 ( ˜ δ 3 1 ) ⊥ , (8)  \nwhere ˜ δ 1 is the shifted linear density field δ 1 , β n are scaledependent transfer functions. Importantly, the forward model (8) is built from orthogonalized operators that satisfy  \n⟨ ˜ O ⊥ m ˜ O ⊥ n ⟩ = 0 , for m = n, (9)  \n̸  \nwhich allows us to remove UV-sensitive two-loop corrections to the transfer functions [27].  \nNote that our model contains an extra cubic operator δ 3 that has a first non-trivial contribution in the bispectrum at the one loop order [23, 90, 91]. We extract the k -dependent transfer functions from each mock and then fit them with the appropriate EFT templates on mildlynonlinear scales. The transfer functions are computed by taking expectation values of the simulated galaxy density field multipled by appropriate shifted operators. Specifically, the transfer function β i of a Fourier space bias operator ˜ O i ( k ) is computed as  \nβ i ( k ) = ⟨ δ HOD g ( k ) ˜ O ⊥ i ∗ ( k ) ⟩ ⟨| ˜ O ⊥ i ( k ) | 2 ⟩ , (10)  \nwhere δ HOD g is the galaxy density field from simulations. The use of expectation values allows us to connect the transfer functions with renormalized EFT parameters that appear, e.g. in separate universe simulations [92] or in correlation functions. In particular, at the formal one-loop order, we have the following expressions in the k → 0 limit 4  \nβ 1 = b 1 + b ∇ 2 δ k 2 + ( b Γ 3 + b 1 6 + 5 2 b G 2 ) ⟨ ˜ δ 1 ˜ Γ 3 ⟩ ⟨ ˜ δ 1 ˜ δ 1 ⟩ + b 2 2 ⟨ ˜ δ 1 ˜ δ 2 1 ⟩ ⟨ ˜ δ 1 ˜ δ 1 ⟩ + ( b G 2 + 2 b 1 7 ) ⟨ ˜ δ 1 ˜ G 2 ⟩ ⟨ ˜ δ 1 ˜ δ 1 ⟩ -b 1 ⟨ ˜ δ 1 ˜ S 3 ⟩ ⟨ ˜ δ 1 ˜ δ 1 ⟩ , β 2 = b 2 2 , β G 2 = b G 2 + 2 7 b 1 , β 3 = b 3 6 . (11)  \nNote that the constant coefficients here are different from those of [65] because we use the Eulerian bias parameters from Eq. (6) (matched at the cubic order), while the parameters used in [65] are more closely related to the Lagrangian bias parameters. The cubic operator ˜ S 3 is the shifted version of  \nS 3 = ψ 2 ( q ) · ∇ δ 1 ( q ) , (12)  \nwhich is produced by the second order displacement ψ 2 . Note that the presence of this operator is enforced by the equivalence principle.  \nIn practice, we use public Hi-Fi mocks 5 to produce the EFT forward model. A typical snapshot generated with the field-level EFT forward model and its residual with the original simulated galaxy field are shown in Fig. 1.  \nWe extract the bias parameters from the k → 0 limit of the transfer functions using Eq. (11). In practice we use k max = 0 . 4 h Mpc -1 , for which the one-loop EFT models are reliable [93, 94]. We have checked that a more conservative choice of k max = 0 . 3 h Mpc -1 gives results consistent with those of k max = 0 . 4 h Mpc -1 , but with a somewhat larger scatter in their distribution. This is the reason we adopt k max = 0 . 4 h Mpc -1 as a baseline choice.  \nTo account for the scatter in transfer functions on large scales, we adopt error weights for k -bins based on the number of Fourier modes. For b 2 , b G 2 , b 3 , we fit the transfer functions with a polynomial c 0 + c 2 k 2 + c 4 k 4 , and then match c 0 to the constant values of bias parameters in Eq. (11). As far as β 1 is concerned, we calculate the power spectra and cross-spectra of shifted operators in Eq. (11), and use them, along with the best-fit values of b 2 and b G 2 from the previous step, to fit the transfer function at low k , which yields b ∇ 2 δ and b Γ 3 . Plots with typical transfer functions fits can be found in Appendix A.  \nThe final ingredient that we need is the distribution of stochasticity parameters, characterizing the power spectrum of the ε field. In practice, we calculate the error power spectrum as  \nP err ( k ) = ⟨| δ HOD g ( k ) -δ EFT g ( k ) | 2 ⟩ . (13)  \nTheoretical consistency dictates that on large scales it  \nshould match the EFT prediction for the stochastic contribution to the galaxy power spectrum [9, 28, 65]  \nP err ( k ) = 1 ¯ n ( 1 + α 0 + α 1 ( k k NL ) 2 ) , (14)  \nwhere ¯ n = V/N gal is the number density of mock galaxies ( V being the simulation box volume), and we chose k NL = 0 . 45 h Mpc -1 following [28]. Note that the EFT fitting pipelines (e.g. [9, 28]) use parameters P shot and a 0 that are similar to our α 0 and α 1 , respectively. However, as we discuss in detail later, these 'standard' EFT models also absorb additional contributions into the stochastic parameters, which make them somewhat different from our α 0 and α 1 .  \nBefore closing this part, let us comment on errors in our measurements of bias parameters. There are two main sources of errors: the residual UV dependence of the transfer functions on large scales and the numerical noise due to the use of a single realization per mock in our fits. The use of orthogonalized shifted operators and k -dependent transfer functions allowed us to remove the bulk of the UV sensitivity of our bias parameters. Within the original field-level method for cosmic variance cancellation, see e.g. [27], the bias parameters stem from a global fit of a constant to many k -bins, which is sensitive to small scale nonlinearity and higher-order loop effects. In order to reduce this sensitivity, one may introduce a relatively aggressive cutoff for the linear fields ([27] used a Gaussian filter with R s = 20 h -1 Mpc). This cutoff generates a mismatch between the field-level and n -point functions results, as the latter are formally extracted in the R s → 0 limit. 6 This is because many EFT codes, e.g. CLASS-PT [10] use dimensional regularization for loop integrals that compute correlation functions, in which all convergent integrals are done without any cutoff, while the divergent pieces are set to zero. In the field level language this means the matching should be done after the orthogonalization (which removes divergent corrections proportional to the mass variance integrals), and without smoothing of the linear fields (i.e. at R s → 0).  \nIn this work, however, we follow the method [65] and extract the bias parameters from the lowk limits of the  \nFIG. 2. The joint distribution of EFT and HOD parameters extracted from 10,500 HOD mocks for BOSS-like galaxies. Density levels correspond to two-dimensional 1σ and 2σ intervals (i.e. 39.3% and 86.5% of samples). Individual samples are also shown as dots. They are especially pronounced in the tails.  \n<!-- image -->  \nFIG. 3. Conditional distributions of EFT parameters for different values of log M cut (see legend). Other HOD parameters are kept fixed. Density levels correspond to two-dimensional 1σ and 2σ intervals.  \n<!-- image -->  \nG  \n∇  \ntransfer functions, in which case the UV-sensitive higherorder corrections are absorbed into the k 2 and k 4 polynomials of our fit, leaving the constant part unaffected provided one uses a reasonably high cutoff. As an explicit check, we demonstrate that the lowk transfer functions used in our fits are stable w.r.t. variations of the grid resolution, which provides an effective cutoff for all fields, see Appendix C for more detail. This test shows that our results have converged w.r.t. the R s → 0 limit relevant  \nfor the matching to the n -point functions in dimensional regularization.  \nAs far as numerical noise is concerned, results of Ref. [65] obtained with a similar box size suggest that the corresponding error is below the level of the scatter induced by variations of HOD models, implying that  \n15  \nFIG. 4. Same as Fig. 3, but B sat is varied (see legend) while other parameters are kept fixed.  \n<!-- image -->  \nthese effects are negligible for parameter constraints. 7", '2.3. Normalizing flow for density estimation and sampling': "The crux of our approach is the ability to effectively model the marginal and HOD-conditional EFT parameter distributions, p ( θ EFT ) and p ( θ EFT | θ HOD ) respectively, given a set of samples { θ EFT , θ HOD } ∼ p ( θ EFT , θ HOD ) as described above. Ideally, we would like to be able to generate new samples from these distributions (i.e., use them as priors), and evaluate the likeli-  \nnew set of samples under the modeled distributions.  \nGenerative models are a class of statistical models that aim to encode potentially complex target distributions. While high-dimensional distribution modeling is an inherently challenging task, machine learning has brought this into the realm of tractability. Normalizing flows [95], a class of deep generative models, are especially suited for our purposes, as they allow for seamless density estimation as well as sampling, including in the conditional regime. Briefly, normalizing flows model the target density ˆ p ( θ ) as a bijective (invertible) learnable function θ = f φ ( u ) from a simple distribution π ( u ), e.g. a multivariate Gaussian,  \nˆ p ( θ ) = π ( u ) ∣ ∣ ∣ ∣ det ( ∂u ∂θ )∣ ∣ ∣ ∣ = π ( f -1 φ ( θ ) ) ∣ ∣ ∣ det J f -1 φ ( θ ) ∣ ∣ ∣ (15)  \nwhere ∣ ∣ ∣ det J f -1 φ ( θ ) ∣ ∣ ∣ is the Jacobian determinant of the inverse transformation and is by construction easy to compute. φ are the parameters of the learnable transformation, usually modeled through an appropriate invertible neural network.  \nGiven a set of samples { θ } ∼ p ( θ ) (e.g, HOD and EFT parameters), we can maximize Eq. (15) over those samples to train the flow, and thus build an approximation for our target density. The optimal parameters φ ∗ of the transformation are obtained, in practice through stochastic gradient descent or a variant thereof, as  \nφ ∗ = arg max φ ⟨ log ˆ p ( θ ) ⟩ θ ∼ p ( θ ) . (16)  \nThe target density can then easily be sampled from, by drawing u ∼ π ( u ) from the simple (Gaussian) base density and running the forward transformation f φ ( u ), as well as evaluated for a new sample θ ' . The flows are implemented using the nflows 8 library, with training and evaluation performed using PyTorch [96]. Further details of the normalizing flow implementation, training, and validation are given in Appendix B.", '3. DISTRIBUTIONS OF GALAXY BIAS AND HOD PARAMETERS': 'The joint samples of EFT parameters for given HOD models are displayed in Fig. 2. Numerical values of b ∇ 2 δ  \nare given in [ h -1 Mpc] 2 units. We do not show the samples of the HOD parameters as they simply scan over their uniform priors. We see that the EFT parameters are strongly correlated among each other. This is the behavior expected on the basis that the (infinite) entire set of EFT parameters must be produced by only 7 HOD parameters. (In the particle physics context, the situation here is analogous to the textbook matching of chiral perturbation theory to the linear sigma model, where low energy constants are either zero, or obey correlations set by a few parameters of the linear sigma model.) The tightest correlation is between b Γ 3 and b G 2 , which follows the linear law  \nb Γ 3 ≈ -3 . 8 b G 2 -0 . 5 . (17)  \nNote that this law is steeper than the local Lagrangian model prediction b Γ 3 = -23 12 b G 2 , and the co-evolution model relation b Γ 3 ∼ -1 . 5 b G 2 [25]. The strong correlation between b G 2 and b Γ 3 may be an artifact of our forward model, which absorbs the contribution of the Γ 3 operator into β 1 , thus making our measurement of Γ 3 sensitive to details of the fitting procedure. It will be interesting to compare our results with independent measurements of b Γ 3 from a forward model with the shifted Γ 3 operator.  \nAnother important observation is that the higherderivative stochastic counterterm α 1 is very close to zero for most of the samples. A similar pattern was discovered before for the clustering of dark matter halos [65]. It would be interesting to see if this pattern is specific to HOD models, i.e. if other galaxy formation prescriptions can produce larger α 1 matching the natural EFT expectation that this parameter should be proportional to the Lagrangian radius of the host halo. This expectation is supported by the positive correlation between α 1 and b 1 . Note that α 1 is mostly correlated with M cut , M 1 and the satellite environmental assembly bias B sat . These correlations can be understood as follows: non-locality in the galaxy stochasticity is sourced by the stochastic non-locality of the underlying halo (set by M cut ) plus the stochasticity from the satellites, whose spatial distribution increases the non-locality scale, thus generating additional M 1 and B sat dependencies.  \nLet us discuss now the correlation between P shot (deviation of white noise amplitude from ¯ n -1 ) and b 1 . Due to halo exclusion effects [97-100], we expect that samples with large enough b 1 should feature sub-Poisson stochastic shot noise, i.e. have negative α 0 . This is in-  \ned the case, on average. However, massive halos also tend to have more satellites, which would increase ¯ n , and hence compensate for the sub-Poisson trend produced by the halos. This is also is supported by correlations between α 0 and satellite properties such as B sat . Indeed, the observation that satellite galaxy assembly bias produces super-Poissonian sampling was pointed out before in [101]. This is the reason why we do not see a strong sub-Poisson constant shot noise contribution for galaxies that reside in massive halos. A similar trend was pointed out in a previous work on the field-level calibration of α 0 from HOD [40], with which we find excellent agreement.  \nNote that other prominent correlations are b 1 -B cen and b ∇ 2 δ -B sat , which suggest that the measurements of these parameters can be used to diagnose the presence of environmental assembly bias.  \nIn order to visualize the response of the bias parameters to variations of the HOD parameters, we study the conditional distribution p ( θ EFT | θ HOD ). It indicates what are the most likely EFT parameters for a given fixed set of HOD parameters. This distribution is extracted by modeling the original samples with normalizing flows as described in Sec. 2.3.  \nThe fiducial conditional distribution of EFT parameters that we study is calculated for the set of HOD parameters fixed to the best-fit values of the CMASS BOSS sample [35],  \nlog M cut = 12 . 66 , log M 1 = 13 . 66 , α = 1 . 34 , log σ = -0 . 5 , κ = 0 . 03 , B cen = -0 . 43 , B sat = -0 . 22 . (18)  \nAs a first example, we plot p ( θ EFT | θ HOD ) for 7 HOD models with different log M cut and all other parameters fixed. The result is shown in Fig. 3. As anticipated from the simple thresholding picture [20], the higher cutoff mass of the halo implies a lower probability of halo formation and hence higher bias of the host halos and galaxies. Note that the higher-derivative bias b ∇ 2 δ is not very sensitive to the cutoff mass, which is at odds with the naive expectation that this parameter should be proportional to the Lagrangian radius of the host halo. It would be interesting to understand the origin of this behavior. This intuition is, however, confirmed by the higher derivative stochastic counterterm α 1 which grows with cutoff mass.  \nHaving confirmed that our conditional distribution reflects the basic intuition about the galaxy bias, let us  \nconsider now the response of EFT coefficients to a less intuitive parameter, the environment-based assembly (secondary) bias of satellites B sat . The corresponding conditional distributions for 7 samples of galaxies with different B sat (other HOD parameters are fixed) is shown in Fig. 4. The first relevant observation is that local-inmatter-density bias parameters b 1 , b 2 , b 3 display only a weak dependence on B sat . This is consistent with the expectation that these parameters are mostly determined by halo mass. The small residual dependence on B sat can be explained as the tendency to have more/less galaxies for positive/negative B sat , which leads to lower/higher values of local biases. The presence of satellite assembly bias is best reflected by non-local bias parameters b G 2 , b Γ 3 , b ∇ 2 δ and higher-derivative stochastic counterterm α 1 , which is consistent with the fact that assembly bias is an intrinsically non-local property of galaxies. Another interesting observation is that a strong satellite assembly bias leads to super-Poisson shot noise (positive P shot ) regardless of the sign of B sat .  \nThe distributions of EFT parameters conditioned to other HOD models have a similar behavior that can be qualitatively understood on the basis of simple physical arguments such as thresholding.', '4. PNG CONSTRAINTS FROM BOSS WITH HOD-INFORMED PRIORS': "In this section, we re-analyze the BOSS data within initial conditions with non-local primordial nonGaussianity. In this analysis, following the standard practice [51, 52, 102], we fix the 'standard' cosmological parameters to the Planck best-fit values, and vary only the EFT and PNG parameters. Single-field PNG is captured by the equilateral and orthogonal templates identical to the ones used in [51]. Note that we use the orthogonal template that has the correct physical scaling in the squeezed limit.  \nOur dataset is BOSS DR12 galaxy clustering in redshift space [103]. The BOSS DR12 data is split in four chunks: low-z ( z eff = 0 . 38) and high-z ( z eff = 0 . 61), South and North Galactic Caps. For each chunk, our datavector is { P 0 , P 2 , P 4 , Q 0 , B 0 , B 2 , B 4 , α ∥ , α ⊥ } , where P ℓ , B ℓ are window-free power spectrum and bispectrum multipoles, ℓ = 0 , 2 , 4 [14, 17, 24], Q 0 is the realspace power spectrum proxy [94], while α ∥ , α ⊥ are the  \nFIG. 5. Corner plots with 2d and 1d marginalized posterior distribution of the PNG parameters from the full BOSS DR12 dataset and the galaxy bias parameters from BOSS NGC high-z sample. HOD priors on bias parameters are shown in gray. Bias parameters for other samples are shown in Appendix. Density levels correspond to 68% and 95% CL.  \n<!-- image -->  \nG  \npost-reconstructed BAO parameters [104]. We use scale cuts k max P ℓ = 0 . 2 h Mpc -1 , k max B ℓ = 0 . 08 h Mpc -1 , k max Q 0 = 0 . 4 h Mpc -1 and k min = 0 . 01 h Mpc -1 . The covariance matrix is estimated from Multidark Patchy mocks [105, 106]. We assume that EFT parameters in four chunks are independent of each other. Other aspects of our analysis are identical to those of [22, 24].  \nWe first run the usual analysis with the conservative priors on galaxy bias parameters similar to [22]. For ref-  \nerence, the priors on the galaxy bias parameter are  \nb 1 ∈ [1 , 4] , b 2 ∼ N (0 , 1) , b G 2 ∼ N (0 , 1) , b Γ 3 ∼ N ( 23 42 ( b 1 -1) , 1 ) , (19)  \nwhere N ( a, σ ) denotes a normal distribution with mean a and variance σ 2 .  \nAs a second step, we replace the conservative priors from the previous analysis with the marginal likelihood p ( θ EFT ), see Eq. (1). We use p ( θ EFT ) to set priors for pa-", 'BOSS PNG with HOD-informed priors': 'TABLE I. Best-fits and 1d marginalized limits for PNG and galaxy bias parameters from BOSS with HOD-informed priors (left panel) and the usual conservative priors (right panel). Upper scripts (1 , 2 , 3 , 4) refer to NGCz3, SGCz3, NGCz1, and SGCz1 samples (z1= 0 . 38, z3= 0 . 61).', 'BOSS PNG with conservative priors': 'rameters θ EFT = { b 1 , b 2 , b G 2 , b Γ 3 } of the BOSS EFT likelihood. The other EFT parameters from our samples, i.e. { b 3 , b ∇ 2 δ , α 0 , α 1 } are not included in the distribution, i.e. are effectively marginalized over, because of the following reasons.  \nParameters P shot and a 0 in the BOSS EFT likelihood are quite different from α 0 , α 1 that we have extracted from the HOD mocks. First, P shot in the EFT models currently implemented in CLASS-PT [10] absorbs the constant deterministic contribution from the one-loop autospectrum of δ 2 , see [107] for a recent discussion. This is the first reason why P shot is different from our α shot , which corresponds to a truly stochastic component. Second, our α shot captures departures from the Poissonian shot noise in a periodic box geometry, which is different from the number density of actual galaxies in a survey that is subject to additional weights, e.g. FKP (see e.g. a discussion in [9]). Third, stochastic counterterms  \nin the BOSS EFT pipeline additionally absorb contributions from fiber collisions [9, 108, 109]. Given these reasons, we keep using uninformative priors on P shot and a 0 in our analysis. We also do not use a prior on b 3 as this parameter does not appear in our current analysis based on the tree-level bispectrum. Finally, we do not use the prior on b ∇ 2 δ in the current analysis because the contribution of the ∇ 2 δ operator to the galaxy power spectrum in redshift space is degenerate with other redshift space counterterms that we have not studied yet. Their precision measurement with the redshift-space field-level EFT of [71] is left for future work. Since the PNG parameters are primarily degenerate with the quadratic bias parameters b 2 and b G 2 , we do not expect our results to depend significantly on the RSD counterterms and parameters { P shot , a 0 } discussed above. For this reason, we adopt standard conservative priors for them in all analyses presented here.  \nAdditionally, we have marginalized over the nonGaussian bias parameter b ζ within conservative priors of [51] in both analyses. In principle, one should extract priors on these parameters from mocks too. However, marginalization over b ζ does not have a noticeable impact on the single-field PNG constraints presented here. We leave a calibration of these parameters from simulations for future work.  \nA comment is in order on the redshift-dependence of the EFT parameters. We have calibrated priors at z = 0 . 5, which is somewhat different from z eff ≈ 0 . 4 and 0 . 6 of BOSS chunks. Within the HOD models, the galaxy distribution is determined by local-in-time properties of halos. Therefore, a redshift dependence of bias parameters can be fully compensated by a change of HOD parameters such as M cut , which are varied in our samples. Physically, we expect some residual sensitivity to the past evolution, e.g. to merger and assembly histories, which is captured within EFT by the expansion along the past fluid trajectory [20, 110]. These effects appear only at the two-loop order in the EFT, and hence are irrelevant for our analysis based on the one-loop calculation. All in all, the redshift dependence of the EFT parameters is adequately captured by variations of HOD parameters at the current precision level.  \nOur main results are displayed in figure 5 and in table I. A corner plot with EFT parameters for each chunk can be found in Appendix A. The first relevant observation is that the posteriors for non-linear galaxy bias parameters shrink significantly after applying the HOD priors. In particular, the 1d marginalized errorbars on b 2 , b G 2 and b Γ 3 shrink by a factor of few. To quantify the improvement in a more rigorous manner, we introduce a figure of merit (FoM) for bias parameters, defined along the lines of the FoM for the dark energy equation of state. Namely, for each individual BOSS data slice we write  \nFoM bias = [det C ( b 1 , b 2 , b G 2 , b Γ 3 )] -1 / 2 , (20)  \nwhere C is the covariance matrix of the bias parameters after marginalizing over other parameters in the chain. In the Gaussian case, our FoM is proportional to the inverse volume of the 4-dimensional ellipsoid enclosing the posterior distribution of the bias coefficients. A relevant parameter in our comparison is the FoM ratio between the old and new bias parameter measurements, e.g. for  \nthe NGCz3 we have  \nFoM bias | HODpriors FoM bias ∣ ∣ cons . priors ∣ ∣ ∣ ∣ ∣ NGCz3 = 60 . 3 , (21)  \nwhich can be interpreted as a factor of ≈ 60 reduction in cumulative inverse variance of the posterior bias parameter distribution. The cumulative figure of merit of bias parameters from all four BOSS data chunks is 1 . 2 · 10 7 .  \nThe second important observation is that our HODbased measurements of bias parameters are in perfect agreement with the EFT results based on the noninformative priors: the HOD-based posteriors are located almost at the centers of posterior densities from the standard EFT analysis. This is a non-trivial consistency check of both our prior calibration technique and the BOSS EFT pipeline since our measurements of nonlinear bias parameters are dominated by the galaxy bispectrum data 9 , while the HOD-informed priors were extracted from the field-level fits to HOD mocks. One can also note that the constraints on the linear bias parameter b 1 do not improve. This is because the data itself is much more informative than the prior for this parameter, cf. the standard BOSS posterior and the HOD prior in Fig. 5. Indeed, our HOD samples are consistent with a wide range of b 1 , while the actual measurements from BOSS are at the level of few percent.  \nThe third important observation is that the 1d marginalized errorbar on f equil NL has narrowed by ≈ 40%, from f equil NL = 806 +490 -510 (conservative priors) to f equil NL = 0 . 323 +0 . 29 -0 . 31 (HOD-informed). The main channel of improvement here is the breaking of degeneracy between f equil NL and b G 2 in the galaxy bispectrum. Since f ortho NL was not significantly correlated with bias parameters to begin with, the improvement for this parameter is only at ≈ 10%, f equil NL = -8 . 2 +130 -150 (conservative priors) vs. f equil NL = 100 +130 -140 (HOD-informed). In line with the previous discussion, let us introduce the FoM for non-Gaussian amplitudes  \nFoM PNG = [ det C ( f equil NL , f ortho NL ) ] -1 / 2 , (22)  \nwhich calculates the inverse area (up to a factor of 1 /π ) of the error ellipse in the ( f equil NL , f ortho NL ) plane. Then we  \nhave  \nFoM PNG | HODpriors FoM PNG ∣ ∣ cons . priors = 1 . 76 , (23)  \nor a ≈ 40% reduction of the total 2-dimensional posterior variance of PNG parameters.  \nFinally, let us mention that the analysis with the conservative priors implies a marginal 2 σ preference for nonzero f equil NL , which cannot be physical as the central value f equil NL = 800 is in strong tension with the CMB measurements [102]. The preference for non-zero f equil NL reduces in our new analysis, which highlights the importance of physically-motivated priors for the full-shape analysis.', '5. COMPARISON WITH SMALL-SCALE MEASUREMENTS': 'An interesting application of the method presented here is the comparison between simulation-based models fitting clustering data on small scales, see e.g. [35-38], and perturbation theory based models constrained by large scales. In particular, we use the conditional model p ( θ EFT | θ HOD ) to convert the HOD posterior chains from [35] into EFT bias parameters. These posteriors were obtained by fitting the two-point correlation function and density split statistics of the CMASS galaxy sample in the scale range 1 < r < 150 h -1 Mpc. The model used for inference is a neural network trained to reproduce the clustering of HOD-based galaxy mocks.  \nIn Figure 6, we compare the estimated bias parameters, denoted as HOD fits, with those found in this paper, denoted as EFT fits. Interestingly, although the range of scales used is widely different, both posteriors agree remarkably. Note that some important differences beyond scale ranges are the parameters that are being fitted, since [35] also fits the cosmological parameters, and the statistics used.', '6. DISCUSSION AND OUTLOOK': "We have presented a new framework for deriving informative priors on EFT 'nuisance' parameters from galaxy formation simulations. For concreteness, we focused here on the HOD models, although our approach can be straightforwardly applied to other galaxy formation models. The central object of our method is the  \nconditional distribution between parameters of galaxy formation models and EFT parameters. We build this distribution by generating a large sample of EFT parameters extracted from HOD mocks with the field-level technique. In this work, we have used the resulting distribution of the perturbative galaxy bias parameters as a prior in the analysis of the BOSS galaxy clustering data in the context of primordial non-Gaussianity. We have found a significant improvement in the PNG constraints, as well as in the measurement of the galaxy bias parameters. Importantly, the values of bias parameters that we have extracted are fully consistent with previous analyses of BOSS based on conservative priors, as well as with the measurements based on the short-scale density-split statistic. This impressive consistency can be interpreted as a sign of convergence of various independent techniques to model the large-scale structure of the Universe.  \nWe note that in this work we have focused on the real space analysis of galaxy clustering that allows us to study perturbative galaxy bias parameters. The informative priors on these parameters are sufficient for our main goal here: to improve limits on single-field primordial non-Gaussianity. A more general analysis, including full redshift-space nuisance parameters and HOD models, is currently in progress.  \nA next step in our program will be to carry out a full re-analysis of the BOSS data including cosmological parameters such as σ 8 , Ω m etc. Another important research direction is to extend our calibration of EFT parameters to more general HOD models relevant, e.g. for eBOSS quasars [113-116] and emission line galaxies [9, 117], as well as to other galaxy formation models, such as abundance matching and hydrodynamical simulations.  \nNote that the approach presented here can be used for a new type of emulator where small scale data are modeled with simulations, and then large scale clustering is reconstructed from the small simulation box using EFT. Similar ideas have been put forward before (see e.g. [84, 118]). A central tool in this emulator is the conditional distribution between EFT and simulation parameters, identical to the one that we derived here for HOD models. The work on this emulator is currently underway.  \nFinally, it would also be interesting to extend our approach to imaging clustering data and intensity mapping, see [84] for recent work in this direction.  \nFIG. 6. Comparison of the EFT biased parameters found by small scale simulation-based fits to density split and two-point statistics based on HOD models (HOD fits), and those found through EFT fits to the power spectrum and the bispectrum on large scales (EFT fits).  \n<!-- image -->  \nCode to reproduce and use the simulation-based priors p ( θ EFT ) and p ( θ EFT | θ HOD ) is available at https:// github.com/smsharma/eft-hod .", 'ACKNOWLEDGMENTS': "We thank Kazuyuki Akitsu, Stephen Chen, Marko Simonovi'c, Fabian Schmidt, and Marcel Schmittfull for useful discussions. The authors also thank the IAIFI Astro-ML Hackathon where the first draft of this work was completed. Monte Carlo Markov Chains with PNG parameters were generated with the Montepython code [119, 120]. This work is supported by the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.  \norg/ ). This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics of U.S. Department of Energy under grant Contract Number DE-SC0012567. AO acknowledges financial support from the Swiss National Science Foundation (grant no CRSII5 193826). MWT acknowledges financial support from the Simons Foundation (Grant Number 929255).", 'Appendix A: Additional plots': 'Fitting the EFT parameters from the field-level model is visualized in Fig. 7 where we show the data transfer functions and their EFT template fits for a typical mock from our sample. The left panel of Fig. 8 shows the extraction of the stochastic α 0 and α 1 parameters by fitting  \nFIG. 7. Transfer functions from our forward model and their EFT template fits for a typical HOD mock.  \n<!-- image -->  \nFIG. 8. Left panel: Error power spectrum from the forward model and its EFT template fit for a typical HOD mock in our sample. The number density of this mock is ¯ n = 3 . 5 · 10 -4 [ h Mpc -1 ] 3 . Right panel : Power spectra of a typical HOD galaxy field, the deterministic part of the EFT forward model, and the residual between the two (the error power spectrum P err ).  \n<!-- image -->  \nthe error power spectrum from the data to the EFT template (14). Typical power spectra of the HOD galaxies and the EFT forward model are shown in the right panel of Fig. 8. The difference between the EFT and the HOD density fields by construction is the error field density, whose power spectrum is also shown in Fig. 8.  \nIn Fig. 9 we show a corner plot with marginalized posterior distribution for PNG parameters plus linear and quadratic bias parameters for all four BOSS data chunks.', 'Appendix B: Normalizing flow training and validation': 'Normalizing flows are implemented using the nflows library, with training and evaluation performed using PyTorch [96]. Masked Autoregressive Flows [121] are used, with 6 flow transformation, each parameterized using a 2-layer masked autoregressive neural network [122] with GELU activations. 10% of the samples are held out for validation, and the flow is trained with batch size 128 for either 30,000 or 50,000 steps (for the EFT and HOD-conditional EFT distributions, respectively) using the Adam optimizer [123] with learning rate 3 × 10 -4 .  \nSince we have a finite number of sampled points for the HOD and EFT parameters, it is imperative to validate that features of our approximated distribution are not artificially sculpted due these points. Figure 10 shows the training (dashed red) and validation (blue) losses over the course of training, with validation being done every 200 steps on a held out set of points. We can see that the loss asymptotes for both the marginal EFT parameter fit (left) and the HOD-conditional fit (right), with no signs of overfitting.  \nFinally, as an additional coverage test, we perform a simulation-based calibration analysis [124] of our normalizing flow model. To that end we generate 100 samples from our trained normalizing flow model and compare them to validation samples. In this process, we calculate the rank of each validation sample among the set of 100 posterior samples drawn from the flow, and plot aggregate statistics, shown in Fig. 11. With perfect calibration, these ranks should be uniformly distributed. No significant deviation from uniformity can be observed, validating the calibration and convergence properties of  \nour learned posterior.', 'Appendix C: Test of UV-sensitivity': "In this section, we provide a convergence test for our transfer function measurements. This analysis also provides an estimate of the residual UV-dependence of the bias parameters extracted with our field-level technique.  \nThe application of the CIC window leads to an implicit smoothing of all fields in our calculations. Note that this smoothing is quite soft due to the particular shape of the CIC window. The cell size of the grid is a proxy of the filtering scale. In our baseline analysis, we use 256 grid points N g , resulting in an effective smoothing scale R s ≃ πk -1 Nyquist = 6 h -1 Mpc.  \nIn Fig. 12 we display transfer functions for one of the HOD mocks in our catalog with four choices of the grid size, N g = 64 , 128 , 256 and 512, which corresponds to the smoothing lengths πk -1 Nyquist = 24 , 12 , 6 , 3 h -1 Mpc, respectively. As expected, we see that using a relatively large cutoff R s affects the transfer functions quite significantly. In particular, the R s = 8 h -1 Mpc transfer functions for the non-linear bias parameters take different values in the k → 0 limit. In contrast, the lowk limit of β 1 is unaffected and yields the renormalized linear bias parameter b 1 [65].  \nThe appearance of the scale dependence can be understood from perturbation theory arguments. For example, at the one-loop order β 2 reads,  \nβ 2 ( k ) = b 2 2 + ⟨ ˜ δ ⊥ 2 G ⊥ 2 ⟩ ⟨ ˜ δ ⊥ 2 δ ⊥ 2 ⟩ . (C1)  \nThe rightmost term above vanishes in the k → 0 limit because the power spectrum of δ 2 is constant on large scales. However, this constant vanishes quickly as one lowers the cutoff, i.e. for R s ≃ 10 h -1 Mpc its value is reduced by a factor of 3 compared to the R s → 0 limit. Hence, the scale dependence of the second term starts playing a more important role as one lowers the cutoff, and ultimately contaminates the lowk limit of β 2 . On top of that, there is an additional scale dependence generated by higher-order loop corrections.  \nIn our particular case, however, Fig. 12 suggests that the cutoff-dependence of the transfer functions, especially their lowk limit, is negligible for the momentum cut k max = 0 . 4 h Mpc -1 used in the analysis.  \nFIG. 9. Corner plots with 2d and 1d marginalized posterior distribution of PNG and galaxy bias parameters b 1 , b 2 , b G 2 from BOSS DR12 power spectra and bispectra for four independent data chunks: NGCz3 (1), SGCz3 (2), NGCz1 (3), SGCz3 (4).  \n<!-- image -->  \n- [1] A. Aghamousa et al. (DESI), (2016), arXiv:1611.00036 [astro-ph.IM].\n- [2] R. Laureijs et al. (EUCLID), (2011), arXiv:1110.3193 [astro-ph.CO].\n- [3] v. Ivezi'c et al. (LSST), Astrophys. J. 873 , 111 (2019), arXiv:0805.2366 [astro-ph].\n- [4] R. Akeson et al. , (2019), arXiv:1902.05569 [astroph.IM].\n- [5] D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, JCAP 1207 , 051 (2012), arXiv:1004.2488 [astroph.CO].\n- [6] J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore,  \nFIG. 10. Training (dashed red line) and validation (blue) line losses over the course of training, for the marginal EFT parameter distribution (left) and HOD-conditional EFT parameter distribution (right).  \n<!-- image -->  \nJHEP 09 , 082 (2012), arXiv:1206.2926 [astro-ph.CO].  \n- [7] M. M. Ivanov, (2022), arXiv:2212.08488 [astro-ph.CO].\n- [8] T. Nishimichi, G. D'Amico, M. M. Ivanov, L. Senatore, M. Simonovi'c, M. Takada, M. Zaldarriaga, and P. Zhang, Phys. Rev. D 102 , 123541 (2020), arXiv:2003.08277 [astro-ph.CO].\n- [9] M. M. Ivanov, (2021), arXiv:2106.12580 [astro-ph.CO].\n- [10] A. Chudaykin, M. M. Ivanov, O. H. E. Philcox, and M. Simonovi'c, Phys. Rev. D 102 , 063533 (2020), arXiv:2004.10607 [astro-ph.CO].\n- [11] S.-F. Chen, Z. Vlah, E. Castorina, and M. White, JCAP 03 , 100 (2021), arXiv:2012.04636 [astro-ph.CO].\n- [12] G. D'Amico, L. Senatore, and P. Zhang, JCAP 01 , 006 (2021), arXiv:2003.07956 [astro-ph.CO].\n- [13] M. M. Ivanov, E. McDonough, J. C. Hill, M. Simonovi'c, M. W. Toomey, S. Alexander, and M. Zaldarriaga, Phys. Rev. D 102 , 103502 (2020), arXiv:2006.11235 [astro-ph.CO].\n- [14] A. Chudaykin, K. Dolgikh, and M. M. Ivanov, Phys. Rev. D 103 , 023507 (2021), arXiv:2009.10106 [astroph.CO].\n- [15] M. M. Ivanov, M. Simonovi'c, and M. Zaldarriaga, JCAP 05 , 042 (2020), arXiv:1909.05277 [astro-ph.CO].\n- [16] G. D'Amico, J. Gleyzes, N. Kokron, D. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Mar'ın, (2019), arXiv:1909.05271 [astro-ph.CO].\n- [17] O. H. E. Philcox and M. M. Ivanov, Phys. Rev. D 105 , 043517 (2022), arXiv:2112.04515 [astro-ph.CO].\n- [18] S.-F. Chen, Z. Vlah, and M. White, JCAP 02 , 008 (2022), arXiv:2110.05530 [astro-ph.CO].\n- [19] S. Alam et al. (BOSS), Mon. Not. Roy. Astron. Soc. 470 , 2617 (2017), arXiv:1607.03155 [astro-ph.CO].\n- [20] V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rept.", '733 , 1 (2018), arXiv:1611.09787 [astro-ph.CO].': "- [21] D. Wadekar, M. M. Ivanov, and R. Scoccimarro, Phys. Rev. D 102 , 123521 (2020), arXiv:2009.00622 [astroph.CO].\n- [22] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonovic, and M. Zaldarriaga, (2022), arXiv:2211.14899 [astro-ph.CO].\n- [23] O. H. E. Philcox, M. M. Ivanov, G. Cabass, M. Simonovi'c, M. Zaldarriaga, and T. Nishimichi, Phys. Rev. D 106 , 043530 (2022), arXiv:2206.02800 [astro-ph.CO].\n- [24] M. M. Ivanov, O. H. E. Philcox, G. Cabass, T. Nishimichi, M. Simonovi'c, and M. Zaldarriaga, Phys. Rev. D 107 , 083515 (2023), arXiv:2302.04414 [astro-ph.CO].\n- [25] A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta, and A. G. S'anchez, (2021), arXiv:2102.06902 [astro-ph.CO].\n- [26] F. Beutler et al. (BOSS), Mon. Not. Roy. Astron. Soc. 466 , 2242 (2017), arXiv:1607.03150 [astro-ph.CO].\n- [27] M. M. Abidi and T. Baldauf, JCAP 1807 , 029 (2018), arXiv:1802.07622 [astro-ph.CO].\n- [28] M. M. Ivanov, O. H. E. Philcox, T. Nishimichi, M. Simonovi'c, M. Takada, and M. Zaldarriaga, Phys. Rev. D 105 , 063512 (2022), arXiv:2110.10161 [astro-ph.CO].\n- [29] A. Barreira, T. Lazeyras, and F. Schmidt, (2021), arXiv:2105.02876 [astro-ph.CO].\n- [30] A. A. Berlind and D. H. Weinberg, Astrophys. J. 575 , 587 (2002), arXiv:astro-ph/0109001.\n- [31] Z. Zheng, A. A. Berlind, D. H. Weinberg, A. J. Benson, C. M. Baugh, S. Cole, R. Dave, C. S. Frenk, N. Katz, and C. G. Lacey, Astrophys. J. 633 , 791 (2005), arXiv:astro-ph/0408564.\n- [32] Z. Zheng, A. L. Coil, and I. Zehavi, Astrophys. J. 667 ,  \nFIG. 11. Simulation-based calibration rank statistics for EFT parameters. Each subplot shows the distribution of ranks for a different parameter: b 1 (top left), b 2 (top right), b G 2 (bottom left), and b Γ3 (bottom right). The horizontal dashed line represents the expected uniform distribution for perfect calibration.  \n<!-- image -->", '760 (2007), arXiv:astro-ph/0703457.': "[33] R. H. Wechsler and J. L. Tinker, Ann. Rev. Astron. Astrophys. 56 , 435 (2018), arXiv:1804.03097 [astroph.GA].  \n[34] Y. Kobayashi, T. Nishimichi, M. Takada, and H. Miyatake, Phys. Rev. D 105 , 083517 (2022), arXiv:2110.06969 [astro-ph.CO].  \n[35] E. Paillas et al. , (2023), arXiv:2309.16541 [astroph.CO].  \n- [36] C. Cuesta-Lazaro et al. , (2023), arXiv:2309.16539 [astro-ph.CO].\n- [37] C. Hahn et al. , (2023), arXiv:2310.15246 [astro-ph.CO].\n- [38] G. Valogiannis, S. Yuan, and C. Dvorkin, (2023), arXiv:2310.16116 [astro-ph.CO].  \nFIG. 12. Dependence of the transfer functions on the Eulerian grid resolution, which acts as a cutoff in the field-level calculation. We show results for the number of grid points N g = 64 , 128 , 256 , 512, which correspond to effective smoothing scales πk -1 Nyquist = 24 , 12 , 6 , 3 h -1 Mpc, respectively.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nP  \n- [39] A. Barreira, G. Cabass, F. Schmidt, A. Pillepich, and D. Nelson, JCAP 12 , 013 (2020), arXiv:2006.09368 [astro-ph.CO].\n- [40] N. Kokron, J. DeRose, S.-F. Chen, M. White, and R. H. Wechsler, Mon. Not. Roy. Astron. Soc. 514 , 2198 (2022), arXiv:2112.00012 [astro-ph.CO].\n- [41] M. Zennaro, R. E. Angulo, S. Contreras, M. PellejeroIb'a˜nez, and F. Maion, Mon. Not. Roy. Astron. Soc. 514 , 5443 (2022), arXiv:2110.05408 [astro-ph.CO].\n- [42] J. M. Sullivan, U. Seljak, and S. Singh, JCAP 11 , 026 (2021), arXiv:2104.10676 [astro-ph.CO].\n- [43] N. Arkani-Hamed, P. Creminelli, S. Mukohyama, and M. Zaldarriaga, JCAP 04 , 001 (2004), arXiv:hepth/0312100.\n- [44] M. Alishahiha, E. Silverstein, and D. Tong, Phys. Rev. D 70 , 123505 (2004), arXiv:hep-th/0404084.\n- [45] L. Senatore, Phys. Rev. D 71 , 043512 (2005), arXiv:astro-ph/0406187.\n- [46] X. Chen, M.-x. Huang, S. Kachru, and G. Shiu, JCAP 01 , 002 (2007), arXiv:hep-th/0605045.\n- [47] P. Creminelli, M. A. Luty, A. Nicolis, and L. Senatore, JHEP 12 , 080 (2006), arXiv:hep-th/0606090.\n- [48] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, JHEP 03 , 014 (2008), arXiv:0709.0293 [hep-th].\n- [49] C. Cheung, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, JCAP 02 , 021 (2008), arXiv:0709.0295 [hep-th].\n- [50] L. Senatore, K. M. Smith, and M. Zaldarriaga, JCAP 01 , 028 (2010), arXiv:0905.3746 [astro-ph.CO].\n- [51] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonovi'c, and M. Zaldarriaga, (2022), arXiv:2201.07238 [astro-ph.CO].\n- [52] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonovi'c, and M. Zaldarriaga, (2022), arXiv:2204.01781 [astro-ph.CO].\n- [53] G. D'Amico, M. Lewandowski, L. Senatore, and P. Zhang, (2022), arXiv:2201.11518 [astro-ph.CO].\n- [54] S.-F. Chen, P. Chakraborty, and C. Dvorkin, (2024),", 'arXiv:2401.13036 [astro-ph.CO].': "- [55] E. Castorina et al. , JCAP 09 , 010 (2019), arXiv:1904.08859 [astro-ph.CO].\n- [56] A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques, and J. Nore˜na, JCAP 05 , 015 (2021), arXiv:2010.14523 [astro-ph.CO].\n- [57] T. Lazeyras, A. Barreira, F. Schmidt, and V. Desjacques, JCAP 01 , 023 (2023), arXiv:2209.07251 [astroph.CO].\n- [58] A. Barreira, JCAP 11 , 013 (2022), arXiv:2205.05673 [astro-ph.CO].\n- [59] A. Barreira and E. Krause, JCAP 10 , 044 (2023), arXiv:2302.09066 [astro-ph.CO].\n- [60] D. Green, Y. Guo, J. Han, and B. Wallisch, (2023), arXiv:2311.04882 [astro-ph.CO].\n- [61] S. Saito, T. Baldauf, Z. Vlah, U. Seljak, T. Okumura, and P. McDonald, Phys. Rev. D 90 , 123522 (2014), arXiv:1405.1447 [astro-ph.CO].\n- [62] M. Schmittfull, T. Baldauf, and U. Seljak, Phys. Rev. D 91 , 043530 (2015), arXiv:1411.6595 [astro-ph.CO].\n- [63] T. Lazeyras and F. Schmidt, JCAP 1809 , 008 (2018), arXiv:1712.07531 [astro-ph.CO].\n- [64] F. Schmidt, F. Elsner, J. Jasche, N. M. Nguyen, and G. Lavaux, JCAP 01 , 042 (2019), arXiv:1808.02002 [astro-ph.CO].\n- [65] M. Schmittfull, M. Simonovi'c, V. Assassi, and M. Zaldarriaga, Phys. Rev. D 100 , 043514 (2019), arXiv:1811.10640 [astro-ph.CO].\n- [66] F. Elsner, F. Schmidt, J. Jasche, G. Lavaux, and N.-M. Nguyen, JCAP 01 , 029 (2020), arXiv:1906.07143 [astroph.CO].\n- [67] G. Cabass and F. Schmidt, JCAP 04 , 042 (2020), arXiv:1909.04022 [astro-ph.CO].\n- [68] C. Modi, S.-F. Chen, and M. White, Mon. Not. Roy. Astron. Soc. 492 , 5754 (2020), arXiv:1910.07097 [astroph.CO].\n- [69] F. Schmidt, JCAP 04 , 032 (2021), arXiv:2009.14176 [astro-ph.CO].  \n- [70] F. Schmidt, G. Cabass, J. Jasche, and G. Lavaux, JCAP 11 , 008 (2020), arXiv:2004.06707 [astro-ph.CO].\n- [71] M. Schmittfull, M. Simonovi'c, M. M. Ivanov, O. H. E. Philcox, and M. Zaldarriaga, JCAP 05 , 059 (2021), arXiv:2012.03334 [astro-ph.CO].\n- [72] T. Lazeyras, A. Barreira, and F. Schmidt, JCAP 10 , 063 (2021), arXiv:2106.14713 [astro-ph.CO].\n- [73] J. Stadler, F. Schmidt, and M. Reinecke, JCAP 10 , 069 (2023), arXiv:2303.09876 [astro-ph.CO].\n- [74] H. Rubira and F. Schmidt, JCAP 01 , 031 (2024), arXiv:2307.15031 [astro-ph.CO].\n- [75] N.-M. Nguyen, F. Schmidt, B. Tucci, M. Reinecke, and A. Kosti'c, (2024), arXiv:2403.03220 [astro-ph.CO].\n- [76] N. A. Maksimova, L. H. Garrison, D. J. Eisenstein, B. Hadzhiyska, S. Bose, and T. P. Satterthwaite, Mon. Not. Roy. Astron. Soc. 508 , 4017 (2021), arXiv:2110.11398 [astro-ph.CO].\n- [77] N. Aghanim et al. (Planck), (2018), arXiv:1807.06209 [astro-ph.CO].\n- [78] S. Yuan, L. H. Garrison, B. Hadzhiyska, S. Bose, and D. J. Eisenstein, Mon. Not. Roy. Astron. Soc. 510 , 3301 (2022), arXiv:2110.11412 [astro-ph.CO].\n- [79] B. Hadzhiyska, S. Bose, D. Eisenstein, L. Hernquist, and D. N. Spergel, MNRAS 493 , 5506 (2020), arXiv:1911.02610 [astro-ph.CO].\n- [80] X. Xu, I. Zehavi, and S. Contreras, MNRAS 502 , 3242 (2021), arXiv:2007.05545 [astro-ph.GA].\n- [81] V. Assassi, D. Baumann, D. Green, and M. Zaldarriaga, JCAP 1408 , 056 (2014), arXiv:1402.5916 [astroph.CO].\n- [82] F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Phys. Rept. 367 , 1 (2002), arXiv:astroph/0112551 [astro-ph].\n- [83] Y. B. Zeldovich, Astron. Astrophys. 5 , 84 (1970).\n- [84] A. Obuljen, M. Simonovi'c, A. Schneider, and R. Feldmann, Phys. Rev. D 108 , 083528 (2023), arXiv:2207.12398 [astro-ph.CO].\n- [85] L. Senatore and M. Zaldarriaga, JCAP 1502 , 013 (2015), arXiv:1404.5954 [astro-ph.CO].\n- [86] T. Baldauf, M. Mirbabayi, M. Simonovi'c, and M. Zaldarriaga, Phys. Rev. D92 , 043514 (2015), arXiv:1504.04366 [astro-ph.CO].\n- [87] D. Blas, M. Garny, M. M. Ivanov, and S. Sibiryakov, JCAP 1607 , 052 (2016), arXiv:1512.05807 [astroph.CO].\n- [88] D. Blas, M. Garny, M. M. Ivanov, and S. Sibiryakov, JCAP 1607 , 028 (2016), arXiv:1605.02149 [astroph.CO].\n- [89] M. M. Ivanov and S. Sibiryakov, JCAP 1807 , 053 (2018), arXiv:1804.05080 [astro-ph.CO].\n- [90] A. Eggemeier, R. Scoccimarro, and R. E. Smith,  \n- (2018), arXiv:1812.03208 [astro-ph.CO].\n- [91] G. D'Amico, Y. Donath, M. Lewandowski, L. Senatore, and P. Zhang, (2022), arXiv:2211.17130 [astro-ph.CO].\n- [92] T. Lazeyras, C. Wagner, T. Baldauf, and F. Schmidt, JCAP 1602 , 018 (2016), arXiv:1511.01096 [astroph.CO].\n- [93] A. Chudaykin, M. M. Ivanov, and M. Simonovi'c, Phys. Rev. D 103 , 043525 (2021), arXiv:2009.10724 [astroph.CO].\n- [94] M. M. Ivanov, O. H. E. Philcox, M. Simonovi'c, M. Zaldarriaga, T. Nischimichi, and M. Takada, Phys. Rev. D 105 , 043531 (2022), arXiv:2110.00006 [astro-ph.CO].\n- [95] D. Rezende and S. Mohamed, in International conference on machine learning (PMLR, 2015) pp. 1530-1538.\n- [96] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. , Advances in neural information processing systems 32 (2019).\n- [97] R. Casas-Miranda, H. J. Mo, R. K. Sheth, and G. Boerner, Mon. Not. Roy. Astron. Soc. 333 , 730 (2002), arXiv:astro-ph/0105008.\n- [98] A. Cooray and R. K. Sheth, Phys. Rept. 372 , 1 (2002), arXiv:astro-ph/0206508.\n- [99] T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Phys. Rev. D 88 , 083507 (2013), arXiv:1305.2917 [astro-ph.CO].\n- [100] T. Baldauf, S. Codis, V. Desjacques, and C. Pichon, Mon. Not. Roy. Astron. Soc. 456 , 3985 (2016), arXiv:1510.09204 [astro-ph.CO].\n- [101] A. P. Hearin, A. R. Zentner, F. C. van den Bosch, D. Campbell, and E. Tollerud, Mon. Not. Roy. Astron. Soc. 460 , 2552 (2016), arXiv:1512.03050 [astro-ph.CO].\n- [102] Y. Akrami et al. (Planck), Astron. Astrophys. 641 , A9 (2020), arXiv:1905.05697 [astro-ph.CO].\n- [103] S. Alam et al. (BOSS), Mon. Not. Roy. Astron. Soc. 470 , 2617 (2017), arXiv:1607.03155 [astro-ph.CO].\n- [104] O. H. E. Philcox, M. M. Ivanov, M. Simonovi'c, and M. Zaldarriaga, JCAP 05 , 032 (2020), arXiv:2002.04035 [astro-ph.CO].\n- [105] F.-S. Kitaura et al. , Mon. Not. Roy. Astron. Soc. 456 , 4156 (2016), arXiv:1509.06400 [astro-ph.CO].\n- [106] S. A. Rodr'ıguez-Torres et al. , Mon. Not. Roy. Astron. Soc. 460 , 1173 (2016), arXiv:1509.06404 [astro-ph.CO].\n- [107] M. M. Ivanov, Phys. Rev. D 109 , 023507 (2024), arXiv:2309.10133 [astro-ph.CO].\n- [108] C. Hahn, R. Scoccimarro, M. R. Blanton, J. L. Tinker, and S. A. Rodr'ıguez-Torres, Mon. Not. Roy. Astron. Soc. 467 , 1940 (2017), arXiv:1609.01714 [astro-ph.CO].\n- [109] M. M. Ivanov, M. Simonovi'c, and M. Zaldarriaga, Phys. Rev. D 101 , 083504 (2020), arXiv:1912.08208 [astroph.CO].\n- [110] M. Mirbabayi, F. Schmidt, and M. Zaldarriaga, JCAP 1507 , 030 (2015), arXiv:1412.5169 [astro-ph.CO].\n- [111] O. H. E. Philcox, Phys. Rev. D 103 , 103504 (2021), arXiv:2012.09389 [astro-ph.CO].\n- [112] O. H. E. Philcox, (2021), arXiv:2107.06287 [astroph.CO].\n- [113] M. Ata et al. , Mon. Not. Roy. Astron. Soc. 473 , 4773 (2018), arXiv:1705.06373 [astro-ph.CO].\n- [114] J. Hou et al. , Mon. Not. Roy. Astron. Soc. 500 , 1201 (2020), arXiv:2007.08998 [astro-ph.CO].\n- [115] R. Neveux et al. , Mon. Not. Roy. Astron. Soc. 499 , 210 (2020), arXiv:2007.08999 [astro-ph.CO].\n- [116] A. Chudaykin and M. M. Ivanov, (2022), arXiv:2210.17044 [astro-ph.CO].\n- [117] A. de Mattia et al. , Mon. Not. Roy. Astron. Soc. 501 , 5616 (2021), arXiv:2007.09008 [astro-ph.CO].\n- [118] C. Modi and O. H. E. Philcox, (2023), arXiv:2309.10270", '[astro-ph.CO].': '- [119] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302 , 001 (2013), arXiv:1210.7183 [astro-ph.CO].\n- [120] T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24 , 100260 (2019), arXiv:1804.07261 [astro-ph.CO].\n- [121] G. Papamakarios, T. Pavlakou, and I. Murray, Advances in neural information processing systems 30 (2017).\n- [122] M. Germain, K. Gregor, I. Murray, and H. Larochelle, in International conference on machine learning (PMLR, 2015) pp. 881-889.\n- [123] D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2014).\n- [124] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman, arXiv e-prints , arXiv:1804.06788 (2018), arXiv:1804.06788 [stat.ME].'}