# Small-scale signatures of primordial non-Gaussianity in k-Nearest Neighbour
cumulative distribution functions

William R. Coulton,1 Tom Abel,2,3,4 and Arka Banerjee5
1 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY 10010, USA
2 Kavli Institute for Particle Astrophysics and Cosmology, Stanford
University, 452 Lomita Mall, Stanford, CA 94305, USA
3 Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford,
CA 94305, USA
4 SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA
94025, USA
5 Department of Physics, Indian Institute of Science Education and Research,
Homi Bhabha Road, Pashan, Pune 411008, India

###### Abstract

Searches for primordial non-Gaussianity in cosmological perturbations are a
key means of revealing novel primordial physics. However, robustly extracting
signatures of primordial non-Gaussianity from non-linear scales of the late-
time Universe is an open problem. In this paper, we apply k-Nearest Neighbor
cumulative distribution functions, kNN-CDFs, to the quijote-png simulations to
explore the sensitivity of kNN-CDFs to primordial non-Gaussianity. An
interesting result is that for halo samples with $M_{h}<10^{14}$ M⊙/h, the
kNN-CDFs respond to equilateral PNG in a manner distinct from the other
parameters. This persists in the galaxy catalogs in redshift space and can be
differentiated from the impact of galaxy modelling, at least within the halo
occupation distribution (HOD) framework considered here. kNN-CDFs are related
to counts-in-cells and, through mapping a subset of the kNN-CDF measurements
into the count-in-cells picture, we show that our results can be modeled
analytically. A caveat of the analysis is that we only consider the HOD
framework, including assembly bias. It will be interesting to validate these
results with other techniques for modeling the galaxy–halo connection, e.g.,
(hybrid) effective field theory or semi-analytical methods.

##  1 Introduction

Many theories of the early Universe predict that the statistical distribution
of primordial potential perturbations is close to Gaussian, but with small
deviations (see e.g., Chen, 2010; Achúcarro et al., 2022; Meerburg et al.,
2019, for recent reviews). The structure of these deviations, known as
primordial non-Gaussianity (PNG), encodes the details of the physical
processes governing the evolution of the Universe at that epoch. Primordial
properties ranging from the number of particles present, the masses and spins
of these particles, the strength of interactions, and primordial symmetries
all leave distinct non-Gaussian signatures (e.g. Maldacena, 2003; Creminelli &
Zaldarriaga, 2004; Alishahiha et al., 2004; Chen et al., 2007; Meerburg et
al., 2009; Arkani-Hamed & Maldacena, 2015; Cabass et al., 2023d). Thus,
characterizing the statistical distribution of primordial perturbations is a
powerful way to reveal new information on the early Universel and probe energy
scales far beyond the reach of terrestrial experiments. In this work, we focus
on three templates of non-Gaussianity that are most relevant for Large Scale
Structure (LSS) - local, equilateral, and orthogonal (Komatsu & Spergel, 2001;
Senatore et al., 2010). Each of these generates a unique signature in the
primordial bispectrum, which characterize the skewness of the primordial
perturbations as a function of scale.

To date, observational studies of primordial non-Gaussianity have been driven
by measurements of the bispectrum of the cosmic microwave background (CMB)
anisotropies (e.g., Komatsu et al., 2003; Planck Collaboration IX, 2020) and
the large-scale distribution of galaxies (D’Amico et al., 2022a, b; Cabass et
al., 2022a, b; Cabass et al., 2023b). Whilst no signatures of primordial non-
Gaussianity have yet been detected, large regions of primordial model space
have already been ruled out. New experiments, such as the Dark Energy
Spectroscopic Instrument, SPHEREX, the Simons Observatory, and CMB-S4, will
provide dramatically expanded and more precise data sets that will
significantly improve upon current bounds (Doré et al., 2014; DESI
Collaboration et al., 2016; Abazajian et al., 2016; Ade et al., 2019).
However, these constraints have not yet reached the regions of greatest
theoretical interest that divide qualitatively different regions, such as
strong and weakly coupled physics (Cabass et al., 2023c).

The bispectrum has been used extensively for two reasons: first, it is the
optimal statistic to constrain local, equilateral and orthogonal non-
Gaussianity in the CMB and in the very large-scale distribution of galaxies
(Babich, 2005; Philcox, 2021). Second, analytical tools have been developed
that can accurately model these observations (e.g. Baumann et al., 2012;
Carrasco et al., 2012; Cabass et al., 2023a, with the latter for a review).
However, for measurements of the small-scale distribution of galaxies, where
the signal-to-noise ratio (SNR) is high, and where the relation to the
primordial anisotropies is non-linear, these statements break down. The non-
linear evolution redistributes information from the primordial bispectrum to
not only the late-time bispectrum, but also to the trispectrum, pentaspectrum,
and beyond (where the trispectrum and pentaspectrum are the kurtosis and 5th
moment as a function scale). This means that analyses based purely on the
late-time bisepctrum are not accessing all of the available information. One
approach is to include these higher-order correlation functions in the
analysis, however it is highly challenging to compute these statistics
(Philcox et al., 2021). Further, the ability to model the bispectrum relies on
a perturbative analysis which is typically only valid on large scales. A final
challenge of this approach is that the non-linear processes governing
structure evolution generate late-time non-Gaussianities, even in the absence
of PNG. These non-Gaussianities can mimic the bispectrum signatures of PNG and
thereby bias inferences. When removing these biases, we need to marginalize
over the uncertainties in our understanding of these processes, for example in
how galaxies occupy halos. This significantly degrades the resulting PNG
constraints, especially for equilateral non-Gaussianity (Baldauf et al., 2016;
Lazanu et al., 2017; Baumann & Green, 2022; Cabass et al., 2023c).

In this work we explore the efficacy of k-nearest neighbour cumulative
distributions functions (kNN-CDFs), an alternative summary statistic to the
hierarchy of $N$-point correlations, to constrain PNG. kNN-CDFs describe the
volume-averaged probability of finding at least k objects, in our case dark
matter halos or galaxies, within a sphere of radius $R$. Recent work (Banerjee
& Abel, 2021a, b; Wang et al., 2022; Banerjee & Abel, 2023) have shown that
kNN-CDFs are a powerful way of analyzing large scale structure data sets. In
particular, kNN-CDFs can break parameter degeneracies that are found with
other statistical probes Banerjee et al. (2022). We investigate whether the
response of kNN-CDFs to PNG is distinct from other processes and therefore
whether kNN-CDFs can separate PNG from late-time non-Gaussianities generated
by nonlinear gravitational evolution and galaxy formation. kNN-CDFs are
closely related to the counts-in-cell (CiC) summary statistic. CiCs have been
extensively studied (Bernardeau & Valageas, 2000; Valageas, 2002a; Bernardeau
et al., 2015; Bernardeau & Reimberg, 2016; Uhlemann et al., 2016), including
for constraining PNG from the dark matter field Uhlemann et al. (2018);
Friedrich et al. (2020), and accurate analytical models have been developed
for them (see e.g., Uhlemann et al., 2020, and references therein). By
exploiting the relationship between CiCs and kNN-CDFs, we can obtain
analytical models that describe our results and thereby replicate one
desirable feature of bispectra analyses.

To examine the impact of primordial non-Gaussianity on kNNs we first use the
quijote-png suite of simulations (Coulton et al., 2022a). These simulations
were designed to test PNG analysis methods and have been used to studying
bispectrum statistics of the matter (Coulton et al., 2022a; Jung et al.,
2022b) and halo fields (Coulton et al., 2022b; Jung et al., 2022a), the halo
mass function (Jung et al., 2023) and machine learning statistics (Jung et
al., 2023; Flöss & Meerburg, 2023). Combined with the original quijote suite
of simulations (Villaescusa-Navarro et al., 2020), we can explore how kNNs
respond to cosmological parameters jointly with primordial non-Gaussianity.

