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).

##### Covariance

We compute jackknife covariances for 2D-$k$NNs using the 10 realizations
provided. We divide each realization into 125 cubic chunks, each of size
($400h^{-1}$Mpc)3. We visualize the resulting correlation matrix in Figure A1
and A2 of Yuan et al. (2024). Both covariances are invertible. The RD-$k$NN
covariance matrix has a condition number of 1$\times 10^{7}$ and the DD-$k$NN
covariance matrix has a condition number of 4$\times 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.

##### Priors

This analysis employs cosmology priors based on the $\Lambda$CDM parameter
range covered by the AbacusSummit simulations,

$\displaystyle\omega_{\mathrm{cdm}}\in[0.099,0.140],\qquad\sigma_{8}\in[0.68,0.94],\qquad
n_{s}\in[0.90,1.02],$ $\displaystyle\omega_{\rm b}\in[0.021,0.024]\,,$ (26)

with the prior distribution specified by the AbacusSummit parameter envelope,
$N_{\rm ur}=2.0328$ and $\alpha_{s}=0$. For the HOD parameters, we adopt the
following flat priors:

$\displaystyle\log_{10}M_{\rm
cut}\sim\mathcal{U}[12.0,14.5],\qquad\log_{10}M_{1}\sim\mathcal{U}[13.0,15.0],\qquad\alpha\sim\mathcal{U}[0.5,1.5],$
$\displaystyle\alpha_{\rm{vel,c}}\sim\mathcal{U}[0.0,1.0],\qquad\alpha_{\rm{vel,s}}\sim\mathcal{U}[0.2,1.8],\qquad\log_{10}\sigma\sim\mathcal{U}[-3.5,1.5],$
$\displaystyle\kappa\sim\mathcal{U}[0.0,2.0],\qquad B_{\rm
cent}\sim\mathcal{U}[-1.0,1.0],\qquad B_{\rm sat}\sim\mathcal{U}[-1.0,1.0]\,.$
(27)

####  5.5.4 Analysis choices

##### Scale 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.

