# Robust cosmological inference from non-linear scales with $k$-th nearest
neighbor statistics

Sihan Yuan, Tom Abel, and Risa H. Wechsler
Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,
452 Lomita Mall, Stanford, CA 94305, USA
Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA
94305, USA
SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA
94025, USA E-mail<EMAIL_ADDRESS>0000-0002-5992-7586 0000-0003-2229-011X

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the methodology for deriving accurate and reliable cosmological
constraints from non-linear scales ($<50h^{-1}$Mpc) with $k$-th nearest
neighbor ($k$NN) statistics. We detail our methods for choosing robust minimum
scale cuts and validating galaxy–halo connection models. Using cross-
validation, we identify the galaxy–halo model that ensures both good fits and
unbiased predictions across diverse summary statistics. We demonstrate that we
can model $k$NNs effectively down to transverse scales of $r_{p}\sim
3h^{-1}$Mpc and achieve precise and unbiased constraints on the matter density
and clustering amplitude, leading to a 2% constraint on $\sigma_{8}$. Our
simulation-based model pipeline is resilient to varied model systematics,
spanning simulation codes, halo finding, and cosmology priors. We demonstrate
the effectiveness of this approach through an application to the Beyond-2p
mock challenge. We propose further explorations to test more complex
galaxy–halo connection models and tackle potential observational systematics.

###### keywords:

cosmology: large-scale structure of Universe – galaxies: haloes – methods:
statistical – methods: numerical

††pubyear: 2022††pagerange: Robust cosmological inference from non-linear
scales with $k$-th nearest neighbor statistics–B

##  1 Introduction

The spatial distribution of galaxies presents one of the most powerful probes
of the fundamental properties of the universe. A new generation of wide-area
spectroscopic surveys, such as the Dark Energy Spectroscopic Instrument (DESI;
Levi et al., 2013), the Subaru Prime Focus Spectrograph (PFS; Takada et al.,
2014), the ESA Euclid satellite mission (Laureijs et al., 2011), and the NASA
Roman Space Telescope (WFIRST; Spergel et al., 2013) will enable us to study
the 3D large-scale structure (LSS) of the universe with unprecedented
precision. This precision promises insights into critical questions, including
the cause of the universe’s accelerated expansion and the theory of gravity,
the nature of dark matter and neutrinos, the physics of the primordial
universe, and the detailed process of galaxy formation.

Over recent decades, LSS analyses have seen significant advancements. On large
scales ($>50\,h^{-1}\,$Mpc), concise analytic models based on perturbation
theory (of the density contrast) are sufficiently accurate. Given that density
field on large scales is well approximated by a Gaussian random field, the
2-point correlation function (2PCF, or its Fourier pair the power spectrum)
captures its complete information content. Consequently, large-scale 2PCF
analyses using perturbative methods have become pivotal in modern cosmology.
However, there is a wealth of additional information on smaller, non-linear
scales. These scales contain complex features in the density field and offer
the highest signal-to-noise clustering measurements, but pose significant
analytical challenges due to the intermingling of cosmology and galaxy
physics.

Addressing these challenges and harnessing the full potential of small-scale
information necessitate alternative modeling frameworks. While gravity-only
$N$-body simulations can predict matter density accurately to smaller scales,
the sheer computational expense and the precision requirements of modern
surveys mean we need large volume, high-resolution simulations. To bridge this
gap, surrogate models or emulators interpolate between a limited set of
simulated cosmologies, providing efficient predictors for arbitrary
cosmologies (e.g. DeRose et al., 2019b; Nishimichi et al., 2019; Maksimova et
al., 2021). This method, termed simulation-based modeling, has been
successfully implemented in various recent studies (e.g. Lange et al., 2022;
Yuan et al., 2022a; Kobayashi et al., 2022; Chapman et al., 2022; Zhai et al.,
2023).

Still, the credibility of these models on smaller scales remains a concern.
Large-scale perturbative models use a sequence of bias parameters, keeping
them physics-agnostic, to marginalise over the connection between galaxies and
dark matter. On smaller scales, however, simulation-based methods require
well-motivated models linking galaxies to dark matter haloes. The most
commonly used approach is Halo Occupation Distribution (HOD; Peacock & Smith,
2000; Scoccimarro et al., 2001; Berlind & Weinberg, 2002; Berlind et al.,
2003; Zheng et al., 2005, 2007) models due to its speed and flexibility. In
its simplest form, only halo mass is used, and the HOD can be summarised as
$P(N_{\mathrm{galaxy}}|M_{\mathrm{halo}})$; a commonly used form of this model
is summarised by Zheng et al. (2007) (also see section 3.2). More recent
studies using additional data and summary statistics have found the need to
include a secondary halo property (e.g. Yuan et al., 2021; Wang et al., 2022;
Beltz-Mohrmann et al., 2023; Contreras et al., 2023), an effect known as
galaxy assembly bias (Gao et al., 2005; Wechsler et al., 2006). Although the
HOD model has undergone extensive refinement, its inherent empirical nature
and simplicity indicates that it may fail in detail when confronted with ever
more precise data. Thus, a pivotal question is ensuring that the model
connecting connect galaxies to dark matter is sufficiently flexible to avoid
systematic biases in the cosmology inference.

To date, most simulation-based studies have adopted a fiducial model for the
galaxy–halo connection without conducting extensive tests to demonstrate its
robustness. Some studies have shown convergence between a few different HOD
models (Yuan et al., 2022a), or for a wider range of models for the
galaxy–halo connection (Reddick et al., 2014), but they have not typically
demonstrated that such models span the necessary model space. This is
essential in the era of precision cosmology given that it has been shown that
different galaxy assembly bias assumptions can impact cosmological parameters
(Lange et al., 2019). Here, we present new tests that can significantly
increase our confidence in the assumed galaxy model.

Beyond the question of robustness, modeling non-linear scales requires
advanced summary statistics. While the large-scale density field is adequately
represented by the 2PCF, non-linear scales require a more complex approach.
The $k$-th nearest neighbor statistics ($k$NN Banerjee & Abel, 2021a; Banerjee
et al., 2022; Yuan et al., 2023a) have emerged as a promising solution due to
their informative nature, computational ease, and interpretative clarity.

In this work, we describe in detail the methodology we use to derive precise
yet dependable cosmology constraints with $k$NNs from non-linear scales in the
Beyond-2p blind mock challenge (Beyond-2p collaboration et al., prep). The
paper is structured as follows. Section 2 describes the blind mock data.
Section 3 details our simulation-based modeling framework, including
simulations, the HOD, our summary statistics, and emulation. Our results,
accompanied by methods for selecting scale cuts and validating HOD modeling,
are presented in section 4. We discuss our methodology’s limitations and
forecast potential advancements in section 5. Conclusions are drawn in section
6.

blackThroughout this paper, we adopt the Planck 2018 $\Lambda$CDM cosmology
($\Omega_{c}h^{2}=0.1200$, $\Omega_{b}h^{2}=0.02237$, $\sigma_{8}=0.811355$,
$n_{s}=0.9649$, $h=0.6736$, $w_{0}=-1$, and $w_{a}=0$) (Planck Collaboration
et al., 2020).

##  2 Blind mock data

To test our cosmology inference framework, we use the blind mock galaxy
catalogs described in Beyond-2p collaboration et al. (prep). blackThe
significance of this mock lies in that the true cosmology deviates
significantly from Planck and is kept secret until all analyses were done and
finalised. The galaxy model is also kept secret and will be hidden in the
future to enable future participation in the blind challenge. In context, this
is the first blind mock challenge specifically designed for analyses pipelines
that target non-linear scales, with sufficient volume to match the new
generation of cosmology surveys.

The mock catalogs are created by populating galaxies in $N$-body simulation
halo catalogs at $z=1$, following an HOD prescription designed by the creators
of the Beyond-2p challenge. While the mocks assume a flat $\Lambda$CDM
cosmology in which the CMB acoustic scale $\theta_{\star}$ is fixed to the
value adopted in AbacusSummit black(see section 3.1), the input values of
other cosmological parameters and the HOD model are kept blind. For creating
the halo catalogs, dark-matter-only simulations are run in
$(2\,h^{-1}\,\mathrm{Gpc})^{3}$ boxes using the GINKAKU $N$-body solver
(Nishimichi et al., prep). The initial conditions are generated using the
second-order Lagrangian perturbation theory (Crocce et al., 2006) at $z=49$,
and haloes are identified using the Rockstar halo finder (Behroozi et al.,
2013). The redshift-space distortions are implemented by modulating each
galaxy’s position along the $z$-axis with the redshift-space displacement
determined by its velocity. blackAn additional ten
$(2\,h^{-1}\,\mathrm{Gpc})^{3}$ boxes are generated with the same cosmology
and HOD, but different initial phases for covariance calculations.

blackIt is important to point out that these mocks are tuned to roughly
produce the correct spatial number density and linear bias of DESI Luminous
Red Galaxies, but we were not given any additional details about the
underlying model beyond the fact that it is an HOD. We do not know the form of
the HOD, the value of any HOD parameters, or the inclusion of any non-vanilla
extensions such as galaxy assembly bias. The fact that the HOD is implemented
on top of a new simulation and halo code means that there will be significant
and complex differences between the mock HOD and any existing implementations.
A key aim of this paper is to quantify to what extent we can marginalise over
these differences and uncertainties about the model.

##  3 Methodology