This paper is structured as follows: in Section 2 we briefly review kNN-CDFs
and in Section 3 we describe the simulations used in this work. In Section 4
we apply kNN-CDFs to catalogs of dark matter halos and characterize the key
features induced by PNG and their similarity to features arising from other
key parameters. In Section 5 repeat this analysis on a set of mock galaxy
catalogs, compare the simulated catalogs to the CiC model and perform a Fisher
forecast of the constraining power. We present our conclusions and outlook in
Section 6. In Appendix A we discuss how different choices in the definition of
our sample impact the results and in Appendix B we discuss the convergence of
our numerical Fisher forecasts.

Figure 1:  The normalized response of five dark matter halo kNN-CDFs to
primordial non-Gaussianity, cosmological and bias parameters. These responses
are measured from the quijote-png simulations for all halos with $M_{h}\geq
3.2\times 10^{13}$ M⊙/h at $z=0$ and, for ease of comparison, each is
normalized by its largest value. The response of the kNN-CDFs to equilateral
non-Gaussianity is different from the other parameters.  Figure 2:  The
signature of primordial non-Gaussianity in the kNN-CDFs displays a strong mass
dependence as is demonstrated by examining three halo mass samples: a high
mass sample ($M_{h}\geq 1\times 10^{14}$ M⊙/h, green ), an intermediate mass
sample ($6\times 10^{13}$M⊙/h $\leq M_{h}<1\times 10^{14}$ M⊙/h, orange) and a
lower mass sample ($3.2\times 10^{13}$M⊙/h $\leq M_{h}<6\times 10^{13}$ M⊙/h,
blue). For comparison we plot the total halo sample as a black dotted line.
Figure 3:  An examination of the kNN-CDFs obtained from halos at $z=0.5$ for
four different halo mass samples. The configuration is otherwise the same as
Fig. 2.  Figure 4:  A comparison of the constraints on the amplitude of
equilateral non-Gaussianity, $f_{\mathrm{NL}}^{\mathrm{equil.}}$, from the
bispectrum and power spectrum up to $k=0.5$h/Mpc (as reported in Coulton et
al., 2022b), with those from the kNN-CDFs. This analysis uses the
$M_{h}>3.2\times 10^{13}$ halo sample from the 1 Gpc3 box at $z=0.0$ and
considers two cases: constraining only $f_{\mathrm{NL}}^{\mathrm{equil.}}$ and
constraining $f_{\mathrm{NL}}^{\mathrm{equil.}}$ whilst marginalizing an
effective bias parameter, the minimum halo mass of the catalog
$M_{\mathrm{min}}$. The distinct impact of primordial non-Gaussianity on the
kNN-CDFs for this sample means that, unlike the bispectrum and power spectrum,
the signature of primordial non-Gaussianity is not strongly degenerate with
the bias parameter.

##  2 Overview of k-Nearest Neighbour cumulative distributions functions

k-Nearest Neighbour cumulative distributions functions simply measure the
volume-averaged probability of finding at least $k$ objects with a sphere of
radius, $R$. In this work we denote these statistics as ${\rm kNN}(R)$. They
provide an alternative means of accessing the information contained within all
the orientation averaged $N$-point correlation functions. There are several
useful features of kNN-CDFs: first they can be computed in a very efficient
manner (see e.g., Banerjee & Abel, 2021a, for details). Second, they naturally
can be applied to catalogs of objects, as is obtained from observations,
without needing the data to be grided (see e.g, Jing, 2005; Sefusatti et al.,
2016, for issues arising from gridding data sets). Third, the kNN-CDFs sample
all regions of the volume equally rather than focusing on over-dense regions,
yielding sensitivity to underdense regions in the volume. In fact, the $1{\rm
NN}$-CDF is directly related to the Void Probability Function (VPF) (White,
1979).

The analysis here largely follows the methods described in Banerjee & Abel
(2021a) and we refer the reader to Banerjee & Abel (2021a) for more details.
To measure the CDFs we use the following procedure:

  1. 1.

We generate a set of distributed volume-filling points distributed on a grid.
We call these the query points. To ensure dense sampling, we typically use
$10$ times as many query points as there are data points.

  2. 2.

We build a $k$d-tree using the data points in a chosen simulation (halo
positions or galaxy positions). For each query point, we use this tree
structure to compute the distances from the query points to the $k$-th nearest
neighbor data points. In this terminology $k=1$ refers the closest data point
to a query point, $k=2$ refers to the second nearest data point from a query
point, and so on.

  3. 3.

For a particular $k$, we sort the list of distances to generate the empirical
$k$NN-CDF for that $k$. We repeat the same step for different values of $k$.

  4. 4.

We then repeat the measurements for each simulation, and compute the average
at a given cosmology and set of galaxy parameters.

To compute the response of the kNN-CDFs to different parameters we use finite
differences as

$\displaystyle\frac{\partial\overline{kNN}}{\partial\theta}=\left<\frac{{kNN}|_{\theta=\theta+\delta\theta}-{kNN}|_{\theta=\theta-\delta\theta}}{2\delta\theta}\right>,$
(1)

where $\delta\theta$ is a small step in the parameter of interest and the
expected is computed as the average over the different simulation
realizations.

The most significant difference between this work and Banerjee & Abel (2021a)
is that we consider all objects in the catalog, rather using a fixed number of
objects. This was used to allow easier comparison to previous quijote-png
analyses (e.g., Coulton et al., 2022b). In Appendix A, we present the results
for samples with fixed number density and find a qualitatively similar
picture.

##  3 Simulations

We use the quiojte and quijote-png simulations for a detailed description we
refer the reader to Villaescusa-Navarro et al. (2020); Coulton et al. (2022a).
These are a suite of gadget-3 (Springel, 2005) dark-matter, N-body simulations
run at a cosmology consistent with Planck (Planck Collaboration VI, 2020):
$\Omega_{m}=$0.3175, $\Omega_{\Lambda}=0.6825$, $\Omega_{b}=0.049$,
$\sigma_{8}=0.834$, $h=0.6711$ and $n_{s}=0.9624$. The initial conditions were
modified to include three shapes of primordial bispectrum: local, equilateral
and orthogonal. The primordial bispectrum,
$B_{\Phi}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$, is defined as

$\displaystyle\langle\Phi(\mathbf{k}_{1})\Phi(\mathbf{k}_{2})\Phi(\mathbf{k}_{3})\rangle=\delta^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})B_{\Phi}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$
(2)

where $\Phi(\mathbf{k})$ is the primordial potential at wavenumber
$\mathbf{k}$. The three primordial bispectra considered here probe a range of
different physical process (see e.g, Chen, 2010; Meerburg et al., 2019;
Achúcarro et al., 2022, for overviews) and have the forms:

$\displaystyle B^{\mathrm{local}}_{\Phi}(k_{1},k_{2},k_{3})=$ $\displaystyle
2f_{\mathrm{NL}}^{\mathrm{local}}P_{\Phi}(k_{1})P_{\Phi}(k_{2})+\text{ 2
perm.},$ (3) $\displaystyle
B^{\mathrm{equil.}}_{\Phi}(k_{1},k_{2},k_{3})=6f_{\mathrm{NL}}^{\mathrm{equil.}}\Big{[}-P_{\Phi}(k_{1})P_{\Phi}(k_{2})+\text{
2 perm.}$
$\displaystyle-2\left(P_{\Phi}(k_{1})P_{\Phi}(k_{2})P_{\Phi}(k_{3})\right)^{\frac{2}{3}}+P_{\Phi}(k_{1})^{\frac{1}{3}}P_{\Phi}(k_{2})^{\frac{2}{3}}P_{\Phi}(k_{3})$
$\displaystyle+\text{5 perm.}\Big{]},$ (4)

