The last few years have seen the emergence of a wide array of novel techniques for analyzing highprecision data from upcoming galaxy surveys which aim to extend the statistical analysis of galaxy clustering data beyond the linear regime and the canonical twopoint 2pt statistics. We test and benchmark some of these new techniques in a community data challenge Beyond2pt initiated during the Aspen 2022 Summer Program LargeScale Structure Cosmology beyond 2Point Statistics whose first round of results we present here. The challenge dataset consists of highprecision mock galaxy catalogs for clustering in real space redshift space and on a light cone. Participants in the challenge have developed endtoend pipelines to analyze mock catalogs and extract unknown masked cosmological parameters of the underlying LambdaCDM models with their methods. The methods represented are densitysplit clustering nearest neighbor statistics BACCO power spectrum emulator void statistics LEFTfield fieldlevel inference using effective field theory EFT and joint power spectrum and bispectrum analyses using both EFT and simulationbased inference. In this work we review the results of the challenge focusing on problems solved lessons learned and future research needed to perfect the emerging beyond2pt approaches. The unbiased parameter recovery demonstrated in this challenge by multiple statistics and the associated modeling and inference frameworks supports the credibility of cosmology constraints from these methods. The challenge data set is publicly available and we welcome future submissions from methods that are not yet represented.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
A ParameterMasked Mock Data Challenge for BeyondTwoPoint Galaxy Clustering Statistics
{'No Header': '25  \n18  \n30', 'The Beyond-2pt Collaboration': 'Elisabeth Krause, 1 Yosuke Kobayashi, 1, 2 Andrés N. Salcedo, 1 Mikhail M. Ivanov, 3 Tom Abel, 4, 5, 6 Kazuyuki Akitsu, 7 Raul E. Angulo, 8, 9 Giovanni Cabass, 10 Sofia Contarini, 11, 12, 13 Carolina Cuesta-Lazaro, 14, 15, 16 ChangHoon Hahn, 17 Nico Hamaus, 18, 19 Donghui Jeong, 20, 21 Chirag Modi, 22, 23 Nhat-Minh Nguyen, 24, 25 Takahiro Nishimichi, 2, 26, 27 Enrique Paillas, 28, 29 Marcos Pellejero Ibañez, 30 Oliver H. E. Philcox, 31, 32 Alice Pisani, 33, 22, 34, 17 Fabian Schmidt, 35 Satoshi Tanaka, 26 Giovanni Verza, 36, 22 Sihan Yuan, 4, 6 Matteo Zennaro, 37  \n1 Department of Astronomy/Steward Observatory, The University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA 2 Department of Astrophysics and Atmospheric Sciences, Faculty of Science, Kyoto Sangyo University, Motoyama, Kamigamo, Kita-ku, Kyoto 603-8555, Japan  \n3 Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n4  \nKavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA  \n5 Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA  \n6 SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA  \n- 7 Theory Center, Institute of Particle and Nuclear Studies, High Energy Accelerator Research Organization(KEK), Tsukuba, Ibaraki 305-0801, Japan  \n8 Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain  \n9 IKERBASQUE, Basque Foundation for Science, E-48013, Bilbao, Spain  \n10 Division of Theoretical Physics, Ruđer Bošković Institute, Zagreb HR-10000, Croatia  \n- 11 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy  \n12 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy  \n13 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, I-40129 Bologna, Italy 14 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n- 15 Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n16 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA  \n17 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA  \nUniversitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 München, Germany  \n19 Excellence Cluster ORIGINS, Bolzmannstr. 2, 85748 Garching, Germany  \n20  \nDepartment of Astronomy and Astrophysics and Institute for Gravitation and the Cosmos, The Pennsylvania State University,  \nUniversity Park, PA 16802, USA  \n21 School of Physics, Korea Institute for Advanced Study, Seoul, South Korea  \n22 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA  \n23 Center for Computational Mathematics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA  \n- 24 Leinweber Center for Theoretical Physics, University of Michigan, 450 Church St, Ann Arbor, MI 48109-1040  \nDepartment of Physics, College of Literature, Science and the Arts, University of Michigan, 450 Church St, Ann Arbor, MI 48109-1040  \n26  \nCenter for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto  \n606-8502, Japan  \n27 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  \n28 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L 3G1, Canada  \n29 Department of Physics and Astronomy, University of Waterloo, ON N2L 3G1, Canada  \nInstitute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ , UK  \n31 Department of Physics, Columbia University, New York, NY 10027, USA  \n32 Simons Society of Fellows, Simons Foundation, New York, NY 10010, USA  \n- 33 Aix-Marseille University, CNRS/IN2P3, CPPM, 163 Av. de Luminy, 13009, Marseille, France  \nkrausee@arizona.edu, yosukekobayashi@arizona.edu, ansalcedo@arizona.edu  \n- ∗ We use to mask (and derived forms) in place of the formerly common to blind throughout this manuscript.  \n34  \nThe Cooper Union for the Advancement of Science and Art, 41 Cooper Square, New York, NY 10003, USA 35 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany  \n36 Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003, USA 37 University of Oxford, Astrophysics, Denys Wilkinson Building, Keble Rd, Oxford OX1 3RH, Regno Unito', 'ABSTRACT': "The last few years have seen the emergence of a wide array of novel techniques for analyzing highprecision data from upcoming galaxy surveys, which aim to extend the statistical analysis of galaxy clustering data beyond the linear regime and the canonical two-point (2pt) statistics. We test and benchmark some of these new techniques in a community data challenge 'Beyond-2pt', initiated during the Aspen 2022 Summer Program 'Large-Scale Structure Cosmology beyond 2-Point Statistics,' whose first round of results we present here. The challenge dataset consists of high-precision mock galaxy catalogs for clustering in real space, redshift space, and on a light cone. Participants in the challenge have developed end-to-end pipelines to analyze mock catalogs and extract unknown ('masked') cosmological parameters of the underlying Λ CDM models with their methods. The methods represented are density-split clustering, nearest neighbor statistics, BACCO power spectrum emulator, void statistics, LEFTfield field-level inference using effective field theory (EFT), and joint power spectrum and bispectrum analyses using both EFT and simulation-based inference. In this work, we review the results of the challenge, focusing on problems solved, lessons learned, and future research needed to perfect the emerging beyond-2pt approaches. The unbiased parameter recovery demonstrated in this challenge by multiple statistics and the associated modeling and inference frameworks supports the credibility of cosmology constraints from these methods. The challenge data set is publicly available at this repository /github , and we welcome future submissions from methods that are not yet represented.  \nKeywords: Cosmological parameters from large-scale structure (340) - Large-scale structure of the universe(902) - Astrostatistics techniques(1886)", '1. INTRODUCTION': "Cosmic large-scale structure-the large-scale distribution of galaxies in the late universe-is shaped by matter density fluctuations in the early universe, the growth rate of the fluctuations, and the expansion rate of the smooth cosmic background. Ongoing and upcoming large-scale structure surveys will collectively map the distribution of galaxies over an extensive fraction of the sky at unprecedented depth (Aghamousa et al. 2016; Laureijs et al. 2011; Takada et al. 2014). These observations will enable precision tests of fundamental physics, placing tight constraints on inflation, neutrino masses, and novel properties of dark matter, gravity, and dark energy.  \nOn large scales, matter and galaxy distributions preserve the almost perfectly Gaussian statistics from the initial conditions according to the standard inflationary models. The two-point (2pt) correlation function or power spectrum captures all information in a Gaussian field. They therefore have been the linchpin of cosmological inference and interpretation.  \nThe situation, however, is beginning to change. Larger and higher-resolution simulations (e.g. Nishimichi et al. 2019; Heitmann et al. 2019; Ishiyama et al. 2021; Maksimova et al. 2021; Schaye et al. 2023; DeRose et al. 2023) and more realistic galaxy mocks (e.g. Behroozi et al. 2019; DeRose et al. 2019; Balaguera-Antolínez et al. 2023; Yuan et al. 2023; Stevens et al. 2023; To et al. 2024) have been developed to extend the inferences to smaller scales (e.g. Reid et al. 2014; Reddick et al. 2014; Lange et al. 2022; Yuan et al. 2022; Lange et al. 2023; Storey-Fisher et al. 2024),. Furthermore, the gain from small-scale information will become less limited by signal-to-noise ratios in galaxy samples from upcoming galaxy surveys.  \nOn smaller scales, 2pt statistics cannot capture all the information available in the data: even with Gaussian initial perturbations, non-Gaussian features naturally emerge in the galaxy distribution as a result of nonlinear gravitational evolution and galaxy formation. Beyond-2pt statistics aim to extract information from non-Gaussian features arising from nonlinear clustering. Examples of the first cosmological constraints from galaxy clustering using beyond-2pt statistics include bispectrum analyses of PSCz (Feldman et al. 2001), 2dF (Verde et al. 2002), VVDS (Marinoni et al. 2008), SDSS DR6-DR7 (Gaztanaga et al. 2009), and SDSS-III BOSS galaxies (Gil-Marín et al. 2017). More recent  \nexamples include analyses of the galaxy bispectrum (Philcox & Ivanov 2022; D'Amico et al. 2022a; Ivanov et al. 2023; Hahn et al. 2023), three-point function (Sugiyama et al. 2023), skew spectrum (Hou et al. 2024), density-split clustering (Paillas et al. 2023), cosmic voids (Hamaus et al. 2020; Contarini et al. 2023), and wavelet scattering transform of the galaxy density field (Valogiannis et al. 2023; Régaldo-Saint Blancard et al. 2023).While the discussion here focuses on galaxy clustering, beyond-2pt statistics are similarly gaining traction in weak lensing analyses (e.g., Petri et al. 2015; Harnois-Déraps et al. 2021; Gatti et al. 2022; Heydenreich et al. 2022; Cheng et al. 2024).  \nIn most cases, these analyses provided competitive or complementary constraints to those derived from standard 2-pt analyses applied to the same data. Beyond-2pt statistics are poised to mature into a prominent role in cosmological inference from forthcoming galaxy clustering data as nonlinear clustering becomes more accessible and statistical power increases, owing to improvements in the inference models and the clustering signal-to-noise, respectively.  \nWith great statistical power comes great systematic responsibility. How (in)sensitive are the beyond-2pt statistics to their modeling choices? Do they respond in known ways to observational systematics? Modeling beyond-2pt statistics and thus physics on the associated nonlinear scales adds considerable complexities in both models of matter density perturbation and matter-galaxy connection, i.e. galaxy bias. Current beyond-2pt statistics adopt a wide array of approaches for each modeling step: the matter field is modeled with a variety of perturbative approaches or N-body simulations, and the matter-galaxy connection is usually encoded via bias expansions, Halo Occupation Distribution (HOD) models, or Sub-Halo Abundance Matching (SHAM) schemes.  \nIn addition to modeling systematics, how well do we understand the likelihoods and covariances of these novel statistics? Can we evaluate and estimate them at the precision required by current and future surveys? Specifically, intractable non-Gaussian likelihoods 1 must be approximated with simulations. Furthermore, analytic covariances are not always available, hence simulations and mocks are again required to estimate the covariances. The difficulties inherent in the generation of numerous high-resolution simulations and mock data add another layer of complexity to such analyses, which could potentially compromise their robustness. Consequently, the first aim of this data challenge is to validate a variety of beyond-2pt statistics and their modeling approaches.  \nOur second aim is to address the following question: How much information can we (robustly) extract from galaxy clustering? Although individual comparisons of information content between standard and beyond-2pt statistics have been made elsewhere, no direct comparison between the beyond-2pt methods presented here was attempted with the same survey volume and with the same galaxy density. This challenge therefore aims to facilitate a community exercise to study the information content in beyond-2pt statistics by providing a standard set of simulated mock data. A critical quantity that controls the amount of beyond-linear information is the maximum cutoff scale in each analysis, k max . As each beyond-2pt analysis might have a different sensitivity to k max , we let each group determine their own k max up to which they can trust their model against model misspecifications. However, to avoid fine-tuning while simultaneously encouraging better uncertainty quantification and systematic control in each analysis, we mask the parameters of the mock catalogs, which include both the cosmological parameters and the galaxy HOD parameterization(s) plus their parameter values.  \nFor galaxy clustering, only a single public, parameter-masked mock challenge exists prior to this work-the PTchallenge , introduced in Nishimichi et al. (2020)-that specifically targets 2pt statistics. PTchallenge was cornerstone for the establishment of EFT as a mature tool to analyze the full shape of the galaxy power spectrum. In particular, it demonstrates the ability of one-loop EFT theories to recover masked cosmological parameters from both amplitude and shape of the (nonlinear) power spectrum at a sub-percent precision.  \nImportant technical lessons from PTchallenge are: (1) learning parameter degeneracies and optimizing analysis choices with respect to them, (2) developing a methodology to select scale cuts, (3) understanding the role of EFT parameters and their priors. Since the original publication, the true cosmology of the PTchallenge has been kept masked and the PTchallenge continues to serve as a testing ground for new galaxy power spectrum models, e.g., Lagrangian EFT (Chen et al. 2020, 2021), a simulation- plus HOD-based emulator (Kobayashi et al. 2020, 2022), and a perturbation-based shape-fitting approach (Brieden et al. 2021, 2022). However, the PTchallenge only made available data in the form of redshift-space power spectrum measurements.  \nOur 'Beyond-2pt' mock challenge extends the PTchallenge in Nishimichi et al. (2020) in several ways. First, we present the data directly at the level of mock galaxy catalogs, rather than summary statistics thereof. Estimator  \nvalidation is therefore part of the challenge for each clustering statistics. Second, we present mock catalogs at three different complexity levels: real-space or redshift-space snapshots and a light cone. These mocks cover a redshift range around z ∼ 1 with different masked flat Λ CDM cosmologies and different mock galaxy populations with masked HOD parameterization(s). The range of complexity levels enables participation by statistics and methods at different maturity levels, while also facilitating robustness tests of different modeling prescriptions. Third, prior to this, many prescriptions were only tested on mock galaxies resembling those found in the BOSS CMASS sample within a Plancklike cosmology. This challenge-extending both galaxy and cosmology models-therefore marks a step forward in this regard as well.  \nFinally, we extend the challenge to beyond-2pt statistics-showing their constraints on cosmological parameters, [Ω m , σ 8 ] , side-by-side with the corresponding constraints from BACCO P, a representative of 2pt-statistics, for all setupsand dissecting the information content they capture. The present challenge presents a first of such benchmark and sets the stage for more detailed comparisons in future. With eight independent analyses submitted from seven international teams, this challenge represents the first community effort toward developing optimal strategies to extract cosmological information from upcoming galaxy surveys.  \nThe mocks and analyses presented in this paper assume a flat ΛCDM cosmology, which can be described by five parameters: (1) the dimensionless Hubble parameter h , (2) the matter fluctuation or primordial power spectrum amplitude, parameterized by either σ 8 or A s , (3) the power spectrum spectral index n s , plus (4-5) the density parameter of cold dark matter Ω cdm and baryons Ω b . Different analyses adopt priors on either the density parameters Ω x or the physical density parameters ω x = Ω x h 2 . Though no analysis varies neutrino mass, some assume massless ( ω ν = 0 ) while others adopt minimal-mass neutrinos ( ω ν = 0 . 0006442 ). The main results of this challenge are summarized and presented in marginalized constraints on σ 8 and the total matter density, Ω m = Ω cdm +Ω b +Ω ν .  \nTo set the stage for the Summary of Results in Sec. 2, we briefly review the information content in galaxy clustering, in the context of the ΛCDM model (with Gaussian initial conditions), as considered throughout this paper. In linear theory, all cosmological information encoded in galaxy clustering is captured by the linear power spectrum, particularly its scale dependence and amplitude. The scale dependence is sensitive to both the initial density perturbations, captured by the spectral index n s , and the stress-energy components that determine the background evolution and growth of perturbations. The physical density parameters ω m , b determine two characteristic scales: the Hubble horizon at matter-radiation equality and the sound horizon at photon-baryon decoupling, which are imprinted on the broadband shape of the matter power spectrum and baryonic acoustic oscillation (BAO) wiggles, respectively. All the challenge mocks include this information.  \nIn real space, i.e., without accounting for observational effects of galaxy peculiar motion, the amplitude of the (linear) galaxy power spectrum at redshift z is proportional to the product of the amplitude of primordial fluctuations, the linear growth factor D ( z ) and linear galaxy bias b 1 . That is, the cosmology parameter σ 8 is completely degenerate with the unknown linear galaxy bias parameter. Observational effects encode additional information: the redshift space distortion (RSD) effect contains information on the growth rate f , which breaks the degeneracy between σ 8 and b 1 , and the Alcock-Paczynski (AP) effect contains geometry information, which constrains Ω m (within Λ CDM). Growth information from RSD is included in redshift-space and light-cone mocks, while geometry information from AP is only available on the light cone.  \nBeyond the linear regime, quasi- and nonlinear clustering and collapsed structures encode additional information with different parameter degeneracies. The results from this challenge underscore the robust extraction of such nonlinear information either through the nonlinear matter power spectrum ( BACCO P) or through statistics that have access to higher-order n -point functions (all other teams). At the power spectrum level, nonlinear evolution introduces additional smearing of the BAO feature and an enhancement of small-scale power (Rimes & Hamilton 2005; Crocce & Scoccimarro 2008; Angulo et al. 2021). The distinct dependencies of beyond-2pt statistics on cosmology and astrophysics break parameter degeneracies and improve constraints on cosmological parameters. We refer to individual analysis sections, and references therein, for detailed discussions on the sources of nonlinear information.  \nThe remainder of the paper is organized as follows. To help readers navigate the following sections, we provide tables summarizing the different analyses (Tab. 1) and the common notation (Tab. 2). In Sec. 2, we summarize and discuss the results of eight parameter-masked analyses. Fig. 1-4 present the key results, as constraints on the parameter combination [Ω m , σ 8 ] from the redshift-space mocks in Fig. 1-2 or from the light-cone mock in Fig. 3, and as constraints on the parameter σ 8 from the real-space mocks in Fig. 4. We introduce the suite of mock catalogs-the main data product of the 'Beyond-2pt' challenge-in Sec. 3, and the parameter (un)masking procedure in Sec. 4. In  \nFigure 1. 2D marginalized constraints on Ω m and σ 8 for parameter-masked analyses of redshift-space mocks (mean of 10 realizations), marginalized over the remaining cosmological parameters of flat ΛCDM and nuisance parameters specific to each method. We list the scale cuts due to nonlinear modeling for different analyses, which is specified in terms of redshift-space separation ( s ), projected radius ( r p ) or Fourier mode ( k ), depending on the analysis method.  \n<!-- image -->  \nm  \nm  \nSec. 5, each team describes their analysis method. We describe post-unmasking reanalyses and resulting lessons for method refinement and optimization of constraining power in Sec. 6. We summarize the results of this challenge and discuss implications for future mock challenges and beyond-2pt analyses in Sec. 7.", '2. SUMMARY OF RESULTS': 'In this section, we summarize the main results of the paper, shown in Fig. 1-4, and provide an overview of the eight participating analyses in Tab. 1.  \nThis challenge is based on a series of mock galaxy catalogs created from N-body simulations with a flat Λ CDM cosmology and HOD galaxy-halo connection models. It contains three levels of increasing realism:  \n- 1. real-space galaxy distribution from a simulation snapshot with periodic boundary conditions (10 realizations with the same cosmology + HOD parameters at z = 1 ),\n- 2. redshift-space galaxy distribution from a simulation snapshot with periodic boundary conditions (10 realizations with the same cosmology + HOD parameters at z = 1 ),  \n- 3. light-cone galaxy mock emulating the observational data most closely (one realization, uniform coverage across 0 . 8 < z < 1 . 3 and with a simple footprint bounded by lines of constant right ascension and declination).  \nInitially, the organizers communicated to the analysis teams only that the challenge catalogs are HOD-based galaxy mocks in flat Λ CDM cosmologies without observational systematics, hiding the true cosmological parameters and the HOD model used to generate the simulations from the analysis teams. Each analysis team then analyzed these parameter-masked mocks while documenting their model, inference, and analysis choices including scale cuts and parameter priors, before submitting their parameter posteriors to the organizers, who subsequently unmasked the fractional distances between posterior means and the ground truth.  \nThe challenge results presented here indicate whether a given analysis returns parameters consistent with the ground truth given their reported error bars. For real-space and redshift-space analyses, analysis teams report results corresponding to the mean of 10 realizations analyzed with the covariance of a single realization volume to minimize parameter biases due to cosmic variance. Assuming a Gaussian data likelihood, in this setup the probability of a ≥ 1 σ fluctuation due to cosmic variance is 0.16%. Hence, biases in the inferred parameters likely indicate model misspecification or incomplete uncertainty modeling. The challenge results do not directly establish whether the reported parameter uncertainties fully reflect the true uncertainties (measurement and model). This is particularly relevant for small-scale analyses and simulation-based methods, based on empirical galaxy-halo connection parameterizations and priors. Hence, the results should not be regarded as a direct quantitative comparison between different analysis methods or summary statistics.  \nWereport results ordered by number of participating analyses for the different mocks. All results on the redshift-space and light-cone mock catalogs are from analyses in flat Λ CDMwith priors that are much broader than current-generation observational constraints (Planck Collaboration et al. 2020; Alam et al. 2021; DESI Collaboration et al. 2024). However, the cosmology priors differ between all analyses, as illustrated in Fig. 5. We note that the void size function (VSF) analysis of the light-cone mock adopts informative priors on three parameters that are weakly constrained by the VSF, with width of three times the current Planck2018 (Planck Collaboration et al. 2020) uncertainties. The analyses on the real-space mock reported here fix all cosmological parameters except σ 8 to their true values.', '2.1. Results on redshift-space mocks': "The redshift-space snapshots were analyzed by most teams, comprising of the following analyses: density-split clustering (Section 5.6), two different flavors of nearest neighbor statistics (kNN; Section 5.5), a joint power spectrum plus bispectrum analysis using effective field theory (EFT P+B; Section 5.2), a hybrid-EFT power spectrum analysis with the BACCO emulator ( BACCO P; Section 5.1), and a simulation-based-inference analysis of the power spectrum plus bispectrum (SBI P+B; Section 5.4). The SBI P+B analysis was trained on simulations with volume (1 h -1 Gpc) 3 , while all other analyses assume a covariance matrix corresponding to the sample variance of a (2 h -1 Gpc) 3 box, resulting in degradation in constraining power. While it would be convenient to assume a simple survey volume rescaling to compare SBI P+B results with other analyses, we refrain from making such a rescaling as constraining power is a function of both survey volume and scale cuts ( k max ), which need to be calibrated anew when changing the survey volume.  \nFigure 1 shows the marginalized constraints in the Ω m -σ 8 plane inferred from the different summary statistics. Remarkably, all analysis teams successfully recover the input cosmology within their 1 -σ confidence region. This result from a parameter-masked challenge further demonstrates the maturity of these reportedly 'novel' statistics and their potential for analyses of near-term data.  \nFigure 2 shows the marginalized posteriors in 1D for Ω m and σ 8 . The nominal distance between the mean and the ground truth ( i.e. bias) in Ω m is the smallest for the EFT P+B, BACCO and density split analyses, while for σ 8 , the unmasking showed the lowest parameter bias for the EFT P+B, kNN RD, and SBI P+B analyses. After unmasking, further analysis by the BACCO team identified a large emulation uncertainty of the hexadecapole as potential source of this parameter bias (c.f. Sect. 5.1.3), which had been unnoticed in previous tests on smaller simulation volumes. We further discuss this post-unmasking reanalysis in Sect. 6. Focusing on the analyses with fixed θ ∗ -the angular size of the sound horizon-we find similar error bars on Ω m from the density-split, EFT P+B and BACCO methods, implying that, for this particular parameter and mock data, the bulk of the information comes from quasi-linear (BAO) scales. The tightest error bars on σ 8 are obtained with the kNN and BACCO methods.  \nBoth density-split clustering and kNN statistics are modeled with emulators built on the AbacusSummit simulations, which fix θ ∗ , while the SBI P+B analysis is built on the Quijote simulations which assume no such constraint. To ease  \nthe comparison, the EFT P+B and BACCO P teams ran additional analyses imposing the same prior. A comparison of the EFT P+B constraints in both parameter spaces (top and bottom part of Fig. 2) indicates that this difference in cosmology priors has limited impact on the constraints shown here. The EFT P+B and SBI P+B results are quite similar in terms of σ 8 constraining power, despite the difference in survey volume. For Ω m however, the SBI P+B constraint is significantly wider than the EFT P+B result. We believe that the wide SBI P+B posterior on Ω m is a result of trimming the largest scales in the SBI P+B analysis, which down weights large scale modes that are important in the Ω m recovery from the turnover of the galaxy power spectrum (cf. Ivanov et al. 2020a; Philcox et al. 2021b).  \nredshift-space snapshots (mean of 10 realizations), analyzed in flat ΛCDM  \nFigure 2. 1D marginalized constraints on Ω m and σ 8 for parameter-masked analyses of redshift-space mocks (mean of 10 realizations, errors of 1 box), marginalized over the remaining cosmological parameters of flat ΛCDM and nuisance parameters specific to each method.  \n<!-- image -->  \nFigure 3. 1D marginalized constraints on Ω m and σ 8 for parameter-masked analyses of light-cone mock (single realization), marginalized over the remaining cosmological parameters of flat ΛCDM and nuisance parameters specific to each method. Void-size cuts are specified in terms of the effective radius R and the cleaned (re-scaled) radius ˜ R .  \n<!-- image -->  \nFigure 4. 1D marginalized constraints on σ 8 for analyses of real-space mocks (mean of 10 realizations), in a reduced parameter space with all cosmological parameters except σ 8 fixed at the fiducial values. Crosses indicate the maximum a posteriori (MAP) estimates, to illustrate projection effects.  \n<!-- image -->  \nt  \nT ab le 1. Ov erview of participating analyses and their analysis ingredien ts. a Cosmology priors are listed for the baseline v alidation and original unmasking submission. T eams BACCO P and EFT P+B kindly reran their analyses with differen cosmology priors (but o th e rwise iden t i c al analysis c hoices) to facilitate comparisons.", '2.2. Results on light-cone mock': 'Unlike the redshift-space data discussed above, there exists only one realization for the light-cone mock. Hence scatter in recovered parameter values is expected due to cosmic variance of the summary statistics measurements. Since different analyses include (substantially) different scales, they are subject to different cosmic variance modes and may be scattered in different directions.  \nThe effective volume of this mock catalog is similar to that of the BOSS DR12 LRG sample. For this mock, we have submissions from the EFT P+B, the void size function (VSF; Section 5.7), and the void-galaxy cross-correlation function (VGCF; Section 5.7) teams; the parameter constraints are shown in Fig. 3. Both void methods and EFT P+B successfully recover the input values of Ω m within 95% CL, and the EFT pipeline has additionally recovered the true input value of σ 8 . We note that the true cosmological parameters are outside the informative prior (on h , ω b and n s ) adopted by the VSF analysis. However, post-unmasking analyses indicate at most minor bias in Ω m due to the miscentered prior. The impact of prior width on the VSF constraining power is discussed in Sec. 5.7.4. The void team has not submitted the σ 8 results, as its accurate inference from the VSF is not possible without an independent calibration of the tracer bias inside voids, e.g. taken from realistic mocks of the sample in question (c.f., Contarini et al. 2023, and Section 5.7.2). For analyzing the VGCF, in this paper only the Alcock-Paczynski test is exploited for cosmological inference, which exclusively constrains parameters describing the background evolution.', '2.3. Results on real-space mock': 'The real-space challenge boxes were analyzed with the EFT-based field-level Bayesian inference (EFT FBI) approach using the LEFTfield code (Section 5.3), EFT P+B, and the BACCO power spectrum emulator. As for the redshift-space analyses, the summary statistics measurements were averaged over the 10 realizations of the Λ CDM real-space galaxy catalog from a periodic box, and all BACCO and EFT P+B inferences were performed assuming a covariance of a (2 h -1 Gpc) 3 box. For EFT FBI, we show the average of the 10 independent posteriors.  \nAll three analysis teams fixed all cosmological parameters to the simulation true values in this challenge, except σ 8 . EFT P+B and BACCO had originally validated and unmasked their analyses of the real-space mocks in the full flatΛ CDM parameter space. After unmasking, both teams kindly reran their analyses with all cosmological parameters except σ 8 set to their truth values to enable comparison with EFT FBI. This allows a direct comparison between σ 8 constraints obtained by the three statistics, as the EFT FBI team and their LEFTfield pipeline currently cannot perform a full inference over the Λ CDM cosmological parameter space. This simplified setup however still offers insights into the information content in (nonlinear) clustering beyond the 2-point function. The results are presented in Figure 4. The constraints of the BACCO analysis are obtained with only the real-space power spectrum, and are thus heavily impacted by the b 1 -σ 8 degeneracy present at leading order in the power spectrum. The degeneracy is broken only by non-linear contributions to the power spectrum. The EFT P+B and FBI teams use extended data vectors to mitigate this b 1 -σ 8 degeneracy.  \nAll analysis teams recover the true value of σ 8 within 95% CL of one single simulation box. At face value, the EFT FBI analysis achieves the tightest constraints on σ 8 . Compared with the EFT P+B analysis, the FBI analysis assumed a more restricted model for galaxy bias and galaxy stochasticity. The two EFT teams subsequently attempted to quantify the impact of their modeling differences in post-unmasking studies; we refer readers to Sec. 5.3 and Sec. 6 for the results and discussions.', '3. MOCK DATA': 'The challenge organizers present three different types of mock catalogs based on N-body simulations: real-space snapshot, redshift-space snapshot, and light cone-all using the HOD prescription.', '3.1. N-body simulations and halo catalogs': '3.1.1. Flat Λ CDM snapshot mocks  \nAfter agreement on the range for the cosmological parameters described in Section 4, the organizers created the flat Λ CDM mock catalogs in z = 1 . Ten realizations have been created for each of real- and redshift-space mocks, each of which is a cubic simulation box with a comoving side length of 2 h -1 Gpc . The N -body simulations have been run by using the cosmological N -body simulation code GINKAKU (Nishimichi, Tanaka & Yoshikawa in prep.). This code employs the Tree Particle-Mesh (TreePM) method to compute the gravitational force in an expanding periodic box in comoving coordinates. The short-range tree force is implemented based on the Framework for Developing Particle  \nSimulators ( FDPS ; Iwasawa et al. 2016; Namekata et al. 2018), a public library for general particle simulations. FDPS facilitates scalable computations on modern supercomputer systems, ensuring optimized workload balance with efficient domain decomposition. The tree force is further accelerated by SIMD instructions implemented in the Phantom-GRAPE library (Nitadori et al. 2006; Tanikawa et al. 2012, 2013). Details of the long-range PM force can be found in Yoshikawa & Fukushige (2005) (see also Ishiyama et al. 2009, 2012, for recent implementations).  \nThe initial conditions are generated using second-order Lagrangian perturbation theory (2LPT; Scoccimarro 1998; Crocce et al. 2006) at z = 49 . Similarly to the existing N -body simulation ensembles like AbacusSummit or Dark Quest , the organizers have run the simulations without the effect of massive neutrino dynamics, and have incorporated massive neutrinos only in the matter transfer functions with the abundance ω ν = 0 . 000644 . The transfer function has been computed at z = 0 and rescaled to the starting redshift z = 49 using the linear growth factor with matter density Ω m = ( ω b + ω cdm + ω ν ) h -2 , but without the scale dependence induced by massive neutrinos.  \nFor each realization, dark-matter halos are populated with mock galaxies as described in Sect. 3.2. We keep the choice of halo finder masked as differences in halo finder can contribute to model misspecification in the highly nonlinear regime (e.g., Hahn et al. 2023; Modi et al. 2023). Neither real-space nor redshift-space mocks include the AlcockPaczynski effect (Alcock & Paczynski 1979). This effect contributes to constraining the cosmological parameters in a more realistic galaxy spectroscopic survey setting and is taken into account in the light-cone mocks described next below.', '3.1.2. Flat Λ CDM light-cone mocks': 'The single-realization Λ CDM light-cone mock is based on one of the publicly available light cone halo catalogs from the AbacusSummit set of simulations (Maksimova et al. 2021; Hadzhiyska et al. 2021) generated with the abacus cosmological N -body code (Garrison et al. 2019). The specific AbacusSummit cosmology realization was chosen to be within the parameter range of all emulators participating in this challenge (c.f. Sect. 4). The organizers shared that the light-cone mock is one of the AbacusSummit cosmologies with participants only after unmasking.  \nTwo teams, ( EFT P& B and Cosmic Voids ), submitted results for the light-cone mock. Both analyses rely on analytic modeling prescriptions, eliminating concerns that the same halo mock catalog might inadvertently be used in training an emulator for this challenge.', '3.2. HOD galaxy mocks': "To generate mock galaxy catalogs for our challenge, we take advantage of the HOD formalism (e.g. Berlind & Weinberg 2002; Bullock et al. 2002; Zheng et al. 2005; Gonzalez-Perez et al. 2018; Salcedo et al. 2022a; Yuan et al. 2023b) In practice, we specify a parameterized form for the mean halo-mass-dependent occupation of galaxies and stochastically populate the dark matter halos with galaxies. The galaxy distribution within a host halo broadly traces the host's internal structure, but may differ in detail in a halo-mass-independent way. We may also in principle include some level of galaxy assembly bias, which refers to the possibility for the occupation of halos of a given mass to depend on properties other than mass (e.g. Hearin et al. 2016; McEwen & Weinberg 2018; Xu et al. 2021; Salcedo et al. 2022b; Wang et al. 2022; Beltz-Mohrmann et al. 2023; Contreras et al. 2023b; Zhai et al. 2023b). In this context we may implement galaxy assembly bias with respect to halo properties present in our simulation catalogs, or environmental properties of the simulation particle distribution. The light-cone mock may include redshift evolution of some of our HOD parameters. In our redshift-space mocks (including the light cone) we also add galaxy peculiar velocities that may include velocity bias (e.g. van den Bosch et al. 2005; Reid et al. 2014; Guo et al. 2016; Yuan et al. 2018; Anbajagane et al. 2022; Lange et al. 2022; Beltz-Mohrmann et al. 2023; Kwan et al. 2023; Zhai et al. 2023a; Kwon & Hahn 2024). The positions of mock galaxies are computed by modulating the real-space coordinates along the z -axis according to their peculiar velocities.  \nWe emphasize that in constructing the mock galaxy catalogs in this challenge, we have remained broadly within the existing literature on the galaxy-halo connection and have by no means attempted to confound the challenge participants.", '4.1. Pre-Unmasking Information': 'Initially, the only information communicated to the analysis teams was the cosmology parameter space (flat w CDM) and that the galaxy assignment is based on a HOD technique, without specifying the detailed HOD implementation.  \nDuring initial test runs on these w CDM mocks, several teams realized that the (still masked) cosmologies were likely outside the parameter support of some emulator-based methods.  \nFigure 5. 1D marginalized cosmology priors in flat ΛCDM for the challenge mocks (solid grey line) and those adopted by individual analysis teams. Analyses that fix θ ∗ are shown with dashed lines, analyses with free θ ∗ are shown as dotted lines. While there are only four/five independent parameters in flat ΛCDM with/without θ ∗ constraint, analysis teams specify their priors in different parameterizations. To facilitate comparison, we project priors in different parameters, noting that not all panels are independent.  \n<!-- image -->  \ns  \nb  \nThe analysis teams then agreed on common range for the cosmological parameters, which correspond to the intersection of the parameter spaces supported by the AbacusSummit , Aemulus and BACCO emulators, which includes a hard constraint on the angular scale of the sound horizon at decoupling θ ∗ imposed in the AbacusSummit simulation suite (Maksimova et al. 2021).  \nA script implementing this parameter restriction was shared with all analysis teams 2 , and the challenge organizers chose new sets of cosmological parameters to generate flat Λ CDM mock data within the specified parameter range. The cosmology prior for the mock catalogs, as well as the priors of individual analyses (documented in Sec. 5) are shown in Fig. 5.  \nEach analysis team documented their analysis choices and unmasking criteria, as well as potential caveats of their analysis (which could be revisited in the event of an unmasking surprise).  \nWe show only parameter offsets to not unmask future participants and enable continued use of the mock datasets as a benchmark for testing novel analysis techniques. The organizers encourage future submissions from other analysis teams and commit to continuous unmasking.', '4.2. Unmasking Process': 'There is only one parameter unmasking stage. Each analysis team chooses when they are ready to unmask and share their parameter chains with the challenge organizers, who produce the plots described in Sec. 2 and share values of the (marginalized) parameter shifts ∆Ω m and ∆ σ 8 with the analysis team. Analysis teams agreed not to share parameter values with other teams but to only show plots with offset parameter values.  \nTable 2. Notation for parameter inference variables.  \nU  \n∼ U', '4.3. Post-Unmasking Analyses': 'In the event that an analysis team decides to adjust their analysis choice(s) after unmasking, they will document all post-unmasking tests in the paper and commit to showing the original unmasking result along with the post-unmasking update.  \nIn case some novel analyses may require expanded discussions and future work to reach the maturity required for precision cosmology analyses, post-unmasking analyses and discussions are encouraged for this paper.', '4.4. Accommodations': 'In practice, the organizers unmasked submissions from all analysis teams but one in a joint video conference. To encourage submissions from analyses from emerging methods with less extensive validation, the unmasking plot shared in the unmasking video conference showed results from all submissions with anonymized labels (e.g., \'Team Blue" and \'Team Yellow"). Each analysis team was informed only of their own label in advance and the organizers did not share with participants the list of participating analyses. Organizers and participants had agreed on the option for analysis teams to withdraw from the challenge after unmasking to provide anonymity to unexpected results. Ultimately, no analysis team made use of the option to withdraw their submission.  \nThe organizers communicated extensively with analysis teams to facilitate direct comparisons between analyses while preserving the parameter masking. Specifically, the organizers coordinated additional submissions from BACCO and EFT P+B with different cosmology priors: i) Λ CDM without hard θ ∗ constraint (c.f. Sec. 4.1) to enable comparison with the SBI P+B analyses based on Quijote simulations without θ ∗ constraint (Fig. 2), ii) fixing all parameters except σ 8 to the true cosmology of the real-space mocks for comparison with EFT FBI results (Fig. 4). For the latter comparison, truth values were shared after unmasking of the BACCO and EFT P+B submissions for the full parameter space. In addition, the organizers shared the normalized initial conditions of real-space mock box1 and truth values of all cosmological parameters except σ 8 with the EFT FBI team.  \nFeedback from participants suggests these accommodations enabled broader participation and further improved the participants\' experience.', '5. ANALYSIS METHODS': 'The following subsections detail the modeling, inference and parameter unmasking choices of the different analysis teams.  \nNotation -We adopt the notation in Tab. 2 for common variables. For example, with this notation, the Gaussian data likelihood adopted in most inferences is given by  \n-2 log L [ ̂ d | Ω ] = ( m ( Ω ) -̂ d ) T C -1 ( m ( Ω ) -̂ d ) + constant . (1)  \nThroughout, we assume a flat Λ CDM cosmology, and use H ( z ) for the Hubble rate at redshift z and D A ( z ) for the angular diameter distance.  \nHOD models -Analyses of galaxy clustering that extend into the highly nonlinear regime require a model for the galaxy-halo connection that relates galaxies to individual halos, for example HOD models (cf. Sec. 3.2 for extended references). This is in contrast to the EFT bias expansion relating the smoothed galaxy density field to the matter density on large scales.  \nThree of the participating analyses (SBI P+B 5.4, k NN statistics 5.5, density-split clustering 5.6) follow the HOD approach, but with somewhat different implementations and model extensions. Specifically, analysis 5.4 uses one parameterization of the HOD model (see Hahn et al. 2023, for details) while analyses 5.5 and 5.6 use a somewhat different parameterization (see Yuan et al. 2021, for details). For brevity, we refer to the two models as HOD1 and HOD2, respectively, in this overview.  \nFirst, we briefly introduce the vanilla HOD formalism. Statistically, the HOD can be summarized as a probabilistic distribution P ( n g | X h ) , where n g is the number of galaxies of the given halo, and X h is some set of halo properties. Typically the galaxy population is divided into central and satellite populations, and assumes the central galaxy occupation follows a Bernoulli distribution whereas the satellites follow a Poisson distribution.  \nThe vanilla HOD model has been extensively used to describe magnitude-limited galaxy samples (Zheng et al. 2007). It parameterizes the mean galaxy occupation as  \n¯ n cent ( M ) = f ic 2 erfc [ log 10 ( M cut /M ) √ 2 σ ] , (2)  \n¯ n sat ( M ) = [ M -M 0 M 1 ] α , (3)  \nwhere the five vanilla parameters characterizing the model are M cut , M 1 , σ, α, M 0 , augmented by a central galaxy incompleteness parameter f ic , defined to be 0 < f ic ≤ 1 . M cut characterizes the minimum halo mass to host a central galaxy and σ describes the steepness of the transition from 0 to f ic in the number central galaxies. M 0 gives the minimum halo mass to host a satellite galaxy, which is also commonly parameterized as κM cut = M 0 , M 1 characterizes the typical halo mass that hosts one satellite galaxy and α is the power law index on the satellite galaxy occupation.  \nBoth HOD models slightly deviate from this vanilla form. HOD1 fixes f ic = 1 , effectively requiring all massive halos to host a central galaxy. HOD2 varies f ic implicitly by rescaling the predicted number density to match the observed number density. HOD2 also adds a modulation term ¯ n cent ( M ) to the satellite occupation function to largely remove satellites from halos without centrals.  \nThe two models also differ in their ways of determining galaxy positions and velocities once the number of galaxies per halo is computed. In both models, the position and velocity of the central galaxy are set to be the same as those the halo center. However, for satellite galaxies, HOD1 assigns the satellite galaxy positions within the virial radius of the halo following a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997) while the velocities are solved from Jeans equations. In HOD2, the satellite galaxies are randomly assigned to halo particles with uniform weights, each satellite inheriting the position and velocity of its host particle.  \nBoth models also augment the vanilla HOD by adding additional flexibilities (c.f. Hahn et al. 2023; Yuan et al. 2021, for implementation details):  \n- · Velocity bias : The galaxy velocities may not perfectly follow the dark-matter halo and particle velocities. HOD1 rescales galaxy velocities with parameter η cen and η sat , which set the velocity dispersions of central and satellite galaxies relative to halo velocity dispersion: σ cen = η cen σ cen and σ sat = η sat σ sat . In HOD2, we define velocity bias parameters α vel , c and α vel , s to modulate the peculiar velocities of the central and satellite galaxies with respect to the host halo center, respectively. In this definition, α vel , c = 0 and α vel , s = 1 indicates no velocity bias.  \nTable 3. Summary of HOD parameterizations used by the SBI (HOD1), k NN and density split (HOD2) analysis teams.  \n- · Galaxy assembly bias : Galaxy occupation can also depend on secondary halo properties beyond halo mass, a phenomenon commonly referred to as galaxy assembly bias or galaxy secondary bias (See Wechsler & Tinker 2018, for a review). HOD1 adds the halo concentration as a secondary dependency in galaxy occupation via the Heaviside Assembly bias model, which is described in detail in Hearin et al. (2016). HOD2 offers the option of using either the halo concentration or the halo environment in a 5 h -1 Mpc filter as the secondary dependency. Both models adopt two assembly bias parameters to modulate the central and satellite occupations separately.\n- · Baryonic effects : Both HOD implementations account for baryonic effects by modulating the radial distribution of satellite galaxies relative to the halo density profile. HOD1 includes a parameter η conc which sets the ratio between the concentration of satellite and halo profile. HOD2 includes a parameter s that modulates the radial satellite galaxy profile with s = 0 indicating no radial bias, s > 0 indicating a more extended (less concentrated) profile of satellites relative to the halo, and vice versa for negative s .  \nNote that the three HOD-based analyses all employ different priors for the HOD parameters, motivated by sensitivity studies specific to each statistic or based on conservative estimates, as detailed in their respective analysis sections.', '5.1. BACCO Hybrid emulator 3': 'In this section, we describe the results of the Beyond-2point challenge obtained by the BACCO -hybrid emulator approach. This emulator was presented in Zennaro et al. (2021) and Pellejero Ibañez et al. (2023) and has been thoroughly tested on SubHalo Abundance Matching extended techniques (SHAMe, Contreras et al. 2021a,b), and survey-based HOD techniques (Euclid Collaboration et al. 2023; Nicola et al. 2024). The model has been further extended to study Intrinsic Alignments in Maion et al. (2023) and to generate field level predictions in Pellejero-Ibañez et al. (2024).  \n5.1.1. Data & Estimators', 'BACCO hybrid emulator': "Figure 6. Left panel: Power spectrum fit of the BACCO hybrid emulator model to the mean of the 10 Λ CDM mock boxes. The error bars correspond to the Gaussian approximation of the covariance matrix at the volume of one of the boxes. Right panel: Same as the left panel but for the 10 Λ CDM mock boxes in redshift space. We show the three multipoles used in this work.  \n<!-- image -->  \nFigure 7. Left panel: Posterior distribution inferred by the BACCO hybrid emulator model to the mean of the 10 Λ CDM mock boxes in real space as a function of k max . Right panel: Same as the left panel but for the 10 Λ CDM mock boxes in redshift space. We have displaced values by subtracting the mean at the k max = 0 . 5[ h/ Mpc] result to keep the inferred cosmology masked.  \n<!-- image -->  \nThis analysis focuses on the power spectrum P ( k ) for the Λ CDMboxes in real space and the first non-zero multipoles ( ℓ = 0 , 2 , 4 ) of the power spectrum P ℓ ( k ) for the redshift space Λ CDM boxes, shown in Fig. 6. We measure these power spectra by mapping the galaxy positions into a mesh of 1024 3 cells and then perform a Fast Fourier transform as usual applying interlacing and 'triangular shape cloud' as the deposit method.", '5.1.2. Model': "The modeling of the clustering statistics in the so-called 'hybrid' approaches (Modi et al. 2020; Pellejero Ibañez et al. 2022) contains two ingredients: i) a map from Lagrangian, qqq , to Eulerian, xxx , space using the displacement field as measured on N -body simulations, ψψψ ( qqq ) , and ii) a functional relation between the matter and galaxy density fields i.e. a bias model, F ( δ L ( qqq )) .  \nN -body displacements -Regarding the first ingredient, we can define completely the Eulerian overdensity field as:  \nFor the BACCO -Hybrid model ψψψ ( qqq ) is measured in N -body simulations (as opposed to PT) by comparing the Lagrangian position of simulation particles with their position in a snapshot at any desired redshift.  \n1 + δ tr ( xxx ) = ∫ d 3 qw( qqq ) δ D ( xxx -qqq -ψψψ ( qqq )) . (4)  \nGalaxy-halo connection -We employ a 2nd-order Lagrangian bias model (Matsubara 2008; Desjacques et al. 2018):  \nw( qqq ) = F ( δ L ( qqq )) =1 + b 1 δ L ( qqq ) + b 2 ( δ 2 L ( qqq ) -⟨ δ 2 L ( qqq ) ⟩ ) + b s ( s 2 ( qqq ) -⟨ s 2 ( qqq ) ⟩ ) + b ∇ ∇ 2 δ L ( qqq ) . (5)  \nHere δ L ( qqq ) stands for the linear field and s 2 is the traceless part of the tidal field, s 2 = s ij s ij = ( ∂ i ∂ j ϕ ( qqq ) -1 / 3 δ K ij δ L ( qqq )) 2 , with ϕ ( qqq ) the linear gravitational potential. The function w( qqq ) weighs the importance of different Lagrangian fields in representing the density of tracers at a given xxx .  \nRedshift space -We base our redshift-space modelling on our work in Pellejero Ibañez et al. (2022). In order to account for this effect, we substitute ψψψ ( qqq ) with ψψψ s ( qqq ) accounting for the shifts in the line-of-sight direction due to the velocity field as follows,  \nψψψ s ( qqq ) = ψψψ ( qqq ) + ˆ qqq z · vvv tr ( qqq ) aH ˆ qqq z , (6)  \nwhere we define the velocity vvv tr based on the velocities of the N -body simulation. Concretely, vvv tr ( xxx ) = vvv ( xxx ) , the velocity of the matter particle, if the tracer is outside of a halo and vvv tr ( xxx ) = vvv halo ( xxx ) if the tracer is inside of a halo. This distortion roughly accounts for the so-called 'Kaiser effect' (Kaiser 1987) but also incorporates additional contributions resulting from the nonlinearity of the halo velocity. we can model the effect of the intra-cluster velocities onto the redshift-space galaxy field by applying a convolution along the line-of-sight direction δ s tr ( sss ) ∗ z [(1 -f sat ) δ D ( s z ) + f sat exp( -λ FoG s z )] .  \nEmulator -The BACCO hybrid emulator employed a suite of high-resolution simulations (first introduced in Angulo et al. 2021) together with cosmology rescaling (Angulo & White 2010) to densely sample a target cosmological parameter space. Then, we compute the Eulerian fields weighted by their corresponding Lagrangian bias fields. We estimate their power and cross-power spectra and, finally, we use this data to train a neural network. As a result, this provides accurate and extremely fast predictions of the 2-point statistics, which makes it possible to use our model in cosmological data analyses, as presented in Pellejero Ibañez et al. (2023). At this point, we further include two noise terms to account for shot-noise to the power spectrum in the form of Noise = 1 / ¯ n ( ϵ 1 + ϵ 2 k 2 ) .", '5.1.3. Inference': 'Priors -Our cosmological parameter priors are defined by the limits of the emulator. They are given by  \nΩ m ∼ U [0 . 23 , 0 . 4] , σ 8 , c ∼ U [0 . 65 , 0 . 9] , Ω b ∼ U [0 . 04 , 0 . 06] , n s ∼ U [0 . 92 , 1 . 01] . (7)  \nAll beyond Λ CDM parameters are set to their default values. Specifically, neutrino mass is set to zero, and evolving Dark Energy parametrizations are not taken into account. For bias parameters, we use uninformative priors, given by:  \nb 1 ∼ U [0 , 2] , b 2 ∼ U [ -2 , 2] , b s ∼ U [ -3 , 3] , b ∇ ∼ U [ -6 , 6] , (8) λ FoG ∼ U [0 , 1] , f sat ∼ U [0 , 1] , ϵ 1 ∼ U [0 , 2] , ϵ 2 ∼ U [ -4 , 4] .  \nCovariance -We make use of a Gaussian covariance. In real space, this is diagonal and proportional to the amplitude of P ( k ) and shot noise. In redshift space, it is block diagonal in the multipoles with correlations given by the Wigner symbols. We are aware that this covariance is not accurate when pushing to scales as small as k ≈ 0 . 5 h/ Mpc since nonlinearities are not taken into account. However, for this challenge, the number density of mocks is such that the shot noise dominates the signal at scales of k ≈ 0 . 3 h/ Mpc. This shot noise contribution will dominate at scales where nonlinearities become relevant, making their contribution negligible.  \nLikelihood -Due to the central limit theorem, it is safe to assume a Gaussian likelihood shape for the power spectrum in this case.  \nSampling & Validation -We make use of the public code MULTINEST 4 Bayesian inference tool for recovering credibility intervals, see Feroz & Hobson (2008), Feroz et al. (2009), and Feroz et al. (2019) for more details. The best-fit values are also extracted from the maximum likelihood values of these chains. We set the number of live points to 1200 and the evidence tolerance to 0.08.', '5.1.4. Analysis Choices': "Scale cuts -We analyze real- and redshift-space power spectra down to nonlinear scales of k max = 0 . 5 h/ Mpc, motivated by our findings in Zennaro et al. (2021), Zennaro et al. (2022), Pellejero Ibañez et al. (2022), and Pellejero Ibañez et al. (2023). Specifically, in Zennaro et al. (2022) we created thousands of SHAMe mocks with different physical parameters to validate the model and to compute physical priors on the bias parameters, including a mock that closely resembles the hydrodynamical simulation Illustris TNG galaxy clustering (Contreras et al. 2023a). These scale cuts pass all the parameters drift tests described in the unmasking criteria subsection and shown in Fig. 7.  \nUnmasking criteria -We follow the following main criteria before unmasking to assert the robustness of our measurements:  \n- · We verify that the reduced χ 2 value of the best-fit value is ≈ 0.1 for the Λ CDM mocks sample mean. The value χ 2 red = 0 . 1 is determined for the real space mocks, while the value χ 2 red = 0 . 12 was determined for the redshift-space mocks.\n- · We determine if any of the parameter posteriors intersect the priors' limits. This is remedied by extending the prior range for the nuisance parameters. In the case of cosmological parameters, however, our emulator provides the priors. For example, we discovered that the σ 8 -b 1 degeneracy broadens the constraints and partially affects the real space emulator priors in the 2 -3 σ region. However, the 1 -2 σ region is well within our emulator's priors, so we believe these results to be reliable.\n- · We investigate the consistency of the inferred parameters with respect to various k max values. In particular, we compute the evolution of the posterior distribution at k max = { 0 . 1 , 0 . 2 , 0 . 3 , 0 . 4 , 0 . 5 } h/ Mpc and find no discernible change, as illustrated in Fig. 7.\n- · We test for redshift-space analysis whether adding or removing multipoles gives consistent cosmological parameter constraints. We find no tensions between results containing only monopole, monopole plus quadrupole, or monopole plus quadrupole and hexadecapole.\n- · We include our theory error budget based on the emulator uncertainties. This is discussed in Appendix A of Pellejero Ibañez et al. (2023). For the test, we add these error estimates in quadrature to the covariance matrix diagonal elements. We find no strong dependence of the recovered cosmological marginalized values on this theory error budget. We note that the theory error should always be included, however, if this error is small compared to the data error, its effect becomes negligible.  \nCaveats -We caution that several assumptions could potentially bias our constraints: First, the choice of a Gaussian covariance with no off-diagonal elements. Second, the possible underestimation of the emulator errors at the cosmology of interest. Third, the priors of the emulator. Even though we expect the true cosmologies to lie within these priors,  \nstrong degeneracies, such as the σ 8 -b 1 direction in real space, might make our priors too informative. Fourth, possible model simplifications. To construct the BACCO emulator hybrid model, we make use of several approximations in the galaxy-to-matter and in the galaxy-velocity assignments. Specifically, the current galaxy velocity model lacks implementation of the mass-dependent Finger-of-God effect and strong velocity bias. These model misspecifications might bias our cosmological parameter estimations.  \nPost-unmasking studies -After unmasking, we found that our results for σ 8 exhibited a 1 σ shift compared to the true values. This is not entirely unexpected, as cosmic variance can introduce shifts in the contours, and locating them within a 1 -3 σ range does not imply a failure of the model. Additionally, our assumption of a diagonal Gaussian covariance matrix exacerbates this issue by ignoring non-diagonal correlation terms that arise on small scales.  \nHowever, to mitigate the impact of cosmic variance, parameter estimation was performed using the mean of 10 independent mocks. Naively, this should reduce the contour size by a factor of approximately √ 10 ∼ 3 . Taking this into account, our contours appear to be at approximately ∼ 3 σ from the true values. To confirm this, it would be necessary to reanalyze the data with reduced error bars corresponding to a volume 10 times larger than each individual mock. Unfortunately, all tests conducted on the BACCO hybrid model in redshift space thus far have been based on the assumption of a BOSS-like sample (Pellejero-Ibanez et al. 2022; Pellejero Ibañez et al. 2023) with an effective volume of V eff = 2 . 8Gpc 3 ∼ 0 . 9Gpc 3 /h 3 and our emulator's noise level remains too high for the galaxy sample and survey volume of this challenge.  \nNevertheless, we can investigate whether incorporating or removing theory errors impacts our predictions. As previously mentioned, we already conducted this test by including errors in quadrature. Based on Pellejero Ibañez et al. (2023), we are aware that typical uncertainties in the emulator account for ∼ 0 . 5% in the monopole, ∼ 1% in the quadrupole, and ∼ 10% in the hexadecapole amplitudes. During this examination, we did not realize the low values of the hexadecapole, as depicted in Figure 7. Consequently, the theory error we included in the noisiest of our estimations was underestimated.  \nTherefore, we decided to present the measurement of cosmological parameters without the hexadecapole, as shown in dashed lines in Figure 7. The shift found is of around half a sigma and the slight decrease in constraining power deems our results unbiased with respect to the true values. It is important to reiterate that this is a post-unmasking finding, and further tests on the model are required to rule out potential failures when studying the HOD models employed in this work.", '5.2. EFT P & B: Analysis Method 5': "In this section, we describe the analysis of the Beyond 2-point mocks with the Effective Field Theory of Large Scale Structure, focusing on the large-scale power spectrum and bispectrum. The key theoretical underpinnings of this are described, e.g. in Baumann et al. (2012); Carrasco et al. (2012); Ivanov (2023) 6 , with recent applications to galaxy surveys and other large-scale structure data shown in e.g., Ivanov et al. (2020a,b); Philcox & Ivanov (2022); D'Amico et al. (2020); Chen et al. (2022). We apply this methodology to all challenge catalogs.  \nInformation Sources -We briefly discuss the various sources of information in the full shape of the galaxy power spectrum, which set the parameter degeneracy directions. More detail on the sources of information can be found in Ivanov et al. (2020a). From the shape of the galaxy power spectrum, one can determine ω b , ω m ≡ Ω m h 2 , and n s regardless of the distance to the sample. These parameters then predict absolute scales such as that of matterradiation-equality and the sound horizon ( k eq and k BAO respectively).  \nIn contrast, distance information is encoded in the angular scales θ ∥ and θ ⊥ , for both BAO and equality. This constrains the Hubble parameter h , and, if multiple redshifts are available, Ω m through the growth rate evolution. For periodic boxes, there is no distance information, since all quantities are measured in h -1 Mpc units, thus all distances are equivalent to H -1 0 .  \nFrom growth, we constrain b 1 σ 8 ( z data ) from the real-space power spectrum or b 3 1 σ 4 8 ( z data ) from the bispectrum. Loop corrections additionally yield b 2 1 σ 4 8 ( z data ) and infrared resummation directly measures σ 8 ( z data ) (see Senatore & Zaldarriaga 2015; Baldauf et al. 2015; Blas et al. 2016a,b; Ivanov & Sibiryakov 2018; Vasudevan et al. 2019), though these effects are comparatively small at the high redshifts of this challenge. In redshift space we instead measure  \nb 1 σ 8 ( z data ) from the monopole and fσ 8 ( z data ) from the quadrupole and hexadecapole power spectrum moments. Loops and the bispectrum monopole help break degeneracies and give directly f ( z data ) and σ 8 ( z data ) .", '5.2.1. Data & Estimators': "In all cases, we analyze the large-scale power spectrum and bispectrum extracted from the mock catalogs. We use P ( k ) and B ( k 1 , k 2 , k 3 ) in real space; in redshift space we consider the power spectrum monopole, quadrupole, and hexadecapole ( P 0 , P 2 , and P 4 ), but restrict to the bispectrum monopole (since higher moments do not add significant signal, see Ivanov et al. 2023). For the redshift-space mocks, we supplement the data vector with the real-space power spectrum proxy, Q 0 , which gives additional information on smaller scales without bias from fingers-of-God (Ivanov et al. 2022b; see also Scoccimarro 2004, D'Amico et al. 2024).  \nFor the periodic box mocks, correlators are estimated using Fast Fourier Transforms, as usual (Scoccimarro 2015). For the light-cone mock catalog, we use the window-free estimators of Philcox (2021a,b), as implemented in the Spectra-Without-Windows code 7 , making use of the mask and random file. In this case, we split the sample into two redshift bins of equal density, with effective redshifts z 1 = 0 . 92 and z 2 = 1 . 17 . Knowledge of redshift evolution is necessary to break the geometric degeneracy.", '5.2.2. Model': "Perturbation theory -We model the power spectrum and bispectrum with the Effective Field Theory of Large Scale Structure, as implemented in CLASS-PT (Chudaykin et al. 2020) 8 The power spectrum multipoles are modeled using the (infrared-resummed) one-loop theory (Ivanov et al. 2020a), and we use tree-level theory for the bispectrum (Ivanov et al. 2022a), noting that higher loops give limited gains (Philcox et al. 2022; D'Amico et al. 2022b).  \nGalaxy-matter connection -We assume the following bias expansion up to third-order (renormalized) operators (Ivanov et al. 2020a):  \nδ g = b 1 δ + 1 2 b 2 δ 2 + b G 2 G 2 + b Γ 3 Γ 3 , (9)  \nwhere δ g and δ are the galaxy and matter overdensities, G 2 and Γ 3 are Galileon tidal operators. The model additionally includes stochastic contributions from shot-noise P shot , A shot , B shot , a 0 , a 2 (including k 2 scale dependence), and counterterms { c 0 , c 2 , c 4 , ˜ c } , encapsulating small-scale physics such as halo formation and velocity effects. This makes no assumptions on the form of the galaxy-halo connection, except that it is statistically isotropic and homogeneous on large scales, and obeys Einstein's equivalence principle (Desjacques et al. 2018).", '5.2.3. Inference': 'Priors -We assume the following priors on cosmological parameters for all analyses:  \nω cdm ∼ U [ -∞ , ∞ ] , 10 9 A s ∼ U [0 . 5 , 5] , n s ∼ U [0 . 87 , 1 . 07] , ω b ∼ U [0 . 01 , 0 . 035] , (10)  \nb 1 ∈ flat [0 , 4] , b 2 ∼ N (0 , 1 2 ) , b G 2 ∼ N (0 , 1 2 ) , b Γ 3 ∼ N ( 23 42 ( b 1 -1) , 1 2 ) , (11)  \nwith the angular size of the sound horizon θ ∗ fixed to the value known to all participants. We fix the neutrino mass to zero and assume a flat Universe in all cases. For bias parameters, we use weakly informative priors, given by Philcox & Ivanov (2022):  \nwhere N ( µ, σ 2 ) indicates a Gaussian distribution with mean µ and variance σ 2 . Note that these priors are consistent with the recent measurements from Ivanov et al. (2024), as well as earlier results of Lazeyras & Schmidt (2018); Abidi & Baldauf (2018). Similarly, we use Gaussian priors for the counterterms c 0 , c 2 , c 4 , c 1 , ˜ c and stochasticity parameters a 0 , a 2 , P shot , A shot , B shot :  \nc 0 [ Mpc /h ] 2 ∼ N (0 , 30 2 ) , c 2 [ Mpc /h ] 2 ∼ N (30 , 30 2 ) , c 4 [ Mpc /h ] 2 ∼ N (0 , 30 2 ) , ˜ c [ Mpc /h ] 4 ∼ N (500 , 500 2 ) , c 1 [ Mpc /h ] 2 ∼ N (0 , 51 2 ) , P shot ∼ N (0 , 1 2 ) , a 0 ∼ N (0 , 1 2 ) , a 2 ∼ N (0 , 1 2 ) , 2 2 (12)  \nB shot ∼ N (1 , 1 ) , A shot ∼ N (0 , 1 ) ,  \nwhere we use the same convention 9 as Ivanov et al. (2022a); Philcox & Ivanov (2022); Philcox et al. (2022), but with ¯ n ≈ 5 × 10 -4 h 3 Mpc -3 of the challenge boxes.  \nCovariance -We assume a Gaussian covariance matrix for the power spectrum and multipoles bispectrum monopole, defined explicitly in Chudaykin & Ivanov (2019); Ivanov et al. (2022a, 2023). This is exact in the linear regime (where modes are uncorrelated), and found to be highly accurate on our scales of interest, due to the high shot-noise, and limited cosmological information available at highk (Wadekar et al. 2020). The covariance is diagonal in k , but includes correlations between different power spectrum multipoles, where necessary. We do not include cross-covariance between the power spectrum and bispectrum, since this is formally of higher order in EFT, and found to be unnecessary in Ivanov et al. (2022a) for a similar choice of scale cuts.  \nThe power spectrum covariance is computed using the specific realizations of the power spectrum multipoles (since these are measured at high significance). For the bispectrum covariance, we adopt an iterative procedure, first computing the covariance using a fiducial cosmology, then updating with the best-fit parameters from a likelihood analysis. This procedure is found to converge quickly.  \nLikelihood -On large-scales, the data is well described by a Gaussian likelihood. This holds by the central limit theorem, and also follows from perturbation theory. The covariance can be computed precisely on large scales, and the number of bins is reasonable.  \nSampling & validation -We sample the various parameters with MontePython (Brinckmann & Lesgourgues 2019), analytically marginalizing over nuisance parameters that enter the model linearly (Philcox et al. 2021a). This model has been validated on several suites of simulations, including MultiDark-Patchy, Las Damas, the Perturbation Theory Challenge, Nseries, and Outer Rim (Ivanov et al. 2023, 2020a; Philcox & Ivanov 2022; Nishimichi et al. 2020; Philcox et al. 2022; Ivanov 2021; Ivanov et al. 2022b; Chudaykin & Ivanov 2023; Ivanov et al. 2023). We do not include a theoretical error covariance (due to our choice of scale cuts), but note that this can be, in principle, incorporated (Baldauf et al. 2016; Chudaykin & Ivanov 2019; Chudaykin et al. 2021b). We consider our MCMC chains converged when they satisfy the Gelman-Rubin criterion R -1 < 0 . 03 .', '5.2.4. Analysis Choices': 'Scale cuts -The EFT P+B analyses adopt the following scale cuts (in h Mpc -1 units):  \n- · real-space mock (periodic box, z = 1 , V = 8 (Gpc/ h ) 3 × 10 ): k P ≤ 0 . 3 , k B ≤ 0 . 15\n- · redshift-space mock (periodic box, z = 1 , V = 8 (Gpc/ h ) 3 × 10 ): k P ℓ ≤ 0 . 20 , 0 . 20 < k Q 0 ≤ 0 . 4 , k B 0 ≤ 0 . 08 .\n- · light-cone mock ( z = [0 . 8 , 1 . 3] , V = 5 . 3 (Gpc/ h ) 3 (for our fiducial cosmology Ω m = 0 . 31 , h = 0 . 676 )): k P ℓ ≤ 0 . 25 , 0 . 25 < k Q 0 ≤ 0 . 4 , k B 0 ≤ 0 . 08  \nThese are motivated by previous works and tests on large-volume simulations (Ivanov et al. 2020a; Philcox & Ivanov 2022; Nishimichi et al. 2020; Ivanov et al. 2022a,b) and validated for this challenge through parameter drift plots for all challenge mocks. As an example, we show parameter posteriors for the redshift-space mock for four different choices of k P ℓ max = 0 . 15 , 0 . 2 , 0 . 25 , 0 . 3 h Mpc -1 in Figure 8. The parameter drifts in this plot suggests that results at k P ℓ max = 0 . 3 h Mpc -1 are biased, because one can clearly see an upward shift of σ 8 by 1 σ , consistent with the effects of two-loop corrections previously observed, e.g. in Nishimichi et al. (2020); Chudaykin et al. (2021b). The posteriors for the smaller scale cuts appear consistent with each other, and we select k P ℓ max = 0 . 2 h Mpc -1 as a baseline choice in order to be conservative. The same methodology is applied for the light-cone and real-space analyses, i.e. studying parameter drifts plus using estimates for the theoretical error based on the two-loop estimates (see also (Chudaykin et al. 2021b; Ivanov et al. 2022a) for studies of real-space clustering in EFT). Using this methodology, all challenge mock catalogs were systematically analysed and tested (see consistency checks below). The particular analysis choices applied here were also cross-validated on other mock catalogs such as Nseries, PT Challenge, and Las Damas. We choose a more aggressive scale cut for the light cones as their volume is small compared to the periodic box simulations, and hence the theory systematic error is smaller than cosmic variance down to smaller scales. In this case, the parameter drift  \nplot is also less conclusive because it is hard to disentangle parameters shifts due to bias from those resulting from a statistical fluctuation due to the addition of new data at a higher k max .  \nWe also note that the bispectrum is restricted to larger scales due to the tree-level modeling. The lower k max for redshift-space snapshots (compared to real space ones) is also supported by expectations of the impact of velocity effects which become non-perturbative at larger scales (noting that the nonlinear scale is k NL ∼ σ -1 FoG in redshift space, as opposed to R -1 halo in real space, if halo formation effects dominate), see Ivanov (2021); Ivanov et al. (2022b) for detailed discussions.  \nConsistency checks -We have run a number of consistency checks that include: variation of the bispectrum covariance matrix, stability w.r.t. scale cut changes, exclusion of Q 0 , B 0 , P 4 from the datavector, and the change of fiducial cosmology in the case of the light-cone mock. We have not found any significant effect when perturbing our baseline choices within a reasonable range.  \nAs far as the choice of fiducial cosmology is concerned, we have explicitly tested that using Ω m = 0 . 31 (our baseline choice), or Ω m from the best-fit to the data has negligible difference on the extracted cosmological parameters, even when Ω m is actually significantly different (by up to ∼ 30% ) from its fiducial value 0.31. In this particular example of a 30% difference (which constitutes ≃ 10 σ Ω m ), the shift in the mean of Ω m was found to be 0 . 1 σ Ω m , while the shift in σ 8 was found to be 0 . 3 σ σ 8 . This implies that the standard Alcock-Paczynski parametrization is quite accurate even for large deviations from the fiducial cosmological parameter values.  \nUnmasking criteria -We allow our results to be unmasked after visual confirmation of the relevant posteriors and their degeneracy directions, as well as simple tests based on the goodness of fit and analysis of variation with k max .  \nFor the real-space and redshift-space mocks we additionally analyzed the mean data vectors (from the 10 realizations) with the covariances that correspond to the full volume of 10 boxes. We found some non-negligible bias for k P ℓ max = 0 . 25 h Mpc -1 for the redshift-space mock in this analysis, which was another reason to stick to our baseline choice k P ℓ max = 0 . 2 h Mpc -1 . For the light-cone mock, however, we only have one realization for the volume which is approximately 40% smaller than that of one snapshot. The two-loop systematic error is smaller than the cosmic variance in this case even at k P ℓ max = 0 . 25 h Mpc -1 , and therefore we could adopt this more optimistic scale cut in our analysis.  \nCaveats -Finally, we caution that failure modes of our analysis could include: overly optimistic scale-cuts (for example due to extremely large Fingers-of-God effects), cosmological or nuisance parameters lying outside the parameter ranges, or galaxy-halo connections violating our symmetry assumptions. The posteriors should be interpreted with caution, as for many cases (e.g. real-space analyses) their shape is highly non-Gaussian. This may lead to parameter projection effects: i.e. 1d marginalized contours may be off due to marginalization over non-Gaussian posteriors, as it was discussed in detail e.g. in Ivanov et al. (2020a); Chudaykin et al. (2021a); Philcox & Ivanov (2022); Ivanov et al. (2023). Note that our choice of k max = 0 . 3 h Mpc -1 for the real-space analyses is somewhat aggressive (done in part in order to facilitate the sampling), so we expect (based on our tests on mock data) that the theory systematic error may be as large as ∼ 1 σ .', '5.3. EFT FBI: Analysis Method 10': "In this section, we review the Effective Field Theory Field-level Bayesian Inference (EFT FBI). EFT FBI directly models and optimally extracts cosmological information from the three-dimensional galaxy density field within the EFT framework. This requires an explicit sampling of all Fourier modes of the initial conditions-in our case, corresponding to 90 3 ≃ 730 , 000 parameters-in addition to bias, stochastic, and cosmological parameters. We refer readers to Schmidt et al. (2019) for the first EFT FBI groundwork and Nguyen et al. (2024) for the most recent advance. We further note parallel efforts on cosmological field-level inference from galaxy clustering 11 , within (Schmittfull et al. 2019; Cabass 2021; Schmittfull et al. 2021; Cabass et al. 2024) and outside of the EFT framework (Seljak et al. 2017; Ramanah et al. 2019; Andrews et al. 2022; Bayer et al. 2023), as well as the most recent advances within the community of machine learning for cosmology (e.g Cuesta-Lazaro & Mishra-Sharma 2023; Doeser et al. 2023; Ho et al. 2024).  \nOur EFT FBI analysis uses LEFTfield , a Lagrangian, EFT-based forward modeling of cosmological density fields. The code was first introduced in Schmidt (2021a) with extensions and validations later shown in Babić et al. (2022); Kostić et al. (2023); Stadler et al. (2023); Tucci & Schmidt (2023); Nguyen et al. (2024).  \nFigure 8. Λ CDM redshift-space mock analysis: parameter constraints from EFT P+B for four choices of the power spectrum scale cut k P ℓ max . All parameters are normalzied to fixed values taken from the best-fit at k P ℓ max = 0 . 2 h Mpc -1 (our fiducial scale cut). Physical baryon density ω b and some nuisance parameters are not displayed.  \n<!-- image -->  \nIn this challenge, we analyze the real-space mocks. For the clustering of galaxies in their comoving rest frame, i.e. real space, both the matter δ and galaxy density field δ g contain advection contributions that involve the Lagrangian displacement field. Since matter and galaxies co-move on large scales following Einstein's equivalence principle, the displacement is the same for galaxies and matter. This direct consequence of the equivalence principle therefore ensures that the advection contribution to δ g can be uniquely predicted (Desjacques et al. 2018), hence the biasσ 8 degeneracy breaking.", '5.3.1. Data and Estimators': 'Our data vector is the entire three-dimensional galaxy density field, namely ̂ d ≡ δ g , at z = 1 . 0 . Similarly to Nguyen et al. (2024), we construct a filtered galaxy density field δ g, Λ from the galaxies in each catalog by a nonuniform-touniform Fourier transform (NUFFT).', '5.3.2. Model': 'The LEFTfield forward model we employ in this challenge involves a Gravity model evolving the initial conditions ˆ s to the Eulerian matter density field and a Galaxy model connecting the Eulerian matter density field to the model prediction for the data vector δ g, Λ .  \nWe note that our pre-unmasking analyses in the challenge had concluded before the analysis in Nguyen et al. (2024) started-where the latter analyzed different datasets than this challenge, consisting of mass-selected main halos. Hence, while both works use the same Gravity model , they differ in the Galaxy model . Specifically, our baseline analyses here use a more restricted model compared to Nguyen et al. (2024). We will highlight the differences below and revisit some of the main differences in our post-unmasking study.  \nGravity model -We predict gravitational evolution using the n -th order Lagrangian Perturbation Theory (LPT) framework (Matsubara 2015), as described in Sec. 2 and App. A-B of Schmidt (2021b). The convergence of this n -LPT scheme was demonstrated by Schmidt (2021b) in their Sec. 6, up to n = 6 . In this challenge, we adopt n = 2 , i.e. the 2LPT model (Buchert 1992; Bouchet et al. 1995).  \nGalaxy model -We describe the matter-galaxy connection with the EFT bias expansion:  \nδ g, Λ = ∑ O b O O [ δ Λ ] (13)  \nwhere O are operators constructed out of the filtered Eulerian matter density field δ Λ , together with their associated galaxy bias coefficients b O . Similar to Nguyen et al. (2024), we expand Eq. (13) in the Eulerian basis. However, unlike Nguyen et al. (2024), here we do so only up to second-order operators:  \nwhere  \nO ∈ [ δ, δ 2 -⟨ δ 2 ⟩ , K 2 -⟨ K 2 ⟩ , ∇ 2 δ ] , (14)  \nK 2 ≡ ( K 2 ij ) = ([ ∂ i ∂ j ∇ 2 -1 3 δ ij ] δ ) 2 (15)  \nis the gravitational tidal field squared, and ∇ 2 δ is the leading higher-derivative operator. Note that these coefficients are equivalent to the subset of coefficients b 1 , b 2 , b G 2 , c 0 in the EFT P+B model (Sec. 5.2). Compared to Eq. (14), Nguyen et al. (2024) add the full set of four third-order bias operators, while the EFT P+B team adds a single thirdorder bias term ( Γ 3 ) that is relevant for the one-loop power spectrum. In our EFT-convergence test [Fig. 10] and post-unmasking study [Fig. 22], we additionally consider the full set of third-order bias terms and quantify the impact of on our σ 8 constraints [Fig. 4].', '5.3.3. Inference': 'We sample and estimate the marginalized posterior of the amplitude rescaling parameter α = σ 8 /σ fiducial 8 given by  \nP ( α ∣ ∣ ∣ ̂ d Λ ) = ∫ D ˆ s ∫ dσ ϵ ∫ d { b O } P ( α, { b O } , σ ϵ , ˆ s ∣ ∣ ∣ ̂ d Λ ) ∝ ∫ D ˆ s ∫ dσ ϵ ∫ d { b O } L ( ̂ d Λ ∣ ∣ ∣ α, { b O } , σ ϵ , ˆ s ) P ( α ) P ( { b O } ) P ( σ ϵ ) P (ˆ s ) . (16)  \nwhere the parameters to be explicitly sampled and numerically marginalized over are Φ ≡ { b O , σ ϵ , ˆ s } . Here, σ ϵ denotes the galaxy stochasticity, further described in the Covariance section.  \nThe normalized initial conditions ˆ s (sometimes referred to as \'phases\'), which are drawn from a unit Gaussian prior, are related to the linear density field via  \nδ (1) ( k , z ) = α [ N 3 grid L 3 P L ( k, z ) ] 1 / 2 ˆ s ( k ) . (17)  \nIn Eq. (17), P L ( k, z ) is the linear matter power spectrum in the fiducial cosmology, while L, N grid are the length and grid size of the ˆ s ( k ) grid, respectively. Note that the finite volume L 3 and finite cutoff Λ imply that a finite grid is sufficient to explicitly represent every relevant Fourier mode of the initial conditions.  \nLikelihood -Following Schmidt et al. (2019); Cabass & Schmidt (2020); Schmidt et al. (2020); Schmidt (2021a); Kostić et al. (2023); Nguyen et al. (2024), we adopt a Gaussian likelihood in Fourier space,  \n̸  \nln L ( ̂ d Λ ∣ ∣ ∣ α, ˆ s, { b O } , σ ϵ ) = -1 2 k max ∑ k =0 [ ln 2 πσ 2 ϵ + 1 σ 2 ϵ ∣ ∣ ∣ ̂ d Λ ( k ) -δ g, Λ [ α, ˆ s, { b O } ]( k ) ∣ ∣ ∣ 2 ] + const. , (18)  \nwhere we have used the fact that the galaxy stochasticity is homogeneous and isotropic in the galaxy comoving rest frame, i.e. in real space, and hence has a diagonal covariance matrix in Fourier space C = σ 2 ϵ . Note that Eq. (18) introduces the analysis cutoff scale k max where k max ≤ Λ . Unlike Nguyen et al. (2024), here we choose a cubic sharpk filter to implement the analysis cutoff k max in Eq. (18) and further set k max = Λ / 1 . 4 . A value of k max < Λ is chosen to reduce the magnitude of higher-derivative terms that are controlled by the length scale Λ -1 . We have neither performed a detailed study of the optimal choice of k max / Λ here nor in Nguyen et al. (2024).  \nEq. (18) does not capture a coupling of galaxy stochasticity to large-scale perturbations (\'density-dependent noise\'), or non-Gaussianity of the noise. While the former can be incorporated via a real-space formulation of the likelihood (Cabass & Schmidt 2020), the latter is technically more difficult to incorporate at the field level. Apart from the bias terms (see Galaxy model above), this constitutes the second main model difference to the EFT P+B analysis.  \nCovariance -Even when limiting to a diagonal covariance in Fourier space, galaxy stochasticity is not necessarily white noise on all scales. In fact,  \nσ ϵ = σ ϵ ( k ) = σ ϵ, 0 [ 1 + σ ϵ,k 2 k 2 + σ ϵ,k 4 k 4 + · · · ] , (19)  \nin full generality (Desjacques et al. 2018). Though physically motivated, σ ϵ,k 2 and σ ϵ,k 4 add significant parameter degeneracies and correlations. Therefore, unlike Nguyen et al. (2024), herein we consider only the scale-independent leading-order contribution 12 in our baseline analyses to minimize sampling cost. 13 Eq. (19) then reduces to σ ϵ ( k ) = σ ϵ, 0 . This still leaves open two possibilities: (1) fixing σ ϵ, 0 = σ ϵ Poisson where the latter follows the expectation for Poisson shot noise, and (2) inferring σ ϵ, 0 from the data itself. We choose option (1) in our baseline analyses, and consider the subleading correction ∝ k 2 , as employed in Nguyen et al. (2024), in a post-unmasking analysis described below.  \nPriors on cosmology and bias -We adopt wide uniform priors for the amplitude rescaling parameter α and the bias parameters b O . To summarize,  \nP ( α ) = U (0 . 5 , 1 . 5) , (20)  \nP ( b δ ) = U (0 . 0 , 4 . 0) , P ( b δ 2 ) = P ( b K 2 ) = U ( -4 . 0 , 4 . 0) , P ( b ∇ 2 δ ) = U ( -20 . 0 , 20 . 0) , P ( σ ϵ ) = δ D ( σ ϵ Poisson ) . (21)  \nPriors on initial conditions -As mentioned above, the cosmological prior on the normalized initial conditions ˆ s , in the absence of primordial non-Gaussianity, is a unit Gaussian prior (free-IC). However, we will also consider a special analysis, on box 1 , where we were given the ground-truth initial conditions used to initialize one of the N-body simulations underlying the challenge data set (fixed-IC). To summarize, we assume the following priors on ˆ s :  \nP (ˆ s ) =    δ D (ˆ s -ˆ s true ) (fixed-IC, box 1 ) , N (ˆ s ; 0 , I ) (free-IC) , (22)  \nwhere, in a fixed-IC analysis, δ D indicates the Dirac delta distribution fixing the initial conditions to the true initial conditions ˆ s true of the simulation data while, in a free-IC analysis, N indicates the multivariate normal distribution on normalized initial conditions with zero mean and identity covariance.  \nFigure 9. Consistency between posteriors obtained in ten EFT FBI free-IC analyses on ten data realizations during stage 2. To avoid unmasking the values of bias parameters, we subtract and divide their values in each chain by the corresponding posterior means obtained from the analysis of box 1 .  \n<!-- image -->  \nSampling -We employ Hamiltonian Monte Carlo sampling to explore the initial conditions ˆ s while adopting slice sampling to sample other parameters. We refer readers to Sec 3.2 of Kostić et al. (2023) and references therein for the motivation behind this choice.  \nMCMC convergence and autocorrelation -Unconverged and correlated MCMC chains lead to biases and uncertainties in the estimation of moments of the parameter posteriors. Below, we focus on MCMC convergence diagnostics and autocorrelation estimates for stage-2 analyses (see Unmasking criteria below), as their results can be directly compared to those obtained by other teams, namely EFT P+B and BACCO P.  \n- · MCMC convergence-We run two MCMC chains with distinct initial parameter values for each analysis. We first identify the warmup and equilibrium phases in each chain through the parameter drifts within individual chains and the Gelman-Rubin (G-R) statistics estimated from each and both chains. 14 After removing the warmup phases, we find posteriors estimated from the two chains to be completely consistent.\n- · MCMC autocorrelation-The univariate autocorrelation function (ACF) provides an estimate for the number of effective, i.e. independent, samples n effective in the MCMC chain for the parameter of interest. In the upper right corner of Fig. 9, we show the ACFs of the cosmological and bias parameters as a function of the lag between two samples in the MCMC chain. 15 A null ACF indicates that two MCMC samples are independent draws from the posteriors.\n- · MCMC Effective Sample Size (ESS)-We target a G-R value of ˆ R ∼ 0 . 01 and an average number of ∼ 255 independent samples of α . 16 The MCMC sampling error on the means of α and σ 8 can be estimated from σ α /n effective α where n effective α is the ESS of α . The ESS of ten stage-2 analyses range from 193 to 311, with a mean\n- 14 We observe neither divergence nor slow mixing in MCMC chains across all pre-unmasking and post-unmasking analyses.  \n<latexit sha1\\_base64="1TDkQcakYvvVOOUEsij8jTezlzw=">AAACX3icdVFdS8MwFE3rx3TqrPokvhSHoOBGKzJ9EURffBEUnBuss6RZuoYlbUluxVH6J30TfPGfmG1V/LwQOJxzz83NSZBypsBxXgxzbn5hsbK0XF1ZXautWxub9yrJJKFtkvBEdgOsKGcxbQMDTruppFgEnHaC0eVE7zxSqVgS38E4pX2BhzELGcGgKd969BQbCuyfepyGsD/yc09giKTIBX4qCk+yYQQHjVnXw6eoAA9pwy2K/4xnTtP1DqMP7jolxUOuDeVA36o7TWda9m/glqCOyrrxrWdvkJBM0BgIx0r1XCeFfo4lMMJpUfUyRVNMRnqtnoYxFlT182k+hb2nmYEdJlKfGOwp+9WRY6HUWAS6c7Kv+qlNyL+0XgbhaT9ncZoBjcnsojDjNiT2JGx7wCQlwMcaYCKZ3tUmEZaYgP6Sqg7B/fnk3+D+qOm2mq3b4/r5RRnHEtpBu2gfuegEnaMrdIPaiKBXwzRWjFXjzayYNdOatZpG6dlC38rcfge9s7gO</latexit>  \nFigure 10. Convergence of σ 8 ( z = 0) inferred in our stage-1 and stage-2 analyses, using the same data realization ( box 1 ). In both panels, the error bars indicate 68% lower and upper limits on α . All σ 8 values are quoted relative to the stage-1 analysis using the fiducial second-order bias model at k max = 0 . 1 . We emphasize that the horizontal dashed line does not indicate the ground truth σ 8 . The horizontal bands indicate the 68% and 95% confidence intervals of the σ 8 posterior in that scenario. Left (Stage 1): Convergence of fixed-IC constraints as a function of the analysis (EFT) cutoff scale k max ( Λ ). All σ 8 posterior means are within 1σ of the corresponding value at k max = 0 . 1 h Mpc -1 for the respective bias model (same color). Magenta and blue points in the left panel indicate second-order and third-order bias model, respectively. The blue points are horizontally displaced to make error bars fully visible. Right (Stage 2): Free-IC constraint from box 1 , again relative to the stage-1 result.  \n<!-- image -->  \nand median of 254.4 and 259.5, respectively. Their G-R statistics range from 0.003 (best) to 0.035 (worst), with a mean of 0.0115 and median of 0.0085.  \nPosterior consistency -A unique feature of the EFT FBI (stage-2) analyses is that we individually analyze all ten data realizations in the real-space suite of mocks. Fig. 9 shows the posteriors obtained in all ten analyses, with their labels hidden to avoid revealing the cosmic variance in each realization. The main takeaway here is that all ten posteriors are consistent with each other within the 68% confidence limit.', '5.3.4. Analysis choices': "Scale cuts -We determine the fiducial scale cut k max in our baseline analyses by choosing the maximum scale that maintains convergence of the amplitude rescaling parameter α inferred with different choices for the Galaxy model , e.g. second- versus third-order galaxy bias expansion (Fig. 10), Lagrangian versus Eulerian basis for the bias expansion. Specifically, we define the relative parameter shift  \n∆ α/σ α ≡ ⟨ α ⟩ A -⟨ α ⟩ B ( σ 2 α A + σ 2 α B ) 1 / 2 , (23)  \nwhere the indices A,B label different analyses with their associated posterior means ⟨ α ⟩ and 68% uncertainties σ α . We require Eq. (23) to be less than 2.0 for final submissions, i.e. the LEFTfield result shown in Fig. 4. The scale cut we end up with for our stage-1 and stage-2 submissions is k max = 0 . 1 h Mpc -1 . This then corresponds to the maximum scale up to which the second-order bias expansion is expected to be reliable. We anticipate that a third-order bias expansion will have higher reach in wavenumbers. For more details, see Nguyen et al. (2024).  \nUnmasking criteria -We unmask the EFT FBI analysis in two stages:  \n- 1. In stage 1, we performed the inference on box 1 , with the initial conditions ˆ s being fixed to ˆ s true (fixed-IC) provided by the challenge organizers (only for box 1 ). During stage 1, we performed several EFT- and numericalconvergence tests because analyses are numerically cheap (at fixed initial conditions). After stage-1 submission,  \n- and before stage-2 analysis started, the organizers provided us a confirmation that stage-1 posterior mean of σ 8 is within 2σ of the true σ 8 (without specifying which side of the true σ 8 it actually is). For transparency, below in Fig. 10, we also highlight the 2σ constraint of the stage-1 submission.\n- 2. In stage 2, without any change in our analysis choices, we extended the analysis to all ten boxes, while also allowing ˆ s to be explicitly sampled and marginalized over (free-IC). We emphasize that stage-2 analyses were performed without the true initial conditions ˆ s true of the simulation boxes.  \nIn summary, the EFT FBI analyses and results are conducted and submitted in two stages, where σ 8 is fully masked, but ˆ s is provided for box 1 in stage 1. The reason behind this choice is we initially agreed with the challenge organizers to conduct a fixed-IC analysis only. After the submission of our fixed-IC analysis (stage 1), encouraged by the organizers, we decided to pursue a full free-IC analysis (stage 2), in coordination with the BACCO P and EFT P+B teams.  \nCaveats -In early tests or parallel studies, e.g. Nguyen et al. (2021, 2024), we have identified an issue of ' σ ϵ collapse' when trying to infer the galaxy stochasticity, parametrized by σ ϵ, 0 , from the data itself with wide uniform prior on σ ϵ, 0 : the inferred σ ϵ, 0 (hence σ ϵ ) tends to an unphysically low value. This is a known issue of inference and sampling using a Gaussian likelihood with unknown covariance, especially on noisy data. That is, in this case, galaxy sample with a low number density. Ongoing work is devoted to an extended noise model, and investigating physically-motivated parametrizations and priors for σ ϵ .  \nPost-unmasking studies -Our baseline analyses opted for efficiency in model evaluation and sampling, specifically the ability to scale up our EFT FBI algorithm and individually analyze all the ten realizations of the real-space Λ CDM data suite during our stage-2 analysis. Therefore, we adopted the second-order bias expansion and assumed galaxy stochasticity to be white noise with Poisson amplitude (Eq. (21)). Our primary motivation for a post-unmasking study is to facilitate a comparison between our σ 8 constraint and that reported by the EFT P+B team. To this end, we have reanalyzed box 1 while extending the Galaxy model in the Model and the Covariance in the Inference as follows:  \n- 1. Galaxy model : going from second-order to third-order bias;\n- 2. Covariance : allowing for the galaxy stochasticity to be scale-dependent with unknown amplitudes, i.e. σ ϵ ( k ) = σ ϵ, 0 [ 1 + σ ϵ,k 2 k 2 ] , where σ ϵ, 0 , σ ϵ,k 2 are jointly inferred from the data.  \nWith these new choices, our Model and Inference correspond to the same adopted in Nguyen et al. (2024). For the additional stochastic parameters σ ϵ, 0 and σ ϵ,k 2 , we assume the following flat priors  \nP ( σ ϵ, 0 ) = U (0 . 9 , 1 . 1) σ ϵ Poisson , P ( σ ϵ,k 2 ) = U ( -10 . 0 , 100 . 0) . (24)  \nThe result of this post-unmasking analysis is shown in Fig. 22 under the label 'FBI extended.' The systematic bias in σ 8 is now completely relieved, while the increased parameter vector leads to a larger error bar compared to the preunmasking analysis. Our constraints are further align with the expectations from the EFT FBI constraints presented in Nguyen et al. (2024). Note that here we implement the analysis cutoff k max with a cubic sharpk filter. Therefore our analyses with k max = 0 . 1 h Mpc -1 here involves roughly the same number of modes as (hence should be compared to) the analyses with a spherical k max = 0 . 12 h Mpc -1 in Nguyen et al. (2024). For more discussion and post-unmasking comparisons, we refer to Sec. 6.", '5.4. SBI : P & B method 17': 'In this section, we outline the analysis of the redshift-space mocks using simulation-based inference (SBI).  \nSBI uses computationally simulated forward models of the observed data that are evaluated using parameter values spanning the prior (Cranmer et al. 2020). It leverages the forward models to learn the likelihood distribution of any measurable data-statistics and can, in principle, be applied to analyze many of the summary statistics presented in this paper. Since the likelihood is learned from the forward models, it relaxes the assumptions on the likelihood (e.g., the emulator approach typically assumes a Gaussian likelihood). Furthermore, the forward models can include observational effects to more robustly treat systematics.SBI has been applied to large-scale structure analysis in Alsing et al. (2019); Jeffrey et al. (2021); Hahn et al. (2023, 2022); Tucci & Schmidt (2023) and our analysis here will most closely follow Hahn et al. (2022).', '5.4.1. Data & Estimators': 'We focus on analyzing the power spectrum multipoles P ℓ ( k ) for ( ℓ = 0 , 2 , 4) and bispectrum monopole B 0 ( k 1 , k 2 , k 3 ) using the mean measurement from the 10 Λ CDM redshift-space boxes. Power spectrum multipoles are measured with Fast Fourier Transforms using Nbodykit (Hand et al. 2018) on a 512 3 mesh, while bispectrum is measured using pySpectrum code (Hahn 2020) on a 760 3 mesh. In addition, we supplement our data vector with the galaxy number density, ¯ n . We use both the statistics between k min = 0 . 009 h/ Mpc and k max = 0 . 5 h/ Mpc. This results in the power-spectrum data vector of size 79 × 3 (for three multipoles) and bispectrum data-vector of size 1898.', '5.4.2. Model': "For SBI, we need to accurately model galaxy clustering statistics over a range of cosmology and galaxy formation parameters. To this end, we use Quijote N -Body simulations suite (Villaescusa-Navarro et al. 2020), specifically the high resolution Λ -CDM Latin-hypercube (LH) which consists of 2000 simulations varying over 5 cosmological parametersΩ m , Ω b , h, n s and σ 8 . Each simulation evolves 1024 3 dark matter particles with TreePM Gadget-III code in a volume of 1 Gpc /h . For each simulation, there are two sets of dark matter halos identified with Rockstar (Behroozi et al. 2013) and Friends-of-friends (FoF) halo finder, respectively. Recent work has shown that SBI models trained on one set can lead to biased results on another (Modi et al. 2023), hence we will use both halo catalogs to evaluate robustness of our analysis. However, the final results of this challenge will be derived using halos from Rockstar , which more accurately identifies the position and velocity of halos.  \nNext, we populate these dark matter halos with galaxies. We use an extended 9-parameter HOD model that includes assembly, concentration and velocity biases, summarized in Table 3. Five of these parameters are the same as the standard Zheng07 HOD (Zheng et al. 2007) modellog M min , σ log M , log M 0 , log M 1 , α . These are supplemented with 'assembly' bias parameters ( A c , A s ) which modify the number of centrals and satellites based on the halo concentration 18 , 'concentration' bias parameter ( η conc ) which allows the concentration of satellite galaxies to deviate from NFW profile, and 'velocity' dispersion parameters ( η c , η s ) which rescale the central and satellites velocities over the halo velocity. For further details on implementation of this decorated HOD, we refer the reader to (Hahn et al. 2023; Hearin et al. 2016).", '5.4.3. Inference': "Methodology -SBI requires a training dataset of ( θ, x ) pairs where θ are the model parameters of interest (here, cosmology and HOD parameters) and x are the corresponding summary statistic (here, power spectrum multipoles, bispectrum and number density).  \nx is generated using the simulations and estimators described above. The dataset of { ( θ, x ) } , therefore, corresponds to samples drawn from the joint distribution p ( θ, x ) . More importantly, { ( θ, x ) } can be used to infer the posterior p ( θ | x ) .  \nWhile different methods can be used to infer p ( θ | x ) using { ( θ, x ) } , we use neural density estimators, which have been used to accurately estimate even high dimensional distributions with limited training data. More specifically, we train a conditional neural density estimator, q ϕ , with parameters ϕ to approximate p ( θ | x ) ≈ q ϕ ( θ | x ) . We do this by minimizing the Kullback-Leibler divergence between p ( θ | x ) and q ϕ ( θ | x ) , estimated by the log-probability of q ϕ evaluated over the training dataset.  \nOnce q ϕ ( θ | x ) is successfully trained, we can infer the posterior for observations, x ' , by simply querying the trained q ϕ ∗ to generate samples from the posterior: θ ' ∼ q ϕ ∗ ( θ | x ' )  \nTraining data and priors -To generate the training dataset for SBI, we use the Quijote LH described in the previous section. For each simulation in LH, we sample 10 different HOD parameter values over a prior range and generate 10 galaxy catalogs. Thus in total, we have 20,000 galaxy catalogs. We use catalogs for 1500 of these simulations (i.e. 15,000 galaxy catalogs) for training, 200 for validation, and held out 300 for testing. The prior range of cosmological parameters is set by the bounds of Quijote LH.  \nΩ m ∼ U [0 . 1 , 0 . 5] , σ 8 ∼ U [0 . 6 , 1 . 0] , Ω b ∼ U [0 . 03 , 0 . 07] , n s ∼ U [0 . 8 , 1 . 2] , h ∼ U [0 . 5 , 0 . 9] (25)  \nFor HOD parameters, we shift the central values of prior on log M min , log M 0 and log M 1 for different cosmologies so as to match the given number density on average. This allows us to more efficiently sample the parameter space by focusing only on the regions that will generate simulations loosely consistent with the data. These central values are estimated as follows: we assume a fiducial value of α = 0 . 7 and satellite fraction of 0.2, then for every cosmology we set log M θ min = log M θ 0 = M h , the halo-mass above which the number of halos is the same as the number of centrals, and then estimate log M θ 1 to match the satellites. The resulting priors are as follows:  \nlog M min ∼ U [log M θ min ± 0 . 15] , log M 0 ∼ U [log M θ 0 ± 0 . 2] , log M 1 ∼ U [log M θ 1 ± 0 . 3]  \nFor remaining HOD parameters, we use the following priors for all cosmologies:  \nα ∼ U [0 . 4 , 1 . 0] , σ log M ∼ U [0 . 3 , 0 . 5] , A c , A s ∼ N (0 , 0 . 2) over [ -1 , 1] , η conc ∼ U [0 . 2 , 2 . 0] , η c ∼ U [0 ., 0 . 7] , η s ∼ U [0 . 2 , 2 . 0]  \nWe estimate the power spectrum multipoles and bispectrum for each of the 20,000 simulated galaxy catalogs.  \nPosterior inference -Weuse sbi 19 package to train masked auto-regressive flows as conditional neural density estimators and learn the posterior, q ϕ ( θ | x ) ∼ p ( θ | x ) . To minimize stochasticity in training, we use Weights-and-Biases 20 package and train 200 networks for each data-statistic by varying hyperparameters such as the width and the number of layers, learning rate and batch size. After training, we collect 10 neural density estimators with best validation loss and use them as an ensemble i.e. we construct a mixture distribution with uniform weighting to approximate the posterior.", '5.4.4. Analysis Choices': 'Scale cuts -The Quijote simulations are only (1Gpc /h ) 3 in volume, while the data-challenge simulations are (2Gpc /h ) 3 in volume. Thus in addition to the typical small-scale cut, we also impose a large scale cut on the powerspectrum and bispectrum measurements. We use only the scales between k min = 0 . 009 h/ Mpc and k max = 0 . 5 h/ Mpc. We point out that the smaller volume of our training simulations than the mock data results in our learnt-likelihood, and hence the posterior, being over-dispersed due to larger cosmic variance. Thus we expect our results to be more conservative than if the simulations had the same volume. We also note that in the future, such limitations on the volume of the training simulations can be overcome using hybrid simulation-based inference techniques to combine EFT analysis on the large scale analysis with SBI analysis only on the small scales (Modi & Philcox 2023).  \nValidation -We use our trained ensemble to predict the cosmological parameters over the held-out test dataset to do coverage tests as described in Talts et al. (2018); Hahn et al. (2022). We verify that all the rank histograms are uniformly distributed within the rank scatter, which suggests (but does not guarantee) that the posterior marginal are accurate. To test for model-misspecification, we repeat the exercise with Quijote -FoF latin hypercube and find consistent results.  \nCaveats -The primary source of error with SBI is model mis-specification (Modi et al. 2023; Cannon et al. 2022). Previous analysis (Hahn et al. 2022; Modi et al. 2023) suggest that inference with power-spectrum up to k max = 0 . 5 h/ Mpc should be robust to these modeling differences. However, more powerful statistics such as bispectrum and combined power-spectrum and bispectrum analysis can exacerbate model mis-specification, both due to increased sensitivity and smaller error bars. This parameter-masked challenge provides a unique opportunity to investigate this aspect.', '5.5. kNNs with AbacusSummit emulator 21': 'In this section, we summarize the methodology of analyzing the 2D k -th nearest neighbor statistics ( k NN) with an extended HOD emulator based on the AbacusSummit simulation suite. The detailed methodology is described in a dedicated supporting paper (Yuan et al. 2024).  \nFigure 11. Visualizations of the RDk NN statistics (left) and DDk NN statistics (right), calculated on the redshift-space mock, averaged over 10 realizations. The dotted lines show the contours of a un-clustered Poisson random sample. We only showcase the first 4 k -orders for brevity. Note that we have defined DDk NN = 1 , thus k = 2 is the first meaningful order.  \n<!-- image -->', '5.5.1. Data & Estimators': 'For this analysis, we use two different 2D k NN statistics: (1) the random-data k NN statistics (RDk NN); (2) the data-data k NN statistics (DDk NN). For a detailed description of these statistics, we refer the readers to Yuan et al. (2023a) and Yuan et al. (2024). Conceptually, the k NN statistics captures the probability distribution of the separation between random query points (or in the case of DDk NN, the query points are galaxies) and their k -th nearest galaxies. Thus, in principle, k NNs are parametrized as a function of k and the separation, which we further decompose in 2D to a transverse component r p and a line-of-sight (LOS) component r π . Figure 11 visualizes the RD and DD k NN statistics up to the first four orders, where each panel shows how the cumulative distribution k NN-CDF depends on the separation ( r p , r π ). While the heatmap and the solid contours showcases the k NN measurement on the provided redshift-space mocks, the dashed contours visualize the measurement on an unclustered Poisson random sample of the same size for comparison. The difference between the solid and dashed contours represent the informative features in the statistics.  \nThe exact compression we use to construct our data vectors for this analysis are as follows. We use k = 1 , 2 , 3 , ..., 9 , and for each k we use 8 logarithmic bins along the r p direction between 0 . 63 h -1 Mpc and 63 h -1 Mpc, and 5 logarithmic bins along the r π direction between 0 . 5 h -1 Mpc and 32 h -1 Mpc. We also remove bins where the CDF is less than 0.05 or greater than 0.95 to increase the overall signal-to-noise of our statistics. As a result, the final statistics is summarized with 114 bins in the RD case and 144 bins in the DD case. Our r p and r π binning choices are designed to expose our analysis to the nonlinear and quasi-nonlinear scales. We are not sensitive to large-scale features such as the BAO by design. Our k choices are fairly arbitrary and are mostly chosen to confine us to nonlinear scales and computational efficiency. We reserve the discussion of optimal k choices for a future paper.', '5.5.2. Model': "To forward model the galaxy k NNs, we employ a neural network based emulator that learns the summary statistics as a function of cosmology and extended HOD parameters. In this sub-section, we introduce the relevant model details and emulation techniques.  \nParameterization -For the galaxy-halo connection, we employ the AbacusHOD package for efficiency and beyond vanilla HOD extensions (Yuan et al. 2021). The model notably assigns satellites to halo particles and includes velocity bias, two types of galaxy assembly bias, and baryonic modulations. The model parameters were summarized as model HOD2 in Table 3. Here we expand on the relevant model extensions.  \nVelocity bias is parametrized with two parameters α vel , c and α vel , s . The central velocity bias parameter α vel , c modulates the peculiar velocity of the central galaxy relative to the halo center along the LoS. α vel , c = 0 corresponds to no central velocity bias, i.e. centrals perfectly track the velocities of halo centers. We also define α vel , c as nonnegative.  \nThe satellite velocity bias parameter α vel , s modulates how the satellite galaxy peculiar velocity deviates from that of the local dark matter particle. α vel , s = 1 indicates no satellite velocity bias, i.e. satellites perfectly track the velocity of their underlying particles. Detailed descriptions of our implementation can be found in Yuan et al. (2021).  \nThe assembly bias parameters can be defined against the halo concentration or the local over-densities over a r env = 5 h -1 Mpc tophat filter. B cent = 0 , B sat = 0 indicate no galaxy assembly bias. In this analysis, we consider assembly bias against both halo concentration and local environment when we compare different models.  \nTo summarize, the extended HOD model employed in the k NN analysis consists of 10 parameters: (1) 5 vanilla HOD parameters M cut , M 1 , σ , α , κ ; (2) an incompleteness parameter f ic ; (3) velocity bias parameters α vel , c and α vel , s ; (4) galaxy assembly bias parameters B cent or B sat ; (5) baryonic modulation parameter s .  \nSimulation details -This analysis employs the AbacusSummit simulation suite (Maksimova et al. 2021), a set of large, high-accuracy cosmological N-body simulations using the Abacus N-body code (Garrison et al. 2019), designed to meet and exceed the Cosmological Simulation Requirements of the Dark Energy Spectroscopic Instrument (DESI) survey (Levi et al. 2013). AbacusSummit consists of over 150 simulations, containing approximately 60 trillion particles at 97 different cosmologies. The AbacusSummit suite also uses a new specialized spherical-overdensity based halo finder known as CompaSO (Hadzhiyska et al. 2021). For this analysis, we use the 'base' configuration boxes for forward modeling, each of which contains 6912 3 particles within a (2 h -1 Gpc ) 3 volume, corresponding to a particle mass of 2 . 1 × 10 9 h -1 M ⊙ . 22  \nTo model cosmology dependencies, we use the 85 emulator boxes. The cosmology parameter basis used for our emulator includes 8 parameters spanning the w CDM+ N eff +running parameter space: the baryon density ω b = Ω b h 2 , the cold dark matter density ω cdm = Ω cdm h 2 , the amplitude of structure σ 8 , the spectral tilt n s , running of the spectral tilt α s , the density of massless relics N eff , and dark energy equation of state parameters w 0 and w a ( w ( a ) = w 0 +(1 -a ) w a ). The different cosmologies are indexed by cXXX , where XXX ranges from 000 to 181. The details of each cosmology are described on the AbacusSummit website 23 .  \nEmulator details -We build a neural-network-based emulator to interpolate between the finite set of cosmologies. Specifically, we follow the approach of Yuan et al. (2022), where we take advantage of the high efficiency of the AbacusHOD code and run MCMC chains in the HOD parameter space against the target data vector at each cosmology. We stop the MCMC chains after 20,000 evaluations in each box and select samples whose likelihood is greater than log L > -9000 . This way, we limit the subsequent emulator training to a compact region in the cosmology+HOD parameter space, improving the emulator precision. For the surrogate model, we adopt a fully connected neural network of 5 layers and 500 nodes per layer with Randomised Leaky Rectified Linear Units (RReLU) activation. We train the network following a mini-batch routine with the Adam optimiser and a mean squared loss function, where we use the diagonal terms of the covariance matrix as bin weights. The performance of the emulator is characterized on 9 hold-out cosmologies ( c001-004 and c171-175 ). We refer the readers to Yuan et al. (2024) for details.", '5.5.3. Inference': "Likelihood Function -For this analysis, we assume a Gaussian likelihood function, see Yuan et al. (2024) for validation of the Gaussian likelihood assumptions. To sample the likelihood function, we employ the nested sampling package dynesty (Speagle & Barbary 2018; Speagle 2020a).  \nCovariance -We compute jackknife covariances for 2Dk NNs using the 10 realizations provided. We divide each realization into 125 cubic chunks, each of size ( 400 h -1 Mpc) 3 . We visualize the resulting correlation matrix in Figure A1 and A2 of Yuan et al. (2024). Both covariances are invertible. The RDk NN covariance matrix has a condition number of 1 × 10 7 and the DDk NN covariance matrix has a condition number of 4 × 10 7 . To account for emulation errors, we compute a separate emulator covariance matrix from the hold-out cosmologies and add that onto the data jackknife covariance matrix.  \nPriors -This analysis employs cosmology priors based on the Λ CDM parameter range covered by the AbacusSummit simulations,  \nω cdm ∈ [0 . 099 , 0 . 140] , σ 8 ∈ [0 . 68 , 0 . 94] , n s ∈ [0 . 90 , 1 . 02] , ω b ∈ [0 . 021 , 0 . 024] , (26)  \nwith the prior distribution specified by the AbacusSummit parameter envelope, N ur = 2 . 0328 and α s = 0 . For the HOD parameters, we adopt the following flat priors:  \nlog 10 M cut ∼ U [12 . 0 , 14 . 5] , log 10 M 1 ∼ U [13 . 0 , 15 . 0] , α ∼ U [0 . 5 , 1 . 5] , α vel , c ∼ U [0 . 0 , 1 . 0] , α vel , s ∼ U [0 . 2 , 1 . 8] , log 10 σ ∼ U [ -3 . 5 , 1 . 5] , κ ∼ U [0 . 0 , 2 . 0] , B cent ∼ U [ -1 . 0 , 1 . 0] , B sat ∼ U [ -1 . 0 , 1 . 0] . (27)  \n5.5.4. Analysis choices  \nScale cut validation -For both k NN configurations, the first cut we apply is removing bins where the CDF value measured on the target mock is either less than 0.05 or greater than 0.95. This removes the noisiest bins and Gaussianizes the likelihood function.  \nWe apply an additional scale cut to the DDk NN statistics in r p > 5 h -1 Mpc. We refer the readers to Yuan et al. (2024) for justification of this cut. Essentially, the DDk NN statistic is highly sensitive to small-scale modeling, and we find our models to not be flexible enough to accurately predict the DDk NN statistic at r p < 5 h -1 Mpc given the precision of a ( 2 h -1 Gpc) 3 volume. We do not apply any scale cut to the RDk NN statistic as it is less sensitive to small-scale modeling. Unfortunately, our scale cut at r p < 5 h -1 Mpc comes at significant cost to our constraining power on cosmology. This highlights an important avenue of future work. We discuss this point further in Sect. 6.  \nConsistency Checks -Our supporting paper (Yuan et al. 2024) details how we test and compare different HOD models invoking different extensions. We leverage both goodness-of-fit metrics and Bayesian evidence. We also conduct cross-validation tests where we fit one subset of the data vectors and predict a different set of data vectors.  \nSpecifically, we set RDk NN and DDk NN as separate data vectors as they summarize different information of the density field. RDk NN captures the density information and DDk NN captures the clustering information. We also include the standard 2PCF as a third statistics that share some information with DDk NN but also captures clustering out to larger scales. We start with three different HOD models that include some combinations of baryonic effect treatment and two different types of galaxy assembly bias in addition to the vanilla HOD. We fit all three different HOD models to one of the three statistics and predict one of the other two statistics.  \nOur tests clearly favor one HOD model, which results in good fits and no tension between different subsets of the data vectors. We show the constraints of the three models on σ 8 in Figure 12, where we see that the favored model B indeed achieves the least biased constraints post-unmasking. We do not reveal the exact details of the favored HOD model to keep the challenge for future participants masked. Note that we submitted only the cosmology constraints for the favored model B for the official unmasking.  \nUnmasking criteria -Our unmasked criteria consist of the series of tests that we have described so far. To summarize, we first test the Gaussianity of our summary statistics. As a result, we removed k NN bins where the CDF value is less than 0.05 or greater than 0.95. Then we construct our covariance matrix and make sure the matrices are invertible and have reasonable condition numbers.  \nThe second set of tests were associated to the emulators. Specifically, we ensure that mean emulator error computed on the hold-out tests is sub-dominant compared to the statistical error computed from covariances. In principle, because we are adding the emulator error to the final covariance matrix, the emulator error should not bias our results. However, because the emulator error is intrinsically a function of model parameter values and our treatment of emulator errors (only looking at mean error) is approximate, we would like to make sure potential issues with the emulator error do not dominate our final error budget and lead to significant biases.  \nFinally, we conduct a suite of consistency tests outlined in the 'consistency checks' section as the final part of our unmasking criteria. We explored different scale cuts and different HOD models until we found a configuration (model B and r p > 5 h -1 Mpc) that yields acceptable goodness-of-fits and self-consistent results in our cross-validation tests. These validation tests are critical in increasing our confidence in our results. Again, all these steps were described in detail in Yuan et al. (2024).  \n0.001  \n0.002  \nFigure 12. (Plot adapted from Figure 9 of Yuan et al. 2024) The constraints of k NN statistics on σ 8 for different HOD models and different scale cuts. Model A, B, C refer to the three HOD candidate HOD models. The blue error bars show the 1 σ marginalized constraints. The orange triangles show the maximum likelihood points. The green line represents the truth. The left, middle, right blocks represent the RDk NN constraints, DDk NN constraints with scale cut r p > 5 h -1 Mpc, and DDk NN constraints without minimum scale cuts ( r min = 0 . 67 h -1 Mpc), respectively. Our pre-unmasking consistency tests favored model B, and we see post-unmasking that model B indeed results in the least biased constraints. Comparing the middle block to the right block shows that the inclusion of small scales drastically tightens the constraints on σ 8 . Unfortunately, our models were not flexible enough to self-consistently model r p < 5 h -1 Mpc. This comparison motivates for more robust modeling of small scales.  \n<!-- image -->  \n0.010  \nmodel Amodel Bmodel Cmodel Amodel Bmodel C  \n0.075  \n0.050  \n0.025  \n0.000  \n0.025  \n0.050  \n0.075  \nmodel Amodel Bmodel Cmodel Amodel Bmodel C  \nThe 68% marginalized cosmology constraints with our final results using our favored model are shown in Figure 1 and Figure 2. The RDk NN fit achieves a χ 2 /d.o.f. = 0.44 whereas the DDk NN fit achieves a χ 2 /d.o.f. = 1.23.  \nCaveats -The main limitation is that we had to employ aggressive scale cuts because our extended HOD models were not flexible enough to reliably model scales smaller than 5 h -1 Mpc. In Yuan et al. (2024), we argue that this is specifically due to inability to model the small-scale velocities correctly. Thus, limitations with our galaxy-halo connection modeling is our main area of concern and improvement.  \nAnother issue that we were aware of before unmasking is the potential mis-characterization of emulator errors since we only accounted for the mean errors calculated on a limited set of hold-out tests. The issue with emulator errors could be mitigated with either running coverage tests or by following a likelihood-free inference approach. Both require a large amount of training data, which requires additional simulations.", '5.6. Density-split clustering 24': 'In this section, we describe the methodology for analyzing the density-split clustering statistics (DSC), based on an emulator trained on the AbacusSummit simulations, which can learn about the dependence of these statistics on cosmology and HOD parameters. More details about the implementation of this emulator can be found in CuestaLazaro et al. (2023)', '5.6.1. Data & Estimators': 'Density-split clustering (Paillas et al. 2021, 2023) is a method that characterizes the environmental dependence of galaxy clustering, exploiting the sensitivity of each environment to cosmology. We focus our analysis on measuring the DSC statistics on the 10 Λ CDM mocks in redshift space. Here we briefly summarize the algorithm to measure the DSC data vector:  \nWe begin by painting the galaxy overdensity field to a rectangular grid of 5 h -1 Mpc resolution. We smooth the field with a Gaussian kernel of width R s = 10 h -1 Mpc , and we then sample it at N query random query points, where N query is equal to five times the number of galaxies. These query points are then split into five quintiles, according to the value of the overdensity at their locations. Figure 13 shows an example of the measured overdensity distribution and its splitting into quintiles, as measured from one of the Λ CDM simulations.  \nWe proceed to measure the cross-correlation function between each quintile and the redshift-space galaxy field, as well as the autocorrelation function of the quintiles themselves. Both of these are measured in bins of s and µ ,  \nRD-  \nk  \nNN  \nDD-  \nk  \nNN  \nrp  \n>5  \ns  \nn  \n<!-- image -->  \n∆ ( R s = 10 h ) Figure 13. Probability distribution function of the galaxy overdensity measured at random query points from one of the Λ CDM simulations, where the galaxy density field has been smoothed with a Gaussian filter of radius R s = 10 h -1 Mpc . The points are classified into five density quintiles, shown by the different colours.  \nusing 241 µ bins from -1 to 1, and radial bins of scale-dependent widths: 1 h -1 Mpc width for s ∈ [0 , 4] h -1 Mpc , 3 h -1 Mpc width for s ∈ [4 , 30] h -1 Mpc , and 5 h -1 Mpc width for s ∈ [30 , 150] h -1 Mpc . The correlation functions are then decomposed into multipoles. Here, we restrict the analysis to the monopole and quadrupole moments. Figure 14 shows the measured multipoles, averaged over the 10 simulations.  \nFinally, the data vector is constructed as a concatenation of the monopole and quadrupole moments of the quintilegalaxy cross-correlation functions, ξ QG 0+2 ( s ) , and the quintile autocorrelation functions , ξ QQ 0+2 ( s ) ,  \n̂ d DSC = ( ξ Q=(0 , 1 , 3 , 4)G 0 ( s ) , ξ Q=(0 , 1 , 3 , 4)G 2 ( s ) , ξ QQ=(00 , 11 , 33 , 44) 0 , ξ QQ=(00 , 11 , 33 , 44) 2 ( s ) ) , (28)  \nusing the same radial binning as describe above. When doing so, we discard the middle quintile, Q 2 , before concatenating the data vector. The reason for this is that the information from all five quintiles is redundant, since the sum of all cross-correlation functions averages to zero by construction. Therefore, all the information available from this method is already contained in the remaining four quintiles.', '5.6.2. Model': 'We model the DSC multipoles using a neural-network based emulator that learns about the dependence of these multipoles on cosmology and HOD parameters. The emulator is trained on the AbacusSummit suite of simulations, using an extended-HOD parametrization, as described in detail in Sect. 5.5.  \nParametrization and simulations -Westart from the dark matter halo catalogues from the base AbacusSummit simulation boxes at 85 different cosmologies, which span an 8-dimensional w 0 w a CDM parameter space, including the physical cold dark matter density ω cdm , the physical baryon density ω b , the RMS of density fluctuations σ 8 , the spectral index n s , the running of the spectral index α s , the effective number of massless relics N eff , and the dark energy equation of state parameters, w 0 and w a . However, when we perform the cosmological inference, we restrict ourselves to a baseΛ CDM parameter space, using the priors on cosmological parameters defined in Eq. (31).  \nWe generate a Latin hypercube with 100,000 samples of an extended-HOD parameter space, including velocity bias and environment-based secondary bias, using the prior range defined in Eq. (32). Each cosmology is assigned 1,000 HOD parameters from the Latin hypercube, and dark matter haloes are populated with HOD galaxies from those  \nFigure 14. A visualization of the density-split clustering data vectors, along with the best-fit model from our emulator. Markers and solid lines show the data vectors and the emulator predictions, respectively. Left: multipoles of the quintile-galaxy crosscorrelation functions. Right: multipoles of the quintile autocorrelation functions. Upper and lower panels show the monopole and quadrupole moments, respectively. We also display the difference between the model and the data, in units of the data error. Each color correspond to a different density quintile, as exemplified in Fig. 13.  \n<!-- image -->  \nparameters, using the AbacusHOD code (Yuan et al. 2021). When the resulting number density of an HOD catalogue is larger than the average number density of the Λ CDM simulations ( n gal ≈ 4 . 5 × 10 -4 ( h/ Mpc) -3 ), we invoke an incompleteness parameter f ic and subsample the catalogue to match the target number density. If the number density is lower than the target, we simply keep it as is. Galaxy positions are converted to redshift space by perturbing their real-space positions with their peculiar velocities along the line of sight (taken to be one of the axes of the simulation box).  \nNeural-network emulator -We split the collection of HOD catalogues into training, validation and test sets. We construct two neural-network emulators, one for each correlation type (quintile-galaxy cross-correlation and quintile autocorrelation). The inputs to the neural network are the cosmological and HOD parameters, whereas the outputs are the concatenated monopole and quadrupole for the four quintiles.  \n<!-- image -->  \n<!-- image -->  \nFigure 15. Error sources in the density-split clustering analysis. Left: correlation matrix associated to the sample variance of the data vector. Right: correlation matrix associated to the imperfect emulation of the multipoles. Middle: relative change between the emulator and data covariance. The error associated to the finite size of the training sample, C train , is the same as the data error in this case, so we do not display it.  \n<!-- image -->  \nThe networks are fully connected, using SiLU activation functions, a mean absolute error loss function, and an AdamW optimizer. To avoisd overfitting due to the limited size of our training set, we also introduce dropout. To improve the performance of the model and reduce training time, we decrease the learning rate by a factor of 10 every 5 epochs over which the validation loss does not improve, until reaching a minimum learning rate of 10 -6 . We use the optuna framework 25 to optimize the hyperparameters of the neural network, aiming at minimizing the validation loss. More details related to the architecture and its optimisation can be found in our public repository 26 .', '5.6.3. Inference': 'Likelihood -We fit our emulator to the mean measurement of the 10 realizations of the Λ CDMredshift-space mocks. We assume a Gaussian data likelihood with covariance C tot , which incorporates the total error budget of our analysis. This covariance includes contributions from sample variance in the data vector ( C data ), sample variance in the simulations used for training ( C train ) and the intrinsic emulator error ( C emu ):  \nC tot = C data + C train + C emu . (29)  \nCovariance -To calculate C data , we measure the DSC multipoles from 1500 realizations of the small AbacusSummit simulations, which are 500 h -1 Mpc on a side and use the baseline AbacusSummit cosmology. We choose a combination of HOD parameters that produce multipoles that minimize the χ 2 with respect to the data vector measured from the Λ CDM simulations. We then divide the covariance by a factor of 64 to compensate for the volume difference between the small AbacusSummit boxes and the 2000 h -1 Mpc Λ CDM boxes. We note that here, we have assumed that the noise in the data vector is that associated to a single box of the Λ CDM simulations. To account for the finite volume of the simulations used for training, we repeat the same procedure to calculate C train , so effectively C train = C data for this analysis. Finally, to estimate C emu , which accounts for the error due to an imperfect emulation of the multipoles, we compare the emulator predictions m ( Ω test ) against the multipoles measured from our test set of simulations, ̂ d test , and compute the covariance of the absolute emulator error ∆ test = m ( Ω test ) -d test as  \nwhere the overline denotes an average across test samples.  \n̂ C emu = 1 n test -1 n test ∑ k =1 ( ∆ test k -∆ test )( ∆ test k -∆ test ) ⊤ , (30)  \nFigure 15 shows the correlation matrices associated to C data and C emu , highlighting how the different quintiles, multipoles and correlation types complement each other. We also display the ratio between the emulator and data covariance, where we can see that, while the data error is the dominant source of error at large scales, the emulator error becomes larger for small separation bins.  \nSampler and priors -We sample the posterior distribution using the Dynesty nested sampler (Speagle 2020b). We assume the following flat priors on cosmological parameters:  \nω cdm ∼ U [0 . 1032 , 0 . 140] , σ 8 ∼ U [0 . 678 , 0 . 938] , n s ∼ U [0 . 9012 , 1 . 025] , ω b ∼ U [0 . 0207 , 0 . 0243] (31)  \nWe fix those parameters beyond the baseΛ CDMmodel to the following default values: w 0 = -1 , w a = 0 , N ur = 2 . 0328 , and α s = 0 . For the HOD parameters, we adopt the following flat priors:  \nlog 10 M cut ∼ U [12 . 0 , 13 . 5] , log 10 M 1 ∼ U [12 . 5 , 15 . 0] , α ∼ U [0 . 3 , 1 . 5] , α v, cen ∼ U [0 . 0 , 1 . 0] , α v, sat ∼ U [0 . 0 , 2 . 0] , log 10 σ ∼ U [ -7 . 0 , 0 . 0] , κ ∼ U [0 . 0 , 1 . 5] , B cen ∼ U [ -0 . 5 , 0 . 5] , B sat ∼ U [ -1 . 0 , 1 . 0] (32)', '5.6.4. Analysis choices': 'To validate the baseline setup of our analysis, we run our inference pipeline adopting different choices of scales, summary statistics, HOD models and treatment of errors. The results are summarized in Fig. 16. Our baseline configuration, shown at the top of each panel, uses scales between 1 . 0 h -1 Mpc < s < 151 h -1 Mpc, including the monopole and quadrupole of the quintile-galaxy CCF and the quintile ACF, using quintiles Q 0 , Q 1 , Q 3 , and Q 0 . The extended HOD model includes both velocity and assembly bias, and the emulator systematic error, as well as the variance associated to our finite training sample, are included in the likelihood calculation. We see that adjusting these choices one by one produces changes in the recovered parameter constraints, but overall, all such variations are consistent within 1 σ of the values inferred from our fiducial configuration, highlighting the robustness of the recovered constraints against the tuning of these settings. Below, we comment on each of these aspects, based on the results from Fig. 16.  \nScale cut validation -We observe that the recovered constraints do not change significantly when excluding the very small scale information at s < 30 h -1 Mpc. This is a direct consequence of the inclusion of the systematic emulator error in our analysis, which starts to dominate the error budget on those scales. We also see that the inclusion of large-scale information adds a significant amount of constraining power. Moreover, we have run independent recovery tests on our test set from the AbacusSummit simulations, which further confirm that our fiducial range 1 h -1 Mpc < s < 151 h -1 Mpc seems robust at recovering unbiased cosmological constraints for the volume of the Λ CDM simulations. Based on this, we choose to adopt this scale range for the baseline analysis.  \nChoice of data vector -Concerning the DSC statistics, the quintile-galaxy CCF amounts for most of the information, with the quintile ACF only adding a small contribution to the budget. In terms of the multipoles, the constraints are mostly driven by the monopole. We also see that the quintiles at the extrema, Q 0 and Q 4 , carry most of the information. This is in line with previous findings based on Fisher forecasts for DSC statistics (Paillas et al. 2023).  \nChoice of HOD model -The choice of the HOD model has little impact on ω cdm , but becomes noticeable for σ 8 and n s . More specifically, not including central and satellite velocity bias in the model shifts the mean of the posterior for σ 8 and n s to smaller values, although well within the 1 σ region of the baseline setup. A similar trend is seen when excluding assembly bias. Given the findings of previous studies, which suggest that velocity and assembly bias are important for accurately fitting the small-scale clustering from observational data and simulations (Guo et al. 2015; Yuan et al. 2021, 2022), and to ensure our emulator is robust against the potential presence of such effects in the Λ CDM simulations fitted in this study, we choose to include them in the default analysis.  \nChoice of error budget -We see that the exclusion of the emulator systematic error has an important effect in the analysis. Removing this error has a dramatic increase on the precision at which we recover the spectral tilt n s , and a smaller effect is observed for ω cdm and σ 8 . A similar trend is spotted for the exclusion of the error associated to the training sample, although in a lesser extent. Based on separate tests with the AbacusSummit simulations, we  \ncdm  \nFigure 16. Masked constraints from density-split clustering, assuming different configurations for the analysis. Blue dots show the mean of the marginalized posterior for each parameter, and the error bars show the associated 68% confidence interval. The blue-shadowed region shows the constraints for the baseline configuration adopted in the paper.  \n<!-- image -->  \n<!-- image -->  \n8  \n<!-- image -->  \ns  \nhave concluded that this increase of constraining power comes at the expense of potential biases in the recovered cosmological parameters (Cuesta-Lazaro et al. 2023). Therefore, to ensure the robustness of our inference analysis, we include both error contributions in the likelihood calculation.  \nChoice of likelihood -The choice of a Gaussian likelihood in our analysis is motivated by previous tests we have carried out with the Quijote and AbacusSummit simulations, which show that the likelihood associated to the DSC data vector is indeed close to Gaussian for the scales and the volume of interest (see e.g., Fig. C1 from Paillas et al. (2023)).  \nUnmasking criteria -Having explored the different choices of settings in our inference pipeline, we have found that the resulting cosmological constraints are robust against such choices. Moreover, we have verified that our best-fit model provides a good χ 2 to the measured data vector. Therefore, we allow our results to be unmasked, adopting the baseline configuration for our analysis described in the previous section.  \nCaveats -Possible caveats that could lead to unexpected results in our analysis after unmasking include:  \n- - HOD parameters lying outside the range of our priors.\n- - An HOD model including effects that were not implemented in our emulator, such as satellite profile flexibilities to account for baryonic effects.\n- - An overly-conservative error budget in our likelihood, mainly related to the systematic error coming from the emulator predictions. We have chosen to add the emulator error calculated across the full prior range of our parameters, which could lead to an over-estimation of the error bars at small scales, and consequently, to an unnecessary degradation of the precision of our constraints.\n- - Previous analyses have shown that DSC is particularly sensitive to Alcock-Paczynski distortions (Paillas et al. 2021), which can increase the relative gain of information with respect to standard two-point measurements. Here, we have worked directly with the true galaxy positions in redshift space, ignoring the Alcock-Paczynski effect, which could potetially degrade the relative improvement in parameter constraints found in previous studies.  \n- - Due to computational constraints, we have not explored changing the number of quantiles that the density PDF is split into, nor the width of the kernel that is used to smooth the density field. We followed choices made in previous analyses that were calibrated for the number density and redshift of BOSS CMASS. These choices could be sub-optimal for the galaxy sample used in this study, leading to weaker constraints than what could be achieved when varying those settings.', '5.7. Cosmic Voids 27': 'Cosmic voids, the extended under-dense regions of the cosmic web, are becoming an increasingly important probe in the galaxy-clustering and cosmology communities (Pisani et al. 2019; Moresco et al. 2022). While originating in the 3D distribution of matter, voids can also be identified in tracers of the latter, such as galaxies. In this work we focus on two summary statistics: the void size function (VSF) and the void-galaxy cross-correlation function (VGCF). The VSF specifies the number of voids as a function of their size and is sometimes referred to as void abundance. The VGCF describes the spatial clustering between the centers of voids and their surrounding galaxies in the form of a cross-correlation function. In the case where all spatial separations are expressed in units of individual void sizes, the VGCF is sometimes simply denoted as a void stack . Both techniques have already been applied to SDSS BOSS data (Dawson et al. 2013) to derive constraints on cosmological parameters (e.g., Hamaus et al. 2016, 2020; Contarini et al. 2023, 2024). In this challenge we consider both the VSF and the VGCF measured in the mock light cone.', '5.7.1. Data vector & estimators': 'We rely on the popular void finding algorithm VIDE 28 (Sutter et al. 2014b), which is based on ZOBOV (Neyrinck 2008) to estimate a tracer density field via Voronoi tessellation and makes use of the watershed technique (Platen et al. 2007) to identify local depressions therein. VIDE has been used extensively for void analyses in both simulations (e.g., Hamaus et al. 2014b; Chan et al. 2014; Pollina et al. 2017; Schuster et al. 2019; Contarini et al. 2019; Verza et al. 2019; Kreisch et al. 2022) and observational data (e.g., Sutter et al. 2012; Hamaus et al. 2016; Pollina et al. 2019; Fang et al. 2019; Contarini et al. 2023, 2024). In order to run on light-cone data, VIDE requires a mask to specify the (mock) survey footprint. For this we generate a healpix mask of resolution n side = 128 from the angular distribution of mockgalaxy coordinates provided in the catalog. This resolution is sufficient to capture the simple rectangular geometry with the given ranges in RA and Dec. The final data product of VIDE consists of various catalogs containing different void properties, such as their sky coordinates, redshifts, effective radii, shape and density information. Because VIDE identifies voids in comoving space, angles and redshifts are transformed to comoving coordinates assuming a fiducial cosmology. Therefore, void properties defined in comoving space, such as effective radii, volumes and densities, depend on this assumption. The same fiducial cosmology is used to transform the comoving coordinates of voids back to sky coordinates, which are therefore less affected by model assumptions. Here we focus our analysis on the Λ CDM light-cone mock catalog, where VIDE identifies more than 14 , 000 voids.  \nVoid size function (VSF) -The data vector for the VSF is obtained by measuring the number of voids found by VIDE in logarithmic bins of radius. In our analysis we split the void sample in three equispaced bins of redshift, i.e. 0 . 812 < z ≤ 0 . 968 , 0 . 968 < z ≤ 1 . 124 and 1 . 124 < z ≤ 1 . 280 . We present in Fig. 17 (first plot) the number of identified voids as a function of their effective radius R , for the three aforementioned redshift bins. R is defined as the radius of a sphere whose volume equals that of the void, which is non-spherical in general. To align the measured VSF with the theoretical model and remove spurious voids, subsequent cleaning is performed. The cleaning procedure is described below. The VSF data vector is given by  \nwhere ˜ R i =1 ,...,N r denotes the N r bins in cleaned void size, and z 0 , 1 , 2 the three tomographic bins described above.  \n̂ d VSF = ( ̂ N ( ˜ R i =1 ,...,N r , z 0 ) , ̂ N ( ˜ R i =1 ,...,N r , z 1 ) , ̂ N ( ˜ R i =1 ,...,N r , z 2 ) ) , (33)  \nVoid-galaxy cross-correlation function (VGCF) -The data vector of the VGCF is constructed by counting the number of galaxies around each void center as a function of their separation in directions along and perpendicular to the line of sight, s ∥ and s ⊥ , measured in units of the effective radius R . A stack over a sample of voids then yields a crosscorrelation function between void centers and galaxies (Hamaus et al. 2015). We use the Landy-Szalay estimator (Landy  \n& Szalay 1993) to obtain a measurement of the VGCF data vector via  \n̂ d VGCF ,i ≡ ̂ ξ vg ,i ( s ⊥ , s ∥ ) , (34) where we use 15 linearly spaced bins in each direction, ranging from the void center to about twice the effective void radius. The randoms are constructed by sampling from the smoothed redshift distributions of galaxies and voids identified in the mock data, and uniformly across the sky, with an oversampling factor of 10 (we have checked that a factor of 50 yields indistinguishable results). In order to impose the identical sky footprint as in the mock data, we apply the same healpix mask to the randoms.', '5.7.2. Model': "VSF Parameterization -The VSF model, developed by Sheth & van de Weygaert (2004) and further modified in Jennings et al. (2013) to account for the volume conservation in the transition to nonlinearity, relies on the excursionset theory. The so-called Vdn ( volume conserving ) model is defined as:  \nd n dln R = f ln σ R ( σ R ) V ( R ) dln σ -1 R dln R L ∣ ∣ ∣ R L = R L ( R ) , (35)  \n∣ where σ R represents the root mean square variance of linear matter perturbations on a scale R L and where, to compute the number density of voids, we rely on the multiplicity function:  \nf ln σ R ( σ R ) = 2 ∞ ∑ j =1 exp ( -( jπx ) 2 2 ) jπx 2 sin ( jπD ) , with D ≡ | δ L v | δ L c + | δ L v | and x ≡ D | δ L v | σ R , (36)  \nConsidering linear theory, the function f ln σ R ( σ R ) represents the volume fraction of the Universe occupied by cosmic voids, with radii in the range ( R,R +d R ) . The additional cosmological dependence of this function lies in the density contrasts defining the formation of dark matter halos and cosmic voids, δ L c and δ L v , respectively. To predict the number of voids identified from observations of biased tracers, a number of effects have to be considered:  \n- (i) the link of δ L v to its nonlinear counterpart δ NL c (corresponding to the same threshold but considering nonlinear theory) in the biased tracer field;\n- (ii) the dynamic distortions caused by the peculiar velocities of the tracers;\n- (iii) the impact of geometric distortions.  \nFirstly, to account for the link of δ L v to its nonlinear counterpart in the biased tracer field, the sample of voids analyzed has to be prepared to align with the definition of voids given by the VSF theory. For this purpose, we apply to the void catalog built with VIDE the cleaning algorithm developed by Ronconi & Marulli (2017), publicly available in the libraries CosmoBolognaLib 29 (Marulli et al. 2016). The output of this procedure is a catalog of non-overlapping spheres with cleaned void radius ˜ R , embedding a fixed density contrast in the tracer density field, δ NL v , tr . We chose to fix the latter to δ NL v , tr = -0 . 7 . We then use a bias relation to convert the density contrast in the tracer density field, δ NL v , tr , to the corresponding one in the matter distribution, δ NL v . Galaxy bias in voids cannot be properly described by the large-scale effective galaxy bias (Contarini et al. 2019; Verza et al. 2022). To model it, we therefore rely on a linear function of the effective bias (Contarini et al. 2019, 2021, 2022, 2023, 2024):  \nδ NL v = δ NL v , tr F ( b eff , σ 8 ) , with F ( b eff σ 8 ) = C slope b eff σ 8 + C offset . (37)  \nSecondly, tracer peculiar velocities cause an enlargement along the observer's line of sight, further increasing the mean void radius (Pisani et al. 2015b; Hamaus et al. 2020; Verza et al. 2023). The above parametrization is used to  \n<!-- image -->  \nFigure 17. The first plot shows the number of voids as a function of their effective radius R (pre-cleaned VSF), identified with VIDE in three equispaced redshift bins (see legend). The second, third and fourth plots show the number counts of voids used for the VSF analysis, after the cleaning procedure, as a function of the cleaned radius ˜ R (for details on the procedure, see section 5.7.2), measured in three equispaced redshift bins. We represent the data with round markers and Poisson error bars. We show in blue the best-fit model derived with the Bayesian analysis described in section 5.7.3. The continuous line indicates the median of the model, while the shaded area covers from the 16 -th to the 84 -th percentile.  \n<!-- image -->  \nalso correct for this enlargement (Contarini et al. 2022), allowing us to statistically align voids observed in the redshiftspace galaxy distribution to theoretical voids, identified in real space with unbiased tracers. In Eq. (37), C slope and C offset are redshift-independent coefficients of the linear function F and can be calibrated by using simulations or, as in this case, marginalized over. Adopting the conservative choice of using wide uniform priors on these parameters, we are able to constrain only the parameter Ω m . We underline, however, that exploiting the information derived from mock catalogs to limit the priors of C slope and C offset would allow us to provide cosmological constraints on σ 8 as well (Contarini et al. 2023).  \nThirdly, we account for the Alcock-Paczyński (AP) effect, i.e. geometric distortions due to mismatch between the fiducial and the true cosmology. Following the prescriptions in Sánchez et al. (2017); Hamaus et al. (2020); Correa et al. (2021), we rescale the observed void radius R ∗ to R = q 1 / 3 ∥ q 2 / 3 ⊥ R ∗ , where q ∥ and q ⊥ are defined via  \nr ∥ = H ∗ ( z ) H ( z ) r ∗ ∥ ≡ q ∥ r ∗ ∥ , r ⊥ = D A ( z ) D ∗ A ( z ) r ∗ ⊥ ≡ q ⊥ r ∗ ⊥ , (38)  \nand the asterisk symbol indicates quantities evaluated in the fiducial cosmology. Here, r ∗ ∥ and r ∗ ⊥ are the observed comoving distances between two points at redshift z , projected along the parallel and perpendicular line-of-sight directions, respectively. H ( z ) is the Hubble parameter and D A ( z ) is the comoving angular-diameter distance.  \nThe final VSF model depends on cosmology through the quantities σ R , q ∥ and q ⊥ , i.e. those parameters determining the amplitude of the matter power spectrum and cosmological distances. Additionally, the model depends on the nuisance parameters C slope and C offset , as well as on the effective bias of tracers. The latter is estimated from their two-point correlation function and marginalized over. Hence, the model space that we consider for the VSF is spanned by 7 parameters: Ω m , σ 8 , h, Ω b h 2 , n s , C slope , C offset . In addition, we assume one massive neutrino species of minimal mass. As already mentioned, in this analysis we will provide constraints on Ω m only, because no calibration of the parameters C slope and C offset is available with the cosmic tracers used in this work. Figure 17 presents measurements of the pre-cleaning VSF, and of the cleaned VSF including best-fit models in the considered redshift bins.  \nVGCF Parameterization -We rely on three fundamental assumptions to model the VGCF. Firstly, we make use of the cosmological principle, implying that stacked voids are spherically symmetric and hence statistically isotropic in real space. This means that they can be used as standard spheres to measure distance ratios via the AP effect. Secondly, we assume that the peculiar velocity field u v around the void center follows linear dynamics according to the linearized  \ncontinuity equation (Peebles 1980),  \nu v ( r ) = -f ( z ) H ( z ) 1 + z r r 3 ∫ r 0 δ ( r ' ) r ' 2 d r ' , (39)  \nwhere f ( z ) is the linear growth rate of the matter-density contrast δ . Peculiar velocities cause voids to appear anisotropic in redshift space, because the Doppler effect contributes an additional component to the cosmological redshift along the r ∥ direction. In order to successfully exploit voids as standard spheres, these redshift-space distortions have to be accurately modeled. For this purpose, the validity of Eq. (39) has been thoroughly tested in simulations and was found to be extremely accurate in void environments (Hamaus et al. 2014a; Schuster et al. 2023), opening up the opportunity to use the AP test with voids as a precision probe of cosmology. However, in this mock challenge we do not have access to the full matter distribution, so our third assumption relies on a linear bias relation between the galaxy- and matter-density contrasts, such that ξ vg ( r ) = bδ ( r ) . This relation has been verified in various simulation studies (Sutter et al. 2014a; Pollina et al. 2017, 2019; Contarini et al. 2019; Ronconi et al. 2019; Schuster et al. 2023), finding that the proportionality constant b asymptotes to the linear large-scale tracer bias b 1 for large voids, while smaller voids tend to yield higher, but scale-independent values.  \nWith these ingredients we can model the coordinate transformation between the real-space vector r and redshift-space vector s for the separation of void centers and galaxies,  \ns = r + 1 + z H ( z ) u ∥ = r -f ( z ) b ( z ) r ∥ r 3 ∫ r 0 ξ vg ( r ' ) r ' 2 . (40)  \nThis requires the VGCF in real space, ξ vg ( r ) , which is not directly observable. We can, however, observe the line-ofsight projected VGCF, ˜ ξ vg ( s ⊥ ) = ∫ ξ vg ( s ) d s ∥ , which is insensitive to redshift-space distortions, as it only depends on separations s ⊥ on the plane of the sky. From it we obtain the real-space VGCF by deprojection via the inverse Abel transform (Pisani et al. 2014; Hawken et al. 2017),  \nξ vg ( r ) = -1 π ∫ ∞ r d ˜ ξ vg ( s ⊥ ) d s ⊥ ( s 2 ⊥ -r 2 ) -1 / 2 d s ⊥ . (41)  \nEqs. (39)-(41) fully specify the redshift-space VGCF at linear order in perturbation theory and a closed-form expression for ξ vg ( s ⊥ , s ∥ ) can be derived with the Jacobian of Eq. (40) (Cai et al. 2016; Hamaus et al. 2017, 2020). In order to account for systematic effects, such as inaccuracies in the deprojection procedure and selection effects in the void sample due to sparse sampling of tracers with nonlinear redshift-space distortions (see caveats below), we augment this model with two additional nuisance parameters M (for monopole-like) and Q (for quadrupole-like), and adopt the semi-empirical expression derived in Ref. Hamaus et al. (2022),  \nξ vg ( s ⊥ , s ∥ ) = M { ξ vg ( r ) + f b ξ vg ( r ) + 2 Q f b µ 2 [ ξ vg ( r ) -ξ vg ( r ) ] } , (42)  \nr ⊥ = q ⊥ s ⊥ , r ∥ = q ∥ s ∥ [ 1 -1 3 f b M ξ vg ( r ) ] -1 , (43)  \nwhere µ = r ∥ /r and ξ vg ( r ) = 3 r -3 ∫ r 0 ξ vg ( r ' ) r ' 2 d r ' . We make use of Eq. (40) to map the coordinates from observed redshift space to real space, where the model is evaluated, and include the AP parameters to account for geometric distortions,  \nwhich can be solved via iteration upon setting an initial value of r = s (Hamaus et al. 2020). Because we measure void-centric distances in units of the effective void radius R , which scales as q 1 / 3 ∥ q 2 / 3 ⊥ with the AP parameters, only ratios of q ⊥ and q ∥ appear in Eq. (43). Hence, the final model space for the VGCF is spanned by four parameters; f/b , q ⊥ /q ∥ , M , and Q . To summarize, our modeling of the VGCF involves the following steps:  \n- (i) estimate the real-space VGCF via deprojection of its line-of-sight projected counterpart using Eq. (41);\n- (ii) account for dynamic (redshift-space) distortions along the line of sight using Eq. (40);\n- (iii) account for geometric (AP) distortions using Eq. (43);  \n0.0  \nFigure 18 depicts a representation of our data vector, the VGCF in separations along and perpendicular to the line of sight, along with the best-fit model based on Eq. (42).  \n<!-- image -->  \n⊥  \nFigure 18. Void-galaxy cross-correlation function in separations along and perpendicular to the line of sight, extracted from the entire Λ CDM light cone (color scale with black contour lines). White contours show the best-fit model with a reduced chi-square value shown on the bottom right.  \n- (iv) account for systematic effects arising via the deprojection procedure, or via selection effects due to sparse sampling and nonlinear RSD, by including nuisance parameters for the monopole and quadrupole terms in Eq. (42).", '5.7.3. Inference': 'VSF Covariance -We assume the off-diagonal terms of the VSF covariance matrix to be negligible on the basis of previous works, such as Bayer et al. (2021); Kreisch et al. (2022); Contarini et al. (2023); Pelliciari et al. (2022). For its diagonal elements, we assume Poisson statistics.  \nVSF Likelihood Function -Since the VSF is based on the number counts of cosmic voids, we assume the likelihood to follow the Poisson statistic (Sahlén et al. 2016; Thiele et al. 2023):  \nL [ ̂ d VSF | Ω ] = ∏ i m VSF ,i ( Ω ) ̂ d VSF ,i exp[ -m VSF ,i ( Ω )] d VSF ,i ! , (44)  \n̂ with model parameters Ω listed in Table 4. We use the functions implemented in the CosmoBolognaLib to sample the likelihood.  \nVSF Priors -We impose wide uniform priors for Ω m and σ 8 , i.e. U [0 . 05 , 0 . 9] and U [0 . 2 , 2] , respectively. We impose instead restricted priors on the remaining cosmological parameters, using the uncertainty provided by Planck2018 (Planck Collaboration et al. 2020) multiplied by a factor of three. This latter prior choice is aimed at reducing the parameterspace volume on those parameters that would be weakly constrained by the VSF. This allows us to increase the speed of our MCMC without introducing biased results. Then we marginalize over the nuisance parameters C slope and C offset , assigning uniform priors as well. For a summary of all our parameters and priors, see Table 4.  \nVGCF Covariance -We estimate the covariance matrix C of the VGCF by means of jackknife resampling the selected void sample, which is spatially non-overlapping. To this end we apply Eq. (34) after removal of one void at a time and calculate the covariance over all jackknife samples.', 'Beyond-2pt mock challenge': 'Table 4. Parameters and prior bounds for the VSF (left) and VGCF (right) analyses.  \nU  \n-  \nVGCF Likelihood Function -We adopt a Gaussian likelihood for the VGCF, which has previously been validated on the QPM (White et al. 2014) and PATCHY (Kitaura et al. 2016) mocks in Hamaus et al. (2017, 2020). We apply the Hartlap correction (Hartlap et al. 2007) to estimate the inverse covariance matrix and maximize the likelihood with respect to the parameter vector Ω = ( f/b, q ⊥ /q ∥ , M , Q ) . To sample the posterior probability distribution of these parameters, we use the affine-invariant MCMC ensemble sampler emcee (Foreman-Mackey et al. 2019). Two parameters of our model explicitly depend on cosmology. One via the growth rate, which can be approximated as f ( z ) ≃ Ω 0 . 6 m ( z ) (Lahav et al. 1991), the other one via the AP parameter q ⊥ /q ∥ ∝ D A ( z ) H ( z ) . The latter is particularly well constrained by the AP effect from voids and is therefore a sensitive probe of the expansion history (Lavaux & Wandelt 2012; Hamaus et al. 2015). In a flat Λ CDM cosmology (neglecting the presence of neutrinos and radiation), the product D A ( z ) H ( z ) is fully determined by the parameter Ω m . We obtain a posterior on the latter by sampling from a Gaussian likelihood containing measurements of q ⊥ /q ∥ and their uncertainty within each redshift bin.  \nVGCF Priors -We impose identical uniform priors for each of our model parameters, Ω ∼ U [ -10 , 10] . These prior boundaries are wide enough to be uninformative, and we have checked that our results are insensitive to this particular choice. For our cosmological parameter of interest in the Λ CDM light-cone mock, we further impose Ω m ∼ U [0 , 1] .', '5.7.4. Analysis choices': 'VSF scale cut validation -We select voids above the cleaned radius ˜ R min = 33 h -1 Mpc , for all the considered redshift bins, i.e. 0 . 812 < z ≤ 0 . 968 , 0 . 968 < z ≤ 1 . 124 , 1 . 124 < z ≤ 1 . 280 . This selection is applied to avoid those spatial scales affected by a loss of void counts and depends on the resolution of the catalog. We measure the void counts as a function of the cleaned radius fixing the value of ˜ R min , and the number of resulting non-empty radius bins are [6 , 7 , 6] , respectively. This is shown in Fig. 17, where we represent the measured void abundances and the VSF best-fit model derived by the MCMC analysis. The agreement between data and theory is such that all measurements are reproduced by the model within 68% uncertainty.  \nVSF consistency checks -We performed a number of tests to check the stability of the results with the choice of the minimum cleaned void radius, i.e. ˜ R min . As expected, when increasing the value of ˜ R min (and thus reducing the number of voids analyzed), the constraining power is reduced. This is shown in the left panel of Fig. 19, where we report the posterior distribution of Ω m derived by imposing ˜ R min = [33 , 34 , 35] h -1 Mpc . Over a certain scale, that is above the region where the void counts incompleteness is stronger, the cosmological constraints on Ω m are all consistent within 68% of uncertainty. However, we note a moderate dependency of the best-fit value of Ω m with the choice of ˜ R min , which could impact the robustness of our results (see caveats below).  \nWe also assessed the impact of assuming different fiducial cosmologies, varying in particular the value of Ω m . The results of this test are reported in the right panel of Fig. 19. We find constraints consistent within 68% of uncertainty in all the analyzed cases, i.e. in a range of ± 0 . 04 from our fiducial value of Ω m . Moreover, we tested different void radius and redshift binning choices, finding consistent results despite the number of cleaned voids being low, thus affected by statistical noise. Finally, we performed consistency checks on the impact of galaxy effective bias modeling and on the bias relation in voids from Eq. (37)), finding a negligible impact on the posterior distribution of Ω m .  \nVGCF scale cut validation -In contrast to most standard probes of large-scale structure, the VGCF analysis is not restricted to the largest scales. In fact, linear theory accurately describes the VGCF on all scales, even arbitrarily close to the void center. This is because the dynamics inside voids remain linear, in accordance with Eq. (39), which  \nFigure 19. Posterior distribution for the value of Ω m based on the modeling of the VSF in the Λ CDM light cone. The mean of the reference posterior (in blue in each panel) has been subtracted to mask the results. We show in the left panel the impact of different cuts on void radius, indicating with shaded regions the 68% confidence interval around the maximum of the posterior distribution. The right panel illustrates the effect of varying the fiducial Ω m value in our analysis pipeline.  \n<!-- image -->  \nhas been tested in simulations for a wide range of void sizes down to only a few Mpc in effective radius (Schuster et al. 2023).  \nHowever, both the void identification method and estimators for the VGCF are affected by sparse sampling of tracers (Sutter et al. 2014a; Cousinou et al. 2019; Schuster et al. 2023). A characteristic scale of a galaxy distribution is its mean galaxy separation (mgs), below which discreteness effects from sparse sampling become important. For example, in this regime random Poisson noise can create or disrupt void detections and counts-in-shell estimators for the VGCF can return biased results due to low or no particle-count statistics. We therefore restrict our VGCF analysis to voids whose effective radius is larger than a multiple of the mean galaxy separation. This can be implemented as a redshift-dependent cut with  \nwhere N mgs is a tuning parameter and ¯ n ( z ) is the mean mock-galaxy density as a function of redshift. The choice of N mgs is a trade-off between maximizing the void sample size (lower N mgs ) and minimizing the number of spurious voids (higher N mgs ). In addition, we restrict our sample to the largest 50% of all voids passing this cut. We repeated our analysis for a range of values with N mgs ∈ [1 , 5] and find consistent results with decreasing uncertainty in our posterior constraints down to a value of N mgs = 3 . Below that the posteriors start shifting and stop shrinking (see the left panel of Fig. 20), indicating the onset of stochastic bias or noise in the void sample. However, we select N mgs = 5 as our more conservative default, which results in a minimum void radius of R ≃ 50 . 6 h -1 Mpc and corresponds to the data vector shown in Fig. 18.  \nR > N mgs ( 4 π 3 ¯ n ( z ) ) -1 / 3 , (45)  \nFurthermore, we have the option to restrict or subdivide the provided redshift range. The choice of redshift binning is a compromise between detecting the redshift evolution of the VGCF on the one hand, and maintaining sufficient statistics for its estimation on the other. For example, in a w CDM cosmology, measuring the AP parameter in multiple redshift bins is necessary to break the degeneracy between w and Ω m entering D A ( z ) H ( z ) . However, in Λ CDM a single redshift bin is sufficient to determine Ω m . Hence, for the Λ CDM light-cone mock we use the entire redshift range available without binning to estimate the VGCF with reduced statistical noise.  \nVGCF consistency checks -A strategy to gain confidence in our analysis is to investigate its dependence on the assumed fiducial cosmology. We began using a flat Λ CDM cosmology with Planck2015 (Planck Collaboration et al. 2016) parameters, both in the void finder and for the construction of the VGCF data vector. The best-fit cosmology obtained from our posterior is then used to update our fiducial cosmology and this procedure is repeated until convergence. In practice, only one such iteration was necessary in updating our void catalog to obtain a stable result within our error margins. The right panel of Fig. 20 shows the impact of varying our fiducial value of Ω m when constructing the VGCF data vector and running our analysis pipeline. It is evident that the posterior distribution moves diametrically  \n<!-- image -->  \nFigure 20. Posterior distributions for the value of Ω m based on the AP test with the VGCF in the Λ CDM light cone. The mean of the first distribution (solid blue line in each panel) has been subtracted to mask the results. The left panel shows the influence of different cuts on void effective radius in units of the mean galaxy separation (mgs), while the right panel illustrates the impact of variations in the assumed fiducial values for Ω m in our analysis pipeline.  \n<!-- image -->  \nto changes in the fiducial value of Ω m . This allows us to identify a fiducial value which agrees most closely with its corresponding posterior mean. We select this posterior as our final result.  \nUnmasking criteria -We calculate the reduced χ 2 statistic for the best-fit model of the VGCF with parameters Ω ∗ as  \nχ 2 red = -2 N dof log L ( ̂ d VGCF | Ω ∗ ) , (46)  \nwith N dof = N bin -N par degrees of freedom, where N bin is the number of bins in the data and N par the number of model parameters. A value close to unity indicates a satisfactory model fit, but does not guarantee unbiased parameter constraints. However, combined with the iteration strategy over fiducial cosmologies described above, it can be used to indicate convergence to an optimal result.  \nCaveats -This section describes failure modes for our analysis that could impact the final result for one or both the considered void statistics. These can entail small shifts for the posterior, or complete failure of predicted values.  \n- · Void catalogs can be contaminated by spurious voids due to sparse sampling (e.g., Neyrinck 2008; Sutter et al. 2014a). Recent work has shown that small voids are more likely to be spurious (Pisani et al. 2015b; Cousinou et al. 2019; Schuster et al. 2023), therefore our analysis cuts based on void size should reduce the probability to include spurious voids in both the VSF and VGCF measurements. An additional purpose of this cut in the VSF analysis is to exclude spatial scales that are affected by incompleteness, i.e. loss of void counts (Pisani et al. 2015a). This occurs due to the sparsity of tracers and is currently not modeled by the VSF theory. The posteriors of both analyses mildly depend on the choice of these cuts, we therefore caution that they can have an impact on our cosmological constraints. This additional source of error is so far neither quantified, nor included as a systematic error in our analysis.\n- · Despite the geometric distortion correction applied, our final results exhibit a mild dependence on the fiducial cosmology assumed. This indicates residual geometric effects that may play a role, but currently this dependence is not quantified as a potential systematic error.\n- · Based on previous studies (Pelliciari et al. 2022; Contarini et al. 2023) the VSF analysis assumes the off-diagonal terms of the covariance matrix to be null. However, we acknowledge the possibility that this assumption could have a minor influence on our results. Also, the assumed Poisson errors associated with the VSF are generally smaller than the true uncertainty related to this measure, so the errors on cosmological parameters may be underestimated. The jackknife approach to determine the covariance of the VGCF can be limited by the finite sample size of voids, which typically leads to an overestimation of uncertainties (Hamaus et al. 2022). Moreover, a super-sample covariance is not included in our analysis, but we expect it to have a negligible impact on the final cosmological constraints (for the VSF, see Bayer et al. (2022)).  \n- · Possible systematic errors related to the void finding and cleaning algorithms can also have an impact. A void sample identified via sparse tracers in redshift space is subject to various selection effects (e.g., Correa et al. 2021, 2022). While we can limit the impact of sparse sampling via Eq. (45), residual selection biases may remain present in our sample. A full exploration of their correlation with various void properties and their mitigation is beyond the scope of this work, however.\n- · The impact of HOD models on void statistics has barely been investigated in the literature so far. Because these models are usually tuned to reproduce the two-point statistics of galaxies, their performance on void statistics is largely unknown and potentially less accurate. For example, a particular choice of HOD model may produce a strong Finger-of-God effect, which results in elongated features along the line of sight. We do observe some elongation beyond the void boundaries in the VGCF shown in Fig. 18. While this may be statistically insignificant, Fingers of God could impact the inner shape of the VGCF resulting in overly flattened contours. In turn, this would bias the AP parameter. Furthermore, the HOD scheme may affect the VSF modeling via the bias-in-void parametrization in Eq. (37).\n- · The deprojection based on Eq. (41) applied to the VGCF is exact for noiseless data, but can amplify statistical fluctuations in noisy data by large amounts (Pisani et al. 2014; Hamaus et al. 2020). The survey characteristics, like tracer density, redshift range, and sky coverage provided in the challenge mocks are more restricted than those available in BOSS, or expected by Euclid, for example. This limits the total number of bins (e.g., in redshift and effective radius) we can divide our void sample by to a single one. We therefore expect stronger systematic errors due to deprojection effects in this analysis.  \nFinally, we emphasize that the current analysis of the presented void statistics does not make use of any external simulation suites that can be used for (nuisance-) parameter calibration, template fitting, and the training of emulators. Such tools are helpful to significantly reduce parameter biases and uncertainties, but require a substantial amount of pre-processing and computing resources. We leave such more extensive approaches to future work.  \nPost-unmasking studies -After unmasking, we realized that the informative prior on h, n s , ω b adopted by the baseline VSF analysis excluded the true cosmology of the light-cone mock. In subsequent analyses, we repeated the VSF analysis with (a) the informative cosmology prior re-centered to include the true cosmology and (b) the broad cosmology prior of the BACCO analysis (Eq. (7)). The central value of Ω m is robust to both prior variations, and the wide prior broadens the marginalized 1 -σ uncertainty on Ω m from 12 . 7% to 15 . 6% .  \nLessons learned -One of the most relevant lessons learned in this challenge concerns the importance of gaining confidence for the treatment of small voids. In particular, performing a masked analysis reveals that when we are uncertain of the outcome, it is tempting to be overly conservative with scale cuts. Our first tests based on less conservative cuts resulted in best-fit cosmologies with a significantly different mean galaxy number density as compared to the one provided by the challenge. Although we were aware that this mismatch could have had multiple origins apart from the background cosmology (e.g. the HOD prescription), we still opted for very cautious scale cuts due to this observation.  \nIn hindsight, however, it is reassuring to realize that the initial results based on larger void samples that extend down to smaller effective radii actually lie closer to the truth with smaller error bars (see the left panels of Figs. 19 and 20). This suggests we could have been more ambitious with including smaller voids in our analyses, but also that our level of confidence in them was too low for doing so. For the future, it demonstrates the importance for increasing our understanding with regards to the systematics of small voids, so that their constraining power can be fully exploited.', '6. DISCUSSION': 'After unmasking, several teams further investigate the accuracy and precision of their original submission. Several of the results of these post-unmasking reanalyses are shown in Fig. 21-22. The discussions in this section underline the importance of accounting for systematics due to model and analysis choices.', '6.1. Cosmology priors': 'Two post-unmasking analyses illustrate the impact of informative cosmology priors on parameter constraints: For the original unmasking submission, all analysis teams chose their own priors (which all differ from the priors from which the mock cosmologies are drawn), as illustrated in Fig. 5. While the posteriors on the target parameters ( Ω m , σ 8 )', 'redshift-space snapshots (mean of 10 realizations), analyzed in flat ΛCDM': 'Figure 21. 1D marginalized constraints on Ω m and σ 8 for analyses of redshift-space mocks, including post-unmasking reanalyses shown with light opacity.  \n<!-- image -->  \nare substantially narrower than the priors in these parameters for all analyses, several analyses adopted informative priors for some of the other ΛCDM parameters. The light green symbols in Fig. 21 show the impact of imposing the BACCO cosmology priors (Eq. (7)) on the EFT P+B analysis, which originally used the widest cosmology priors among all participating analysis teams (Eq. (10)). The approximately 15% improvement in 1D-marginalized constraints is primarily caused by the informative prior on n s . On the other hand, the VSF analysis chose informative priors on h , ω b and n s that excluded the light-cone mock cosmology at several σ . Post-unmasking analyses indicate only a minor bias in Ω m due to the miscentered prior and a 23% degradation in constraining power with the wider cosmology prior of the BACCO analysis. These examples highlight the importance of homogenizing cosmology priors for future detailed comparisons of different analysis methods.', '6.2. Post-unmasking BACCO analysis refinement': 'The BACCO team identified a large emulation uncertainty of the hexadecapole previously unnoticed that became significant due to the increased simulation volume of the challenge mocks compared to their previous validation (c.f. Sect. 5.1.3). The light blue symbols show a post-unmasking reanalysis of monopole and quadrupole only, with all other analysis settings held fixed.', '6.3. Post-unmasking EFT analyses comparisons': 'The EFT P+B and EFT FBI teams adopted different modeling choices in their pre-unmasking analyses, which are partly motivated by the complexities of numerical implementation and computation. To investigate the impact of these differences on their respective constraints, after unmasking, both teams have run new analyses with assumptions that aim to draw closer to the baseline analysis of the other team. In Fig. 22, we show results of these post-unmasking studies (light opacity), next to the original pre-unmasking results (full opacity).  \nTo first understand the difference between the two EFT pre-unmasking analyses, note that the EFT FBI preunmasking analysis described the matter-galaxy connection with a second-order galaxy bias expansion and assumed that the stochastic contribution to the galaxy density field, ε , is Gaussian distributed with a white power spectrum σ Poisson , effectively neglecting the coupling between the stochastic and deterministic fields of the form δ g ⊃ δ m ε (e.g. Desjacques et al. 2018). In contrast, the EFT P+B baseline analysis included (1) all second-order plus the most relevant third-order bias term, (2) the most general galaxy stochasticity model at the given order including a stochastic bispectrum and stochastic-deterministic cross bispectrum contributions (stemming from the δ m ε term in the bias expansion), and (3) O ( k 2 ) corrections to the stochasticity power spectrum.', 'Beyond-2pt Collaboration': "Figure 22. 1D marginalized constraints on σ 8 for analyses of real-space mocks, including post-unmasking reanalyses shown with light opacity.  \n<!-- image -->  \nTo see how these differences between the EFT FBI and EFT P+B analysis affect their result, the EFT P+B team have run a 'P+B restricted' analysis on the real-space mocks, at the same (pre-unmasking) scale cuts k max = { 0 . 3 , 0 . 15 } h Mpc -1 , and fixing all cosmological parameters to their true values-except σ 8 -and with A shot = B shot = b Γ 3 = a 0 = R 2 ∗ = 0 (the higher derivative bias defined as in Chudaykin et al. (2021b)). The first important observation is that EFT P+B restricted analysis recovered σ 8 with a significant bias. This is because all components described above, such as non-Gaussian stochasticity in the bispectrum, third-order operators in the galaxy bias expansion, and stochastic-deterministic coupling are non-negligible at the scales involved in the EFT P+B analyses. 30 In order to ease the comparison in the presence of this bias, the 'P+B restricted' posterior mean in Fig. 22 is re-centered to the 'EFT P+B' posterior mean in the EFT P+B baseline analysis. Thus, the 'P+B restricted' results in Fig. 22 should only be considered in terms of the error bars. A direct comparison to the 'EFT FBI' error bars shows that the nominal variances on σ 8 in both analyzes are relatively similar, suggesting that the difference between these two constraints likely comes from different analysis assumptions.  \nConversely, in their post-unmasking study, the EFT FBI team have run an 'FBI extended' analysis that includes the full set of third-order bias operators in the galaxy bias expansion (i.e., three additional bias terms compared to the 'EFT P+B' pre-unmasking analysis), as well as higher-derivative stochasticity (but still without mode coupling in the stochastic part), at their pre-unmasking scale cut k max = 0 . 1 h Mpc -1 . The 'FBI extended' posterior in Fig. 22 is completely consistent with the ground truth σ 8 , with a noticeably increased width (in between those from the 'EFT P+B' and 'P+B restricted' analyses, pre- and post-unmasking).  \nThe upshot of these post-unmasking analyses and comparisons by the two EFT teams is that we find a broad consistency between the EFT P+B and FBI results, but future work is required to extend this comparison to the full cosmology parameter space, and to redshift-space clustering.", '6.4. Towards correct cosmology constraints from the highly nonlinear regime': "While there is undoubtedly significant signal-to-noise in galaxy clustering in the highly nonlinear regime, the conversion from signal-to-noise to cosmology parameter constraints crucially relies on parameterizations of the galaxy-halo connection in the nonlinear regime being sufficiently flexible to marginalize over all modeling uncertainties that could bias cosmological parameter estimation. Furthermore, models for clustering statistics in the highly nonlinear regime are typically evaluated using simulations and emulators. This type of model evaluation is subject to uncertainties from cosmic variance, emulation errors, and finite training sample size, which may further reduce constraining power.  \nFour participating analyses utilize scales within the highly nonlinear regime, of which three rely on HOD-based models. Here, we summarize pre-unmasking validation tests and post-unmasking reanalyses that illustrate the impact of highly nonlinear scales and discuss limitations to error quantification in the highly nonlinear regime.  \nThe BACCO P analyses with k max = 0 . 2 h Mpc -1 and k max = 0 . 5 h Mpc -1 (Fig. 7) indicate a substantial gain in constraining power from strongly nonlinear scales. The BACCO analysis relies on a 'hybrid' Lagrangian galaxy bias expansion, which has been validated against HOD and SubHalo Abundance Matching extended (SHAMe, Contreras et al. 2021a,b) techniques.  \nThe SBI P+B analysis recovers the input cosmology to better than 1 -σ with a scale cut k max = 0 . 5 h Mpc -1 . While the SBI P+B team did not perform scale cut variations for the challenge mock catalogs, other analyses using the same model for galaxy samples with similar number density show substantial gains in constraining power from scales 0 . 25 h Mpc -1 < k ≤ 0 . 5 h Mpc -1 for the power spectrum (Hahn et al. 2022) and 0 . 3 h Mpc -1 < k ≤ 0 . 5 h Mpc -1 for the bispectrum (Hahn et al. 2023).  \nThe density-split clustering analysis includes measurements down to 1 h -1 Mpc . However, the gain in constraining power from 30 h -1 Mpc to 1 h -1 Mpc scales is limited (c.f. Fig. 16), as the density-split correlation functions are not sensitive to variations in the clustering below the smoothing kernel scale ( R s = 10 h -1 Mpc ) and the small-scale measurement only contains information about the smoothed density PDF. Additionally, the emulator and training-set errors contribute significantly to the total covariance on small scales, further suppressing the information gain. Finally, previous studies (Paillas et al. 2023) showed that density split statistics can efficiently extract information from AP distortions, that were not included in the current redshift-space mocks.  \nThe DDk NN baseline analysis employs the scale cut r p > 5 h -1 Mpc . Post-unmasking reanalyses with even more aggressive scale cuts shown in Fig. 12 do not pass robustness tests, indicating insufficient model flexibility or uncertainty modeling. Hence, further model refinements are required to quantify the cosmological constraining power of aggressive analyses in the 1-halo regime.", '6.4.1. Galaxy-halo connection models': 'We reiterate that the organizers communicated to all participants that the challenge mock catalogs are based on the HOD formalism, and it is important to acknowledge that HODs, like any other model for the galaxy-halo connection, come with their own set of limitations due to assumptions about galaxy formation. While the specific HOD parameterization is not revealed, this information gives HOD-based modeling approaches an inherent advantage. Among the participating analyses, the modeling of density-split clustering, k NN statistics, and SBI P+B are explicitly based on HODs. We refer to Cuesta-Lazaro et al. (2023) for parameter recovery validation of the HOD-based densitysplit statistics model on a SHAM-based galaxy mock, and Hahn et al. (2023) for parameter recovery validation of the HOD-based SBI P+B model against different HOD models and different N-body halo catalogs. The k NN postunmasking reanalyses with different HOD models (Fig. 12) show that k NN-DD analyses with aggressive 1-halo scale cuts are currently subject to either error miscalibration or model misspecification, even within the HOD framework. Future research on more flexible parameterizations will be required to pass comprehensive parameter recovery tests (see Reddick et al. 2014, for an early example of comprehensive non-HOD recovery tests) and to obtain well-calibrated error bars. Such increased modeling flexibility will likely be accompanied by a degradation in constraining power, especially when accounting for non-HOD galaxy-halo connection models (see Wechsler & Tinker 2018, for a review). Hence, the advance of novel analysis methods that exploit the highly nonlinear regime requires developing and validating more realistic mocks to establish a more accurate challenge framework that can encapsulate a broader set of models.  \nValidation of the uncertainties reported by individual methods will ultimately require a suite of (ideally, parametermasked) mock catalogs with variations in cosmology and galaxy-halo connection model. However, not just any suite will do, and the design of such a suite is essential for enabling Bayesian quantification of model performance and uncertainty calibration. Ideally, such a suite would be generated by drawing from posteriors of galaxy-halo connection model(s) that satisfy the following:  \n- 1. The mock-generating model(s) should be consistent with observations within some measurement uncertainty.\n- 2. The mock-generating model(s) should be predicated upon different assumptions from the models in the analysis being tested.  \nFor this program to be successful, realism matching simulations to observations is essential. This is a challenging task because we do not yet have a suite of simulations available that fully matches all observational data. Still, paying attention to developing accurate and flexible mock sky surveys that have the right set of included systematics should be a high priority so that the field is not misled either by missing important systematics or trying to account for unphysical systematics. Furthermore, any Bayesian quantification of model sufficiency and uncertainty validation critically depends on the (hyper) parameter priors for the test suite, which should be determined in data space rather than through arbitrary parameterization choices. This development will benefit from continued improvements in cosmological hydro-dynamical simulations (for a review see e.g. Crain & van de Voort 2023) as well as from the  \nincreasing diversity of semi-empirical methods for generating realistic mock galaxy populations (e.g. Behroozi et al. 2019; Hearin et al. 2021; Wechsler et al. 2022; Hearin et al. 2023; Kwon et al. 2023).', '7. CONCLUSION': "It is well-established that the galaxy density field contains valuable information beyond the power spectrum, and many 'novel' statistics and analysis methods have emerged and matured over the last decade to exploit this information. This makes it timely to survey the state of these analyses. In this paper, we present a parameter-masked mock challenge for beyond-2pt galaxy clustering statistics. The challenge data set consists of mock catalogs created from N-body simulations with a flat Λ CDM cosmology and HOD galaxy-halo connection models, with parameter values known only to the organizers. While all parameter values are masked, analysis teams optimize scale cuts and determine other analysis choices (e.g., nuisance parameters and their priors) for each analysis method and submit one result per method for unmasking. Upon unmasking, the organizers share plots of the relative parameter biases and uncertainties of the target parameters ( Ω m , σ 8 ) but do not share information on other cosmological parameters, HOD parameterization, or HOD parameters. Post-unmasking analyses are encouraged to enable continued method development but need to be clearly labeled as such. The main results of this mock challenge can be summarized in three themes:  \n- · Design of pre-unmasking analysis strategies and validation studies. The priors for the mock catalogs (e.g., σ 8 ∼ U [0 . 68 , 0 . 9] and any HOD parameters; c.f. Sec. 3) are significantly broader than implicit priors in (nonmasked) data analyses from previous observations. The inability to iterate on Sec. 2 provided an incentive for the analysis teams to develop consistency checks and unmasking criteria, summarized in Sec. 5, even if this robustness comes at some cost in constraining power. These validation studies are an essential ingredient for future analyses to meet the accuracy requirements of next-generation data sets.\n- · Constraints from parameter-masked mock challenge. The unmasking results presented in Sec. 2 showcase the competitive constraining power of multiple beyond-2pt statistics and novel analysis methods in a parametermasked mock challenge. This performance of multiple statistics, as well as the associated modeling and inference frameworks, lends credibility to obtaining accurate and precise cosmology constraints via these methods. Further, the consistency across different analysis approaches enables a level of cross-validation on real data that it impossible to achieve for a single method on its own. The combination of multiple statistics and modeling, in principle, will enable the most precise constraints (e.g., P+B; see Banerjee & Abel 2021; Kreisch et al. 2022; Hou et al. 2023; Valogiannis et al. 2023; Bayer et al. 2023; Massara et al. 2023; Storey-Fisher et al. 2024, for other combinations of beyond-2pt and 2pt statistics). Accurate joint covariances and consistent models of the galaxy-halo connection across different methods will however require further research and pipeline developments.\n- · Post-unmasking method refinements. The unmasking results assess the performance of an analysis based on previous method development and masked validation tests. After unmasking, several teams further investigate the accuracy and precision of their original submission (c.f. Sec. 6). With caution against excessive fine-tuning, post-unmasking reanalyses identify directions for model refinements and methodological improvements for future analyses. Additionally, post-unmasking comparisons between different analysis methods provide a starting point for future studies contrasting different approaches to understanding the source of cosmological information beyond the conventional, linear and quasi-linear galaxy 2pt analyses.  \nWhile most of these findings can in principle be obtained by each team on individually generated mock catalogs, an externally organized parameter-masked challenge provides a common benchmark and naturally separates the validation of analysis choices and post-unmasking refinements, enabling a clearer assessment of constraining power. Hence we invite future submissions from other analysis teams and offer to update summary results plots in the beyond-2pt challenge repository with new submissions /github .  \nFrom a participant's perspective, the set-up of the challenge and unmasking procedure, as well as the exchanges with the organizers and other teams, allowed a beneficial constructive atmosphere and strongly encouraged scientific interactions. Overall, the positive work environment throughout the challenge encouraged the collaborative development of validation tests and contributed to a better understanding of the tools used by the various teams.  \nA single parameter-masked mock challenge offers no panacea: models and analysis methods for different statistics evolve (e.g., due to post-unmasking refinements motivated by this study), scale cuts and priors must be calibrated  \n- anew for different survey volumes and galaxy samples. Furthermore, this particular mock challenge was clearly  \nlabeled as consisting of HOD-based galaxies free of observational systematics. Therefore, future parameter-masked mock challenges with more realistic astrophysical and observational complications as well as with statistical precision mirroring the increasing survey volume and galaxy density of future surveys will be required to further method validation. While this challenge provides a vital stress test and performance benchmark at one specific point in cosmology and HOD parameter space, the current setup leaves the uncertainty calibration and validation against model misspecification to the individual analysis teams, using individually generated mock catalogs with parameter and model variations. As discussed in Sec. 6.4, developing a suite of mocks suitable for the validation of parameter constraints, including their uncertainties, from highly nonlinear scales at the accuracy of Stage-IV surveys will be a major research and computing project that is beyond the resources of individual groups. Hence it would be desirable for future community-wide challenges to provide suites of mock catalogs across parameter space and galaxy-halo connection models to facilitate a rigorous assessment of uncertainty quantification.  \nTo conclude, we emphasize the maturity of multiple 'novel' statistics and analysis methods that participated in this parameter-masked mock challenge. The individual constraining power of a particular statistic depends on the specific parameter space considered, and we caution against extrapolating relative constraining power in this challenge to other scenarios. The main strength of this emerging field is the complementarity of approaches, which will enable extensive cross-validations to yield reliable and competitive results.", 'ACKNOWLEDGMENTS': "This mock challenge was initiated during the Aspen Center for Physics 2022 Summer Program 'Large-Scale Structure Cosmology beyond 2-Point Statistics,' co-organized by DJ, EK, Hiranya Peiris, and FS. We are grateful to Hiranya Peiris for co-organizing this workshop and input on the initial design of this challenge, and to the Aspen Center for Physics, supported by the National Science Foundation grant PHY-1607611. We thank Steward Observatory, University of Arizona for hosting a second wokshop for participating analysis teams in spring 2023, with support from the David and Lucile Packard Foundation.  \nThe light-cone mock galaxy catalog is based on the AbacusSummit simulation light cones and we thank the AbacusSummit team for making their data products publicly available. The N -body simulations and subsequent halo catalog creation for producing the Λ CDM snapshot mocks were carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. We further acknowledge High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and RDI and maintained by the UA Research Technologies department.  \nWe are grateful to Boryana Hadzhiyska, Andrew Hearin, Johannes Lange, Ariel Sánchez, and Risa Wechsler for their valuable comments on the manuscript. We further thank Camille Avestruz, Humna Awan, Andrew Hearin, Dragan Huterer, Nick Kokron, Martin Reinecke, Marko Simonović, Julia Stadler, Masahiro Takada, Kuan Wang, Risa Wechsler, and Martin White for helpful discussions.  \nEK, YK and ANS were supported in part by David and Lucile Packard Foundation and a research fellowship from the Alfred P. Sloan foundation. CH was supported by the AI Accelerator program of the Schmidt Futures Foundation. NMNacknowledges support from the Leinweber Foundation. OHEP is a Junior Fellow of the Simons Society of Fellows. The work of TA and SY was supported by the U.S. Department of Energy SLAC Contract DE-AC02-76SF00515. KA acknowledges the support from Fostering Joint International Research (B) under Contract No. 21KK0050. MPI is supported by STFC consolidated grant no. RA5496. AP acknowledges support from the Simons Foundation to the Center for Computational Astrophysics at the Flatiron Institute, as well as support from the European Research Council (ERC) under the European Union's Horizon programme (COSMOBEST ERC funded project, grant agreement 101078174). CCL is supported by the National Science Foundation under Cooperative Agreement PHY2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions).  \nThe EFT FBI analyses were conducted on the COBRA and FREYA HPC clusters at the Max Planck Computing and Data Facility.", 'AUTHOR CONTRIBUTIONS': 'All authors contributed to the interpretation of results and reviewed the manuscript. Author contributions to individual analyses and challenge organization are listed below.  \nElisabeth Krause Challenge : challenge conceptualization and organization, workshop organization, design of mock catalogs, coordination with analysis teams. Pre-unmasking review of analysis sections, general writing and end-to-end  \nediting.  \nYosuke Kobayashi Challenge : design of mock catalogs, N-body simulations, creation of mock galaxy catalogs, communication with analysis teams.  \nAndrés N. Salcedo Challenge : design of mock catalogs, creation of mock galaxy catalogs, communication with analysis teams, writing and editing of general sections.  \nMikhail Ivanov Challenge : customized EFT P+B analysis runs for comparisons with other methods, writing and editing of general sections. EFT P+B analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nTom Abel k NN analysis : design of validation tests, interpretation of results, post-unmasking studies.  \nKazuyuki Akitsu - EFT P+B analysis : analysis runs, interpretation of results, and writing of analysis section.  \nRaul Angulo BACCO P analysis : design of validation tests, interpretation of results, post-unmasking studies.  \nGiovanni Cabass - EFT P+B analysis : design of validation tests and interpretation of results.  \nSofia Contarini - Void analysis : development of analysis pipeline, design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nCarolina Cuesta-Lazaro - Density-Split analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nChangHoon Hahn - SBI P+B analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nNico Hamaus - Void analysis : development of analysis pipeline, design of validation tests, analysis runs, interpretation of results, writing of analysis section, and analysis team coordination.  \nDonghui Jeong Challenge : challenge conceptualization and workshop organization.  \nChirag Modi - SBI P+B analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nNhat-Minh Nguyen Challenge : General writing and end-to-end editing. EFT FBI analysis : development of analysis pipeline, design of validation tests, analysis runs, interpretation of results, post-unmasking studies, and writing of analysis section.  \nTakahiro Nishimishi Challenge : development of N-body simulation code GINKAKU , optimization of accuracy parameters for the use in the challenge, writing of simulation section.  \nEnrique Paillas - Density-Split analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nMarcos Pellejero Ibañez BACCO P analysis : analysis team coordination, design of validation tests, analysis runs, interpretation of results, post-unmasking studies, and writing of analysis section.  \nOliver H. E. Philcox - EFT P+B analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nAlice Pisani - Void analysis : design of validation tests, analyses interpretation, interpretation of results, writing of analysis section, and analysis team coordination.  \nFabian Schmidt Challenge : challenge conceptualization and coordination, workshop organization. EFT FBI analysis : development of analysis pipeline, design of validation test, interpretation of results, post-unmasking studies.  \nSatoshi Tanaka Challenge : development and tuning of N-body simulation code GINKAKU .  \nGiovanni Verza - Void analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nSihan Yuan Challenge : writing of HOD overview. k NN analysis : design of validation tests, analysis runs, interpretation of results, and writing of analysis section.  \nMatteo Zennaro BACCO P analysis : design of validation tests, interpretation of results, post-unmasking studies.', 'REFERENCES': 'Nunes, R. C., Vagnozzi, S., Kumar, S., Di Valentino, E., & Mena, O. 2022, Phys. Rev. D, 105, 123506,  \nYuan, S., Garrison, L. H., Eisenstein, D. J., & Wechsler, R. H. 2022, MNRAS, doi: 10.1093/mnras/stac1830 Yuan, S., Garrison, L. H., Hadzhiyska, B., Bose, S., & Eisenstein, D. J. 2021, MNRAS, 510, 3301, doi: 10.1093/mnras/stab3355  \nYuan, S., Hadzhiyska, B., & Abel, T. 2023, Mon. Not. Roy. Astron. Soc., 520, 6283, doi: 10.1093/mnras/stad550 Yuan, S., Zamora, A., & Abel, T. 2023a, MNRAS, 522, 3935, doi: 10.1093/mnras/stad1275  \n- Yuan, S., Zhang, H., Ross, A. J., et al. 2023b, arXiv e-prints, arXiv:2306.06314,  \ndoi: 10.48550/arXiv.2306.06314  \n- Zennaro, M., Angulo, R. E., Contreras, S., Pellejero-Ibáñez, M., & Maion, F. 2022, MNRAS, 514, 5443,  \ndoi: 10.1093/mnras/stac1673  \nZennaro, M., Angulo, R. E., Pellejero-Ibáñez, M., et al. 2021, arXiv e-prints, arXiv:2101.12187. https://arxiv.org/abs/2101.12187 Zhai, Z., Percival, W. J., & Guo, H. 2023a, MNRAS, 523, 5538, doi: 10.1093/mnras/stad1793 Zhai, Z., Tinker, J. L., Banerjee, A., et al. 2023b, ApJ, 948, 99, doi: 10.3847/1538-4357/acc65b Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760, doi: 10.1086/521074 Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005,  \nApJ, 633, 791, doi: 10.1086/466510'}