We use a simulation-based model for this analysis. In this section, we
introduce the various layers of the model, including the simulation suite, the
HOD model, the $k$NN and 2PCF summary statistics, and our neural-net-based
emulator.

###  3.1 AbacusSummit simulations

The AbacusSummit simulation suite (Maksimova et al., 2021) is a set of large,
high-accuracy cosmological N-body simulations using the Abacus N-body code
(Garrison et al., 2019, 2021), designed to 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. For this
analysis, we use exclusively the “base” configuration boxes in the simulation
suite, each of which contains $6912^{3}$ particles within a
$(2h^{-1}$Gpc$)^{3}$ volume, corresponding to a particle mass of $2.1\times
10^{9}h^{-1}M_{\odot}$. 111For more details, see
https://abacussummit.readthedocs.io/en/latest/abacussummit.html The
AbacusSummit suite also uses a specialised spherical-overdensity based halo
finder known as CompaSO (Hadzhiyska et al., 2022c).

Figure 1: Training and test cosmologies in $\Lambda$CDM parameter space. Top
panels in each row show the distribution of each parameter. Note that we omit
$h$ as it is fixed to the acoustic scale.

Here we mainly use the AbacusSummit cosmology grid, a set of 85 base boxes run
at 85 distinct cosmologies, at fixed initial phase. The 85 cosmologies are
tagged c000-181 non-consecutively and visualised in Figure 1 (see Maksimova et
al. 2021 for details). The details of each cosmology are described on the
AbacusSummit website.
222https://abacussummit.readthedocs.io/en/latest/cosmologies.html

While we only conduct our analysis in $\Lambda$CDM cosmology space in this
paper, we train the emulator in a larger $w$CDM+$N_{\mathrm{eff}}$+running
space with 8 parameters: the baryon density $\omega_{b}=\Omega_{b}h^{2}$, the
cold dark matter density $\omega_{\mathrm{cdm}}=\Omega_{\mathrm{cdm}}h^{2}$,
the amplitude of structure $\sigma_{8}$, the spectral tilt $n_{s}$, running of
the spectral tilt $\alpha_{s}$, the density of massless relics
$N_{\mathrm{eff}}$, and dark energy equation of state parameters $w_{0}$ and
$w_{a}$ ($w(a)=w_{0}+(1-a)w_{a}$).

A set of observational constraints were followed in the design of the
parameter grid. For example, $w_{a}$ is varied while holding the equation of
state at $z=0.333$ constant, so that the low-redshift cosmic distance scale is
minimally changed. Similarly, changes in $N_{\mathrm{eff}}$ are compensated by
changes in $\omega_{\mathrm{cdm}}$ and $n_{s}$, to minimise changes in the
CMB. A flat spatial curvature is assumed for all cosmologies, and the Hubble
constant $H_{0}$ is chosen to match the CMB acoustic scale $\theta_{*}$ as
measured by Planck 2018 (P18) cosmology (Planck Collaboration et al., 2020).
We briefly summarise the cosmology choices here.

c000 corresponds to the fiducial P18 cosmology, specifically the mean
estimates of the Planck TT,TE,EE+lowE+lensing likelihood chains. The 25 boxes
with different phases are at c000.

c001-004 correspond to four secondary cosmologies. Specifically, a WMAP9 + ACT
+ SPT (Calabrese et al., 2017), a thawing dark energy model
($w_{0}=-0.7,w_{a}=-0.5$), a model with extra relativistic density
($N_{\mathrm{eff}}=3.7$), and a model with lower amplitude of clustering (c000
but with $\sigma_{8}$ = 0.75).

c100-126 is a linear derivative grid set up around c000, with symmetric pairs
along all 8 parameter axes. The grid also includes additional pairs along
$\omega_{\mathrm{cdm}}$, $\sigma_{8}$, $n_{s}$, $w_{0}$, and $w_{a}$ with
smaller step sizes.

c130-181 forms an emulator grid that provides a wider coverage of the
8-dimensional parameter space. This was done by placing electrostatic points
on the surface of an 8-dimensional ellipsoid, whose extent was chosen to be 3
to 8 standard deviations beyond current constraints from the combination of
CMB and large-scale structure data. The grid excludes anti-podal reflected
points and includes extra excursion along $\sigma_{8}$.

These 85 cosmologies span the $w$CDM+$N_{\mathrm{eff}}$+running parameter
space and form the basis of our simulation-based model. It is worth pointing
out that while we train our emulators over this broader cosmology model space,
the analysis presented in the rest of this paper only considers the
$\Lambda$CDM space, fixing the other parameters to their Planck2018 values.
Additionally, the neutrino mass is fixed at 60meV and the Hubble parameter $h$
is fully degenerate with matter density by fixing the acoustic scale across
all 85 cosmologies. We refer the readers to Maksimova et al. 2021 for full
description and justification of the cosmology design.

###  3.2 The Halo Occupation Distribution (HOD)

The galaxy–halo connection model we use for generating the realistic mocks and
for the forward model is known as the Halo Occupation Distribution (HOD; e.g.
Zheng et al., 2005, 2007), which probabilistically populates dark matter
haloes with galaxies according to a set of halo properties. For a Luminous Red
Galaxy (LRG) sample, the HOD is well approximated by a vanilla model given by:

$\displaystyle\bar{n}_{\mathrm{cent}}^{\mathrm{LRG}}(M)$
$\displaystyle=\frac{f_{\mathrm{ic}}}{2}\mathrm{erfc}\left[\frac{\log_{10}(M_{\mathrm{cut}}/M)}{\sqrt{2}\sigma}\right],$
(1) $\displaystyle\bar{n}_{\mathrm{sat}}^{\mathrm{LRG}}(M)$
$\displaystyle=\left[\frac{M-\kappa
M_{\mathrm{cut}}}{M_{1}}\right]^{\alpha}\bar{n}_{\mathrm{cent}}^{\mathrm{LRG}}(M),$
(2)

where the five vanilla parameters characterizing the model are
$M_{\mathrm{cut}},M_{1},\sigma,\alpha,\kappa$. The parameter
black$M_{\mathrm{cut}}$ sets the halo mass at which the mean central
occupation is 0.5, while $M_{1}$ characterises the typical halo mass that
hosts one satellite galaxy. $\sigma$ describes the steepness of the transition
from 0 to 1 in the number of central galaxies, $\alpha$ is the power law index
on the number of satellite galaxies, and $\kappa M_{\mathrm{cut}}$ gives the
minimum halo mass to host a satellite galaxy. We have added a modulation term
$\bar{n}_{\mathrm{cent}}^{\mathrm{LRG}}(M)$ to the satellite occupation
function to remove satellites from haloes without centrals333There is evidence
that such central-less satellites may exist in a realistic stellar-mass
selected catalogue (Jiménez et al., 2019). We include this term for
consistency with previous works, but it should have minimal impact on
clustering.. We have also included an incompleteness parameter
$f_{\mathrm{ic}}$, which is a downsampling factor controlling the overall
number density of the mock galaxies. This parameter is relevant when trying to
match the observed mean density of the galaxies in addition to clustering
measurements. By definition, $0<f_{\mathrm{ic}}\leq 1$.

In addition to determining the number of galaxies per halo, the standard HOD
model also dictates the position of velocity of the galaxies. For the central
galaxy, its position and velocity are set to be the same as those of the halo
center, specifically the L2 subhalo center-of-mass for the CompaSO haloes. For
the satellite galaxies, they are randomly assigned to halo particles with
uniform weights, each satellite inheriting the position and velocity of its
host particle.

To model redshift-space distortion, we also include an additional HOD
extension known as velocity bias, which biases the velocities of the central
and satellite galaxies. This is shown to to be a necessary ingredient in
modeling BOSS LRG redshift-space clustering on small scales (e.g. Guo et al.,
2015; Yuan et al., 2021). Velocity bias has also been identified in
hydrodynamical simulations and measured to be consistent with observational
constraints (e.g. Ye et al., 2017; Yuan et al., 2022b).

blackWe parametrise velocity bias through two additional parameters:
$\alpha_{\mathrm{vel,c}}$ is the central velocity bias parameter, which
modulates the peculiar velocity of the central galaxy relative to the halo
center along the line-of-sight (LoS). Specifically in this model, the central
galaxy velocity along the LoS is thus given by

$v_{\mathrm{cent,z}}=v_{\mathrm{L2,z}}+\alpha_{\mathrm{vel,c}}\delta
v(\sigma_{\mathrm{LoS}}),$ (3)

where $v_{\mathrm{L2,z}}$ denotes the LoS component of the central subhalo
velocity, $\delta v(\sigma_{\mathrm{LoS}})$ denotes the Gaussian scatter, and
$\alpha_{\mathrm{vel,c}}$ is the central velocity bias parameter. By
definition, $\alpha_{\mathrm{vel,c}}=0$ corresponds to no central velocity
bias. We also define $\alpha_{\mathrm{vel,c}}$ as non-negative, as negative
and positive $\alpha_{c}$ are fully degenerate observationally.

blackThe second parameter is $\alpha_{\mathrm{vel,s}}$, the satellite velocity
bias parameter, which modulates how the satellite galaxy peculiar velocity
deviates from that of the local dark matter particle. Specifically, the
satellite velocity is given by

$v_{\mathrm{sat,z}}=v_{\mathrm{L2,z}}+\alpha_{\mathrm{vel,s}}(v_{\mathrm{p,z}}-v_{\mathrm{L2,z}}),$
(4)