and

$\displaystyle B^{\mathrm{ortho-LSS}}_{\Phi}(k_{1},k_{2},k_{3})=$
$\displaystyle 6f_{\mathrm{NL}}^{\mathrm{ortho-
LSS}}\left(P_{\Phi}(k_{1})P_{\Phi}(k_{2})P_{\Phi}(k_{3})\right)^{\frac{2}{3}}\Bigg{[}$
$\displaystyle-\left(1+\frac{9p}{27}\right)\frac{k_{3}^{2}}{k_{1}k_{2}}+\textrm{2
perms}+\left(1+\frac{15p}{27}\right)\frac{k_{1}}{k_{3}}$
$\displaystyle+\textrm{5 perms}-\left(2+\frac{60p}{27}\right)$
$\displaystyle+\frac{p}{27}\frac{k_{1}^{4}}{k_{2}^{2}k_{3}^{2}}+\textrm{2
perms}-\frac{20p}{27}\frac{k_{1}k_{2}}{k_{3}^{2}}+\textrm{2 perms}$
$\displaystyle-\frac{6p}{27}\frac{k_{1}^{3}}{k_{2}k_{3}^{2}}+\textrm{5
perms}+\frac{15p}{27}\frac{k_{1}^{2}}{k_{3}^{2}}+\textrm{5 perms}\Bigg{]},$
(5)

where $P_{\Phi}(k)$ is the primordial potential power spectrum,
$f_{\mathrm{NL}}^{\mathrm{X}}$ is the amplitude of each bispectrum and

$\displaystyle p=\frac{27}{-21+\frac{743}{7(20\pi^{2}-193)}}\,.$ (6)

We refer the reader to Coulton et al. (2022a) for a detailed description of
the implementation of these bispectra.

For each shape 500 simulations are run with an amplitudes of the primordial
bispectrum, $f_{\mathrm{NL}}^{X}=100$, where $X$ denotes the shape, and 500
with $f_{\mathrm{NL}}^{X}=-100$. The seeds of the $f_{\mathrm{NL}}^{X}=100$
and $f_{\mathrm{NL}}^{X}=-100$ simulations are matched to reduce cosmic
variance. The quiojte simulations varied a set of cosmological parameters
above and below the fiducial value, for use in Fisher forecasts. We used
simulations that varied the amplitude of the linear matter fluctuations on
smoothed on $8{\rm Mpc}/h$ scales, $\sigma_{8}$, the Hubble constant, $h$, the
fractional density of matter, $\Omega_{m}$, and the primordial spectral tilt,
$n_{s}$. For each parameter there are 500 simulations with the parameter
perturbed above and 500 perturbed below the fiducial value, again with matched
seeds. We also use the 15,000 simulations run at the fiducial cosmology to
compute covariance matrices.

For the analysis of the dark matter halos, we analyze the same samples as used
Coulton et al. (2022b); Jung et al. (2022a). Specially we use the friends-of-
friends (FoF, Davis et al., 1985) halo catalog at redshifts, z, $z=0.0$ and
$z=0.5$, and only include halos with mass, $M_{h}$, $M_{h}\geq 3.2\times
10^{13}$ M⊙/h. We work in redshift space by displacing the halos along the
line of sight ($\hat{\mathbf{z}}$ axis) according to their velocity.

We use a halo occupation distribution (HOD) to generate mock galaxy catalogs
from the simulations. Within the HOD framework used here, galaxies are
assigned to halos probabilistically based solely on the halo mass, i.e.
$P(N_{\mathrm{gal}}|M_{h})$. In this work we use the Zheng et al. (2007)
formulation that decomposes the total number of galaxies in a halo into
central and satellite contributions as
$N_{\mathrm{gal}}=N_{\mathrm{central}}+N_{\mathrm{satellite}}$. The central
galaxies follow a Bernoulli distribution with mean

$\displaystyle\langle
N_{\mathrm{central}}\rangle=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{\log
M_{h}-\log M_{\mathrm{min}}}{\sigma_{\log M}}\right)\right]$ (7)

and the satellite galaxies are Poissonian distribution with rate

$\displaystyle\langle N_{\mathrm{satellite}}\rangle=\langle
N_{\mathrm{central}}\rangle\left(\frac{M_{h}-M_{0}}{M_{1}}\right)^{\alpha}.$
(8)

The parameters $M_{\mathrm{min}}$ and $\sigma_{\log M}$ set the minimum mass
of halos that host galaxies and the width of the transition to hosting a
central galaxy. The parameters $M_{0}$, $M_{1}$ and $\alpha$ control the power
law distribution of the satellite galaxies. The central galaxies are placed at
the center of the halo, with the halo’s velocity, whilst the satellite
galaxies are distributed according to a NFW profile with velocities set
acording to the isotropic Jeans equations (Navarro et al., 1996; Łokas &
Mamon, 2001). We use these velocities to displace the galaxies along the line
of sight to produce catalogs in redshift space.

Biagetti & et al. (prep) derived a set of best fit HOD parameters for the
quijote-png simulations, such that the galaxy catalogs matched the CMASS BOSS
galaxy survey at $z=0.5$. We used a set of parameters motivated by those fits:
$M_{\mathrm{min}}=2.2\times 10^{13}$ M⊙/h, $M_{0}=2.8\times 10^{13}$ M⊙/h,
$M_{1}=1.78\times 10^{14}$ M⊙/h, $\sigma_{\log M}=0.15$ and $\alpha=0.5$. We
use these parameters to generate catalogs at $z=0.0$. The purpose of this
choice is to demonstrate properties of the kNN-CDFs and provide a comparison
to previous works at $z=0.0$. Thus this galaxy catalog is not designed to
match any specific experiment. Note that minimum dark matter halo mass used
for the galaxy catalogs is $M_{h}=1.3\times 10^{1}3M_{\odot}/h$, so lower than
that used in analyses of the dark matter halos.

##  4 The impact of PNG on dark matter halo k-Nearest Neighbour cumulative
distribution functions

In Fig. 1 we show how kNN-CDFs distributions for five different numbers of
neighbours respond to PNG, variations in cosmological parameters and a simple
bias parameter ($M_{\mathrm{min}}$). Interestingly the kNN-CDFs statistics,
when compared across different numbers of nearest neighbours, respond
differently to equilateral non-Gaussianity than to all other parameters.

To understand this further we break the halo catalog into three subsets: a
high mass sample, $M_{h}\geq 1\times 10^{14}$ M⊙/h, an intermediate sample,
$6\times 10^{13}$M⊙/h $\leq M_{h}\leq 1\times 10^{14}$ M⊙/h , and a low mass
sample $3.2\times 10^{13}$M⊙/h $\leq M_{h}\leq 6\times 10^{13}$ M⊙/h. The
results are shown in Fig. 2. There is a complex mass dependence of these
responses and for some mass bins the local PNG shows a similar response. These
results show similarities to effects seen in the halo mass function (see e.g.,
LoVerde et al., 2008; Wagner et al., 2010; Jung et al., 2023); the response of
the halo mass function to equilateral and local PNG changes sign at
$~{}7\times 10^{13}M_{\odot}$/h and $1\times 10^{14}M_{\odot}$/h. These
similarities suggest a common underlying cause. Changes in the number density
will impact the kNN-CDFs; the simplest case of a rescaling of the mass
function leads to a horizontal shift of the kNN-CDF. The response of the halo
mass function to PNG is complex (see e.g, Fig 1 of Jung et al., 2023) with the
number of high mass halos being boosted, whilst the number of low mass halos
is reduced. As different mass halos are clustered to differing extents we
expect a complex signature in the kNN-CDFs. A second suggestive piece of
evidence to support this hypothesis is that the redshift evolution of the kNN
signature mirrors the effects seen in the mass function. This can be seen in
Fig. 3 where we show the response of the kNN-CDFs at redshift $z=0.5$, again
split in four mass samples. The distinctive feature moves to lower mass, as
does the signature in the halo mass function.