We apply an additional scale cut to the DD-$k$NN statistics in
$r_{p}>5h^{-1}$Mpc. We refer the readers to Yuan et al. (2024) for
justification of this cut. Essentially, the DD-$k$NN statistic is highly
sensitive to small-scale modeling, and we find our models to not be flexible
enough to accurately predict the DD-$k$NN statistic at $r_{p}<5h^{-1}$Mpc
given the precision of a ($2\,h^{-1}$Gpc)3 volume. We do not apply any scale
cut to the RD-$k$NN statistic as it is less sensitive to small-scale modeling.
Unfortunately, our scale cut at $r_{p}<5h^{-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.

##### Consistency 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.

Specifically, we set RD-$k$NN and DD-$k$NN as separate data vectors as they
summarize different information of the density field. RD-$k$NN captures the
density information and DD-$k$NN captures the clustering information. We also
include the standard 2PCF as a third statistics that share some information
with DD-$k$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.

Our 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 $\sigma_{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.

##### Unmasking 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.

The 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.

Finally, 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}>5h^{-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).

The $68\%$ marginalized cosmology constraints with our final results using our
favored model are shown in Figure 1 and Figure 4. The RD-$k$NN fit achieves a
$\chi^{2}$/d.o.f. = 0.44 whereas the DD-$k$NN fit achieves a $\chi^{2}$/d.o.f.
= 1.23.

##### Caveats

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.

Another 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.

Figure 12: (Plot adapted from Figure 9 of Yuan et al. 2024) The constraints of
$k$NN statistics on $\sigma_{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\sigma$ marginalized constraints. The orange triangles
show the maximum likelihood points. The green line represents the truth. The
left, middle, right blocks represent the RD-$k$NN constraints, DD-$k$NN
constraints with scale cut $r_{p}>5h^{-1}$Mpc, and DD-$k$NN constraints
without minimum scale cuts ($r_{\mathrm{min}}=0.67h^{-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 $\sigma_{8}$. Unfortunately, our
models were not flexible enough to self-consistently model $r_{p}<5h^{-1}$Mpc.
This comparison motivates for more robust modeling of small scales.

###  5.6 Density-split clustering303030Authors: Enrique Paillas, Carolina
Cuesta-Lazaro.

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 Cuesta-Lazaro 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 $\Lambda$CDM mocks in redshift space.
Here we briefly summarize the algorithm to measure the DSC data vector:

We begin by painting the galaxy overdensity field to a rectangular grid of
$5\,h^{-1}{\rm Mpc}$ resolution. We smooth the field with a Gaussian kernel of
width $R_{s}=10\,h^{-1}{\rm Mpc}$, and we then sample it at $N_{\rm query}$
random query points, where $N_{\rm 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 $\Lambda$CDM simulations.

We 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 $\mu$,
using 241 $\mu$ bins from -1 to 1, and radial bins of scale-dependent widths:
$1\,h^{-1}{\rm Mpc}$ width for $s\in[0,4]\,h^{-1}{\rm Mpc}$, $3\,h^{-1}{\rm
Mpc}$ width for $s\in[4,30]\,h^{-1}{\rm Mpc}$, and $5\,h^{-1}{\rm Mpc}$ width
for $s\in[30,150]\,h^{-1}{\rm 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.

Finally, the data vector is constructed as a concatenation of the monopole and
quadrupole moments of the quintile-galaxy cross-correlation functions,
$\xi_{0+2}^{\rm QG}(s)$, and the quintile autocorrelation functions ,
$\xi_{0+2}^{\rm QQ}(s)$,

$\widehat{\mathbf{d}}_{\rm{DSC}}=\left(\xi_{0}^{\rm
Q=(0,1,3,4)G}(s),\xi_{2}^{\rm Q=(0,1,3,4)G}(s),\xi_{0}^{\rm
QQ=(00,11,33,44)},\xi_{2}^{\rm QQ=(00,11,33,44)}(s)\right)\,,$ (28)

using the same radial binning as describe above. When doing so, we discard the
middle quintile, ${\rm 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.

Figure 13: Probability distribution function of the galaxy overdensity measured at random query points from one of the $\Lambda$CDM simulations, where the galaxy density field has been smoothed with a Gaussian filter of radius $R_{s}=10\,h^{-1}{\rm Mpc}$. The points are classified into five density quintiles, shown by the different colours. |
---|---
|
Figure 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 cross-correlation 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. Figure 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, ${\bf C}_{\rm train}$, is the same
as the data error in this case, so we do not display it.

####  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.

##### Parametrization and simulations

We start 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 $\omega_{\rm cdm}$, the physical baryon density $\omega_{b}$, the RMS
of density fluctuations $\sigma_{8}$, the spectral index $n_{s}$, the running
of the spectral index $\alpha_{\rm s}$, the effective number of massless
relics $N_{\rm 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-$\Lambda$CDM parameter space, using the priors on
cosmological parameters defined in Eq. (5.6.3).

We 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. (5.6.3). Each cosmology is assigned 1,000
HOD parameters from the Latin hypercube, and dark matter haloes are populated
with HOD galaxies from those parameters, 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 $\Lambda$CDM simulations ($n_{\rm
gal}\approx 4.5\times 10^{-4}\,(h/{\rm Mpc})^{-3}$), we invoke an
incompleteness parameter $f_{\rm 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).

##### Neural-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.

The 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 framework323232Available at Github.com/optuna/optuna 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 repository333333Available at
Github.com/florpi/sunbird.

####  5.6.3 Inference

##### Likelihood

We fit our emulator to the mean measurement of the 10 realizations of the
$\Lambda$CDM redshift-space mocks. We assume a Gaussian data likelihood with
covariance $\mathsf{C}_{\rm{tot}}$, which incorporates the total error budget
of our analysis. This covariance includes contributions from sample variance
in the data vector ($\mathsf{C}_{\rm data}$), sample variance in the
simulations used for training ($\mathsf{C}_{\rm train}$) and the intrinsic
emulator error ($\mathsf{C}_{\rm emu}$):

$\mathsf{C}_{\rm tot}=\mathsf{C}_{\rm data}+\mathsf{C}_{\rm
train}+\mathsf{C}_{\rm emu}\,.$ (29)

##### Covariance

To calculate $\mathsf{C}_{\rm data}$, we measure the DSC multipoles from 1500
realizations of the small AbacusSummit simulations, which are
$500\,h^{-1}\text{Mpc}$ on a side and use the baseline AbacusSummit cosmology.
We choose a combination of HOD parameters that produce multipoles that
minimize the $\chi^{2}$ with respect to the data vector measured from the
$\Lambda$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}\text{Mpc}$ $\Lambda$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 $\Lambda$CDM simulations. To account for the finite volume of the
simulations used for training, we repeat the same procedure to calculate
$\mathsf{C}_{\rm train}$, so effectively $\mathsf{C}_{\rm
train}=\mathsf{C}_{\rm data}$ for this analysis. Finally, to estimate
$\mathsf{C}_{\rm emu}$, which accounts for the error due to an imperfect
emulation of the multipoles, we compare the emulator predictions
$\mathbf{m}(\mathbf{\Omega}^{\rm test})$ against the multipoles measured from
our test set of simulations, $\widehat{\mathbf{d}}^{\rm test}$, and compute
the covariance of the absolute emulator error $\Delta^{{\rm
test}}=\mathbf{m}(\mathbf{\Omega}^{\rm test})-\widehat{\mathbf{d}}^{\rm test}$
as

${\bf C}_{\rm emu}=\frac{1}{n_{\rm test-1}}\sum_{k=1}^{n_{\rm
test}}\left(\Delta^{{\rm test}_{k}}-\overline{\Delta^{{\rm
test}}}\right)\left(\Delta^{{\rm test}_{k}}-\overline{\Delta^{{\rm
test}}}\right)^{\top}\,,$ (30)

where the overline denotes an average across test samples.

Figure 15 shows the correlation matrices associated to $\mathsf{C}_{\rm data}$
and $\mathsf{C}_{\rm 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.

##### Sampler and priors

We sample the posterior distribution using the Dynesty nested sampler
(Speagle, 2020b). We assume the following flat priors on cosmological
parameters:

$\displaystyle\omega_{\rm
cdm}\sim\mathcal{U}[0.1032,0.140],\qquad\sigma_{8}\sim\mathcal{U}[0.678,0.938],\qquad
n_{s}\sim\mathcal{U}[0.9012,1.025],$ $\displaystyle\omega_{\rm
b}\sim\mathcal{U}[0.0207,0.0243]$ (31)

We fix those parameters beyond the base-$\Lambda$CDM model to the following
default values: $w_{0}=-1$, $w_{a}=0$, $N_{\rm ur}=2.0328$, and
$\alpha_{s}=0$. For the HOD parameters, we adopt the following flat priors:

$\displaystyle\log_{10}M_{\rm
cut}\sim\mathcal{U}[12.0,13.5],\qquad\log_{10}M_{1}\sim\mathcal{U}[12.5,15.0],\qquad\alpha\sim\mathcal{U}[0.3,1.5],$
$\displaystyle\alpha_{v,{\rm
cen}}\sim\mathcal{U}[0.0,1.0],\qquad\alpha_{v,{\rm
sat}}\sim\mathcal{U}[0.0,2.0],\qquad\log_{10}\sigma\sim\mathcal{U}[-7.0,0.0],$
$\displaystyle\kappa\sim\mathcal{U}[0.0,1.5],\qquad B_{\rm
cen}\sim\mathcal{U}[-0.5,0.5],\qquad B_{\rm sat}\sim\mathcal{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}\text{Mpc}<s<151\,h^{-1}\text{Mpc}$, including the monopole and
quadrupole of the quintile-galaxy CCF and the quintile ACF, using quintiles
${\rm Q}_{0}$, ${\rm Q}_{1}$, ${\rm Q}_{3}$, and ${\rm 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\sigma$
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.

Figure 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.

##### Scale cut validation

We observe that the recovered constraints do not change significantly when
excluding the very small scale information at $s<30h^{-1}\text{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}\text{Mpc}<s<151\,h^{-1}\text{Mpc}$ seems robust at
recovering unbiased cosmological constraints for the volume of the
$\Lambda$CDM simulations. Based on this, we choose to adopt this scale range
for the baseline analysis.

##### Choice 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, ${\rm Q}_{0}$ and
${\rm 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).

##### Choice of HOD model

The choice of the HOD model has little impact on $\omega_{\rm cdm}$, but
becomes noticeable for $\sigma_{8}$ and $n_{s}$. More specifically, not
including central and satellite velocity bias in the model shifts the mean of
the posterior for $\sigma_{8}$ and $n_{s}$ to smaller values, although well
within the $1\sigma$ 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 $\Lambda$CDM simulations
fitted in this study, we choose to include them in the default analysis.

##### Choice 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 $\omega_{\rm cdm}$ and $\sigma_{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 have 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.

##### Choice 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)).

##### Unmasking 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 $\chi^{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.

##### Caveats

Possible caveats that could lead to unexpected results in our analysis after
unmasking include:

  * -

HOD parameters lying outside the range of our priors.

  * -

An HOD model including effects that were not implemented in our emulator, such
as satellite profile flexibilities to account for baryonic effects.

  * -

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.

  * -

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.

  * -

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 Voids343434Authors: Sofia Contarini, Giovanni Verza, Nico
Hamaus, Alice Pisani.

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 VIDE363636Available at
Bitbucket.org/cosmicvoids/vide_public/wiki/Home (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_{\mathrm{side}}=128$ from the angular
distribution of mock-galaxy 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
$\Lambda$CDM light-cone mock catalog, where VIDE identifies more than $14,000$
voids.

##### Void 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\leq 0.968$,
$0.968<z\leq 1.124$ and $1.124<z\leq 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

$\widehat{\mathbf{d}}_{\rm{VSF}}=\left(\widehat{N}(\tilde{R}_{i=1,\ldots,N_{r}},z_{0}),\widehat{N}(\tilde{R}_{i=1,\ldots,N_{r}},z_{1}),\widehat{N}(\tilde{R}_{i=1,\ldots,N_{r}},z_{2})\right)\,,$
(33)

where $\tilde{R}_{i=1,\ldots,N_{r}}$ denotes the $N_{r}$ bins in cleaned void
size, and $z_{0,1,2}$ the three tomographic bins described above.

##### Void-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_{\parallel}$ and $s_{\perp}$,
measured in units of the effective radius $R$. A stack over a sample of voids
then yields a cross-correlation function between void centers and galaxies
(Hamaus et al., 2015). We use the Landy-Szalay estimator (Landy & Szalay,
1993) to obtain a measurement of the VGCF data vector via

$\widehat{d}_{{\mathrm{VGCF},i}}\equiv\widehat{\xi}_{\mathrm{v}\mathrm{g},i}(s_{\perp},s_{\parallel})\;,$
(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 excursion-set theory. The so-
called Vdn (volume conserving) model is defined as:

$\frac{{\rm d}n}{{\rm d}\ln R}=\frac{f_{\ln\sigma_{\rm R}}(\sigma_{\rm
R})}{V(R)}\,\frac{{\rm d}\ln\sigma_{\rm R}^{-1}}{{\rm d}\ln
R_{\mathrm{L}}}\biggr{\rvert}_{R_{\mathrm{L}}=R_{\mathrm{L}}(R)}\,,$ (35)

where $\sigma_{\rm R}$ represents the root mean square variance of linear
matter perturbations on a scale $R_{\mathrm{L}}$ and where, to compute the
number density of voids, we rely on the multiplicity function:

$\begin{gathered}f_{\ln\sigma_{\rm R}}(\sigma_{\rm
R})=2\sum_{j=1}^{\infty}\,\exp{\bigg{(}-\frac{(j\pi x)^{2}}{2}\bigg{)}}\,j\pi
x^{2}\,\sin{\left(j\pi D\right)}\,,\\\ \text{with}\quad
D\equiv\frac{|\delta_{\mathrm{v}}^{\mathrm{L}}|}{\delta_{\mathrm{c}}^{\mathrm{L}}+|\delta_{\mathrm{v}}^{\mathrm{L}}|}\,\quad\text{and}\quad
x\equiv\frac{D}{|\delta_{\mathrm{v}}^{\mathrm{L}}|}\sigma_{\rm
R}\,,\end{gathered}$ (36)

Considering linear theory, the function $f_{\ln\sigma_{\rm R}}(\sigma_{\rm
R})$ represents the volume fraction of the Universe occupied by cosmic voids,
with radii in the range $(R,R+{\rm d}R)$. The additional cosmological
dependence of this function lies in the density contrasts defining the
formation of dark matter halos and cosmic voids,
$\delta_{\mathrm{c}}^{\mathrm{L}}$ and $\delta_{\mathrm{v}}^{\mathrm{L}}$,
respectively. To predict the number of voids identified from observations of
biased tracers, a number of effects have to be considered:

  1. (i)

the link of $\delta_{\mathrm{v}}^{\mathrm{L}}$ to its nonlinear counterpart
$\delta_{\mathrm{c}}^{\mathrm{NL}}$ (corresponding to the same threshold but
considering nonlinear theory) in the biased tracer field;

  2. (ii)

the dynamic distortions caused by the peculiar velocities of the tracers;

  3. (iii)

the impact of geometric distortions.

Firstly, to account for the link of $\delta_{\mathrm{v}}^{\mathrm{L}}$ 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 CosmoBolognaLib373737Available at
Gitlab.com/federicomarulli/CosmoBolognaLib (Marulli et al., 2016). The output
of this procedure is a catalog of non-overlapping spheres with cleaned void
radius $\tilde{R}$, embedding a fixed density contrast in the tracer density
field, $\delta_{\mathrm{v,tr}}^{\mathrm{NL}}$. We chose to fix the latter to
$\delta_{\mathrm{v,tr}}^{\mathrm{NL}}=-0.7$. We then use a bias relation to
convert the density contrast in the tracer density field,
$\delta_{\mathrm{v,tr}}^{\mathrm{NL}}$, to the corresponding one in the matter
distribution, $\delta_{\mathrm{v}}^{\mathrm{NL}}$. 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):

$\begin{gathered}\delta_{\mathrm{v}}^{\mathrm{NL}}=\frac{\delta_{\mathrm{v,tr}}^{\mathrm{NL}}}{\mathcal{F}(b_{\mathrm{eff}},\sigma_{8})}\,,\text{
with}\\\
\mathcal{F}(b_{\mathrm{eff}}\,\sigma_{8})=C_{\mathrm{slope}}\,b_{\mathrm{eff}}\,\sigma_{8}+C_{\mathrm{offset}}\,.\end{gathered}$
(37)

Secondly, 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
also correct for this enlargement (Contarini et al., 2022), allowing us to
statistically align voids observed in the redshift-space galaxy distribution
to theoretical voids, identified in real space with unbiased tracers. In Eq.
(37), $C_{\rm slope}$ and $C_{\rm offset}$ are redshift-independent
coefficients of the linear function $\mathcal{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 $\Omega_{\rm m}$. We underline, however,
that exploiting the information derived from mock catalogs to limit the priors
of $C_{\rm slope}$ and $C_{\rm offset}$ would allow us to provide cosmological
constraints on $\sigma_{8}$ as well (Contarini et al., 2023).

Thirdly, 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_{\parallel}^{1/3}q_{\perp}^{2/3}R^{*}$, where $q_{\parallel}$ and
$q_{\perp}$ are defined via

$\begin{gathered}r_{\parallel}=\frac{H^{*}(z)}{H(z)}\,r^{*}_{\parallel}\equiv
q_{\parallel}\,r^{*}_{\parallel}\,,\\\ r_{\perp}=\frac{D_{\rm
A}(z)}{D^{*}_{\rm A}(z)}\,r^{*}_{\perp}\equiv
q_{\perp}\,r^{*}_{\perp}\,,\end{gathered}$ (38)

and the asterisk symbol indicates quantities evaluated in the fiducial
cosmology. Here, $r^{*}_{\parallel}$ and $r^{*}_{\perp}$ 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_{\rm A}(z)$ is the comoving angular-diameter
distance.

The final VSF model depends on cosmology through the quantities $\sigma_{\rm
R}$, $q_{\parallel}$ and $q_{\perp}$, i.e. those parameters determining the
amplitude of the matter power spectrum and cosmological distances.
Additionally, the model depends on the nuisance parameters $C_{\rm slope}$ and
$C_{\rm 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: $\Omega_{\rm m},\sigma_{8},h,\Omega_{\rm b}h^{2},n_{\rm s},C_{\rm
slope},C_{\rm offset}$. In addition, we assume one massive neutrino species of
minimal mass. As already mentioned, in this analysis we will provide
constraints on $\Omega_{\rm m}$ only, because no calibration of the parameters
$C_{\rm slope}$ and $C_{\rm 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.

Figure 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 $\tilde{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.

##### VGCF 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
$\mathbf{u}_{\mathrm{v}}$ around the void center follows linear dynamics
according to the linearized continuity equation (Peebles, 1980),

$\mathbf{u}_{\mathrm{v}}(\mathbf{r})=-f(z)\frac{H(z)}{1+z}\,\frac{\mathbf{r}}{r^{3}}\int_{0}^{r}\delta(r^{\prime})\,r^{\prime\,2}\,\mathrm{d}r^{\prime}\;,$
(39)

where $f(z)$ is the linear growth rate of the matter-density contrast
$\delta$. Peculiar velocities cause voids to appear anisotropic in redshift
space, because the Doppler effect contributes an additional component to the
cosmological redshift along the $\mathbf{r}_{\parallel}$ 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 $\xi_{\mathrm{v}\mathrm{g}}(r)=b\delta(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.

With these ingredients we can model the coordinate transformation between the
real-space vector $\mathbf{r}$ and redshift-space vector $\mathbf{s}$ for the
separation of void centers and galaxies,

$\mathbf{s}=\mathbf{r}+\frac{1+z}{H(z)}\,\mathbf{u}_{\parallel}=\mathbf{r}-\frac{f(z)}{b(z)}\,\frac{\mathbf{r}_{\parallel}}{r^{3}}\int_{0}^{r}\xi_{\mathrm{v}\mathrm{g}}(r^{\prime})\,r^{\prime\,2}\;.$
(40)

This requires the VGCF in real space, $\xi_{\mathrm{v}\mathrm{g}}(r)$, which
is not directly observable. We can, however, observe the line-of-sight
projected VGCF,
$\tilde{\xi}_{\mathrm{v}\mathrm{g}}(s_{\perp})=\int\xi_{\mathrm{v}\mathrm{g}}(\mathbf{s})\,\mathrm{d}s_{\parallel}$,
which is insensitive to redshift-space distortions, as it only depends on
separations $s_{\perp}$ 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),

$\xi_{\mathrm{v}\mathrm{g}}(r)=-\frac{1}{\pi}\int_{r}^{\infty}\frac{\mathrm{d}\tilde{\xi}_{\mathrm{v}\mathrm{g}}(s_{\perp})}{\mathrm{d}s_{\perp}}\left(s_{\perp}^{2}-r^{2}\right)^{-1/2}\mathrm{d}s_{\perp}\;.$
(41)

Eqs. (39)–(41) fully specify the redshift-space VGCF at linear order in
perturbation theory and a closed-form expression for
$\xi_{\mathrm{v}\mathrm{g}}(s_{\perp},s_{\parallel})$ 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 $\mathcal{M}$ (for
monopole-like) and $\mathcal{Q}$ (for quadrupole-like), and adopt the semi-
empirical expression derived in Ref. Hamaus et al. (2022),

$\xi_{\mathrm{v}\mathrm{g}}(s_{\perp},s_{\parallel})=\mathcal{M}\left\\{\xi_{\mathrm{v}\mathrm{g}}(r)+\frac{f}{b}\overline{\xi}_{\mathrm{v}\mathrm{g}}(r)+2\mathcal{Q}\,\frac{f}{b}\mu^{2}\left[\xi_{\mathrm{v}\mathrm{g}}(r)-\overline{\xi}_{\mathrm{v}\mathrm{g}}(r)\right]\right\\}\;,$
(42)

where $\mu=r_{\parallel}/r$ and
$\overline{\xi}_{\mathrm{v}\mathrm{g}}(r)=3r^{-3}\\!\int_{0}^{r}\xi_{\mathrm{v}\mathrm{g}}(r^{\prime})\,r^{\prime
2}\,\mathrm{d}r^{\prime}$. 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,

$r_{\perp}=q_{\perp}s_{\perp}\;,\qquad
r_{\parallel}=q_{\parallel}s_{\parallel}\left[1-\frac{1}{3}\frac{f}{b}\mathcal{M}\,\overline{\xi}_{\mathrm{v}\mathrm{g}}(r)\right]^{-1}\\!\;,$
(43)

which 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_{\parallel}^{1/3}q_{\perp}^{2/3}$ with the AP parameters, only ratios of
$q_{\perp}$ and $q_{\parallel}$ appear in Eq. (43). Hence, the final model
space for the VGCF is spanned by four parameters; $f/b$,
$q_{\perp}/q_{\parallel}$, $\mathcal{M}$, and $\mathcal{Q}$. To summarize, our
modeling of the VGCF involves the following steps:

  1. (i)

estimate the real-space VGCF via deprojection of its line-of-sight projected
counterpart using Eq. (41);

  2. (ii)

account for dynamic (redshift-space) distortions along the line of sight using
Eq. (40);

  3. (iii)

account for geometric (AP) distortions using Eq. (43);

  4. (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).

Figure 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).

Figure 18: Void-galaxy cross-correlation function in separations along and
perpendicular to the line of sight, extracted from the entire $\Lambda$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.

####  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.

##### VSF 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):

$\mathcal{L}[\widehat{\mathbf{d}}_{\rm{VSF}}|\mathbf{\Omega}]=\prod_{i}\frac{m_{{\rm{VSF}},i}(\mathbf{\Omega})^{\widehat{d}_{{\rm{VSF}},i}}\,\exp{\left[-m_{{\rm{VSF}},i}(\mathbf{\Omega})\right]}}{\widehat{d}_{{\rm{VSF}},i}!}\,,$
(44)

with model parameters $\mathbf{\Omega}$ listed in Table 4. We use the
functions implemented in the CosmoBolognaLib to sample the likelihood.

##### VSF Priors

We impose wide uniform priors for $\Omega_{\rm m}$ and $\sigma_{8}$, i.e.
$\mathcal{U}[0.05,0.9]$ and $\mathcal{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 parameter-space 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_{\rm slope}$ and $C_{\rm offset}$, assigning uniform priors as
well. For a summary of all our parameters and priors, see Table 4.

VSF

Parameter | Prior
---|---
$\Omega_{\mathrm{m}}$ | $\mathcal{U}[0.05,0.9]$
$\sigma_{8}$ | $\mathcal{U}[0.2,2]$
$h$ | $\mathcal{U}[0.657,0.683]$
$\Omega_{\mathrm{b}}\ h^{2}$ | $\mathcal{U}[0.0216,0.0224]$
$n_{\mathrm{s}}$ | $\mathcal{U}[0.9535,0.9763]$
$C_{\mathrm{slope}}$ | $\mathcal{U}[-50,10]$
$C_{\mathrm{offset}}$ | $\mathcal{U}[-10,50]$

VGCF

Parameter | Prior
---|---
$\Omega_{\mathrm{m}}$ | $\mathcal{U}[0,1]$
$f/b$ | $\mathcal{U}[-10,10]$
$q_{\perp}/q_{\parallel}$ | $\mathcal{U}[-10,10]$
$\mathcal{M}$ | $\mathcal{U}[-10,10]$
$\mathcal{Q}$ | $\mathcal{U}[-10,10]$

Table 4: Parameters and prior bounds for the VSF (left) and VGCF (right)
analyses.

##### VGCF Covariance

We estimate the covariance matrix $\mathsf{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.

##### VGCF 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
$\mathbf{\Omega}=(f/b,q_{\perp}/q_{\parallel},\mathcal{M},\mathcal{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)\simeq\Omega_{\mathrm{m}}^{0.6}(z)$
(Lahav et al., 1991), the other one via the AP parameter
$q_{\perp}/q_{\parallel}\propto D_{\mathrm{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 $\Lambda$CDM cosmology (neglecting the presence of
neutrinos and radiation), the product $D_{\mathrm{A}}(z)H(z)$ is fully
determined by the parameter $\Omega_{\mathrm{m}}$. We obtain a posterior on
the latter by sampling from a Gaussian likelihood containing measurements of
$q_{\perp}/q_{\parallel}$ and their uncertainty within each redshift bin.

##### VGCF Priors

We impose identical uniform priors for each of our model parameters,
$\mathbf{\Omega}\sim\mathcal{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 $\Lambda$CDM light-cone mock, we further impose
$\Omega_{\mathrm{m}}\sim\mathcal{U}[0,1]$.

####  5.7.4 Analysis choices

##### VSF scale cut validation

We select voids above the cleaned radius $\tilde{R}_{\mathrm{min}}=33\ h^{-1}\
\mathrm{Mpc}$, for all the considered redshift bins, i.e. $0.812<z\leq 0.968$,
$0.968<z\leq 1.124$, $1.124<z\leq 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 $\tilde{R}_{\mathrm{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.

Figure 19: Posterior distribution for the value of $\Omega_{\rm m}$ based on
the modeling of the VSF in the $\Lambda$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 $\Omega_{\mathrm{m}}$ value in our analysis
pipeline.

##### VSF 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. $\tilde{R}_{\rm min}$. As
expected, when increasing the value of $\tilde{R}_{\rm 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 $\Omega_{\rm m}$ derived by imposing $\tilde{R}_{\rm
min}=[33,34,35]\ h^{-1}\ \mathrm{Mpc}$. Over a certain scale, that is above
the region where the void counts incompleteness is stronger, the cosmological
constraints on $\Omega_{\rm m}$ are all consistent within $68\%$ of
uncertainty. However, we note a moderate dependency of the best-fit value of
$\Omega_{\rm m}$ with the choice of $\tilde{R}_{\rm min}$, which could impact
the robustness of our results (see caveats below).

We also assessed the impact of assuming different fiducial cosmologies,
varying in particular the value of $\Omega_{\rm 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
$\pm 0.04$ from our fiducial value of $\Omega_{\rm 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 $\Omega_{\rm m}$.

##### VGCF 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 has been tested in simulations for a wide
range of void sizes down to only a few Mpc in effective radius (Schuster et
al., 2023).

However, 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

$R>N_{\mathrm{mgs}}\left(\frac{4\pi}{3}\bar{n}(z)\right)^{-1/3}\;,$ (45)

where $N_{\mathrm{mgs}}$ is a tuning parameter and $\bar{n}(z)$ is the mean
mock-galaxy density as a function of redshift. The choice of
$N_{\mathrm{mgs}}$ is a trade-off between maximizing the void sample size
(lower $N_{\mathrm{mgs}}$) and minimizing the number of spurious voids (higher
$N_{\mathrm{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_{\mathrm{mgs}}\in[1,5]$ and find consistent results with decreasing
uncertainty in our posterior constraints down to a value of
$N_{\mathrm{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_{\mathrm{mgs}}=5$ as
our more conservative default, which results in a minimum void radius of
$R\simeq 50.6h^{-1}\mathrm{Mpc}$ and corresponds to the data vector shown in
Fig. 18.

Furthermore, 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 $\Omega_{\mathrm{m}}$ entering
$D_{\mathrm{A}}(z)H(z)$. However, in $\Lambda$CDM a single redshift bin is
sufficient to determine $\Omega_{\mathrm{m}}$. Hence, for the $\Lambda$CDM
light-cone mock we use the entire redshift range available without binning to
estimate the VGCF with reduced statistical noise.

Figure 20: Posterior distributions for the value of $\Omega_{\mathrm{m}}$
based on the AP test with the VGCF in the $\Lambda$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
$\Omega_{\mathrm{m}}$ in our analysis pipeline.

##### VGCF 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 $\Lambda$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 $\Omega_{\mathrm{m}}$ when
constructing the VGCF data vector and running our analysis pipeline. It is
evident that the posterior distribution moves diametrically to changes in the
fiducial value of $\Omega_{\mathrm{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.

##### Unmasking criteria

We calculate the reduced $\chi^{2}$ statistic for the best-fit model of the
VGCF with parameters $\mathbf{\Omega}^{*}$ as

$\chi^{2}_{\mathrm{red}}=-\frac{2}{N_{\mathrm{dof}}}\log\mathcal{L}(\widehat{\mathbf{d}}_{\rm{VGCF}}|\mathbf{\Omega}^{*})\;,$
(46)

with $N_{\mathrm{dof}}=N_{\mathrm{bin}}-N_{\mathrm{par}}$ degrees of freedom,
where $N_{\mathrm{bin}}$ is the number of bins in the data and
$N_{\mathrm{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.

##### Caveats

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.

  * •

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.

  * •

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.

  * •

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)).

  * •

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.

  * •

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).

  * •

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.

Finally, 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.

##### Post-unmasking studies

After unmasking, we realized that the informative prior on $h,\,n_{\rm
s},\,\omega_{\rm 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 $\Omega_{\rm{m}}$ is robust to both prior
variations, and the wide prior broadens the marginalized $1-\sigma$
uncertainty on $\Omega_{\rm{m}}$ from $12.7\%$ to $15.6\%$.

##### Lessons 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.

In 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.

Figure 21: 1D marginalized constraints on $\Omega_{\rm{m}}$ and $\sigma_{8}$
for analyses of redshift-space mocks, including post-unmasking reanalyses
shown with light opacity.

###  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 ($\Omega_{\rm{m}},\,\sigma_{8}$) are
substantially narrower than the priors in these parameters for all analyses,
several analyses adopted informative priors for some of the other
$\Lambda\mathrm{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_{\rm{s}}$. On the other hand, the VSF analysis chose informative priors on
$h$, $\omega_{\rm{b}}$ and $n_{\rm{s}}$ that excluded the light-cone mock
cosmology at several $\sigma$. Post-unmasking analyses indicate only a minor
bias in $\Omega_{\rm{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).

To first understand the difference between the two EFT pre-unmasking analyses,
note that the EFT FBI pre-unmasking analysis described the matter-galaxy
connection with a second-order galaxy bias expansion and assumed that the
stochastic contribution to the galaxy density field, $\varepsilon$, is
Gaussian distributed with a white power spectrum $\sigma_{\rm Poisson}$,
effectively neglecting the coupling between the stochastic and deterministic
fields of the form $\delta_{g}\supset\delta_{m}\varepsilon$ (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 $\delta_{m}\varepsilon$ term in the bias expansion), and
(3) $\mathcal{O}(k^{2})$ corrections to the stochasticity power spectrum.

To 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_{\mathrm{max}}=\\{0.3,0.15\\}\,h\,\text{Mpc}^{-1}$, and fixing all
cosmological parameters to their true values—except $\sigma_{8}$—and with
$A_{\rm shot}=B_{\rm shot}=b_{\Gamma_{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 $\sigma_{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.383838It is possible that these
terms are less important on large scales, e.g. at $k\sim
0.1\,h\,\text{Mpc}^{-1}$ corresponding to $k_{\mathrm{max}}$ in the EFT FBI
analyses. 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 $\sigma_{8}$ in both analyzes are relatively similar, suggesting that the
difference between these two constraints likely comes from different analysis
assumptions.

Conversely, 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_{\mathrm{max}}=0.1\,h\,\text{Mpc}^{-1}$. The “FBI
extended” posterior in Fig. 22 is completely consistent with the ground truth
$\sigma_{8}$, with a noticeably increased width (in between those from the
“EFT P+B” and “P+B restricted” analyses, pre- and post-unmasking).

The 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.

Figure 22: 1D marginalized constraints on $\sigma_{8}$ for analyses of real-
space mocks, including post-unmasking reanalyses shown with light opacity.

###  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.

Four 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.

The BACCO P analyses with $k_{\rm{max}}=0.2h\,\text{Mpc}^{-1}$ and
$k_{\rm{max}}=0.5h\,\text{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.

The SBI P+B analysis recovers the input cosmology to better than $1{-}\sigma$
with a scale cut $k_{\rm{max}}=0.5h\,\text{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.25h\,\text{Mpc}^{-1}<k\leq 0.5h\,\text{Mpc}^{-1}$ for the power spectrum
(Hahn et al., 2022) and $0.3h\,\text{Mpc}^{-1}<k\leq 0.5h\,\text{Mpc}^{-1}$
for the bispectrum (Hahn et al., 2023).

The density-split clustering analysis includes measurements down to
$1\,h^{-1}\,\mathrm{Mpc}$. However, the gain in constraining power from
$30\,h^{-1}\,\mathrm{Mpc}$ to $1\,h^{-1}\,\mathrm{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}{\rm 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.

The DD-$k$NN baseline analysis employs the scale cut $r_{p}>5h^{-1}\rm{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
density-split 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
post-unmasking 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.

Validation of the uncertainties reported by individual methods will ultimately
require a suite of (ideally, parameter-masked) 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:

  1. 1.

The mock-generating model(s) should be consistent with observations within
some measurement uncertainty.

  2. 2.

The mock-generating model(s) should be predicated upon different assumptions
from the models in the analysis being tested.

For 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 increasing 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 $\Lambda$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
($\Omega_{\mathrm{m}}$, $\sigma_{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:

  * •

Design of pre-unmasking analysis strategies and validation studies. The priors
for the mock catalogs (e.g., $\sigma_{8}\sim\mathcal{U}[0.68,0.9]$ and any HOD
parameters; c.f. Sec. 3) are significantly broader than implicit priors in
(non-masked) 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.

  * •

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 parameter-masked 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.

  * •

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.

While 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 .

From 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.

A 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
anew for different survey volumes and galaxy samples. Furthermore, this
particular mock challenge was clearly labeled 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.

To 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.

The 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 $\Lambda$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.

We 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.

EK, 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. NMN
acknowledges 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).

The 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.

Elisabeth 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 editing.
Yosuke Kobayashi – Challenge: design of mock catalogs, N-body simulations,
creation of mock galaxy catalogs, communication with analysis teams.
Andrés N. Salcedo – Challenge: design of mock catalogs, creation of mock
galaxy catalogs, communication with analysis teams, writing and editing of
general sections.
Mikhail 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.

Tom Abel – $k$NN analysis: design of validation tests, interpretation of
results, post-unmasking studies.
Kazuyuki Akitsu – EFT P+B analysis: analysis runs, interpretation of results,
and writing of analysis section.
Raul Angulo – BACCO P analysis: design of validation tests, interpretation of
results, post-unmasking studies.
Giovanni Cabass – EFT P+B analysis: design of validation tests and
interpretation of results.
Sofia Contarini – Void analysis: development of analysis pipeline, design of
validation tests, analysis runs, interpretation of results, and writing of
analysis section.
Carolina Cuesta-Lazaro – Density-Split analysis: design of validation tests,
analysis runs, interpretation of results, and writing of analysis section.
ChangHoon Hahn – SBI P+B analysis: design of validation tests, analysis runs,
interpretation of results, and writing of analysis section.
Nico Hamaus – Void analysis: development of analysis pipeline, design of
validation tests, analysis runs, interpretation of results, writing of
analysis section, and analysis team coordination.
Donghui Jeong – Challenge: challenge conceptualization and workshop
organization.
Chirag Modi – SBI P+B analysis: design of validation tests, analysis runs,
interpretation of results, and writing of analysis section.
Nhat-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.
Takahiro Nishimishi – Challenge: development of N-body simulation code
GINKAKU, optimization of accuracy parameters for the use in the challenge,
writing of simulation section.
Enrique Paillas – Density-Split analysis: design of validation tests, analysis
runs, interpretation of results, and writing of analysis section.
Marcos 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.
Oliver H. E. Philcox – EFT P+B analysis: design of validation tests, analysis
runs, interpretation of results, and writing of analysis section.
Alice Pisani – Void analysis: design of validation tests, analyses
interpretation, interpretation of results, writing of analysis section, and
analysis team coordination.
Fabian 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.
Satoshi Tanaka – Challenge: development and tuning of N-body simulation code
GINKAKU.
Giovanni Verza – Void analysis: design of validation tests, analysis runs,
interpretation of results, and writing of analysis section.
Sihan Yuan – Challenge: writing of HOD overview. $k$NN analysis: design of
validation tests, analysis runs, interpretation of results, and writing of
analysis section.
Matteo Zennaro – BACCO P analysis: design of validation tests, interpretation
of results, post-unmasking studies.

## References

  * Abidi & Baldauf (2018) Abidi, M. M., & Baldauf, T. 2018, JCAP, 07, 029, doi: 10.1088/1475-7516/2018/07/029
  * Aghamousa et al. (2016) Aghamousa, A., et al. 2016.  https://arxiv.org/abs/1611.00036
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533, doi: 10.1103/PhysRevD.103.083533
  * Alcock & Paczynski (1979) Alcock, C., & Paczynski, B. 1979, Nature, 281, 358, doi: 10.1038/281358a0
  * Alsing et al. (2019) Alsing, J., Charnock, T., Feeney, S., & Wandelt, B. 2019, MNRAS, 488, 4440, doi: 10.1093/mnras/stz1960
  * Anbajagane et al. (2022) Anbajagane, D., Aung, H., Evrard, A. E., et al. 2022, MNRAS, 510, 2980, doi: 10.1093/mnras/stab3587
  * Andrews et al. (2022) Andrews, A., Jasche, J., Lavaux, G., & Schmidt, F. 2022.  https://arxiv.org/abs/2203.08838
  * Angulo & White (2010) Angulo, R. E., & White, S. D. M. 2010, MNRAS, 405, 143, doi: 10.1111/j.1365-2966.2010.16459.x
  * Angulo et al. (2021) Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869, doi: 10.1093/mnras/stab2018
  * Babić et al. (2022) Babić, I., Schmidt, F., & Tucci, B. 2022, JCAP, 08, 007, doi: 10.1088/1475-7516/2022/08/007
  * Balaguera-Antolínez et al. (2023) Balaguera-Antolínez, A., et al. 2023, Astron. Astrophys., 673, A130, doi: 10.5281/zenodo.7330776
  * Baldauf et al. (2015) Baldauf, T., Mirbabayi, M., Simonović, M., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 043514, doi: 10.1103/PhysRevD.92.043514
  * Baldauf et al. (2016) —. 2016.  https://arxiv.org/abs/1602.00674
  * Banerjee & Abel (2021) Banerjee, A., & Abel, T. 2021, MNRAS, 500, 5479, doi: 10.1093/mnras/staa3604
  * Baumann et al. (2012) Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, JCAP, 07, 051, doi: 10.1088/1475-7516/2012/07/051
  * Bayer et al. (2022) Bayer, A. E., Liu, J., Terasawa, R., et al. 2022, arXiv e-prints, arXiv:2210.15647.  https://arxiv.org/abs/2210.15647
  * Bayer et al. (2023) —. 2023, Phys. Rev. D, 108, 043521, doi: 10.1103/PhysRevD.108.043521
  * Bayer et al. (2023) Bayer, A. E., Seljak, U., & Modi, C. 2023, in 40th International Conference on Machine Learning.  https://arxiv.org/abs/2307.09504
  * Bayer et al. (2021) Bayer, A. E., Villaescusa-Navarro, F., Massara, E., et al. 2021, ApJ, 919, 24, doi: 10.3847/1538-4357/ac0e91
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, Astrophys. J., 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Beltz-Mohrmann et al. (2023) Beltz-Mohrmann, G. D., Szewciw, A. O., Berlind, A. A., & Sinha, M. 2023, ApJ, 948, 100, doi: 10.3847/1538-4357/acc576
  * Berlind & Weinberg (2002) Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587, doi: 10.1086/341469
  * Blas et al. (2016a) Blas, D., Garny, M., Ivanov, M. M., & Sibiryakov, S. 2016a, JCAP, 07, 052, doi: 10.1088/1475-7516/2016/07/052
  * Blas et al. (2016b) —. 2016b, JCAP, 07, 028, doi: 10.1088/1475-7516/2016/07/028
  * Bouchet et al. (1995) Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, Astron. Astrophys., 296, 575.  https://arxiv.org/abs/astro-ph/9406013
  * Brieden et al. (2021) Brieden, S., Gil-Marín, H., & Verde, L. 2021, JCAP, 12, 054, doi: 10.1088/1475-7516/2021/12/054
  * Brieden et al. (2022) —. 2022, JCAP, 06, 005, doi: 10.1088/1475-7516/2022/06/005
  * Brinckmann & Lesgourgues (2019) Brinckmann, T., & Lesgourgues, J. 2019, Phys. Dark Univ., 24, 100260, doi: 10.1016/j.dark.2018.100260
  * Brooks & Gelman (1998) Brooks, S. P., & Gelman, A. 1998, Journal of Computational and Graphical Statistics, 7, 434, doi: 10.1080/10618600.1998.10474787
  * Buchert (1992) Buchert, T. 1992, Mon. Not. Roy. Astron. Soc., 254, 729
  * Bullock et al. (2002) Bullock, J. S., Wechsler, R. H., & Somerville, R. S. 2002, MNRAS, 329, 246, doi: 10.1046/j.1365-8711.2002.04959.x
  * Cabass (2021) Cabass, G. 2021, JCAP, 01, 067, doi: 10.1088/1475-7516/2021/01/067
  * Cabass & Schmidt (2020) Cabass, G., & Schmidt, F. 2020, JCAP, 04, 042, doi: 10.1088/1475-7516/2020/04/042
  * Cabass & Schmidt (2020) Cabass, G., & Schmidt, F. 2020, J. Cosmology Astropart. Phys, 2020, 051, doi: 10.1088/1475-7516/2020/07/051
  * Cabass et al. (2024) Cabass, G., Simonović, M., & Zaldarriaga, M. 2024, Phys. Rev. D, 109, 043526, doi: 10.1103/PhysRevD.109.043526
  * Cai et al. (2016) Cai, Y.-C., Taylor, A., Peacock, J. A., & Padilla, N. 2016, MNRAS, 462, 2465, doi: 10.1093/mnras/stw1809
  * Cannon et al. (2022) Cannon, P., Ward, D., & Schmon, S. M. 2022, arXiv e-prints, arXiv:2209.01845, doi: 10.48550/arXiv.2209.01845
  * Carrasco et al. (2012) Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, JHEP, 09, 082, doi: 10.1007/JHEP09(2012)082
  * Chan et al. (2014) Chan, K. C., Hamaus, N., & Desjacques, V. 2014, Phys. Rev. D, 90, 103521, doi: 10.1103/PhysRevD.90.103521
  * Chen et al. (2021) Chen, S.-F., Vlah, Z., Castorina, E., & White, M. 2021, JCAP, 03, 100, doi: 10.1088/1475-7516/2021/03/100
  * Chen et al. (2020) Chen, S.-F., Vlah, Z., & White, M. 2020, JCAP, 07, 062, doi: 10.1088/1475-7516/2020/07/062
  * Chen et al. (2022) —. 2022, JCAP, 02, 008, doi: 10.1088/1475-7516/2022/02/008
  * Cheng et al. (2024) Cheng, S., Marques, G. A., Grandón, D., et al. 2024, arXiv e-prints, arXiv:2404.16085, doi: 10.48550/arXiv.2404.16085
  * Chudaykin et al. (2021a) Chudaykin, A., Dolgikh, K., & Ivanov, M. M. 2021a, Phys. Rev. D, 103, 023507, doi: 10.1103/PhysRevD.103.023507
  * Chudaykin & Ivanov (2019) Chudaykin, A., & Ivanov, M. M. 2019, JCAP, 11, 034, doi: 10.1088/1475-7516/2019/11/034
  * Chudaykin & Ivanov (2023) Chudaykin, A., & Ivanov, M. M. 2023, Phys. Rev. D, 107, 043518, doi: 10.1103/PhysRevD.107.043518
  * Chudaykin et al. (2020) Chudaykin, A., Ivanov, M. M., Philcox, O. H. E., & Simonović, M. 2020, Phys. Rev. D, 102, 063533, doi: 10.1103/PhysRevD.102.063533
  * Chudaykin et al. (2021b) Chudaykin, A., Ivanov, M. M., & Simonović, M. 2021b, Phys. Rev. D, 103, 043525, doi: 10.1103/PhysRevD.103.043525
  * Contarini et al. (2021) Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021, doi: 10.1093/mnras/stab1112
  * Contarini et al. (2023) Contarini, S., Pisani, A., Hamaus, N., et al. 2023, ApJ, 953, 46, doi: 10.3847/1538-4357/acde54
  * Contarini et al. (2024) —. 2024, A&A, 682, A20, doi: 10.1051/0004-6361/202347572
  * Contarini et al. (2019) Contarini, S., Ronconi, T., Marulli, F., et al. 2019, MNRAS, 488, 3526, doi: 10.1093/mnras/stz1989
  * Contarini et al. (2022) Contarini, S., Verza, G., Pisani, A., et al. 2022, A&A, 667, A162, doi: 10.1051/0004-6361/202244095
  * Contreras et al. (2021a) Contreras, S., Angulo, R. E., & Zennaro, M. 2021a, MNRAS, 504, 5205, doi: 10.1093/mnras/stab1170
  * Contreras et al. (2021b) —. 2021b, MNRAS, 508, 175, doi: 10.1093/mnras/stab2560
  * Contreras et al. (2023a) Contreras, S., Chaves-Montero, J., & Angulo, R. E. 2023a, MNRAS, 525, 3149, doi: 10.1093/mnras/stad2434
  * Contreras et al. (2023b) Contreras, S., Angulo, R. E., Springel, V., et al. 2023b, MNRAS, 524, 2489, doi: 10.1093/mnras/stac3699
  * Correa et al. (2022) Correa, C. M., Paz, D. J., Padilla, N. D., et al. 2022, MNRAS, 509, 1871, doi: 10.1093/mnras/stab3070
  * Correa et al. (2021) Correa, C. M., Paz, D. J., Sánchez, A. G., et al. 2021, MNRAS, 500, 911, doi: 10.1093/mnras/staa3252
  * Cousinou et al. (2019) Cousinou, M. C., Pisani, A., Tilquin, A., et al. 2019, Astronomy and Computing, 27, 53, doi: 10.1016/j.ascom.2019.03.001
  * Crain & van de Voort (2023) Crain, R. A., & van de Voort, F. 2023, ARA&A, 61, 473, doi: 10.1146/annurev-astro-041923-043618
  * Cranmer et al. (2020) Cranmer, K., Brehmer, J., & Louppe, G. 2020, Proceedings of the National Academy of Science, 117, 30055, doi: 10.1073/pnas.1912789117
  * Crocce et al. (2006) Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, Mon. Not. Roy. Astron. Soc., 373, 369, doi: 10.1111/j.1365-2966.2006.11040.x
  * Crocce & Scoccimarro (2008) Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533, doi: 10.1103/PhysRevD.77.023533
  * Cuesta-Lazaro & Mishra-Sharma (2023) Cuesta-Lazaro, C., & Mishra-Sharma, S. 2023.  https://arxiv.org/abs/2311.17141
  * Cuesta-Lazaro et al. (2023) Cuesta-Lazaro, C., Paillas, E., Yuan, S., et al. 2023.  https://arxiv.org/abs/2309.16539
  * D’Amico et al. (2022a) D’Amico, G., Donath, Y., Lewandowski, M., Senatore, L., & Zhang, P. 2022a.  https://arxiv.org/abs/2206.08327
  * D’Amico et al. (2022b) —. 2022b.  https://arxiv.org/abs/2211.17130
  * D’Amico et al. (2020) D’Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 05, 005, doi: 10.1088/1475-7516/2020/05/005
  * D’Amico et al. (2024) D’Amico, G., Senatore, L., Zhang, P., & Nishimichi, T. 2024, J. Cosmology Astropart. Phys, 2024, 037, doi: 10.1088/1475-7516/2024/01/037
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * DeRose et al. (2019) DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2019, arXiv e-prints, arXiv:1901.02401, doi: 10.48550/arXiv.1901.02401
  * DeRose et al. (2023) DeRose, J., Kokron, N., Banerjee, A., et al. 2023, J. Cosmology Astropart. Phys, 2023, 054, doi: 10.1088/1475-7516/2023/07/054
  * DESI Collaboration et al. (2024) DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2024, arXiv e-prints, arXiv:2404.03002, doi: 10.48550/arXiv.2404.03002
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1, doi: 10.1016/j.physrep.2017.12.002
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rept., 733, 1, doi: 10.1016/j.physrep.2017.12.002
  * Doeser et al. (2023) Doeser, L., Jamieson, D., Stopyra, S., et al. 2023.  https://arxiv.org/abs/2312.09271
  * Dutton & Macciò (2014) Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359, doi: 10.1093/mnras/stu742
  * Elsner et al. (2020) Elsner, F., Schmidt, F., Jasche, J., Lavaux, G., & Nguyen, N.-M. 2020, JCAP, 01, 029, doi: 10.1088/1475-7516/2020/01/029
  * Euclid Collaboration et al. (2023) Euclid Collaboration, Pezzotta, A., Moretti, C., et al. 2023, arXiv e-prints, arXiv:2312.00679, doi: 10.48550/arXiv.2312.00679
  * Fang et al. (2019) Fang, Y., Hamaus, N., Jain, B., et al. 2019, MNRAS, 490, 3573, doi: 10.1093/mnras/stz2805
  * Feldman et al. (2001) Feldman, H. A., Frieman, J. A., Fry, J. N., & Scoccimarro, R. 2001, Phys. Rev. Lett., 86, 1434, doi: 10.1103/PhysRevLett.86.1434
  * Feroz & Hobson (2008) Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449, doi: 10.1111/j.1365-2966.2007.12353.x
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10, doi: 10.21105/astro.1306.2144
  * Foreman-Mackey et al. (2019) Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, The Journal of Open Source Software, 4, 1864, doi: 10.21105/joss.01864
  * Garrison et al. (2019) Garrison, L. H., Eisenstein, D. J., & Pinto, P. A. 2019, MNRAS, 485, 3370, doi: 10.1093/mnras/stz634
  * Gatti et al. (2022) Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509, doi: 10.1103/PhysRevD.106.083509
  * Gaztanaga et al. (2009) Gaztanaga, E., Cabre, A., Castander, F., Crocce, M., & Fosalba, P. 2009, Mon. Not. Roy. Astron. Soc., 399, 801, doi: 10.1111/j.1365-2966.2009.15313.x
  * Gil-Marín et al. (2017) Gil-Marín, H., Percival, W. J., Verde, L., et al. 2017, Mon. Not. Roy. Astron. Soc., 465, 1757, doi: 10.1093/mnras/stw2679
  * Gonzalez-Perez et al. (2018) Gonzalez-Perez, V., Comparat, J., Norberg, P., et al. 2018, MNRAS, 474, 4024, doi: 10.1093/mnras/stx2807
  * Guo et al. (2015) Guo, H., Zheng, Z., Zehavi, I., et al. 2015, MNRAS, 446, 578, doi: 10.1093/mnras/stu2120
  * Guo et al. (2016) Guo, H., Zheng, Z., Behroozi, P. S., et al. 2016, MNRAS, 459, 3040, doi: 10.1093/mnras/stw845
  * Hadzhiyska et al. (2021) Hadzhiyska, B., Eisenstein, D., Bose, S., Garrison, L. H., & Maksimova, N. 2021, arXiv e-prints, arXiv:2110.11408.  https://arxiv.org/abs/2110.11408
  * Hahn (2020) Hahn, C. 2020, pySpectrum: Power spectrum and bispectrum calculator, Astrophysics Source Code Library, record ascl:2009.014.  http://ascl.net/2009.014
  * Hahn et al. (2019) Hahn, C., Beutler, F., Sinha, M., et al. 2019, Mon. Not. Roy. Astron. Soc., 485, 2956, doi: 10.1093/mnras/stz558
  * Hahn et al. (2022) Hahn, C., Eickenberg, M., Ho, S., et al. 2022, arXiv e-prints, arXiv:2211.00723, doi: 10.48550/arXiv.2211.00723
  * Hahn et al. (2023) Hahn, C., Eickenberg, M., Ho, S., et al. 2023.  https://arxiv.org/abs/2310.15243
  * Hahn et al. (2023) Hahn, C., Eickenberg, M., Ho, S., et al. 2023, J. Cosmology Astropart. Phys, 2023, 010, doi: 10.1088/1475-7516/2023/04/010
  * Hamaus et al. (2017) Hamaus, N., Cousinou, M.-C., Pisani, A., et al. 2017, J. Cosmology Astropart. Phys, 7, 014, doi: 10.1088/1475-7516/2017/07/014
  * Hamaus et al. (2020) Hamaus, N., Pisani, A., Choi, J.-A., et al. 2020, J. Cosmology Astropart. Phys, 2020, 023, doi: 10.1088/1475-7516/2020/12/023
  * Hamaus et al. (2016) Hamaus, N., Pisani, A., Sutter, P. M., et al. 2016, Physical Review Letters, 117, 091302, doi: 10.1103/PhysRevLett.117.091302
  * Hamaus et al. (2015) Hamaus, N., Sutter, P. M., Lavaux, G., & Wandelt, B. D. 2015, J. Cosmology Astropart. Phys, 11, 036, doi: 10.1088/1475-7516/2015/11/036
  * Hamaus et al. (2014a) Hamaus, N., Sutter, P. M., & Wandelt, B. D. 2014a, Physical Review Letters, 112, 251302, doi: 10.1103/PhysRevLett.112.251302
  * Hamaus et al. (2014b) Hamaus, N., Wandelt, B. D., Sutter, P. M., Lavaux, G., & Warren, M. S. 2014b, Physical Review Letters, 112, 041304, doi: 10.1103/PhysRevLett.112.041304
  * Hamaus et al. (2022) Hamaus, N., Aubert, M., Pisani, A., et al. 2022, A&A, 658, A20, doi: 10.1051/0004-6361/202142073
  * Hand et al. (2018) Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160, doi: 10.3847/1538-3881/aadae0
  * Harnois-Déraps et al. (2021) Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623, doi: 10.1093/mnras/stab1623
  * Hartlap et al. (2007) Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399, doi: 10.1051/0004-6361:20066170
  * Hawken et al. (2017) Hawken, A. J., Granett, B. R., Iovino, A., et al. 2017, A&A, 607, A54, doi: 10.1051/0004-6361/201629678
  * Hearin et al. (2023) Hearin, A. P., Chaves-Montero, J., Alarcon, A., Becker, M. R., & Benson, A. 2023, MNRAS, 521, 1741, doi: 10.1093/mnras/stad456
  * Hearin et al. (2021) Hearin, A. P., Chaves-Montero, J., Becker, M. R., & Alarcon, A. 2021, The Open Journal of Astrophysics, 4, 7, doi: 10.21105/astro.2105.05859
  * Hearin et al. (2016) Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016, MNRAS, 460, 2552, doi: 10.1093/mnras/stw840
  * Heitmann et al. (2019) Heitmann, K., et al. 2019, Astrophys. J. Suppl., 245, 16, doi: 10.3847/1538-4365/ab4da1
  * Heydenreich et al. (2022) Heydenreich, S., Brück, B., Burger, P., et al. 2022, A&A, 667, A125, doi: 10.1051/0004-6361/202243868
  * Ho et al. (2024) Ho, M., et al. 2024.  https://arxiv.org/abs/2402.05137
  * Hou et al. (2023) Hou, J., Moradinezhad Dizgah, A., Hahn, C., & Massara, E. 2023, J. Cosmology Astropart. Phys, 2023, 045, doi: 10.1088/1475-7516/2023/03/045
  * Hou et al. (2024) Hou, J., Moradinezhad Dizgah, A., Hahn, C., et al. 2024, arXiv e-prints, arXiv:2401.15074, doi: 10.48550/arXiv.2401.15074
  * Ishiyama et al. (2009) Ishiyama, T., Fukushige, T., & Makino, J. 2009, PASJ, 61, 1319, doi: 10.1093/pasj/61.6.1319
  * Ishiyama et al. (2012) Ishiyama, T., Nitadori, K., & Makino, J. 2012, arXiv e-prints, arXiv:1211.4406.  https://arxiv.org/abs/1211.4406
  * Ishiyama et al. (2021) Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210, doi: 10.1093/mnras/stab1755
  * Ivanov (2021) Ivanov, M. M. 2021, Phys. Rev. D, 104, 103514, doi: 10.1103/PhysRevD.104.103514
  * Ivanov (2023) —. 2023, Effective Field Theory for Large-Scale Structure, doi: 10.1007/978-981-19-3079-9_5-1
  * Ivanov et al. (2024) Ivanov, M. M., Cuesta-Lazaro, C., Mishra-Sharma, S., Obuljen, A., & Toomey, M. W. 2024.  https://arxiv.org/abs/2402.13310
  * Ivanov et al. (2023) Ivanov, M. M., Philcox, O. H. E., Cabass, G., et al. 2023, Phys. Rev. D, 107, 083515, doi: 10.1103/PhysRevD.107.083515
  * Ivanov et al. (2022a) Ivanov, M. M., Philcox, O. H. E., Nishimichi, T., et al. 2022a, Phys. Rev. D, 105, 063512, doi: 10.1103/PhysRevD.105.063512
  * Ivanov et al. (2022b) Ivanov, M. M., Philcox, O. H. E., Simonović, M., et al. 2022b, Phys. Rev. D, 105, 043531, doi: 10.1103/PhysRevD.105.043531
  * Ivanov & Sibiryakov (2018) Ivanov, M. M., & Sibiryakov, S. 2018, JCAP, 07, 053, doi: 10.1088/1475-7516/2018/07/053
  * Ivanov et al. (2020a) Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020a, JCAP, 05, 042, doi: 10.1088/1475-7516/2020/05/042
  * Ivanov et al. (2020b) —. 2020b, Phys. Rev. D, 101, 083504, doi: 10.1103/PhysRevD.101.083504
  * Iwasawa et al. (2016) Iwasawa, M., Tanikawa, A., Hosono, N., et al. 2016, PASJ, 68, 54, doi: 10.1093/pasj/psw053
  * Jeffrey et al. (2021) Jeffrey, N., Alsing, J., & Lanusse, F. 2021, MNRAS, 501, 954, doi: 10.1093/mnras/staa3594
  * Jennings et al. (2013) Jennings, E., Li, Y., & Hu, W. 2013, MNRAS, 434, 2167, doi: 10.1093/mnras/stt1169
  * Kaiser (1987) Kaiser, N. 1987, MNRAS, 227, 1, doi: 10.1093/mnras/227.1.1
  * Kitaura et al. (2016) Kitaura, F.-S., Rodríguez-Torres, S., Chuang, C.-H., et al. 2016, MNRAS, 456, 4156, doi: 10.1093/mnras/stv2826
  * Kobayashi et al. (2022) Kobayashi, Y., Nishimichi, T., Takada, M., & Miyatake, H. 2022, Phys. Rev. D, 105, 083517, doi: 10.1103/PhysRevD.105.083517
  * Kobayashi et al. (2020) Kobayashi, Y., Nishimichi, T., Takada, M., Takahashi, R., & Osato, K. 2020, Phys. Rev. D, 102, 063504, doi: 10.1103/PhysRevD.102.063504
  * Kostić et al. (2023) Kostić, A., Nguyen, N.-M., Schmidt, F., & Reinecke, M. 2023, J. Cosmology Astropart. Phys, 2023, 063, doi: 10.1088/1475-7516/2023/07/063
  * Kreisch et al. (2022) Kreisch, C. D., Pisani, A., Villaescusa-Navarro, F., et al. 2022, ApJ, 935, 100, doi: 10.3847/1538-4357/ac7d4b
  * Kwan et al. (2023) Kwan, J., Saito, S., Leauthaud, A., et al. 2023, ApJ, 952, 80, doi: 10.3847/1538-4357/acd92f
  * Kwon & Hahn (2024) Kwon, K. J., & Hahn, C. 2024, arXiv e-prints, arXiv:2401.12318, doi: 10.48550/arXiv.2401.12318
  * Kwon et al. (2023) Kwon, K. J., Hahn, C., & Alsing, J. 2023, ApJS, 265, 23, doi: 10.3847/1538-4365/acba14
  * Lahav et al. (1991) Lahav, O., Lilje, P. B., Primack, J. R., & Rees, M. J. 1991, MNRAS, 251, 128, doi: 10.1093/mnras/251.1.128
  * Landy & Szalay (1993) Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64, doi: 10.1086/172900
  * Lange et al. (2022) Lange, J. U., Hearin, A. P., Leauthaud, A., et al. 2022, MNRAS, 509, 1779, doi: 10.1093/mnras/stab3111
  * Lange et al. (2023) Lange, J. U., Hearin, A. P., Leauthaud, A., et al. 2023, doi: 10.1093/mnras/stad473
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193, doi: 10.48550/arXiv.1110.3193
  * Lavaux & Wandelt (2012) Lavaux, G., & Wandelt, B. D. 2012, ApJ, 754, 109, doi: 10.1088/0004-637X/754/2/109
  * Lazeyras & Schmidt (2018) Lazeyras, T., & Schmidt, F. 2018, JCAP, 09, 008, doi: 10.1088/1475-7516/2018/09/008
  * Levi et al. (2013) Levi, M., Bebek, C., Beers, T., et al. 2013, arXiv e-prints, arXiv:1308.0847.  https://arxiv.org/abs/1308.0847
  * Lewis (2019) Lewis, A. 2019.  https://arxiv.org/abs/1910.13970
  * Maion et al. (2023) Maion, F., Angulo, R. E., Bakx, T., et al. 2023, arXiv e-prints, arXiv:2307.13754, doi: 10.48550/arXiv.2307.13754
  * Maksimova et al. (2021) Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, Mon. Not. Roy. Astron. Soc., 508, 4017, doi: 10.1093/mnras/stab2484
  * Maksimova et al. (2021) Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, MNRAS, doi: 10.1093/mnras/stab2484
  * Marinoni et al. (2008) Marinoni, C., Guzzo, L., Cappi, A., et al. 2008, Astron. Astrophys., 487, 7, doi: 10.1051/0004-6361:20078891
  * Marulli et al. (2016) Marulli, F., Veropalumbo, A., & Moresco, M. 2016, Astronomy and Computing, 14, 35, doi: 10.1016/j.ascom.2016.01.005
  * Massara et al. (2023) Massara, E., Villaescusa-Navarro, F., Hahn, C., et al. 2023, ApJ, 951, 70, doi: 10.3847/1538-4357/acd44d
  * Matsubara (2008) Matsubara, T. 2008, Phys. Rev. D, 78, 083519, doi: 10.1103/PhysRevD.78.083519
  * Matsubara (2015) —. 2015, Phys. Rev. D, 92, 023534, doi: 10.1103/PhysRevD.92.023534
  * McEwen & Weinberg (2018) McEwen, J. E., & Weinberg, D. H. 2018, MNRAS, 477, 4348, doi: 10.1093/mnras/sty882
  * Modi et al. (2020) Modi, C., Chen, S.-F., & White, M. 2020, Monthly Notices of the Royal Astronomical Society, 492, 5754, doi: 10.1093/mnras/staa251
  * Modi et al. (2023) Modi, C., Pandey, S., Ho, M., et al. 2023, arXiv e-prints, arXiv:2309.15071, doi: 10.48550/arXiv.2309.15071
  * Modi & Philcox (2023) Modi, C., & Philcox, O. H. E. 2023, arXiv e-prints, arXiv:2309.10270, doi: 10.48550/arXiv.2309.10270
  * Moresco et al. (2022) Moresco, M., Amati, L., Amendola, L., et al. 2022, Living Reviews in Relativity, 25, 6, doi: 10.1007/s41114-022-00040-z
  * Namekata et al. (2018) Namekata, D., Iwasawa, M., Nitadori, K., et al. 2018, PASJ, 70, 70, doi: 10.1093/pasj/psy062
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, Astrophys. J., 490, 493, doi: 10.1086/304888
  * Neyrinck (2008) Neyrinck, M. C. 2008, MNRAS, 386, 2101, doi: 10.1111/j.1365-2966.2008.13180.x
  * Nguyen et al. (2021) Nguyen, N.-M., Schmidt, F., Lavaux, G., & Jasche, J. 2021, JCAP, 03, 058, doi: 10.1088/1475-7516/2021/03/058
  * Nguyen et al. (2024) Nguyen, N.-M., Schmidt, F., Tucci, B., Reinecke, M., & Kostić, A. 2024.  https://arxiv.org/abs/2403.03220
  * Nicola et al. (2024) Nicola, A., Hadzhiyska, B., Findlay, N., et al. 2024, J. Cosmology Astropart. Phys, 2024, 015, doi: 10.1088/1475-7516/2024/02/015
  * Nishimichi et al. (2020) Nishimichi, T., D’Amico, G., Ivanov, M. M., et al. 2020, Phys. Rev. D, 102, 123541, doi: 10.1103/PhysRevD.102.123541
  * Nishimichi et al. (2019) Nishimichi, T., et al. 2019, Astrophys. J., 884, 29, doi: 10.3847/1538-4357/ab3719
  * Nitadori et al. (2006) Nitadori, K., Makino, J., & Hut, P. 2006, New Astron., 12, 169, doi: 10.1016/j.newast.2006.07.007
  * Nunes et al. (2022) Nunes, R. C., Vagnozzi, S., Kumar, S., Di Valentino, E., & Mena, O. 2022, Phys. Rev. D, 105, 123506, doi: 10.1103/PhysRevD.105.123506
  * Paillas et al. (2021) Paillas, E., Cai, Y.-C., Padilla, N., & Sánchez, A. G. 2021, MNRAS, 505, 5731, doi: 10.1093/mnras/stab1654
  * Paillas et al. (2023) Paillas, E., et al. 2023.  https://arxiv.org/abs/2309.16541
  * Paillas et al. (2023) Paillas, E., Cuesta-Lazaro, C., Zarrouk, P., et al. 2023, MNRAS, 522, 606, doi: 10.1093/mnras/stad1017
  * Peebles (1980) Peebles, P. J. E. 1980, The large-scale structure of the universe (Princeton University Press, Princeton U.S.A.)
  * Pellejero-Ibañez et al. (2024) Pellejero-Ibañez, M., Angulo, R. E., Jamieson, D., & Li, Y. 2024, MNRAS, doi: 10.1093/mnras/stae489
  * Pellejero Ibañez et al. (2023) Pellejero Ibañez, M., Angulo, R. E., Zennaro, M., et al. 2023, MNRAS, 520, 3725, doi: 10.1093/mnras/stad368
  * Pellejero Ibañez et al. (2022) Pellejero Ibañez, M., Stücker, J., Angulo, R. E., et al. 2022, MNRAS, 514, 3993, doi: 10.1093/mnras/stac1602
  * Pellejero-Ibanez et al. (2022) Pellejero-Ibanez, M., Angulo, R. E., Zennaro, M., et al. 2022, arXiv e-prints, arXiv:2207.06437, doi: 10.48550/arXiv.2207.06437
  * Pelliciari et al. (2022) Pelliciari, D., Contarini, S., Marulli, F., et al. 2022, arXiv e-prints, arXiv:2210.07248, doi: 10.48550/arXiv.2210.07248
  * Petri et al. (2015) Petri, A., Liu, J., Haiman, Z., et al. 2015, Phys. Rev. D, 91, 103511, doi: 10.1103/PhysRevD.91.103511
  * Philcox (2021a) Philcox, O. H. E. 2021a, Phys. Rev. D, 103, 103504, doi: 10.1103/PhysRevD.103.103504
  * Philcox (2021b) —. 2021b, Phys. Rev. D, 104, 123529, doi: 10.1103/PhysRevD.104.123529
  * Philcox & Ivanov (2022) Philcox, O. H. E., & Ivanov, M. M. 2022, Phys. Rev. D, 105, 043517, doi: 10.1103/PhysRevD.105.043517
  * Philcox et al. (2022) Philcox, O. H. E., Ivanov, M. M., Cabass, G., et al. 2022, Phys. Rev. D, 106, 043530, doi: 10.1103/PhysRevD.106.043530
  * Philcox et al. (2021a) Philcox, O. H. E., Ivanov, M. M., Zaldarriaga, M., Simonovic, M., & Schmittfull, M. 2021a, Phys. Rev. D, 103, 043508, doi: 10.1103/PhysRevD.103.043508
  * Philcox et al. (2021b) Philcox, O. H. E., Sherwin, B. D., Farren, G. S., & Baxter, E. J. 2021b, Phys. Rev. D, 103, 023538, doi: 10.1103/PhysRevD.103.023538
  * Pisani et al. (2014) Pisani, A., Lavaux, G., Sutter, P. M., & Wandelt, B. D. 2014, MNRAS, 443, 3238, doi: 10.1093/mnras/stu1399
  * Pisani et al. (2015a) Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015a, Phys. Rev. D, 92, 083531, doi: 10.1103/PhysRevD.92.083531
  * Pisani et al. (2015b) Pisani, A., Sutter, P. M., & Wandelt, B. D. 2015b, ArXiv e-prints.  https://arxiv.org/abs/1506.07982
  * Pisani et al. (2019) Pisani, A., Massara, E., Spergel, D. N., et al. 2019, BAAS, 51, 40.  https://arxiv.org/abs/1903.05161
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Platen et al. (2007) Platen, E., van de Weygaert, R., & Jones, B. J. T. 2007, MNRAS, 380, 551, doi: 10.1111/j.1365-2966.2007.12125.x
  * Pollina et al. (2017) Pollina, G., Hamaus, N., Dolag, K., et al. 2017, MNRAS, 469, 787, doi: 10.1093/mnras/stx785
  * Pollina et al. (2019) Pollina, G., Hamaus, N., Paech, K., et al. 2019, MNRAS, 487, 2836, doi: 10.1093/mnras/stz1470
  * Press et al. (1986) Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Numerical recipes. The art of scientific computing
  * Ramanah et al. (2019) Ramanah, D. K., Lavaux, G., Jasche, J., & Wandelt, B. D. 2019, Astron. Astrophys., 621, A69, doi: 10.1051/0004-6361/201834117
  * Reddick et al. (2014) Reddick, R. M., Tinker, J. L., Wechsler, R. H., & Lu, Y. 2014, ApJ, 783, 118, doi: 10.1088/0004-637X/783/2/118
  * Régaldo-Saint Blancard et al. (2023) Régaldo-Saint Blancard, B., Hahn, C., Ho, S., et al. 2023, arXiv e-prints, arXiv:2310.15250, doi: 10.48550/arXiv.2310.15250
  * Reid et al. (2014) Reid, B. A., Seo, H.-J., Leauthaud, A., Tinker, J. L., & White, M. 2014, MNRAS, 444, 476, doi: 10.1093/mnras/stu1391
  * Rimes & Hamilton (2005) Rimes, C. D., & Hamilton, A. J. S. 2005, Mon. Not. Roy. Astron. Soc., 360, L82, doi: 10.1111/j.1745-3933.2005.00051.x
  * Ronconi et al. (2019) Ronconi, T., Contarini, S., Marulli, F., Baldi, M., & Moscardini, L. 2019, MNRAS, 488, 5075, doi: 10.1093/mnras/stz2115
  * Ronconi & Marulli (2017) Ronconi, T., & Marulli, F. 2017, A&A, 607, A24, doi: 10.1051/0004-6361/201730852
  * Sahlén et al. (2016) Sahlén, M., Zubeldía, Í., & Silk, J. 2016, ApJ, 820, L7, doi: 10.3847/2041-8205/820/1/L7
  * Salcedo et al. (2022a) Salcedo, A. N., Weinberg, D. H., Wu, H.-Y., & Wibking, B. D. 2022a, MNRAS, 510, 5376, doi: 10.1093/mnras/stab3793
  * Salcedo et al. (2022b) Salcedo, A. N., Zu, Y., Zhang, Y., et al. 2022b, Science China Physics, Mechanics, and Astronomy, 65, 109811, doi: 10.1007/s11433-022-1955-7
  * Sánchez et al. (2017) Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640, doi: 10.1093/mnras/stw2443
  * Schaye et al. (2023) Schaye, J., et al. 2023, Mon. Not. Roy. Astron. Soc., 526, 4978, doi: 10.1093/mnras/stad2419
  * Schmidt (2021a) Schmidt, F. 2021a, JCAP, 04, 032, doi: 10.1088/1475-7516/2021/04/032
  * Schmidt (2021b) —. 2021b, JCAP, 04, 033, doi: 10.1088/1475-7516/2021/04/033
  * Schmidt et al. (2020) Schmidt, F., Cabass, G., Jasche, J., & Lavaux, G. 2020, JCAP, 11, 008, doi: 10.1088/1475-7516/2020/11/008
  * Schmidt et al. (2019) Schmidt, F., Elsner, F., Jasche, J., Nguyen, N. M., & Lavaux, G. 2019, JCAP, 01, 042, doi: 10.1088/1475-7516/2019/01/042
  * Schmittfull et al. (2019) Schmittfull, M., Simonović, M., Assassi, V., & Zaldarriaga, M. 2019, Phys. Rev. D, 100, 043514, doi: 10.1103/PhysRevD.100.043514
  * Schmittfull et al. (2021) Schmittfull, M., Simonović, M., Ivanov, M. M., Philcox, O. H. E., & Zaldarriaga, M. 2021, JCAP, 05, 059, doi: 10.1088/1475-7516/2021/05/059
  * Schuster et al. (2023) Schuster, N., Hamaus, N., Dolag, K., & Weller, J. 2023, J. Cosmology Astropart. Phys, 2023, 031, doi: 10.1088/1475-7516/2023/05/031
  * Schuster et al. (2019) Schuster, N., Hamaus, N., Pisani, A., et al. 2019, J. Cosmology Astropart. Phys, 2019, 055, doi: 10.1088/1475-7516/2019/12/055
  * Scoccimarro (1998) Scoccimarro, R. 1998, Mon. Not. Roy. Astron. Soc., 299, 1097, doi: 10.1046/j.1365-8711.1998.01845.x
  * Scoccimarro (2004) Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007, doi: 10.1103/PhysRevD.70.083007
  * Scoccimarro (2015) —. 2015, Phys. Rev. D, 92, 083532, doi: 10.1103/PhysRevD.92.083532
  * Seljak et al. (2017) Seljak, U., Aslanyan, G., Feng, Y., & Modi, C. 2017, JCAP, 12, 009, doi: 10.1088/1475-7516/2017/12/009
  * Senatore & Zaldarriaga (2015) Senatore, L., & Zaldarriaga, M. 2015, JCAP, 02, 013, doi: 10.1088/1475-7516/2015/02/013
  * Sheth & van de Weygaert (2004) Sheth, R. K., & van de Weygaert, R. 2004, MNRAS, 350, 517, doi: 10.1111/j.1365-2966.2004.07661.x
  * Speagle & Barbary (2018) Speagle, J., & Barbary, K. 2018, dynesty: Dynamic Nested Sampling package, Astrophysics Source Code Library.  http://ascl.net/1809.013
  * Speagle (2020a) Speagle, J. S. 2020a, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Speagle (2020b) —. 2020b, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stadler et al. (2023) Stadler, J., Schmidt, F., & Reinecke, M. 2023, J. Cosmology Astropart. Phys, 2023, 069, doi: 10.1088/1475-7516/2023/10/069
  * Stevens et al. (2023) Stevens, A. R. H., Sinha, M., Rohl, A., et al. 2023.  https://arxiv.org/abs/2312.04137
  * Storey-Fisher et al. (2024) Storey-Fisher, K., Tinker, J. L., Zhai, Z., et al. 2024, ApJ, 961, 208, doi: 10.3847/1538-4357/ad0ce8
  * Sugiyama et al. (2023) Sugiyama, N. S., Yamauchi, D., Kobayashi, T., et al. 2023, MNRAS, 523, 3133, doi: 10.1093/mnras/stad1505
  * Sutter et al. (2014a) Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2014a, MNRAS, 442, 462, doi: 10.1093/mnras/stu893
  * Sutter et al. (2012) Sutter, P. M., Lavaux, G., Wandelt, B. D., & Weinberg, D. H. 2012, ApJ, 761, 187, doi: 10.1088/0004-637X/761/2/187
  * Sutter et al. (2014b) Sutter, P. M., Lavaux, G., Wandelt, B. D., et al. 2014b, MNRAS, 442, 3127, doi: 10.1093/mnras/stu1094
  * Takada et al. (2014) Takada, M., et al. 2014, PASJ, 66, R1, doi: 10.1093/pasj/pst019
  * Talts et al. (2018) Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. 2018, arXiv e-prints, arXiv:1804.06788, doi: 10.48550/arXiv.1804.06788
  * Tanikawa et al. (2013) Tanikawa, A., Yoshikawa, K., Nitadori, K., & Okamoto, T. 2013, New A, 19, 74, doi: 10.1016/j.newast.2012.08.009
  * Tanikawa et al. (2012) Tanikawa, A., Yoshikawa, K., Okamoto, T., & Nitadori, K. 2012, New A, 17, 82, doi: 10.1016/j.newast.2011.07.001
  * Thiele et al. (2023) Thiele, L., Massara, E., Pisani, A., et al. 2023, arXiv e-prints, arXiv:2307.07555, doi: 10.48550/arXiv.2307.07555
  * To et al. (2024) To, C.-H., DeRose, J., Wechsler, R. H., et al. 2024, ApJ, 961, 59, doi: 10.3847/1538-4357/ad0e61
  * Tucci & Schmidt (2023) Tucci, B., & Schmidt, F. 2023.  https://arxiv.org/abs/2310.03741
  * Valogiannis et al. (2023) Valogiannis, G., Yuan, S., & Dvorkin, C. 2023.  https://arxiv.org/abs/2310.16116
  * Valogiannis et al. (2023) Valogiannis, G., Yuan, S., & Dvorkin, C. 2023, arXiv e-prints, arXiv:2310.16116, doi: 10.48550/arXiv.2310.16116
  * van den Bosch et al. (2005) van den Bosch, F. C., Weinmann, S. M., Yang, X., et al. 2005, MNRAS, 361, 1203, doi: 10.1111/j.1365-2966.2005.09260.x
  * Vasudevan et al. (2019) Vasudevan, A., Ivanov, M. M., Sibiryakov, S., & Lesgourgues, J. 2019, JCAP, 09, 037, doi: 10.1088/1475-7516/2019/09/037
  * Verde et al. (2002) Verde, L., et al. 2002, Mon. Not. Roy. Astron. Soc., 335, 432, doi: 10.1046/j.1365-8711.2002.05620.x
  * Verza et al. (2023) Verza, G., Carbone, C., Pisani, A., & Renzi, A. 2023, J. Cosmology Astropart. Phys, 2023, 044, doi: 10.1088/1475-7516/2023/12/044
  * Verza et al. (2022) Verza, G., Carbone, C., & Renzi, A. 2022, ApJ, 940, L16, doi: 10.3847/2041-8213/ac9d98
  * Verza et al. (2019) Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019, J. Cosmology Astropart. Phys, 2019, 040, doi: 10.1088/1475-7516/2019/12/040
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro, F., et al. 2020, Astrophys. J. Suppl., 250, 2, doi: 10.3847/1538-4365/ab9d82
  * Wadekar et al. (2020) Wadekar, D., Ivanov, M. M., & Scoccimarro, R. 2020, Phys. Rev. D, 102, 123521, doi: 10.1103/PhysRevD.102.123521
  * Wang et al. (2022) Wang, K., Mao, Y.-Y., Zentner, A. R., et al. 2022, MNRAS, 516, 4003, doi: 10.1093/mnras/stac2465
  * Wechsler et al. (2022) Wechsler, R. H., DeRose, J., Busha, M. T., et al. 2022, ApJ, 931, 145, doi: 10.3847/1538-4357/ac5b0a
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435, doi: 10.1146/annurev-astro-081817-051756
  * White et al. (2014) White, M., Tinker, J. L., & McBride, C. K. 2014, MNRAS, 437, 2594, doi: 10.1093/mnras/stt2071
  * Xu et al. (2021) Xu, X., Zehavi, I., & Contreras, S. 2021, MNRAS, 502, 3242, doi: 10.1093/mnras/stab100
  * Yoshikawa & Fukushige (2005) Yoshikawa, K., & Fukushige, T. 2005, Publications of the Astronomical Society of Japan, 57, 849–860, doi: 10.1093/pasj/57.6.849
  * Yuan et al. (2024) Yuan, S., Abel, T., & Wechsler, R. H. 2024, MNRAS, 527, 1993, doi: 10.1093/mnras/stad3359
  * Yuan et al. (2018) Yuan, S., Eisenstein, D. J., & Garrison, L. H. 2018, MNRAS, 478, 2019, doi: 10.1093/mnras/sty1089
  * Yuan et al. (2022) Yuan, S., Garrison, L. H., Eisenstein, D. J., & Wechsler, R. H. 2022, MNRAS, doi: 10.1093/mnras/stac1830
  * Yuan et al. (2021) Yuan, S., Garrison, L. H., Hadzhiyska, B., Bose, S., & Eisenstein, D. J. 2021, MNRAS, 510, 3301, doi: 10.1093/mnras/stab3355
  * Yuan et al. (2023) Yuan, S., Hadzhiyska, B., & Abel, T. 2023, Mon. Not. Roy. Astron. Soc., 520, 6283, doi: 10.1093/mnras/stad550