where $v_{\mathrm{p,z}}$ denotes the line-of-sight component of particle
velocity, and $\alpha_{\mathrm{vel,s}}$ is the satellite velocity bias
parameter. $\alpha_{\mathrm{vel,s}}=1$ indicates no satellite velocity bias,
i.e. satellites perfectly track the velocity of their underlying particles.

So far, the baseline HOD is parameterised by 8 parameters,
$M_{\mathrm{cut}},M_{1},\sigma,\alpha,\kappa$, $\alpha_{\mathrm{vel,c}}$,
$\alpha_{\mathrm{vel,s}}$, and $f_{\mathrm{ic}}$. We now introduce three
physically motivated extensions in addition to the baseline HOD:

  * •

$A_{\mathrm{cent}}$ or $A_{\mathrm{sat}}$ are the concentration-based
secondary bias parameters for centrals and satellites, respectively. Also
known as galaxy assembly bias parameters. $A_{\mathrm{cent}}=0$ and
$A_{\mathrm{sat}}=0$ indicate no concentration-based secondary bias in the
centrals and satellites occupation, respectively. A positive $A$ indicates a
preference for lower-concentration haloes, and vice versa.

  * •

$B_{\mathrm{cent}}$ or $B_{\mathrm{sat}}$ are the environment-based secondary
bias parameters for centrals and satellites, respectively. The environment is
defined as the mass density within a $r_{\mathrm{env}}=5h^{-1}$Mpc tophat of
the halo center, excluding the halo itself. $B_{\mathrm{cent}}=0$ and
$B_{\mathrm{sat}}=0$ indicate no environment-based secondary bias. A positive
$B$ indicates a preference for haloes in less dense environments, and vice
versa.

  * •

$s$ is the baryon feedback parameter, which modulates how the radial
distribution of satellite galaxies within haloes deviates from the radial
profile of the halo (mimicking baryonic effects). $s=0$ indicates no radial
bias, i.e. satellites follow the dark matter distribution. $s>0$ indicates a
more extended (less concentrated) profile of satellites relative to the dark
matter, and vice versa. It has also been shown that galaxy selection based on
luminosity and SFR can bias the satellite profile (Orsi & Angulo, 2018). Thus,
the $s$ parameter can additionally marginalise over some selection effects.

Throughout the rest of this paper, we conduct fits with three different HOD
models, each including a subset of the three extensions. To keep future
participants blind, we refer to these models as models A, B, and C. A key
objective of this paper is to compare these three different models and
demonstrate that one model is validated by the data while the other two are
not (section 4.2). blackFor completeness, we also tested a fourth model that
includes all aforementioned extensions, but we do not see any improvements in
goodness-of-fits and model evidence so we omit that from this paper. It is
important to keep the HOD model compact so as to not dilute the cosmological
information in the data.

For computational efficiency, we adopt the highly optimised AbacusHOD
implementation, which significantly speeds up the HOD calculation per HOD
parameter combination (Yuan et al., 2021). The code is publicly available as a
part of the abacusutils package at https://github.com/abacusorg/abacusutils.
Example usage can be found at
https://abacusutils.readthedocs.io/en/latest/hod.html.

###  3.3 2D-$k$NN

blackThe $k$-th nearest neighbor ($k$NN) are highly informative statistics
that summarise the spatial distribution of data points by tabulating their
distances to a set of random or uniform volume-filling query points. Banerjee
& Abel (2021a, b) showed through Fisher forecasts that $k$NNs are highly
informative on cosmological parameters and are sensitive to all orders of
$n$-point correlation functions. Operationally, for each $k$, we identify each
query point’s $k$-th nearest data points and record their distances
$r_{i}^{k}$, where $i$ loops through all the query points and $k$ represents
the order. These $r_{i}^{k}$s can then be summarised with a cumulative
distribution function (CDF) as a function of length. The ensemble of these
CDFs for all $k$s forms our $k$NN summary statistics.

Figure 2: RD-$k$NN statistics calculated on the redshift-space Beyond-2p mock,
averaged over 10 realizations. Each panel corresponds to a $k$ order. In each
panel, the colour gradient represents the CDF as a function of projected
separation and line-of-sight separation. The solid lines show the contours of
the CDF. The dotted lines show the contours of a CDF of an un-clustered
Poisson random sample. The difference between the solid and dotted contours
indicates the signatures of clustering. We show only the first 4 $k$ orders
for brevity. Figure 3: Visualizations of the DD-$k$NN statistics calculated on
the redshift-space I mock, averaged over 10 realizations. The lines are
colours are definely similarly to Figure 2. Note that we have defined DD-$1$NN
$=1$, thus $k=2$ is the first meaningful order.

Yuan et al. (2023a) generalised the standard $k$NN to 2D to disentangle the
redshift-space distortion features from the projected galaxy clustering.
Specifically, we decompose the distance between each query–data pair $r$ into
a $r_{p}$ and a $r_{\pi}$ component, and we bin both projections into a 2D
histogram. Then we calculate a 2D CDF where each bin accumulates the counts
from all bins with smaller $r_{p}$ and $r_{\pi}$. Finally, we normalise the
cumulative counts by the total number of query points. Conceptually, the 2D
$k$NN-CDF is exactly analogous to the default 1D $k$NN-CDF, except we tabulate
distances and the cumulative statistics in 2D. Yuan et al. (2023a) showed that
the 2D $k$NNs are richly informative and fully derive other summary statistics
such as the 2-point correlation function (2PCF), counts-in-cells, and the void
probability function (VPF). The study also showed that the 2D $k$NNs place
tight constraints on HOD parameters in a realistic setting.

For this analysis, we adopt the 2D $k$NN statistics to comprehensively extract
information about galaxy density and clustering on non-linear scales.
Specifically, we separately analyze both the query–data $k$NN (RD-$k$NN) and
the data–data $k$NN (DD-$k$NN). blackYuan et al. (2023a) described two flavors
of $k$NNs in detail, but we briefly describe them here. RD-$k$NN represents
the standard setup, where we put down a volume-filling set of uniform or
random query points and we tabulate their distances to the data points. With
DD-$k$NN, we use the data points as query points, thus we tabulate the
distance between data points and data points.

Although the two statistics might seem similar at first glance, they reveal
distinct facets of the galaxy field. The RD-$k$NN is primarily a density
statistic, capturing density counts in specific volumes (counts-in-cells) and
detailing the size and distribution of voids. In contrast, the DD-$k$NN
enumerates the galaxy pairs in group and cluster regions. It also serves as a
generating function for the 2PCF and closely relates to high-order correlation
functions. The unique information provided by the two statistics becomes
apparent in their respective constraints. We also conduct cross-validation
tests between the two statistics to assess the robustness of our results. For
computational efficiency, we adopt the parallel $k$NN implementation written
in Rust444https://crates.io/crates/fnntw/ and wrapped for
Python555https://pypi.org/project/pyfnntw/.

We set up our 2D $k$NN data vector as follows: we use 8 logarithmic bins along
the $r_{p}$ direction between $0.32h^{-1}$Mpc and $63h^{-1}$Mpc, and 5
logarithmic bins along the $r_{\pi}$ direction between $1h^{-1}$Mpc and
$32h^{-1}$Mpc. We include $k=1,2,3,...,9$. For RD-$k$NN, we set up the query
points as a uniformly spaced grid with cell length $10h^{-1}$Mpc. We also
remove bins where the CDF values are less than 0.05 or greater than 0.95 as
those bins are noisy and contain little physical information. As a result, the
DD-$k$NN data vector is of length 143 whereas the RD-$k$NN data vector is of
length 114. The minimum transverse scale after removing these bins is
approximately 3$h^{-1}$Mpc for RD-$k$NN and 0.5$h^{-1}$Mpc for DD-$k$NN.

We justify our scale choices as follows. Our methodology is focused on non-
linear scales. The maximum scale in projected separation is chosen to obtain
some power from linear scales but not rely on well-known large-scale features
such as the BAO. The modest maximum projected separation also reduces any
potential biases in the jackknife covariance calculation. The maximum
separation along the line of sight is designed to capture the finger-of-god
and large-scale Kaiser effects. The fiducial minimum scale reaches into the
1-halo regime, which is well measured given the volume but which we expect to
be significantly affected by halo definition and baryon feedback. After
removing noisy bins, the RD-$k$NN should be largely insensitive to 1-halo
scale physics, but the DD-$k$NN remains sensitive to very small scales. Thus,
it is essential that we validate our modeling of the smallest scales for this
statistic.

Figure 2 and Figure 3 visualise the 2D $k$NN statistics up to the first four
orders. We show the measurement on the redshift-space mocks as well as the
measurement on an unclustered Poisson random sample of the same size for
comparison. The difference between the solid and dashed contours denotes the
informative features in the statistics. The corresponding covariance matrices
are calculated from 1250 jackknife volumes cut from the provided boxes. Each
box is of length $400h^{-1}$Mpc, sufficient for the scales considered. The
covariance matrices are visualised in Figure 13 and Figure 14, respectively.

Figure 4: A qualitative assessment of the Gaussianity of the likelihoods for
RD-$k$NN (top panel) and DD-$k$NN (bottom panel). The histograms show the
distribution of $\chi^{2}$ values as measured from 1883 AbacusSummit small
boxes. The solid line shows an analytic $\chi^{2}$ distribution with degrees
of freedom set to the number of bins, as we expect from Gaussian statistics.
Both statistics are consistent with the analytic distributions and with our
assumption of a Gaussian likelihood.