A key challenge in constraining equilateral non-Gaussianity is disentangling
it from non-Gaussianity introduced by the non-linear evolution of the LSS. The
distinct impact of PNG on the kNN-CDFs, for certain samples, means that the
degeneracy with these late time effects will be significantly reduced. To
explore this we perform a simple Fisher forecast for constraints using kNN-
CDFs measurements. In this forecast we use kNN-CDFs with the following number
of neighbours: 1, 2, 4, 8, 16, 32, 64 and 128. We cut the kNN-CDFs at a
minimum scale of 10 Mpc/h and cut the tails of the distributions, where kNN-
CDF$<0.005$ or kNN-CDF$>0.995$. We use the halo catalog with $M_{h}>3.2\times
10^{13}$M⊙/h at z=0.0.

An interesting question is what likelihood describes the kNN-CDFs.
Characterizing the distribution of the kNN-CDFs is complex and in this work we
consider an alternative avenue. We choose to compress the statistics and then
assume a Gaussian likelihood for the compressed statistics. From other studies
(Anbajagane et al., 2023) , we know that the likelihoods of the CDFs are very
close to Gaussian, as long as we stay away from the tails and do not sample
the CDF too densely. If the likelihood of the kNN-CDFs was known, the data
could be compressed lossesly into a set of summary statistics that, as quasi
maximum likelihood estimators, are Gaussian distributed (see e.g., Lehmann &
Casella, 2006; Alsing & Wandelt, 2018). Here compress the kNN-CDFs
distribution functions using the moped compression (Heavens et al., 2000) as

$\displaystyle\hat{\mathbf{\theta}}_{i}=\frac{\partial\overline{kNN}}{\partial\theta_{i}}\mathcal{C}^{-1}\left(kNN-\overline{kNN}\right),$
(9)

where $\hat{\mathbf{\theta}}_{i}$ are the compressed statistics and
$\overline{kNN}$ and $\mathcal{C}$ are the mean and covariance of the kNN-CDFs
measurements. As the kNN-CDFs are not Gaussian, this compression losses
information. However, the compressed statistics are well approximated by a
Gaussian distribution - this can be understood through the central limit
theorem. We then compute forecast parameter constraints as

$\displaystyle\sigma(\theta_{i})^{2}=F^{-1}_{ii}$ (10)

where the Fisher information is given by

$\displaystyle
F_{ij}=\frac{\partial\hat{\mathbf{\theta_{I}}}}{\partial{\theta}_{i}}\Sigma_{IJ}\frac{\partial\hat{\mathbf{\theta_{J}}}}{\partial{\theta}_{j}},$
(11)

and $\Sigma_{IJ}$ is the covariance of the summary statistics. We compute the
derivatives for compression and the Fisher forecast numerically as described
in Coulton et al. (2022a). We split the simulations into two disjoint sets:
the first set is used for the compression and the second for the Fisher
forecast. A more detailed discussion of Fisher forecasts using compressed
statistics is given in Coulton & Wandelt (2023).

As seen in Fig. 4, marginalizing over a simple bias parameter dramatically
degrades the power-spectrum and bispectrum prediction, whilst leaves the kNN-
CDF result largely unchanged. It is well known that the impact equilateral
non-Gaussianity of the galaxy power-spectrum and bispectrum is highly
degenerate with galaxy bias (Baldauf et al., 2016; Lazanu et al., 2017;
Baumann & Green, 2022); our results, even with a simple bias model, reproduce
this. This degeneracy significantly degrades our constraints and limits the
scales that can be used. Effective field theory (EFT) approaches provide a
systematic way to marginalize over these effects (Baumann et al., 2012),
however no systematic method exists for scales beyond the validity of EFT and
the large degeneracy means obtaining a robust result from non-linear scales
will be difficult. On the other hand, this result demonstrates that kNN-CDFs
can effectively disentangle this bias parameter and PNG. This suggests that
they many provide a path to robust, small-scale PNG measurements and in the
next section we test this hypothesis more stringently with an extended, and
more realistic, bias model. It is also interesting that the unmarginalization
bispectrum and kNN constraints are very similar.

##  5 Extracting signatures of PNG from galaxy k-Nearest Neighbour
distributions

Next we apply the kNN-CDFs to the mock galaxy sample described in Section 3.
The choice of this galaxy sample is motivated by the results of Section 4: it
is interesting to see whether, for a galaxy sample not dominated by
contributions from halos with $M_{h}>10^{14}$ M⊙/h, the equilateral PNG
signature remains.

In this section, we primarily focus on exploring what can be learnt about
equilateral non-Gaussianity with kNN-CDFs.

Figure 5:  The signature of primordial non-Gaussianity, cosmological and HOD
parameters in the kNN-CDFs of the galaxy sample. The response of this galaxy
sample to primordial non-Gaussianity is distinct from the response to
cosmological and HOD parameters.

###  5.1 The impact of PNG on galaxy kNN-CDFs

Using the HOD galaxy sample we compare how cosmological parameters, HOD
parameters, and primordial non-Gaussianity impact the properties of the galaxy
sample. As seen in Fig. 5, the impact of equilateral PNG in this sample
remains distinct from all other contributions. Note that several of the HOD
parameters are highly degenerate. In particular $\sigma_{\log M}$ and $\log
M_{\mathrm{min}}$ are essentially perfectly degenerate.

The fiducial HOD used in this work, as described in Section 3, assumes that
the number of galaxies in each dark matter halo is only a function of the
halo’s mass. However, numerous studies of hydrodynamical simulations and semi-
analytic models have shown that other properties, such as concentration and
environment, are important (e.g., Gao et al., 2005; Gao & White, 2007;
Wechsler et al., 2006; Croton et al., 2007; Li et al., 2008; Bose et al.,
2019; Hadzhiyska et al., 2020; Xu et al., 2021; Hadzhiyska et al., 2021). The
dependence on additional halo properties is called ‘assembly bias’.
Observational evidence that assembly bias is important for current two-point
analyses is mixed: some works find assembly bias to be important (Zentner et
al., 2019; Yuan et al., 2021, 2022; Contreras et al., 2023) and others finding
no strong evidence (Lin et al., 2016; Niemiec et al., 2018; Salcedo et al.,
2022; Yuan et al., 2023a; Rocher et al., 2023). The reconciliation of these
mixed results lies in the scales used in the analyses and the properties of
the samples. Given the increasing sensitivity of upcoming surveys, the
importance of assembly bias in simulations and small scale observations, and
the hints that assembly bias may be more important for beyond two-point
statistics (Yuan et al., 2023c); we explore how one model of assembly bias
impacts our kNN-CDF analysis.

We explore assembly bias using the Luminous Red Galaxy (LRG) extended HOD from
Yuan et al. (2022); Yuan et al. (2023c). This model uses the local density
around dark matter halos as the secondary halo parameter that modulates the
number of halos within a dark matter halo. This is implemented by the
following modifications

$\displaystyle\log M^{\mathrm{new}}_{\mathrm{min}}=\log
M_{\mathrm{min}}+B_{\mathrm{cent}}\left(\delta_{m}^{\mathrm{rank}}-0.5\right),$
$\displaystyle
M^{\mathrm{new}}_{0}=\frac{M_{0}}{M_{\mathrm{min}}}M^{\mathrm{new}}_{\mathrm{min}},\text{\,
and}$ $\displaystyle\log M_{1}^{\mathrm{new}}=\log
M_{1}+B_{\mathrm{sat}}\left(\delta_{m}^{\mathrm{rank}}-0.5\right),$ (12)

where $B_{\mathrm{cent}}$ and $B_{\mathrm{sat}}$ characterize the level of
assembly bias and $\delta_{m}^{\mathrm{rank}}$ is the rank of the local over-
density, smoothed on a scale of $5$Mpc/h.

In Fig. 5, we show how response of the kNN-CDF measurements to variations of
the two assembly bias parameters, around the fiducial values of 0. We see
that, for this galaxy sample, the impact on the kNN-CDFs measurements is not
similar to the signature of equilateral non-Gaussianity. However, further work
is required to test whether this holds for other galaxy samples and other
secondary parameters (such as concentration or local tidal field).

###  5.2 Analytical Modeling

kNN-CDFs are closely related to counts-in-cell statistics, which characterize
the distribution of number counts, or the over density, smoothed over some
scale. More precisely, measuring all the kNN-CDFs for all numbers of nearest
neighbours at a single radius is exactly equivalent to the pdf of counts-in-
cell (CiC) statistic, smoothed with a spherical top hat, evaluated a single
smoothing scale. Thus, a measurement of all nearest neighbours of the kNN-CDFs
is exactly equivalent to measuring the CiC pdf at all smoothing scales.
However, typical kNN-CDFs analyses only consider a subset of all the kNN-CDFs
nearest neighbours distributions (here we have used nearest neighbours
1-8,16,32,64 and 128), whilst typical CiC analyses consider only a handful of
smoothing scales. From this view, the two statistics are highly complementary.

There has been extensive work on CiC statistics (e.g., Fosalba & Gaztanaga,
1998; Valageas, 2001, 2002b; Bernardeau et al., 2014; Uhlemann et al., 2016;
Ivanov et al., 2019; Uhlemann et al., 2020) including on using CiCs
measurements for constraining PNG (e.g, Gaztanaga & Fosalba, 1998; Valageas,
2002c; Uhlemann et al., 2018; Friedrich et al., 2020). Most recently Friedrich
et al. (2020) found that the counts-in-cell matter field pdf is a powerful
probe of primordial non-Gaussianity that can be modelled analytically. This
model for the matter field can be combined with work on the galaxy-matter
connection to obtain predictions for the galaxy CiC pdf. In this work we use
the cosmomentum package developed in Friedrich et al. (2020) with the galaxy-
matter connections from Friedrich et al. (2018, 2022).

These CiC analytical advancements can be equally applied to understand kNN-
CDFs measurements, thanks to the intimate connection between the two
statistics. As analytical modeling of the kNN-CDFs has been considered in past
work (Banerjee et al., 2022), here we simply use the CiC models to validate
our results. To do so we map the kNN-CDF measurements onto the CiC framework
and compare our measurements to the CiC predictions. We perform this
comparison, rather than mapping the analytical predictions on the kNN-CDF
frame, primarily as the model requires a rescaling of the non-linear variance
at each smoothing scale, see Friedrich et al. (2020) for a detailed discussion
of this. This rescaling is performed once for the CiC pdf, as it is a single
smoothing scale, but must be performed for each radius of the kNN-CDFs. As
this equates to a rescaling of each point for a single kNN-CDFs, it is only
through consistency across kNN-CDFs, with many different numbers of nearest
neighbours, can the theory model be rigorously tested. Tests across different
numbers of nearest neighbours of kNN-CDFs effectively equates to the CiC
picture! A more minor reason is that the analytical method is only valid for
the bulk on the CiC pdf and these cuts have a simpler expression in the CiC
framework.

To generate our analytical predictions we use the Lagrangian bias method from
Friedrich et al. (2022), matching the galaxy linear bias, $b_{1}$, and number
density, $\bar{n}$ to that of our simulated galaxy catalogs. Further we use
our simulations to compute the non-linear variance for each smoothing scale
and use this to calibrate our analytical predictions as described in Friedrich
et al. (2020). We refit the parameters of the non-Poissonian shot noise (c.f.
Eq. 24 of Friedrich et al., 2022) to match our simulations. We find the
following parameterization of the non-Poissonian shot noise describes the
bulk-pdf behavior of our simulations

$\displaystyle\frac{\langle N^{2}|\delta_{m}\rangle-\langle
N|\delta_{m}\rangle^{2}}{\langle
N|\delta_{m}\rangle}=\begin{cases}\alpha_{0}+\alpha_{1}\delta_{m},&\text{if
$\delta_{m}<0$}\\\ \alpha_{0},&\text{otherwise},\end{cases}$ (13)

where $\alpha_{0}=0.95$ and $\alpha_{1}=-0.35$ control the deviation from
Poissonian noise. An alternative could be a quadratic bias model, however for
simplicity we use piece-wise linear model.

In Fig. 6 we compare measurements of the derivatives of the first 250 kNN-CDFs
of our galaxy sample to the analytical predictions. For both the smoothing
scales shown, $15$ Mpc/h and $28.25$ Mpc/h, and for all three parameters
($\sigma_{8}$, $\Omega_{m}$ and $f_{\mathrm{NL}}^{\mathrm{equil.}}$) we find
generally good agreement between the simulations and theoretical model. This
agreement provides strong evidence that the signatures of PNG in the kNN-CDFs
are physical, and not an artifact of our simulation methods. This eases any
potential concerns on the limitations of the simulations, such as on the
impact of resolution, finite volumes, the approximation generation of
primordial non-Gaussianity or missing components (such as baryons). The
analytical method also helps develop intuition behind our observed features.

The theoretical model from Friedrich et al. (2020) and Friedrich et al. (2022)
focused on modelling the bulk of the pdf and thus the model is not expected to
work in the tails of the distribution. Further, Friedrich et al. (2020) noted
that resolution effects can lead to small discrepancies between the model and
simulations. Combined these effects are thought to explain the observed
differences seen in our measurements.

Figure 6:  A comparison between the analytical model presented in Section 5.2,
solid lines, and the kNN-CDF measurements converted to counts-in-cells for two
smoothing scales, dotted lines. The error bars denote the error on the mean.

###  5.3 Parameter Constraints

Using these results we consider a Fisher forecast using the galaxy sample. As
in Section 4 we perform the Fisher forecast in the compressed data space. As
$\sigma_{\log M}$ and $\log M_{\mathrm{min}}$ are perfectly degenerate for our
sample, we only consider one of these as a free parameter, $\log
M_{\mathrm{min}}$ . The results are shown in Fig. 7. From 1Gpc3 volume the
kNN-CDFs would be able to constrain equilateral non-Gaussianity to
$\sigma(f_{\mathrm{NL}}^{\mathrm{equil.}})=299$. We note that the resolution
limitations of the simulations mean that a realistic sample, which much higher
number densities, would be able to obtain tighter constraints. The value of
this approach can be seen by comparing to results in (Coulton et al., 2022b;
Jung et al., 2022a, 2023), where the sample halo sample is analyzed with
bispectrum and halo mass function statistics.

Whilst $f_{\mathrm{NL}}^{\mathrm{equil.}}$ is not degenerate with any single
parameter, see Fig. 5, there are some mild degeneracies in the 9 dimensional
space. Quantitatively, the marginalized constraint is 2.4 times larger than
the case when the eight other parameters are fixed. Given the large degree of
flexibility within the 9 parameter space, it is not surprising that
degeneracies appear. The degradation seen here is similar to that found when
marginalizing over the same cosmological parameters and only one, simple bias
parameter in a bispectrum halo analysis (Coulton et al., 2022b; Jung et al.,
2022a), implying that the size of kNN degeneracies is significantly smaller.