For comparison, we also make several references to the 2-point correlation
function (2PCF) in this analysis. Specifically, we make use of the redshift-
space 2PCF $\xi(r_{p},r_{\pi})$, which can be computed using the Landy &
Szalay (1993) estimator:

$\xi(r_{p},r_{\pi})=\frac{DD-2DR+RR}{RR},$ (5)

where $DD$, $DR$, and $RR$ are the normalised numbers of data–data,
data–random, and random–random pair counts in each bin of $(r_{p},r_{\pi})$;
$r_{p}$ and $r_{\pi}$ are the transverse and line-of-sight separations in
comoving units. The redshift space 2PCF $\xi(r_{p},r_{\pi})$ in principle
represents the full information content of the 2PCF.

For the rest of the paper, we adopt a Gaussian likelihood function for both
$k$NN statistics, specifically

$\log
L=-\frac{1}{2}(x_{\mathrm{proposed}}-x_{\mathrm{target}})^{T}\boldsymbol{C}^{-1}(x_{\mathrm{proposed}}-x_{\mathrm{target}}),$
(6)

where $x$ is the $k$NN data vector and $\boldsymbol{C}$ is the covariance
matrix. We do not include a mean density term in the likelihood as the mean
density information is naturally captured by the $k$NNs.

We test the assumption of a Gaussian likelihood as follows. If the likelihood
of the summary statistic is Gaussian distributed, the $\chi^{2}$ values should
also follow a $\chi^{2}$ distribution with degrees of freedom determined by
the number of bins. Figure 4 shows the comparison of 1883 independent
realizations of mock $k$NN statistics calculated on the AbacusSummit small
boxes with analytic $\chi^{2}$ distributions. Both our statistics are
consistent with a $\chi^{2}$ distribution and thus our Gaussian assumptions
are valid. There is a high-$\chi^{2}$ tail that exceeds the analytic
prediction in the DD-$k$NN case. Additional scale cuts, employed in the
following sections, remove this excess.

###  3.4 Forward model and emulator

blackHaving defined the summary statistics and likelihood function, we now
describe our modeling methodology that predicts the relevant summary
statistics given arbitrary input parameters in cosmology and HOD.
Specifically, we use a forward model as follows. Starting from dark matter
only simulations parameterised with cosmological parameters, we populate the
simulated haloes with galaxies using parameterised HOD models. Then we compute
the $k$NN statistics on the resulting mock galaxy density. We compute the
likelihood function by comparing the mock predicted $k$NN with the measured
$k$NN on the target sample, incorporating the jackknife covariance matrix and
assuming a Gaussian likelihood function.

Due to the high computational cost of running large, high-resolution
simulations, we cannot run arbitrary cosmologies on the fly, but instead, we
have to rely on the 85 cosmologies in the AbacusSummit suite to build an
approximate emulator model. To achieve this, we populate each cosmology box
with a few thousand of HOD models and record the resulting $k$NNs to form a
large training set. Specifically, we follow the approach of Yuan et al.
(2022a), 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
approach limits the subsequent emulator training to a compact region in the
cosmology+HOD parameter space, reducing the effects of outliers and improving
the emulator precision.

Having selected the training sample, we now build an emulator or surrogate
model that takes the cosmology and HOD parameters as inputs and outputs the
desired $k$NN statistics as outputs. For the emulator 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.

To test the performance of the emulator, we remove 9 cosmologies (c001-004 and
c171-175) from the training set and reserve them as outsample tests. Figure 1
visualises the distribution of the training and test sets in $\Lambda$CDM
parameter space. At each of the 9 cosmologies, we randomly sample 400 HODs
with $\log L>-9000$ and combine them to form a test set of size 3600. We
compute the mean absolute error of the emulator on all $9\times 400$ test
samples and find a mean error that is sub-dominant to the data error,
approximately 30–60$\%$ of the jackknife error. We further conduct tests to
ensure the emulator is not suffering from systematic biases towards the mean
(see Figure 10 of Yuan et al. 2023b). Finally, we summarise the emulator
errors by computing an emulator covariance matrix from the hold-out test
samples and add it to our jackknife covariance matrices. We provide the
emulator correlation matrices in Figure 16 and Figure 15. The final covariance
matrices going into our likelihood function is a combination of data sample
variance ($\mathbf{C^{\mathrm{jackknife}}}$, accounting for the phase
difference between the data and the model) and emulator errors
($\mathbf{C^{\mathrm{emulator}}}$):

$\mathbf{C^{\mathrm{comb}}}=\mathbf{C^{\mathrm{jackknife}}}+\mathbf{C^{\mathrm{emulator}}}$
(7)

To summarise, our model starts from the AbacusSummit simulations, and we model
galaxy distribution with AbacusHOD, invoking several extensions to the
standard HOD model. We then compute the desired summary statistics and emulate
that as a function of cosmology and HOD using a neural net model. We finally
confront our forward model with the mock data by sampling the parameter
posteriors and deriving constraints.

##  4 Results

In this section, we present the ensemble results with different analysis
choices and then describe the validation tests and model selection tests
carried out to identify our final results presented in Beyond-2p collaboration
et al. (prep). We adopt flat priors for each parameter, and Table 1 summarises
the prior bounds. For posterior sampling, we employ the efficient nested
sampling package dynesty (Speagle & Barbary, 2018; Speagle, 2020). We initiate
each nested sampling chain with 2000 live points and a stopping criterion of
$d\log\mathcal{Z}=0.01$, where $\mathcal{Z}$ is the evidence.

Parameter | Bounds
---|---
$\ln\omega_{b}$ | [-4.5, -2]
$\ln\omega_{\mathrm{cdm}}$ | [-3, -1]
$\sigma_{8}$ | [0.5, 1.1]
$n_{s}$ | [0.8, 1.2]
$\log_{10}M_{\mathrm{cut}}$ | $[12.0,14.5]$
$\log_{10}M_{1}$ | $[13.0,16.0]$
$\log_{10}\sigma$ | $[-3.5,1.5]$
$\alpha$ | $[0.5,1.5]$
$\kappa$ | $[0.0,2.0]$
$\alpha_{c}$ | [0.0, 1.0]
$\alpha_{s}$ | [ 0.2, 1.8]
$B_{\mathrm{cent}}$ | $[-1.0,1.0]$
$B_{\mathrm{sat}}$ | $[-1.0,1.0]$
$s$ | $[-1.0,1.0]$
Table 1: Prior bounds for the $k$NN analyses, for the cosmology parameters,
HOD parameters, and assembly bias parameters.

We note that all analyses in this section were conducted blind and all the
results were finalised before unblinding. We do not present or make reference
to any analysis done post-unblinding. We also do not omit blind analyses that
might be deemed unsatisfactory post unblinding. However, the texts presented
in this section were written post-unblinding and make reference to the true
values for the sake of discussion.

###  4.1 Cosmology posteriors

Figure 5 presents the 2D marginalised posteriors when fitting the RD-$k$NN in
the $\Lambda$CDM cosmology space. Note that for this analysis, we fix
$w_{0}=-1$ and $w_{a}=0$. The Hubble parameter $h$ is fully degenerate with
matter density by fixing the acoustic scale. The three colours correspond to
constraints from the different HOD models. The legend summarises the goodness-
of-fits, and the contours show the 1-2$\sigma$ constraints.

Figure 5: The $\Lambda$CDM cosmology posteriors inferred with the RD-$k$NN
data vector using three different HOD models. The axes show the difference
between the inferred and true values. The three HOD models are summarised in
section 3.2. For the RD-$k$NN data vector, we use $k=1,2,3,...,9$, covering
scales $3h^{-1}\mathrm{Mpc}<r_{p}<63h^{-1}\mathrm{Mpc}$ and
$0.5h^{-1}\mathrm{Mpc}<r_{\pi}<31.6h^{-1}\mathrm{Mpc}$. The fainter solid
lines correspond to maximum likelihood points. We omit $h$ as it is fixed to
the acoustic scale. The three contours show good consistency with each other
and the input values. blackThe contours correspond to $68\%$ and $95\%$
confidence intervals.

We see that all three models correctly recover the true cosmology and derive
mutually consistent cosmology. This shows that the RD-$k$NN is a powerful
statistic that captures non-linear information without being necessarily
sensitive to the details of HOD modeling. This is perhaps not surprising as we
have configured it with a query spacing of $10h^{-1}$Mpc, thus making it
mostly sensitive to densities on scales of a few Megaparsecs or larger and
removing much of the sensitivity to 1-halo physics. It is also worth reminding
that the RD-$k$NN is a density statistic that measures the configuration of
mass around volume-sampled queries, meaning that it is more sensitive to the
mass distribution around voids. This distinguishes RD-$k$NN from clustering
statistics (DD-$k$NN and 2PCF) that measure galaxy–galaxy distances, which
sample mostly highly clustered regions that are strongly affected by 1-halo
physics.

Figure 6: The $\Lambda$CDM cosmology posteriors inferred with the DD-$k$NN
data vector using three different HOD models. The axes show the difference
between the inferred and true values. For the DD-$k$NN data vector, we use
$k=2,3,...,9$, covering scales
$0.67h^{-1}\mathrm{Mpc}<r_{p}<63h^{-1}\mathrm{Mpc}$ and
$0.5h^{-1}\mathrm{Mpc}<r_{\pi}<31.6h^{-1}\mathrm{Mpc}$. The three contours
show mild disagreement at the 2-3$\sigma$ level. The goodness-of-fit as shown
in terms of $\chi^{2}$/d.o.f. is poor. blackThe contours correspond to $68\%$
and $95\%$ confidence intervals.