Given the large parameter space, it is unsurprising that parameter
degeneracies would develop. The information in the kNN-CDFs primarily comes
from small-scale clustering. Therefore, a kNN-CDF analysis could be combined
with other cosmological probes such as baryon acoustic oscillation
measurements, supernova, CMB power spectra measurements or even large-scale
galaxy power spectrum and bispectrum measurements. These probes are all highly
complementary and could be used to break degeneracies with the kNN-CDFs. In
Fig. 7 we show two examples of how this could be beneficial. In the first case
we add a prior on the Hubble parameter, this prior is from Riess et al. (2022)
and in the second case we include independent priors on the cosmological
parameters based on Planck 2018 cosmology (Planck Collaboration et al., 2020).
In both of these cases we see significantly reduced degeneracies
$f_{\mathrm{NL}}^{\mathrm{equil.}}$. With the Planck priors the
$f_{\mathrm{NL}}^{\mathrm{equil.}}$ constraint is only $40\%$ larger than the
case where all other parameters are fixed. This demonstrates that kNN-CDFs can
separate late-time physics, such as how galaxies occupy halos, from the
primordial signatures. Note that these priors are purely demonstrative, future
experiments would offer significantly improved priors.

Figure 7:  Constraints on equilateral non-Gaussianity and a set of
cosmological and HOD parameters. Whilst the impact of equilateral non-
Gaussianity is not degenerate with any single parameter, there are more
complex degeneracies in the large 9 dimensional parameter space. As the kNN-
CDF information primarily comes from small-scale, $<50$ Mpc/h, clustering
these measurements can easily be combined with probes of the background
expansion, large-scale clustering and CMB measurements. Thus, we compare
constraints from only the kNN-CDFs (blue) with constraints including a prior
on the Hubble constant consistent with local distance measurements (Riess et
al., 2022) or a Planck 2018 prior on the cosmological parameters
($\sigma_{8}$, $\Omega_{m}$, $h$ and $n_{s}$). Combining just with the Hubble
prior allows most of the equilateral degeneracies to be removed. These
constraints are from our galaxy sample at $z=0.0$ with a volume of $1$ Gpc3.

##  6 Discussions

In this work we have simulated the impact of three types of primordial non-
Gaussianity on dark matter halo and galaxy kNN-CDFs from the quiojte-png
simulation suite. The signatures of primordial non-Gaussianity in kNN-CDFs are
distinct from changes in cosmological parameters and also from bias parameters
(a simple $M_{\mathrm{min}}$ for dark matter halos and five halo occupation
distribution parameters for the galaxies). The signature of PNG was
characterized by examining multiple mass samples and two redshifts. The mass
and redshift evolution show features similar to the halo mass function and
therefore suggest the kNN-CDFs are accessing related physical effects. The
scales used in this analysis ($10-70$ Mpc/h) are smaller than typical
primordial non-Gaussianity analyses (Nishimichi et al., 2020; D’Amico et al.,
2022b; Philcox et al., 2022) and so present a highly complementary analysis
approach. An interesting future topic would be to combine these approaches.
This would utilize both the optimal analysis of the large scales with the
bispectrum and small scale information probed by the kNN measurements. Further
kNN-CDFs have already been shown to powerfully constrain bias parameters (Yuan
et al., 2023c) and this could further enhance large-scale bispectrum
measurements.

To further validate and build intuition for our results, we compared them to
theoretical predictions. This was done by exploiting the close relationship
between kNN-CDFs and counts-in-cells. The impact of PNG on counts in cells for
dark matter has been studied before (Uhlemann et al., 2018; Friedrich et al.,
2020) and accurate analytic tools have been developed for those results.
Mapping our kNN-CDFs results into the counts-in-cells (CiC) frame we found
reasonable good agreement with the theoretical prediction. This provides a
stringent validation of simulations and helps demonstrate that spurious
artifacts of the initial conditions, which can arise from higher order effects
of the IC generation method, are negligible. Further, this provides a
theoretical framework to further understand the signatures found in the
simulations. The close relationship between kNN-CDFs and CiCs means that our
results can equally be thought of as both an extension of the results of
Friedrich et al. (2020) to dark matter halos and galaxies and an examination
of the CiC space (counts at fixed scales) from a different perspective
(thresholded counts as a function of scale).

In the primordial space, the primordial non-Gaussianities are bispectra as a
function of both scale and configuration – the relative magnitudes of the
three wavevectors. The theoretical model that accurately describes our kNN-
CDFs is only a function of the skewness at different scales of the linear
density field, as demonstrated in Eq. 52 of Friedrich et al. (2020). This has
two consequences for the kNN-CDFs: first, this states that distinct signature
of equilateral PNG, shown in Fig. 1 and Fig. 5, arises due to the scale
dependence of the equilateral PNG skewness. This is connected to the halo mass
function, which shows similar mass dependence and is discussed in Section 4,
as the skewness on different scales controls the response of the halo mass
function to PNG (Chongchitnan & Silk, 2010; LoVerde & Smith, 2011). Second,
this suggests that kNN-CDFs are not able to fully access all the information
contained in the primordial bispectrum, as the skewness averages over the
configuration of the bispectrum and may “wash out” some of the signal. An
interesting question is how much of the primordial signal is lost. To answer
that we compare the information in the skenwness of the linear density field
to the bispectrum of the same field. The results, shown in Table 1, indicate
that for equilateral non-Gaussianity almost all of the information can be
accessed via measurements of the skewness. For the other types of non-
Gaussianity, especially orthogonal, more information is lost.

There are several interesting future directions. First ,the galaxy sample
considered here does not represent an observational sample. It was constructed
to explore the HOD degeneracies for a sample populating the ‘low mass’ halos
in our simulation. Our simulation mass resolution leads to a minimum halo mass
is very close to the mean mass of observational samples (e.g., BOSS, unWISE,
DESI LRGs and More et al., 2015; Krolewski et al., 2020; Yuan et al., 2023a).
Considering a more realistic sample, which requires higher resolution
simulations, would be a valuable next step. Second, whilst the catalogs used
in this analysis were in redshift space, the kNN-CDFs used here were
isotropic. Recent work by Yuan et al. (2023c) has shown that using 2
dimensional kNN-CDFs allows more information to be extracted from cosmological
data sets. A third interesting extension, also proposed in Yuan et al. (2023c)
, is to compute the kNN-CDFs not from a set of random points, but from the
data-points themselves. This statistic was shown to be highly effective in
constraining HOD parameters and thus could be effective for further reducing
degeneracies between PNG and late-time physics. The results in the main text
used all the halos/galaxies in the catalog. This induced a dependence on the
total number of objects in the sample. In Appendix A we explore one means of
obviated this: using a fixed number density. However, for a realistic survey
we would desire to use all observed objects (rather than downsampling to a
desired number density). An alternative would be to normalize the radius by
the mean separation. Similarly this analysis worked in the simplified geometry
of a periodic box at a single redshift. Exploring the impacts of non-trivial
geometries, light-cone effects (e.g. Yuan et al., 2023b), observational masks
and combining multiple redshifts would set the stage for an analysis of
observations. A final direction would be to consider in more detail the
effects of assembly bias – for example including velocity bias or
investigating alternative secondary halo parameters. These topics will be the
subject of future work.

In conclusion, these results suggest that kNN-CDFs could be a powerful
statistic for separating out primordial non-Gaussianity and late-time physical
processes.

## Acknowledgements

The authors are very grateful to Sandy Yuan, Oliver Philcox, Francisco
Villaescusa-Navarro, David Spergel, Oliver Friedrich, Cora Uhlemann and the
Quijote-PNG collaboration for useful discussions. This work was supported by
collaborative visits funded by the Cosmology and Astroparticle Student and
Postdoc Exchange Network (CASPEN). This work was also supported by U.S.
Department of Energy grant DE-AC02-76SF00515 to SLAC National Accelerator
Laboratory managed by Stanford University.