In Figure 6, we show the cosmology posteriors obtained from fitting the
DD-$k$NN. The three models are mutually in up to $2\sigma$ tension with each
other. The large goodness-of-fit values indicate that our model is not
flexible enough to produce the features in the DD-$k$NN. Given these
discrepancies, we either need more flexible models or must apply further scale
cuts to achieve reliable results. Note that this decision was made before
unblinding, and was based primarily on observed tensions among the models and
the poor $\chi^{2}$ value. Without access to additional information about the
target galaxy sample — such as spectral energy distributions, selection
criteria, or imaging data — that would typically be available in a real
survey, we do not have straightforward ways to decouple the galaxy formation
model from cosmology. Hence, we choose to use additional scale cuts to improve
the fit.

Figure 7: The best-fit $\chi^{2}$/d.o.f. when using different minimum scale
cuts. Different colours refer to the three different HOD models. The $x$-axis
iterates through different scale cuts (in units of $h^{-1}$Mpc), showing the
fits grouped vertically by data vector. DD-$k$NN only achieves a good fit at
$r_{p}>5h^{-1}$Mpc; RD-$k$NN achieves a good fit at smaller scales. This
highlights the DD-$k$NN’s high sensitivity to galaxy bias and the details of
the HOD model on small scales.

Figure 7 shows the best-fit $\chi^{2}$ as a function of minimum scale cuts.
The first two columns show the fiducial cuts used for the RD-$k$NN and
DD-$k$NN data vectors. The last two columns show best-fit $\chi^{2}$ when
progressively larger minimum scale cuts are applied to the DD-$k$NNs. We see
that we only achieve a good fit $\chi^{2}$/d.o.f. $\approx 1$ at
$r_{p}>5h^{-1}$Mpc. Figure 8 shows the DD-$k$NN cosmology posterior when an
additional cut $r_{p}>5h^{-1}$Mpc is applied. We see that indeed, the
posterior constraints of the three models are now consistent with each other.
Thus, we conclude that our current HODs can only reliably model DD-$k$NN at
scales greater than $5h^{-1}$Mpc, and we only report our DD-$k$NN constraints
at $r_{p}>5h^{-1}$Mpc for the blind challenge. However, we still include the
full-scale constraints in our discussions to motivate future work.

Figure 8: The $\Lambda$CDM cosmology posteriors inferred with the DD-$k$NN
data vector at $r_{p}>5h^{-1}$Mpc using three different HOD models. Again, we
use $k=2,3,...,9$, covering scales
$5h^{-1}\mathrm{Mpc}<r_{p}<63h^{-1}\mathrm{Mpc}$ and
$0.5h^{-1}\mathrm{Mpc}<r_{\pi}<31.6h^{-1}\mathrm{Mpc}$. The three contours
show good agreement, with significantly improved goodness-of-fit compared to
Figure 6. blackThe contours correspond to $68\%$ and $95\%$ confidence
intervals.

The fact that the current HOD models failed to fit DD-$k$NNs on very small
scales suggests that the DD-$k$NNs are highly informative on HOD modeling. In
Yuan et al. (2023a), we demonstrated that the DD-$k$NN is significantly more
constraining on HOD and assembly bias than the standard redshift-space 2PCF.
We showed that the 2PCF is a strict marginalization of the DD-$k$NN.
Conceptually, DD-$k$NN captures extra information because while the 2PCF
encodes the average number of neighbors of a galaxy at a certain distance, the
DD-$k$NN additionally encodes the $k$ ordering of neighbors at any distance.
Thus, DD-$k$NN additionally captures the phase space and topological
configuration information.

There are several potential reasons why our model cannot produce the DD-$k$NN
below $5h^{-1}$Mpc. First of all, without additional information on the galaxy
sample that would otherwise be available in a realistic survey, we cannot
meaningfully improve our galaxy–halo modeling. For example, the photometric
selection allows us to assess the completeness of the sample at different
magnitudes. The colours and morphology would inform us of the galaxy type. The
spectra can additionally give us statistical descriptions of key galaxy
properties such as stellar mass, star formation rate, and metallicity. All
these sources of information are useful in building up priors on the galaxy
model. For example, Wang et al. (2022) showed that a magnitude-limited sample
in SDSS can be robustly described by an HOD model that includes assembly bias.
A series of recent studies have also combined realistic selections with state-
of-the-art hydrodynamical simulations and semi-analytic models of galaxy
formation to build up priors on the appropriate galaxy–halo connection model
(e.g. Xu et al., 2021; Hadzhiyska et al., 2022b, a; Yuan et al., 2022b).

An additional systematic effect is that the target mock is generated with a
different simulation code and different halo finder, which can all produce
small-scale features that the DD-$k$NN is potentially sensitive to. This is an
issue that we would need to be able to marginalise over in our forward model
if we want to exploit 1-halo scales. We propose to assess this point in a
future post-unblinding analysis where we can disentangle our uncertainty about
the true HOD and the systematic effects associated with simulations and halo
finding.

To summarise, the DD-$k$NNs are potentially highly sensitive to small-scale
effects such as the details of galaxy–halo connection modeling, halo finding,
and simulation codes. For the sake of this analysis, we can protect ourselves
from these systematics by applying a relatively conservative scale cut of
$r_{p}>5h^{-1}$Mpc. In contrast, the RD-$k$NNs appear to be significantly less
sensitive to these effects as long as we choose appropriate query spacing.

Figure 9 summarises the cosmology constraints from the different analysis
choices. These include all blind analyses we carried out for the challenge.
The $y$-axis reports the difference between the inferred and true values, with
the green lines denoting 0. The blue lines show the posterior mean and
$1\sigma$ constraints. The orange triangles show the maximum likelihood
points. Immediately, all our results are remarkably consistent with each
other. We largely derive unbiased cosmology constraints, except in $n_{s}$
when we use the full-scale range in DD-$k$NN. The occasional disagreement
between the maximum likelihood points and the posterior constraints indicates
potential prior volume / projection effects. This warrants more careful
investigation in future application of our methods.

Figure 9: The marginalised $\Lambda$CDM cosmology posteriors inferred with the
two $k$NN data vectors and using three different HOD models. The blue error
bars indicate the (16, 50, 84)$\%$ posterior constraints, and orange triangles
show the maximum likelihood point. The $x$-axis iterates through the HOD
models, with the fits vertically grouped by data vector. The parameter values
are blinded. There is excellent agreement between the RD-$k$NN and DD-$k$NN
($r_{p}>5h^{-1}$Mpc). The full-scale DD-$k$NN ($r_{p}>0.67h^{-1}$Mpc) also
yields reasonably consistent constraints, albeit with significantly tighter
error bars, and moderately biased results in some parameters.

Finally, we highlight the stringent constraints on $\sigma_{8}$ from both
$k$NN statistics. The full DD-$k$NN constrains $\sigma_{8}$ error bars to less
than 0.01. Even after the scale cuts, the $1\sigma$ constraints are still
below $0.02$ for both RD-$k$NN and DD-$k$NN, significantly stronger than
existing $\sim 5\%$ constraints from galaxy clustering in BOSS/SDSS (e.g. Alam
et al., 2017; Lange et al., 2022; Kobayashi et al., 2022; Yuan et al., 2022a;
Zhai et al., 2023). While this extra precision comes directly from the
increased volume in this mock challenge relative to BOSS, what we have
demonstrated is that we can effectively take advantage of the extra volume on
non-linear scales, which is essential for taking full advantage of the data
that will be imminently available with DESI. These constraints also
demonstrate the ability of our novel statistics to capture small-scale
information. Assuming we are able to construct more informed galaxy–halo
connection models with DESI, this should enable us to push our analysis down
to smaller scales and is expected to yield even tighter constraints.

###  4.2 Validation of galaxy–halo connection modeling

To derive reliable cosmology constraints, it is essential to demonstrate that
the galaxy–halo connection model is sufficiently flexible to describe the
relevant features in the galaxy–halo physics, but not so flexible that it
dilutes the cosmology constraints. There are several ways one can inform and
validate the galaxy–halo connection model. We summarise them here:

  * •

Simulated models: Simulated galaxy models such as hydrodynamical simulations
and semi-analytic models are excellent sandboxes to gain physical intuitions
about the galaxy sample and identify necessary ingredients in their
galaxy–halo connection models. blackHowever, hydrodynamical simulations are
currently too small to calibrate our models to the precision necessary for
next-generation surveys. Semi-analytical models are cheaper and can be
generalised to larger volumes, but they are still too expensive right now to
fairly sample the model space while maintaining cosmological volumes (see
Perez et al. 2023 for the current state-of-the-art).

  * •

Galaxy observations: Leveraging existing observations and information about
the galaxy sample is a powerful and data-driven way to understand and
constrain the galaxy–halo connection. For example, one can directly constrain
the halo mass given observed galaxies via gravitational lensing. There are
also patches of the sky that are observed by a variety of facilities (such as
the COSMOS field), yielding deep images across a wide wavelength domain and
spectra. One can leverage these datasets to understand the completeness, mass
function, and other galaxy properties. These observations can be richly
informative but require careful modeling to connect to theory models of dark
matter.

  * •

Mock challenges: Showing unbiased cosmology constraints on simulated galaxy
mocks is perhaps the most direct way to show robustness. This was the primary
motivation for the Beyond-2p challenge. However, the mocks are still generated
relying on some assumptions about galaxy formation physics or the galaxy–halo
connection. Thus, while these tests help develop confidence (the larger the
variety of galaxy models, the better), it does not rule out the possibility
that unforeseen galaxy features in the real Universe can bias the inference.

  * •

Statistical validations: Statistical tests such as goodness-of-fit and cross-
validation are model-agnostic ways to show consistency with data. The idea is
to confront the model with all relevant aspects of the data and show that the
model can describe the data without internal tensions. blackIn this section,
we will first test the flexibility of our HOD models with goodness-of-fit
metrics. Then we select the “correct” HOD model via cross-validations, where
we break the data into multiple parts that contain different information and
demonstrate that the model calibrated on one part can consistently predict
another part. A sufficiently flexible model with correct physical assumptions
should achieve good fits and self-consistently predict all relevant aspects of
data. In the current work, because we have access to multiple summary
statistics that access different subsets of the information content of the
galaxy field, we can conduct this cross-validation test by fitting the model
on one summary statistic and checking if the constrained model can predict
another statistic.

In the context of this mock challenge, we cannot leverage simulated galaxy
models or galaxy observations to learn the galaxy–halo connection model, but
we can conduct statistical tests to assess model validity. To demonstrate this
point, we employed three different HOD models in our analysis, adding
different combinations of assembly bias and baryon feedback prescriptions to
the vanilla HOD. The objective of this section is to compare the three models
with goodness-of-fit metrics and cross-validation tests to show that one of
the three models is preferred by the data and dependable on the relevant
scales.

First, we compare the three models on the basis of their $\chi^{2}$/d.o.f.
(see Figure 7 and the legends of Figure 5–8). For RD-$k$NN, all three models
achieve good fits to the mock data, with models A and B showing slightly
better $\chi^{2}$. For DD-$k$NN, models B and C yield significantly better
fits on small scales. With the $r_{p}>5h^{-1}$Mpc cut, the models A and B
result in better fits. These comparisons are not conclusive but show that
model B performs consistently well in configurations. We do not show joint
fits of the DD-$k$NN and RD-$k$NN due to limitations in modeling the joint
covariance matrix.

The reduced $\chi^{2}$ metric gives a rough characterisation of the goodness-
of-fit but offers a poor assessment of possible over-fitting. Instead, we
perform cross-validation tests where we take the fits on one summary statistic
to predict another summary statistic, relying on the notion that the model
with the right physical assumptions should extrapolate to predict additional
data. Our summary statistics are particularly well suited for this test
because the different summary statistics capture disjoint information in the
density field. Specifically, while the DD-$k$NN and the 2PCF capture the
clustering information and emphasise the structure of density peaks on small
scales, the RD-$k$NN capture the density information and transition between
voids and density peaks (Yuan et al., 2023a). It is thus particularly useful
to conduct cross-validations between a clustering statistic such as (DD-$k$NN
or the 2PCF) and a density statistic (such as RD-$k$NN). blackThis idea is not
completely new. A similar cross-validation was done in Wang et al. (2022)
between clustering and counts-in-cell statistics to validate galaxy assembly
bias models. A similar cross-validation scheme was also used in recent
assessments of potential tensions between galaxy clustering and galaxy-galaxy
lensing statistics (e.g. Leauthaud et al., 2017; Yuan et al., 2020; Lange et
al., 2020; Contreras et al., 2023).

Figure 10: The predictive distribution of the RD-$k$NNs when the models are
fitted to the DD-$k$NNs. The $y$-axis shows the difference between the
RD-$k$NN prediction and the true RD-$k$NN on the blind mock, normalised by the
diagonal jackknife errors. The three shaded lines show the $1\sigma$
predictive distributions of the three different HOD models. The $x$-axis shows
the bin indices in RD-$k$NN. Each large block separated the solid vertical
dashed lines correspond to different $k$s. The sub-blocks separated by dashed
vertical lines denote different $r_{p}$s; each sub-block cycles through the
$r_{\pi}$s. Model B makes a consistent prediction for RD-$k$NNs given
DD-$k$NNs, whereas the other two models do not.

In Figure 10, we show the RD-$k$NN predictive distribution derived when the
models are fit to the DD-$k$NNs. In particular, the $y$-axis shows the
difference between the RD-$k$NN prediction with the three models and the true
RD-$k$NN on the blind mock, normalised by the jackknife errors. The three
shaded lines represent the $1\sigma$ predictions for the three HOD models.
Model B successfully predicts the correct RD-$k$NN, while the two other models
fail. blackThis indicates that model B is both sufficiently flexible to
describe the relevant aspects of the data (good $\chi^{2}$/d.o.f.) and carries
the right physical assumptions such that it can predict additional statistics
it is not fit on. It is also worth noting that all three models make reliable
predictions at the beginning of each sub-block (shown by vertical dashed
lines), corresponding to small $r_{\pi}$ separation, but fail at larger
$r_{\pi}$ values. Thus, it appears that model B might be better at capturing
the correct redshift-space distortion compared to the other models. All three
models make consistent predictions at large $k$ (as in $k$NN, not spatial
frequency), which makes sense because large $k$ orders probe larger scales and
contain less information about the galaxy–halo connection.

Figure 11: Redshift-space 2PCF $\xi(r_{p},r_{\pi})$ predictive distribution in
the case where the models are fitted to the RD-$k$NNs. The $y$-axis on the
left side shows the difference between the predicted $\xi(r_{p},r_{\pi})$ and
the true $\xi(r_{p},r_{\pi})$, normalised by the jackknife errors. We have
broken the 2D $\xi(r_{p},r_{\pi})$ data vector into different rows, each row
corresponding to a $r_{\pi}$, where the $r_{\pi}$ ranges for each row are
shown on the $y$-axis to the right. Model B successfully predicts the 2PCF
down to smaller $r_{p}$ scales than the other two models.

Figure 11 shows the redshift-space 2PCF $\xi(r_{p},r_{\pi})$ predictive
distribution when the models are fitted to the RD-$k$NNs. The $y$-axis on the
left side shows the difference between the predicted $\xi(r_{p},r_{\pi})$ and
the true $\xi(r_{p},r_{\pi})$, normalised by the jackknife errors. We have
broken the 2D $\xi(r_{p},r_{\pi})$ data vector into different rows, each row
corresponding to a $r_{\pi}$. We see that all three models can predict the
2PCF at large $r_{p}$, but none of the three models correctly predict the 2PCF
all the way down to $r_{p}=1h^{-1}$Mpc. However, comparing the three models,
model B predicts the 2PCF consistently down to $r_{p}=2h^{-1}$Mpc whereas the
other two models deviate from the true values at larger scales.

blackTo summarise, model B is strongly preferred by the mock data, both
because of good $\chi^{2}$/d.o.f. and its superior ability to predict
statistics that it was not fit on. Model B performed better in both the
RD-$k$NN prediction and the 2PCF prediction. It was particularly telling when
we fitted the three different HOD models on DD-$k$NN and only model B gave a
fully consistent prediction of the RD-$k$NN.

The trend in the disagreement shown in Figure 11 also shows that the current
failure to model scales smaller than $r_{p}<3h^{-1}$Mpc might be due to
insufficient flexibility in modeling the small-scale finger-of-god. While we
have included velocity bias in both central and satellite galaxies in our
velocity model, the parameterisation only allows for a constant fractional
shift relative to the underlying dark matter. blackThere exists evidence that
a more sophisticated velocity bias model is needed. For example, Ye et al.
(2017) reports significant mass dependencies in hydrodynamical models.
Additionally, the differences between simulation codes and halo finders can
introduce significant complexities in the galaxy velocities. This highlights
the need for more flexible velocity modeling in our pipeline. blackSeveral
other possible reasons for poor fits at very small scales include unrealistic
models used in the mock production and underestimated covariances. It is
particularly important to further validate our covariance matrices due to
their large size and significant off-diagonal contribution. We reserve that
for a future paper.

##  5 Discussion

Figure 12: The marginalised $\Lambda$CDM cosmology posteriors when using
different minimum scale cuts for the DD-$k$NN statistics. For this test, we
have fixed the HOD model to model B. The constraining power of DD-$k$NN
strongly depends on minimum scale cut. We do not vary the minimum scale of
RD-$k$NN as it would require extensive re-computation.

So far, we have presented the methodology and results of applying our $k$NN
analysis pipeline to the Beyond-2p mocks. While we derive strong and unbiased
constraints on cosmology, we highlight a few caveats and provide outlook for
future work in this section.

###  5.1 Addressed and unaddressed systematics

The fact that we derive unbiased cosmology from deeply non-linear scales is
significant for a couple of reasons. We are among the first to show that a
simulation-based method for non-linear scales can recover unbiased cosmology
results when confronted with mock data generated with a completely different
simulation. This shows that our modeling framework can be resilient against
several simulation-related systematics, including gravity codes, simulation
resolution, and halo finding (e.g. Heitmann et al., 2008; Knebe et al., 2011;
Grove et al., 2022). We attribute this resilience to the careful choice of HOD
models and scale cuts. This conclusion is made stronger by the fact that our
mock data was blind and the true cosmology was far from Planck, which shows
that we are not biased by priors.