| 1$\sigma$ Fisher constraint
---|---
PNG Type | Bispectrum | Skewness
Local | 4.4 | 6.1
Equilateral | 18 | 20
Orthogonal | 8.0 | 15
Table 1: A comparison of the constraining power obtainable from measurements
of the skewness, as a function of scale, and the full bispectrum of the linear
field at $z=0.0$. These results are obtained from an analytic Fisher forecast,
with a $k_{\mathrm{max}}=0.5$h/Mpc, for both analyses.

##  Appendix A Fixed Number Density

In this appendix we investigate how a different analysis choice, computing the
kNN-CDFs from a fixed number of data points, alters our conclusions. Using a
fixed number of objects, rather than all objects, removes the sensitivity to
the total number of objects. In Fig. 8 we recompute the dark matter halo kNN-
CDFs using samples of 150000 randomly chosen dark matter halos at $z=0.0$ with
$M_{h}>3.2\times 10^{13}$ $M_{\odot}/h$. The response of all kNN-CDFs to
different parameters is different to the case where all the halos are used.
The smooth bell curves are no longer seen as the shifts in the number of halos
has been removed. Despite the difference in shape, the response to equilateral
non-Gaussianity across the different numbers of neighbours is still distinct
in shape from the other cosmological parameters. In Section 4, we discussed
that the kNN-CDFs may be responding to changes in the halo mass function.
Using a fixed total number of halos does not remove this sensitivity. The
hypothesis was that the kNN-CDF signature is related to the different effects
on the number of high and low mass halos, rather than the total number, and
the different clustering of high and low mass halos. This differential effect
would not be removed by using a fixed number of objects, as considered in this
appendix.

Thus while the quantitative details are different, the results of this
appendix qualitatively match the case where we use all the objects in the
catalog.

Figure 8:  The kNN-CDFs of the halo field, as computed as in Fig. 1, however
these are evaluated at fixed number density of $n_{h}=0.00015$ (h/Mpc)3. Note
that to aid comparison of the derivatives they have all been normalized by the
peak value. The signature of primordial non-Gaussianity is still distinct from
the other parameters. This shows that effect does not arise solely from a
change in the total number of halos.

##  Appendix B Convergence Tests

For simulation-based Fisher forecasts it is vital to test that the results are
converged. Unconverged results arise from the noise associated with using a
finite number of Monte Carlo simulations. In Fig. 9 we test the convergence
our our galaxy Fisher forecast, Section 5.3. The slow change in the forecast
constraints with the number of simulations used implies that the forecast is
sufficiently converged and so reliable.

Figure 9:  The compressed Fisher forecast parameter errors as a function of
the fraction of the total number simulations used in the analysis. The error
bars are estimated by the measured variances across different data spilts. See
Coulton & Wandelt (2023) for a more detailed discussion of compressed Fisher
forecast convergences.