However, one key systematic that this exercise has not fully addressed is
galaxy–halo connection modeling. Both the models and the mock data are
generated from HOD models. We first note that the HOD itself is actually a
broad framework and there is a large amount of potential systematics that can
arise from different implementations of HODs (e.g. Hearin et al., 2016). In
its most general interpretation, the HOD model simply implies a probabilistic
occupation based on halo properties (see Wechsler & Tinker, 2018 for a
review). Thus, there is a large amount of freedom in terms of the specific
function form of the HOD, the probability distribution, and the halo
properties used. Moreover, the fact that the mock data and the models used
different simulations and halo finders means that the behavior of the same HOD
can be different in complex ways. All this is to say that even within the HOD
framework, it is non-trivial to recover unbiased cosmologies from two
different HOD implementations on top of two different halo finders.

Nevertheless, the HOD is an empirical model that is necessarily wrong in
detail when compared to the real Universe (e.g. Hadzhiyska et al., 2020; Xu et
al., 2021; Beltz-Mohrmann et al., 2023). Thus, it is important to repeat this
exercise with blind mocks constructed from more sophisticated galaxy–halo
connection models. blackOne good and available option is to use subhalo
abundance matching (SHAM Vale & Ostriker, 2006; Lehmann et al., 2017), which
relies on the assumed correlation between some galaxy property with some
subhalo property to assign galaxies. SHAM is meaningfully different from the
HOD as it is based on subhalos instead of halos, and it does not rely on
analytic forms. However, SHAM has two drawbacks. First, it is similar to HODs
in the sense that it is still commonly based on halo mass or other proxies of
mass. Second, it also requires large high-resolution simulations with
effective subhalo finding, which can be expensive. Another option is to use
machine learning models to paint galaxies over from hydrodynamical simulations
or semi-analytic models (e.g. Delgado et al., 2022; Xu et al., 2021; Lovell et
al., 2022; McGibbon & Khochfar, 2022; Chittenden & Tojeiro, 2023). These
approaches automatically incorporate a variety of subtleties in galaxy
modeling and can be a powerful tool to generate complex and realistic mocks. A
third option is to use physically motivated empirical models such as
UniverseMachine and Diffstar (Behroozi et al., 2019; Alarcon et al., 2023).
These models do not rely on assumptions built into external simulations and
provide a self-consistent way of modeling galaxy assembly. Both the second and
third options should be highly differentiated from the HOD and would serve
well to test against galaxy–halo modeling systematics.

Finally, this analysis does not address observational systematics such as
survey incompleteness due to survey windows and fiber collisions, or evolution
effects over finite redshift ranges. These effects are non-trivial to correct
for when calculating novel summary statistics such as the $k$NNs. However, in
principle we can forward model these effects using simulation lightcones (Yuan
et al., 2023b; Hahn et al., 2022). Specifically, in Yuan et al. (2023b), we
present both a flexible and performant framework to explore HOD and redshift
evolution models on lightcones, and techniques to correct for survey boundary
and fiber collision effects. Similar realistic lightcone mocks were also
constructed and tested for large photometric surveys in DeRose et al. (2019a),
MacCrann et al. (2018), and To et al. (2021). Combining the forward modeling
framework and the modeling techniques laid out in this paper will be powerful
in exploiting non-linear scale information in upcoming spectroscopic surveys.

###  5.2 The impact of scale cuts

In section 4.1, we chose a fiducial scale cut of $r_{p}>5h^{-1}$Mpc for the
DD-$k$NN statistics. Such a scale cut is indeed necessary for this specific
analysis as our galaxy–halo connection model is not sufficiently flexible to
self-consistently predict features on smaller scales. However, we also pointed
out that in a realistic survey, we would have considerably more information on
the target galaxy sample to better inform the galaxy–halo connection model.
Thus, it is interesting to discuss how our constraints would improve as we
push down to smaller scales with improved models.

Figure 12 shows the 1D marginalised constraints on cosmology as we include
smaller $r_{p}$ in DD-$k$NNs. Clearly, the constraining power of the DD-$k$NN
statistics depends strongly on the minimum scale cuts. At
$r_{p}>0.67h^{-1}$Mpc, we see a 4 times improvement in the $\sigma_{8}$ and
$n_{s}$ constraints and a 2-3 times improvement in the constraints on the
density parameters. This demonstrates that while we can already derive highly
competitive constraints with our fiducial scale cut of $r_{p}>5h^{-1}$Mpc, the
true power of the DD-$k$NN statistics lies in the smaller scales. This
strongly motivates the need for a more informed yet flexible model of
galaxy–halo connection in the advent of next-generation surveys such as DESI.
We note that these constraints are derived with a fixed galaxy–halo connection
model that we already showed to be insufficient to describe the smaller
scales. Thus, while the improvement in constraining power is informative of
the next priority in further developing simulation-based models, it should not
be interpreted as a detailed forecast for future analysis. The exact
constraining power will depend on the modeling specifics.

It is also important to briefly discuss the impact of maximum scale cuts. In
this analysis, we are exclusively interested in non-linear scales, so we chose
a modest maximum scale cut of $r_{p}<63h^{-1}$Mpc, which reaches into the
linear regime but does not capture key large-scale observables such as the BAO
or the matter-radiation equality scale. In principle, we can gain sensitivity
to those scales by extending the $r_{p}$ range and including much higher $k$
orders. We expect this expansion to significantly improve our constraints on
the mass density parameters but also to have minimal effect on our
$\sigma_{8}$ and $n_{s}$ constraints, as those come predominantly from small
scales and RSD. We reserve a full-scale $k$NN analysis for a future paper.

##  6 Conclusions

The key to unlocking the cosmological information on non-linear scales lies in
employing simulation-based models and high-order summary statistics. In this
paper, we present our methodology for recovering robust cosmology constraints
from non-linear scales, leveraging two flavors of the $k$NN statistics, one
capturing the density field (RD-$k$NN) and the other capturing the clustering
(DD-$k$NN). We confront a blind mock with three different HOD models,
including varied prescriptions of assembly bias and baryonic effects built on
top of the AbacusSummit cosmology suite. We demonstrate several key points in
this paper:

  * •

This is the first time a simulation-based model has been shown to derive
strong and unbiased cosmology constraints in a blind mock challenge of
comparable volume to a modern spectroscopic survey. The power of this result
is strengthened by the fact that the mock uses different simulation and halo
codes than the input data, and by the fact that the underlying cosmology is
significantly different than Planck.

  * •

We demonstrated the flexibility and resilience of our approach using multiple
HOD models by employing a mix of goodness-of-fit metrics and cross-validation
tests. These tests give a clear preference for one of the models, which is
able to reproduce a second statistic after being fit to a first. On real data,
more information may be available to further inform the HOD model tests.

  * •

We use goodness-of-fit metrics to establish minimum scale cuts for our
analysis. The RD-$k$NN statistic emerged as less prone to non-linear scale
modeling errors. However, the DD-$k$NN proved to be considerably more
constraining when we employ improved models of small scales.

This work underscores the potential of using beyond-2p statistics on non-
linear scales for cosmology, especially in the context of the next generation
of spectroscopic surveys. We argue that integrating direct insights on the
galaxy–halo connection from data and simulations would likely enable analyses
on even smaller scales.

Moving forward, it would be valuable to re-do this analysis with a variety of
more sophisticated galaxy–halo connection models that are not based on HODs.
That would more convincingly demonstrate the credibility of simulation-based
modeling approaches for cosmology. It is also important to test against
observational systematics such as survey incompleteness and redshift
evolution, for example, using a lightcone-based forward modeling approach like
in Yuan et al. (2023b) and applying these systematics directly. One can also
resort to realistic non-HOD mocks such as in DeRose et al. (2019a), MacCrann
et al. (2018), and To et al. (2021). We look forward to applications of such
approaches in future analyses.

## Acknowledgements

We would like to thank Elisabeth Krause, Yosuke Kobayashi, Andres Salcedo,
Enrique Paillas, and others for useful feedback and suggestions in various
stages of this analysis.

This work was supported by the U.S. Department of Energy through grant DE-
SC0013718 and under DE-AC02-76SF00515 to SLAC National Accelerator Laboratory,
and by the Kavli Institute for Particle Astrophysics and Cosmology. This work
was performed in part at the Aspen Center for Physics, which is supported by
National Science Foundation grant PHY-2210452.

This work used resources of the National Energy Research Scientific Computing
Center (NERSC), a U.S. Department of Energy Office of Science User Facility
located at Lawrence Berkeley National Laboratory, operated under Contract No.
DE-AC02-05CH11231.

The AbacusSummit simulations were conducted at the Oak Ridge Leadership
Computing Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725, through support from projects AST135 and
AST145, the latter through the Department of Energy ALCC program.

## Data Availability

The simulation data are available at
https://abacussummit.readthedocs.io/en/latest/. The AbacusHOD code package is
publicly available as a part of the abacusutils package at
https://github.com/abacusorg/abacusutils. Example usage can be found at
https://abacusutils.readthedocs.io/en/latest/hod.html.

## References

  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
  * Alarcon et al. (2023) Alarcon A., Hearin A. P., Becker M. R., Chaves-Montero J., 2023, MNRAS, 518, 562
  * Banerjee & Abel (2021a) Banerjee A., Abel T., 2021a, MNRAS, 500, 5479
  * Banerjee & Abel (2021b) Banerjee A., Abel T., 2021b, MNRAS, 504, 2911
  * Banerjee et al. (2022) Banerjee A., Kokron N., Abel T., 2022, MNRAS, 511, 2765
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Beltz-Mohrmann et al. (2023) Beltz-Mohrmann G. D., Szewciw A. O., Berlind A. A., Sinha M., 2023, ApJ, 948, 100
  * Berlind & Weinberg (2002) Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587
  * Berlind et al. (2003) Berlind A. A., et al., 2003, ApJ, 593, 1
  * Beyond-2p collaboration et al. (prep) Beyond-2p collaboration et al. in prep, XXX
  * Calabrese et al. (2017) Calabrese E., et al., 2017, Phys. Rev. D, 95, 063525
  * Chapman et al. (2022) Chapman M. J., et al., 2022, MNRAS, 516, 617
  * Chittenden & Tojeiro (2023) Chittenden H. G., Tojeiro R., 2023, MNRAS, 518, 5670
  * Contreras et al. (2023) Contreras S., Chaves-Montero J., Angulo R. E., 2023, MNRAS,
  * Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, Mon. Not. Roy. Astron. Soc., 373, 369
  * DeRose et al. (2019a) DeRose J., et al., 2019a, arXiv e-prints, p. arXiv:1901.02401
  * DeRose et al. (2019b) DeRose J., et al., 2019b, ApJ, 875, 69
  * Delgado et al. (2022) Delgado A. M., Wadekar D., Hadzhiyska B., Bose S., Hernquist L., Ho S., 2022, MNRAS, 515, 2733
  * Gao et al. (2005) Gao L., Springel V., White S. D. M., 2005, MNRAS, 363, L66
  * Garrison et al. (2019) Garrison L. H., Eisenstein D. J., Pinto P. A., 2019, MNRAS, 485, 3370
  * Garrison et al. (2021) Garrison L. H., Eisenstein D. J., Ferrer D., Maksimova N. A., Pinto P. A., 2021, MNRAS, 508, 575
  * Grove et al. (2022) Grove C., et al., 2022, MNRAS, 515, 1854
  * Guo et al. (2015) Guo H., et al., 2015, MNRAS, 446, 578
  * Hadzhiyska et al. (2020) Hadzhiyska B., Bose S., Eisenstein D., Hernquist L., Spergel D. N., 2020, MNRAS, 493, 5506
  * Hadzhiyska et al. (2022a) Hadzhiyska B., et al., 2022a, arXiv e-prints, p. arXiv:2210.10068
  * Hadzhiyska et al. (2022b) Hadzhiyska B., et al., 2022b, arXiv e-prints, p. arXiv:2210.10072
  * Hadzhiyska et al. (2022c) Hadzhiyska B., Eisenstein D., Bose S., Garrison L. H., Maksimova N., 2022c, MNRAS, 509, 501
  * Hahn et al. (2022) Hahn C., et al., 2022, arXiv e-prints, p. arXiv:2211.00723
  * Hearin et al. (2016) Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552
  * Heitmann et al. (2008) Heitmann K., et al., 2008, Computational Science and Discovery, 1, 015003
  * Jiménez et al. (2019) Jiménez E., Contreras S., Padilla N., Zehavi I., Baugh C. M., Gonzalez-Perez V., 2019, MNRAS, 490, 3532
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Kobayashi et al. (2022) Kobayashi Y., Nishimichi T., Takada M., Miyatake H., 2022, Phys. Rev. D, 105, 083517
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64
  * Lange et al. (2019) Lange J. U., van den Bosch F. C., Zentner A. R., Wang K., Hearin A. P., Guo H., 2019, MNRAS, 490, 1870
  * Lange et al. (2020) Lange J. U., Leauthaud A., Singh S., Guo H., Zhou R., Smith T. L., Cyr-Racine F.-Y., 2020, arXiv e-prints, p. arXiv:2011.02377
  * Lange et al. (2022) Lange J. U., Hearin A. P., Leauthaud A., van den Bosch F. C., Guo H., DeRose J., 2022, MNRAS, 509, 1779
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Leauthaud et al. (2017) Leauthaud A., et al., 2017, MNRAS, 467, 3024
  * Lehmann et al. (2017) Lehmann B. V., Mao Y.-Y., Becker M. R., Skillman S. W., Wechsler R. H., 2017, ApJ, 834, 37
  * Levi et al. (2013) Levi M., et al., 2013, arXiv e-prints, p. arXiv:1308.0847
  * Lovell et al. (2022) Lovell C. C., Wilkins S. M., Thomas P. A., Schaller M., Baugh C. M., Fabbian G., Bahé Y., 2022, MNRAS, 509, 5046
  * MacCrann et al. (2018) MacCrann N., et al., 2018, MNRAS, 480, 4614
  * Maksimova et al. (2021) Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS,
  * McGibbon & Khochfar (2022) McGibbon R. J., Khochfar S., 2022, MNRAS, 513, 5423
  * Nishimichi et al. (2019) Nishimichi T., et al., 2019, ApJ, 884, 29
  * Nishimichi et al. (prep) Nishimichi et al. in prep., XXX
  * Orsi & Angulo (2018) Orsi Á. A., Angulo R. E., 2018, MNRAS, 475, 2530
  * Peacock & Smith (2000) Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144
  * Perez et al. (2023) Perez L. A., Genel S., Villaescusa-Navarro F., Somerville R. S., Gabrielpillai A., Anglés-Alcázar D., Wandelt B. D., Yung L. Y. A., 2023, ApJ, 954, 11
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Reddick et al. (2014) Reddick R. M., Tinker J. L., Wechsler R. H., Lu Y., 2014, ApJ, 783, 118
  * Scoccimarro et al. (2001) Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, ApJ, 546, 20
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Speagle & Barbary (2018) Speagle J., Barbary K., 2018, dynesty: Dynamic Nested Sampling package, Astrophysics Source Code Library (ascl:1809.013)
  * Spergel et al. (2013) Spergel D., et al., 2013, arXiv e-prints, p. arXiv:1305.5422
  * Takada et al. (2014) Takada M., et al., 2014, PASJ, 66, R1
  * To et al. (2021) To C.-H., et al., 2021, MNRAS, 502, 4093
  * Vale & Ostriker (2006) Vale A., Ostriker J. P., 2006, MNRAS, 371, 1173
  * Wang et al. (2022) Wang K., Mao Y.-Y., Zentner A. R., Guo H., Lange J. U., van den Bosch F. C., Mezini L., 2022, MNRAS, 516, 4003
  * Wechsler & Tinker (2018) Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435
  * Wechsler et al. (2006) Wechsler R. H., Zentner A. R., Bullock J. S., Kravtsov A. V., Allgood B., 2006, ApJ, 652, 71
  * Xu et al. (2021) Xu X., Kumar S., Zehavi I., Contreras S., 2021, MNRAS, 507, 4879
  * Ye et al. (2017) Ye J.-N., Guo H., Zheng Z., Zehavi I., 2017, ApJ, 841, 45
  * Yuan et al. (2020) Yuan S., Eisenstein D. J., Leauthaud A., 2020, MNRAS, 493, 5551
  * Yuan et al. (2021) Yuan S., Garrison L. H., Hadzhiyska B., Bose S., Eisenstein D. J., 2021, MNRAS, 510, 3301
  * Yuan et al. (2022a) Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022a, MNRAS,
  * Yuan et al. (2022b) Yuan S., Hadzhiyska B., Bose S., Eisenstein D. J., 2022b, MNRAS, 512, 5793
  * Yuan et al. (2023a) Yuan S., Zamora A., Abel T., 2023a, MNRAS,
  * Yuan et al. (2023b) Yuan S., Hadzhiyska B., Abel T., 2023b, MNRAS, 520, 6283
  * Zhai et al. (2023) Zhai Z., et al., 2023, ApJ, 948, 99
  * Zheng et al. (2005) Zheng Z., et al., 2005, ApJ, 633, 791
  * Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760

##  Appendix A Covariance matrices

Figure 13 and Figure 14 show the covariance matrices of the RD-$k$NN and
DD-$k$NN data vectors calculated from 1250 jackknife regions on the Beyond-2p
redshift-space mock.

Figure 13: Visualizations of the RD-$k$NN correlation matrix. The covariance
matrix is computed from 1250 jackknife regions. The “bins” represents the
flattened indices of $k$NN bins along $k$, $r_{p}$, and $r_{\pi}$ axes, in
orders of outer to inner cycles. For ease of visualization, we separate the
bins into blocks via the vertical and horizontal white lines, where each block
corresponds to a $k$. Within each $k$ block, the increasing bin number cycles
through $r_{\pi}$ values at each fixed $r_{p}$.  Figure 14: Visualizations of
the DD-$k$NN correlation matrix. The covariance matrix is computed from 1250
jackknife regions. The indices are arranged in the same way as Figure 13. Note
that for DD-$k$NN, the lowest meaningful value is $k=2$, because the first
neighbor is the galaxy itself.

##  Appendix B Emulator performance

In this section, we present metrics indicating the performance of the $k$NN
emulators. Figure 15 and Figure 16 show the correlation matrices corresponding
to the RD-$k$NN and DD-$k$NN emulator errors. Contrary to the sample variance
covariance matrices shown in Figure 13 and Figure 14, the emulator errors
exhibit larger off-diagonal contributions. The mean absolute error of the
emulators is approximately 0.3 and 0.6 that of the sample variance for
RD-$k$NN and DD-$k$NN, respectively. That translates to a $10$-40$\%$ increase
to the diagonal amplitude of the covariance matrices.

Figure 15: The correlation matrix of the RD-$k$NN emulator covariance matrix.
This plot is constructed similarly to Figure 13, where we use a white grid to
denote different $k$s. The bin indices also assume the same order.  Figure 16:
The correlation matrix of the DD-$k$NN emulator covariance matrix.