## References

  * Abazajian et al. (2016) Abazajian K. N., et al., 2016, arXiv e-prints, p. arXiv:1610.02743
  * Achúcarro et al. (2022) Achúcarro A., et al., 2022, arXiv e-prints, p. arXiv:2203.08128
  * Ade et al. (2019) Ade P., et al., 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Alishahiha et al. (2004) Alishahiha M., Silverstein E., Tong D., 2004, Phys. Rev. D, 70, 123505
  * Alsing & Wandelt (2018) Alsing J., Wandelt B., 2018, MNRAS, 476, L60
  * Anbajagane et al. (2023) Anbajagane D., et al., 2023, arXiv e-prints, p. arXiv:2308.03863
  * Arkani-Hamed & Maldacena (2015) Arkani-Hamed N., Maldacena J., 2015, arXiv e-prints, p. arXiv:1503.08043
  * Babich (2005) Babich D., 2005, Phys. Rev. D, 72, 043003
  * Baldauf et al. (2016) Baldauf T., Mirbabayi M., Simonović M., Zaldarriaga M., 2016, arXiv e-prints, p. arXiv:1602.00674
  * Banerjee & Abel (2021a) Banerjee A., Abel T., 2021a, MNRAS, 500, 5479
  * Banerjee & Abel (2021b) Banerjee A., Abel T., 2021b, MNRAS, 504, 2911
  * Banerjee & Abel (2023) Banerjee A., Abel T., 2023, MNRAS, 519, 4856
  * Banerjee et al. (2022) Banerjee A., Kokron N., Abel T., 2022, MNRAS, 511, 2765
  * Baumann & Green (2022) Baumann D., Green D., 2022, J. Cosmology Astropart. Phys., 2022, 061
  * Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmology Astropart. Phys., 2012, 051
  * Bernardeau & Reimberg (2016) Bernardeau F., Reimberg P., 2016, Phys. Rev. D, 94, 063520
  * Bernardeau & Valageas (2000) Bernardeau F., Valageas P., 2000, A&A, 364, 1
  * Bernardeau et al. (2014) Bernardeau F., Pichon C., Codis S., 2014, Phys. Rev. D, 90, 103519
  * Bernardeau et al. (2015) Bernardeau F., Codis S., Pichon C., 2015, MNRAS, 449, L105
  * Biagetti & et al. (prep) Biagetti M., et al. in prep.
  * Bose et al. (2019) Bose S., Eisenstein D. J., Hernquist L., Pillepich A., Nelson D., Marinacci F., Springel V., Vogelsberger M., 2019, MNRAS, 490, 5693
  * Cabass et al. (2022a) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022a, arXiv e-prints, p. arXiv:2201.07238
  * Cabass et al. (2022b) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022b, arXiv e-prints, p. arXiv:2204.01781
  * Cabass et al. (2023a) Cabass G., Ivanov M. M., Lewandowski M., Mirbabayi M., Simonović M., 2023a, Physics of the Dark Universe, 40, 101193
  * Cabass et al. (2023b) Cabass G., Ivanov M. M., Philcox O. H. E., 2023b, Phys. Rev. D, 107, 023523
  * Cabass et al. (2023c) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2023c, Physics Letters B, 841, 137912
  * Cabass et al. (2023d) Cabass G., Jazayeri S., Pajer E., Stefanyszyn D., 2023d, Journal of High Energy Physics, 2023, 21
  * Carrasco et al. (2012) Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, Journal of High Energy Physics, 2012, 82
  * Chen (2010) Chen X., 2010, Advances in Astronomy, 2010, 638979
  * Chen et al. (2007) Chen X., Huang M.-x., Kachru S., Shiu G., 2007, J. Cosmology Astropart. Phys., 2007, 002
  * Chongchitnan & Silk (2010) Chongchitnan S., Silk J., 2010, ApJ, 724, 285
  * Contreras et al. (2023) Contreras S., Chaves-Montero J., Angulo R. E., 2023, MNRAS,
  * Coulton & Wandelt (2023) Coulton W. R., Wandelt B. D., 2023, arXiv e-prints, p. arXiv:2305.08994
  * Coulton et al. (2022a) Coulton W. R., et al., 2022a, arXiv e-prints, p. arXiv:2206.01619
  * Coulton et al. (2022b) Coulton W. R., et al., 2022b, arXiv e-prints, p. arXiv:2206.15450
  * Creminelli & Zaldarriaga (2004) Creminelli P., Zaldarriaga M., 2004, J. Cosmology Astropart. Phys., 2004, 006
  * Croton et al. (2007) Croton D. J., Gao L., White S. D. M., 2007, MNRAS, 374, 1303
  * D’Amico et al. (2022a) D’Amico G., Lewandowski M., Senatore L., Zhang P., 2022a, arXiv e-prints, p. arXiv:2201.11518
  * D’Amico et al. (2022b) D’Amico G., Donath Y., Lewandowski M., Senatore L., Zhang P., 2022b, arXiv e-prints, p. arXiv:2206.08327
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Doré et al. (2014) Doré O., et al., 2014, arXiv e-prints, p. arXiv:1412.4872
  * Flöss & Meerburg (2023) Flöss T., Meerburg P. D., 2023, arXiv e-prints, p. arXiv:2305.07018
  * Fosalba & Gaztanaga (1998) Fosalba P., Gaztanaga E., 1998, MNRAS, 301, 503
  * Friedrich et al. (2018) Friedrich O., et al., 2018, Phys. Rev. D, 98, 023508
  * Friedrich et al. (2020) Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2020, MNRAS, 498, 464
  * Friedrich et al. (2022) Friedrich O., Halder A., Boyle A., Uhlemann C., Britt D., Codis S., Gruen D., Hahn C., 2022, MNRAS, 510, 5069
  * Gao & White (2007) Gao L., White S. D. M., 2007, MNRAS, 377, L5
  * Gao et al. (2005) Gao L., Springel V., White S. D. M., 2005, MNRAS, 363, L66
  * Gaztanaga & Fosalba (1998) Gaztanaga E., Fosalba P., 1998, MNRAS, 301, 524
  * Hadzhiyska et al. (2020) Hadzhiyska B., Bose S., Eisenstein D., Hernquist L., Spergel D. N., 2020, MNRAS, 493, 5506
  * Hadzhiyska et al. (2021) Hadzhiyska B., Bose S., Eisenstein D., Hernquist L., 2021, MNRAS, 501, 1603
  * Heavens et al. (2000) Heavens A. F., Jimenez R., Lahav O., 2000, MNRAS, 317, 965
  * Ivanov et al. (2019) Ivanov M. M., Kaurov A. A., Sibiryakov S., 2019, J. Cosmology Astropart. Phys., 2019, 009
  * Jing (2005) Jing Y. P., 2005, ApJ, 620, 559
  * Jung et al. (2022a) Jung G., et al., 2022a, arXiv e-prints, p. arXiv:2211.07565
  * Jung et al. (2022b) Jung G., et al., 2022b, ApJ, 940, 71
  * Jung et al. (2023) Jung G., et al., 2023, arXiv e-prints, p. arXiv:2305.10597
  * Komatsu & Spergel (2001) Komatsu E., Spergel D. N., 2001, Phys. Rev. D, 63, 063002
  * Komatsu et al. (2003) Komatsu E., et al., 2003, ApJS, 148, 119
  * Krolewski et al. (2020) Krolewski A., Ferraro S., Schlafly E. F., White M., 2020, J. Cosmology Astropart. Phys., 2020, 047
  * Lazanu et al. (2017) Lazanu A., Giannantonio T., Schmittfull M., Shellard E. P. S., 2017, Phys. Rev. D, 95, 083511
  * Lehmann & Casella (2006) Lehmann E. L., Casella G., 2006, Theory of point estimation. Springer Science
  * Li et al. (2008) Li Y., Mo H. J., Gao L., 2008, MNRAS, 389, 1419
  * Lin et al. (2016) Lin Y.-T., Mandelbaum R., Huang Y.-H., Huang H.-J., Dalal N., Diemer B., Jian H.-Y., Kravtsov A., 2016, ApJ, 819, 119
  * LoVerde & Smith (2011) LoVerde M., Smith K. M., 2011, J. Cosmology Astropart. Phys., 2011, 003
  * LoVerde et al. (2008) LoVerde M., Miller A., Shandera S., Verde L., 2008, J. Cosmology Astropart. Phys., 2008, 014
  * Łokas & Mamon (2001) Łokas E. L., Mamon G. A., 2001, MNRAS, 321, 155
  * Maldacena (2003) Maldacena J., 2003, Journal of High Energy Physics, 2003, 013
  * Meerburg et al. (2009) Meerburg P. D., van der Schaar J. P., Corasaniti P. S., 2009, J. Cosmology Astropart. Phys., 2009, 018
  * Meerburg et al. (2019) Meerburg P. D., et al., 2019, BAAS, 51, 107
  * More et al. (2015) More S., Miyatake H., Mandelbaum R., Takada M., Spergel D. N., Brownstein J. R., Schneider D. P., 2015, ApJ, 806, 2
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Niemiec et al. (2018) Niemiec A., et al., 2018, MNRAS, 477, L1
  * Nishimichi et al. (2020) Nishimichi T., D’Amico G., Ivanov M. M., Senatore L., Simonović M., Takada M., Zaldarriaga M., Zhang P., 2020, Phys. Rev. D, 102, 123541
  * Philcox (2021) Philcox O. H. E., 2021, Phys. Rev. D, 104, 123529
  * Philcox et al. (2021) Philcox O. H. E., Hou J., Slepian Z., 2021, arXiv e-prints, p. arXiv:2108.01670
  * Philcox et al. (2022) Philcox O. H. E., Ivanov M. M., Cabass G., Simonović M., Zaldarriaga M., Nishimichi T., 2022, Phys. Rev. D, 106, 043530
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Planck Collaboration VI (2020) Planck Collaboration VI 2020, A&A, 641, A6
  * Planck Collaboration IX (2020) Planck Collaboration IX 2020, A&A, 641, A9
  * Riess et al. (2022) Riess A. G., et al., 2022, ApJ, 934, L7
  * Rocher et al. (2023) Rocher A., et al., 2023, arXiv e-prints, p. arXiv:2306.06319
  * Salcedo et al. (2022) Salcedo A. N., et al., 2022, Science China Physics, Mechanics, and Astronomy, 65, 109811
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, MNRAS, 460, 3624
  * Senatore et al. (2010) Senatore L., Smith K. M., Zaldarriaga M., 2010, J. Cosmology Astropart. Phys., 2010, 028
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Uhlemann et al. (2016) Uhlemann C., Codis S., Pichon C., Bernardeau F., Reimberg P., 2016, MNRAS, 460, 1529
  * Uhlemann et al. (2018) Uhlemann C., Pajer E., Pichon C., Nishimichi T., Codis S., Bernardeau F., 2018, MNRAS, 474, 2853
  * Uhlemann et al. (2020) Uhlemann C., Friedrich O., Villaescusa-Navarro F., Banerjee A., Codis S., 2020, MNRAS, 495, 4006
  * Valageas (2001) Valageas P., 2001, A&A, 379, 8
  * Valageas (2002a) Valageas P., 2002a, A&A, 382, 412
  * Valageas (2002b) Valageas P., 2002b, A&A, 382, 412
  * Valageas (2002c) Valageas P., 2002c, A&A, 382, 431
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, ApJS, 250, 2
  * Wagner et al. (2010) Wagner C., Verde L., Boubekeur L., 2010, J. Cosmology Astropart. Phys., 2010, 022
  * Wang et al. (2022) Wang Y., Banerjee A., Abel T., 2022, MNRAS, 514, 3828
  * Wechsler et al. (2006) Wechsler R. H., Zentner A. R., Bullock J. S., Kravtsov A. V., Allgood B., 2006, ApJ, 652, 71
  * White (1979) White S. D. M., 1979, MNRAS, 186, 145
  * Xu et al. (2021) Xu X., Zehavi I., Contreras S., 2021, MNRAS, 502, 3242
  * Yuan et al. (2021) Yuan S., Hadzhiyska B., Bose S., Eisenstein D. J., Guo H., 2021, MNRAS, 502, 3582
  * Yuan et al. (2022) Yuan S., Garrison L. H., Hadzhiyska B., Bose S., Eisenstein D. J., 2022, MNRAS, 510, 3301
  * Yuan et al. (2023a) Yuan S., et al., 2023a, arXiv e-prints, p. arXiv:2306.06314
  * Yuan et al. (2023b) Yuan S., Hadzhiyska B., Abel T., 2023b, MNRAS, 520, 6283
  * Yuan et al. (2023c) Yuan S., Zamora A., Abel T., 2023c, MNRAS, 522, 3935
  * Zentner et al. (2019) Zentner A. R., Hearin A., van den Bosch F. C., Lange J. U., Villarreal A. S., 2019, MNRAS, 485, 1196
  * Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760