# The Disordered Heterogeneous Universe:
Galaxy Distribution and Clustering Across Length Scales

Oliver H. E. Philcox<EMAIL_ADDRESS>Department of Astrophysical Sciences,
Princeton University, Princeton, NJ 08540, USA  School of Natural Sciences,
Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA
Center for Theoretical Physics, Department of Physics, Columbia University,
New York, NY 10027, USA  Simons Society of Fellows, Simons Foundation, New
York, NY 10010, USA  Salvatore Torquato<EMAIL_ADDRESS>School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive,
Princeton, NJ 08540, USA  Department of Chemistry, Department of Physics,
Princeton Institute for the Science and Technology of Materials, and Program
in Applied and Computational Mathematics, Princeton University, Princeton, New
Jersey 08540, USA

###### Abstract

The studies of disordered heterogeneous media and galaxy cosmology share a
common goal: analyzing the disordered distribution of particles and/or
building blocks at ‘microscales’ to predict physical properties of the medium
at ‘macroscales’, whether it be a liquid, colloidal suspension, composite
material, galaxy cluster, or entire Universe. The theory of disordered
heterogeneous media provides an array of theoretical and computational
techniques to characterize a wide class of complex material microstructures.
In this work, we apply them to describe the disordered distributions of
galaxies obtained from recent suites of dark matter simulations. We focus on
the determination of lower-order correlation functions, ‘void’ and ‘particle’
nearest-neighbor functions, certain cluster statistics, pair-connectedness
functions, percolation properties, and a scalar order metric to quantify the
degree of order. Compared to analogous homogeneous Poisson and typical
disordered systems, the cosmological simulations exhibit enhanced large-scale
clustering and longer tails in the void and particle nearest-neighbor
functions, due to the presence of quasi-long-range correlations imprinted by
early Universe physics, with a minimum particle separation far below the mean
nearest-neighbor distance. On large scales, the system appears ‘hyperuniform’,
as a result of primordial density fluctuations, whilst on the smallest scales,
the system becomes almost ‘antihyperuniform’, as evidenced by its number
variance. Additionally, via a finite scaling analysis, we compute the
percolation threshold of the galaxy catalogs, finding this to be significantly
lower than for Poisson realizations (at reduced density $\eta_{c}=0.25$ in our
fiducial analysis compared to $\eta_{c}=0.34$), with strong dependence on the
mean density; this is consistent with the observation that the galaxy
distribution contains voids of up to $50\%$ larger radius. However, the two
sets of simulations appear to share the same fractal dimension on scales much
larger than the average inter-galaxy separation, implying that they lie in the
same universality class. We also show that the distribution of galaxies are a
highly correlated disordered system (relative to the uncorrelated Poisson
distribution), as measured by the $\tau$ order metric. Finally, we consider
the ability of large-scale clustering statistics to constrain cosmological
parameters, such as the Universe’s expansion rate, using simulation-based
inference. Both the nearest-neighbor distribution and pair-connectedness
function (which includes contributions from correlation functions of all
order) are found to considerably tighten bounds on the amplitude of quantum-
mechanical fluctuations from inflation at a level equivalent to observing
twenty-five times more galaxies. The pair-connectedness function in particular
provides a useful alternative to the standard three-particle correlation,
since it contains similar large-scale information to the three-point function,
can be computed highly efficiently, and straightforwardly extended to small
scales (though likely requires simulation-based modeling). This work provides
the first application of such techniques to cosmology, providing both a novel
system to test heterogeneous media descriptors, and a tranche of new tools for
cosmological analyses. A range of extensions are possible, including
implementation on observational data; this will require further study on
various observational effects, necessitating high-resolution simulations.

##  I Introduction

From condensed phases of matter to ecological systems to the primordial
distribution of matter in the universe, Nature abounds with examples of
disordered arrangements of interacting entitites that form structures with
diverse geometries and topologies. To understand the collective behavior of
such phenomena, it is vital to have a mathematical formalism that enables a
stochastic description of the constituent objects, particularly with regards
to their spatial distribution and clustering, whether they be carbon atoms,
concrete conglomerates or individual galaxies. The theory of disordered
heterogeneous media [1, 2], which includes techniques from statistical
mechanics [3], provides a natural and powerful machinery with which to equip
ourselves in this venture. In particular, its primary objective is to connect
the properties of the interacting constituents to their large-scale
attributes, such as a material’s bulk transport, mechanical and
electromagnetic properties. This is rigorously done by generally relating the
bulk properties to an infinite set of diverse types of statistical correlation
functions that characterize the microstructures [1], including those that
contain topological information, such as phase connectivity and percolation
characteristics. While this theoretical machinery has been primarily applied
to earth-bound materials, their applicability is far from terrestrial: these
techniques work similarly for any phenomenon that can be treated as complex
disordered heterogeneous media [1], including the spatial distribution and
clustering of galaxies.

Whilst the bulk of cosmological research in the past two decades has focused
on analysis of the ‘cosmic microwave background’ (the radiation signature of
physics in the first $\sim 300\,000$ years, which provides a snapshot of the
early Universe), that concerning the distribution of galaxies using
statistical descriptors has become progressively more important [4, 5, 6],
particularly with the advent of large-scale surveys, such as the forthcoming
Dark Energy Spectroscopic Instrument (DESI) [7] and the Euclid satellite [8].
The distribution of galaxies traces the distribution of matter in the early
Universe [e.g., 4]; as such, it encodes information on a wealth of
cosmological parameters, such as the density of matter. An open question is
how best to analyze the data: most works focus on measuring the correlation
functions of the galaxy distribution, and comparing them to physical models
(either explicitly derived, or numerically simulated) [e.g., 9], though this
is known to be suboptimal in terms of information content. Whilst a number of
alternative statistics have been proposed (including void statistics [e.g.,
10, 11], marked density fields [e.g., 12, 13, 14, 15], Gaussianized fields
[e.g., 16, 17, 18, 19], reconstructed density fields [e.g, 20] field-level
inference [e.g., 21, 22, 23], Minkowski functionals and other topological
descriptors [e.g., 24, 25, 26, 27, 28, 29, 30, 31, 32, 33], and beyond), there
is little consensus on which have practical utility (with most having been
applied only to the dark matter distribution), and few are natural from a
theoretical standpoint. An important insight is that the galaxy distribution
is simply a set of irregularly arranged point-like particles in three-
dimensions; this is mathematically identical to the structure of many
terrestrial materials, including atomic systems, colloids and sphere packing.
As such, both scenarios can be treated with the same mathematical formalisms;
i.e. heterogeneous media and statistical mechanical techniques designed to
quantify the clustering of particles in materials can be used to provide a
practical and well-motivated manner in which to understand the galaxy
distribution.

This work considers the application of a number of statistical descriptors
from the theory of disordered heterogeneous media to characterize structurally
the distribution of galaxies, which we treat as discrete point configuration.
We ask two main questions: (1) what can we learn about cosmology through the
lens of disordered heterogenous media and statistical mechanics? (2) what
condensed matter physics lessons, more broadly, can we learn from cosmological
structures ? As a proof-of-concept, we will consider a number of descriptors
[e.g., 4, 1], including the two- and three-particle correlation functions,
‘void’ and ‘particle’ nearest-neighbor functions, certain cluster statistics,
pair-connectedness functions, percolation properties, and scalar order metrics
to quantify the degree of order. We show how their behaviors in the cosmic
landscape, probed through cosmological dark matter simulations, differs
substantially from that expected from a simple Poissonian distribution of
points as well as well-known homogeneous models of correlated disordered point
patterns, showing that cosmology challenges general expectations of standard
heterogeneous media and statistical mechanical models. In particular, we will
find stronger clustering on large-scales, giving an enhancement in the pair
correlation and pair-connectedness functions and an excess of large scale
voids; these effects arise due to early Universe physics, which source quasi-
long-range correlations in the galaxy distribution and create a hyperuniform
system. Particular interest will be paid to the question of clustering and
phase ‘percolation’; this is a well-understood phenomenon for many models in
condensed matter but can be similarly extended to galaxy distributions, and
yields interesting results, The cosmological case will be found to percolate
faster, but asymptote to the Poisson case if the density is low, and both
scenarios share the same set of critical exponents. Furthermore, we will
consider the utility of descriptors from the theory of disordered
heterogeneous media in cosmological settings, quantifying how they can add
additional information regarding the early and late Universe, finding that the
pair-connectedness function adds significant cosmological information at
minimal additional computational cost. We caution that further work will be
required before the statistics can be applied to observational data: this must
include discussion of redshift-space effects (arising from the conversion of
galaxy redshifts to distance, creating anisotropy with respect to the sample
line-of-sight), and the dependence of descriptors on galaxy properties, such
as luminosities or masses.

Our conclusions to the above questions asked will be the following: (1) the
theory of heterogeneous media and statistical mechanics provides an array of
useful tools that can enhance the utility of galaxy survey datasets,
strengthening the constraints on cosmological parameters and probing novel
features of the distribution, (2) due to its quasi-long-ranged correlations,
galaxy samples exhibit very different behavior to most terrestrial media, and
thus provides an important sandbox for applying and understanding condensed
matter techniques. Whilst we restrict ourselves to galaxy surveys in this
work, they are by no means the only cosmological application of such
statistics: a number of other phenomena could be described by such approaches.
These include two-phase media such as the distribution of cosmic voids (empty
regions of gargantuan extent) and the statistics of ionized hydrogen bubbles
during the ‘reionization’ phase of the Universe. Such areas provide a
bountiful mine from which to derive future work.

The remainder of this work is structured as follows. In §II, we provide an
overview of the statistics used in this work, before presenting the proposed
testing ground (simulated galaxy samples) in §III. Comparison of the
statistics on galactic and Poisson data is shown in §II, with §V providing a
discussion of percolation physics in the two systems. Finally, §VI considers
the utility of a specific statistic, the pair-connectedness function, in
cosmological contexts, before we conclude in §VII.

##  II Statistical Descriptors of Point Configurations

In the section, we define the various statistical descriptors of point
configurations that will be used in the remainder of this work. We principally
adopt notation from the statistical mechanics community (particularly
following [1]), though connect this to the cosmological terminology, when
relevant. Although we principally work in $\mathbb{R}^{3}$, most of the
following discussion remains relevant in other metric spaces. A schematic
illustrating the various statistical descriptors considered in this work is
shown in Fig. 2.

Figure 1: Schematic depicting the various statistical descriptors used in this
work. The black points indicate the positions of random particles (here
visualized in $\mathbb{R}^{2}$), with the groups of colored circles (each of
diameter $D$) demonstrating the clusters. The pair correlation function,
$g_{2}$, counts pairs of particles belonging both to the same cluster (pink
lines), and to different clusters (grey lines), whereas the pair-connectedness
function, $P_{2}$, contains only particles within the same cluster (pink
lines). We show also the nearest-neighbor functions: $H_{V}$ encodes the
distance between a randomly positioned point and the nearest galaxy (red
arrows), whilst $H_{P}$ gives the separation between a galaxy and its closest
neighbor (blue lines). Further details on these statistics are given in §II.

Figure 2: Example of the cosmological simulations used in this work. This
shows a $(250\times 250\times 50)h^{-3}\mathrm{Mpc}^{3}$ slice of the dark
matter distribution from a single Quijote simulation, with the colorbar
indicating the fractional density of simulation particles in each pixel. The
colored points show the positions of dark matter halos; these are used as a
proxy for galaxies. To facilitate analysis of percolation and connectedness,
galaxies are assigned to clusters (indicated by the various colors), via a
burning (or ‘friends-of-friends’) algorithm with a separation corresponding to
reduced density, $\eta$, here set to $0.2$. Throughout this work, length is
given in cosmologists’ units of Megaparsecs divided by the reduced expansion
rate $h\approx 0.7$), in which the mean pairwise particle separation is
$\approx 8h^{-1}\mathrm{Mpc}$.

###  II.1 Correlation Functions

The fundamental quantity describing a discrete set of $N$ points in some large
region in $\mathbb{R}^{d}$ of volume $V$ is the $N$-particle probability
density functions, $P_{N}$, which is defined such that
$P_{N}(\bm{r}_{1},\cdots,\bm{r}_{N})\mathrm{d}\bm{r}_{1}\cdots\bm{r}_{N}$ is
the probability of finding the first particle within $\mathrm{d}\bm{r}_{1}$ of
$\bm{r}_{1}$, the second within $\mathrm{d}\bm{r}_{2}$ of $\bm{r}_{2}$ et
cetera [e.g., 4, 1]. Of more practical use is the $n$-particle probability
density function, marginalized over the positions of the other $(N-n)\geq 0$
particles: this is defined as

$\displaystyle\rho_{n}(\bm{r}_{1},\cdots,\bm{r}_{n})=\frac{N!}{(N-n)!}\int\mathrm{d}\bm{r}_{n+1}\cdots\mathrm{d}\bm{r}_{N}\,P_{N}(\bm{r}_{1},\cdots,\bm{r}_{n},\bm{r}_{n+1},\cdots,\bm{r}_{N}),$
(1)

with
$\rho_{n}(\bm{r}_{1},\cdots,\bm{r}_{n})\mathrm{d}\bm{r}_{1}\cdots\bm{r}_{n}$
being proportional to the probability of finding one indistinguishable
particles within $\mathrm{d}\bm{r}_{1}$ of $\bm{r}_{1}$ et cetera.

For a statistically homogeneous medium, the one-particle density function is a
constant, i.e. $\rho_{1}(\bm{r})\equiv\bar{\rho}$, which is the mean number
density (number of points per unit volume in the infinite-volume or
‘thermodynamic’ limit), commonly labelled $\bar{n}$ in cosmological contexts
(with $\bar{\rho}$ often used to refer to the mean energy density of the
Universe). More generally, for statistically homogeneous systems,
$\rho_{n}(\bm{r}_{1},\cdots,\bm{r}_{n})$ is translationally invariant,
enabling us to re-express it as follows:

$\displaystyle\rho_{n}(\bm{r}_{1},\cdots,\bm{r}_{n})={\bar{\rho}}^{n}g_{n}(\bm{r}_{12},\cdots,\bm{r}_{1n}),$
(2)

where $g_{n}(\bm{r}_{12},\ldots,\bm{r}_{1n})$ is the $n$-particle correlation
function (closely related to the cosmologists’ $n$-point correlation
function), which depends on the relative positions $\bm{r}_{12},\cdots$, where
$\bm{r}_{ij}\equiv\bm{r}_{j}-\bm{r}_{i}$. The two-particle or pair correlation
function $g_{2}$ is particularly important in applications, and is
schematically illustrated in Fig. 2. For translationally invariant point
configurations without long-range order,
$g_{n}(\bm{r}_{12},\cdots,\bm{r}_{1n})\rightarrow 1$ when the points (or
‘particles’) are mutually far from one another, i.e., as
$|\bm{r}_{ij}|\rightarrow\infty$ ($1\leq i<j<n$),
$\rho_{n}(\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{n})\rightarrow\rho^{n}$. Thus,
the deviation of $g_{n}$ from unity provides a measure of the degree of
spatial correlations (positive or negative) between the particles. Note that
for a translationally invariant Poisson (spatially uncorrelated) point
configurations, $g_{n}=1$ is unity for all values of its argument. If the
point configuration is in addition statistically isotropic, the functions
$g_{n}$ are invariant under joint rotations of $\bm{r}_{ij}$, such that
$g_{2}$ is a function only of $|\bm{r}_{12}|\equiv r_{12}$, and $g_{3}$
depends only on $r_{12}$, $r_{13}$ and $r_{23}$.111Due to the conversion of
cosmological redshifts into distances, the observed galaxy density is not
isotropic, but distorted along the line-of-sight to the galaxy survey,
$\hat{\bm{n}}$. As such, the pair correlation function depends on an
additional angle, i.e. $g_{2}(r)\to g_{2}(r,\hat{\bm{r}}\cdot\hat{\bm{n}})$.
We neglect this dependence in this work, but note that it will be important
when the statistics discussed herein are applied to observational data.

In cosmology it is commonplace to consider not the probability density
function of the full ensemble, but a set of realizations of the microscopic
density (often known as the ‘density field’), each denoted by
$\hat{\rho}(\bm{r})$. This gives the probability that there is a particle
within $\mathrm{d}\bm{r}$ of $\bm{r}$ for a specific field (for sufficiently
small $\mathrm{d}\bm{r}$), and can be represented as a sum of $N$ Dirac
deltas: $\hat{\rho}(\bm{r})=\sum_{i=1}^{N}\delta_{\rm D}(\bm{r}-\bm{x}_{i})$.
Averaging over realizations (denoted by the expectation operator
$\mathbb{E}$), we can relate $\hat{\rho}(\bm{r})$ to the $n$-particle
probability density functions

$\displaystyle\mathbb{E}\left[\hat{\rho}(\bm{r}_{1})\cdots\hat{\rho}(\bm{r}_{n})\right]\equiv\rho_{n}(\bm{r}_{1},\cdots,\bm{r}_{n}).$
(3)

If the field is statistically homogeneous, and the volume (i.e.
$\int\mathrm{d}\bm{r}$) sufficiently large, this is equivalent to a spatial
average via the ergodic principle. In cosmological contexts, (3) is usually
adopted, with $\hat{\rho}(\bm{r})$ often referred to as
$\bar{\rho}\left[1+\delta(\bm{r})\right]$ for overdensity field $\delta$.
Additionally, it is conventional to work with disconnected correlation
functions, $\xi^{(n)}$, often known as ‘$n$-point correlation functions’: the
first few satisfy

$\displaystyle\xi^{(2)}(\bm{r}_{12})$ $\displaystyle\equiv$
$\displaystyle\mathbb{E}\left[\delta(\bm{r}_{1})\delta(\bm{r}_{2})\right]\equiv
g_{2}(\bm{r}_{12})-1$ (4)
$\displaystyle\xi^{(3)}(\bm{r}_{12},\bm{r}_{13},\bm{r}_{23})$
$\displaystyle\equiv$
$\displaystyle\mathbb{E}\left[\delta(\bm{r}_{1})\delta(\bm{r}_{2})\delta(\bm{r}_{3})\right]\equiv
g_{3}(\bm{r}_{12},\bm{r}_{13},\bm{r}_{23})-g_{2}(\bm{r}_{12})-g_{2}(\bm{r}_{13})-g_{2}(\bm{r}_{23})+2$

[e.g., 4], and are all zero under Poisson statistics. We will principally work
with the full $g_{n}$ functions in this work, adopting statistical mechanics
conventions.

A particularly important descriptor is the structure factor, $\mathcal{S}(k)$,
which is related to the Fourier transform of the total correlation function
$h(r)\equiv g_{2}(r)-1$:

$\displaystyle\mathcal{S}(\bm{k})\equiv 1+\bar{\rho}\,\tilde{h}(\bm{k})\equiv
1+\bar{\rho}\,\int\mathrm{d}\bm{r}\,e^{i\bm{k}\cdot\bm{r}}\left[h(\bm{r})\right].$
(5)

For a Poisson point distribution, $\mathcal{S}(\bm{k})=1$ for all $\bm{k}$.
The structure factor and the cosmologists’ power spectrum [4], $P(\bm{k})$,
are trivially related to one another via
$\mathcal{S}(\bm{k})={\bar{\rho}}P(\bm{k})$.

The structure factor provides a useful way to quantify large-scale (low
wavenumber $k\equiv|\bm{k}|$) correlation and fluctuation properties of a
point configuration and plays a central role in the hyperuniformity concept.
Hyperuniform states of matter are correlated systems that are characterized by
an anomalous suppression of long-wavelength (i.e., large-length-scale) density
fluctuations compared to those found in garden-variety disordered systems,
such as ordinary fluids and amorphous solids [34, 35]. A hyperuniform (or
superhomogeneous [36]) many-particle system in $d$-dimensional Euclidean space
$\mathbb{R}^{d}$ is one in which (normalized) density fluctuations are
completely suppressed at very large length scales, implying that the structure
factor $\mathcal{S}(\bm{k})$ tends to zero in the infinite-wavelength limit,
i.e.,

$\lim_{|\bm{k}|\to 0}\mathcal{S}(\bm{k})=0.$ (6)

Equivalently, a hyperuniform system is one in which the number variance
$\sigma^{2}_{{}_{N}}(R)\equiv\langle N(R)^{2}\rangle-\langle N(R)\rangle^{2}$
of particles within a spherical observation window of radius $R$ grows more
slowly than the window volume in the large-$R$ limit, i.e., slower than
$R^{d}$. Typical disordered systems, such as liquids and structural glasses,
have the standard asymptotic volume scaling $\sigma^{2}_{N}(R)\sim R^{d}$ and
hence are not hyperuniform. For general translationally invariant point
configation in $\mathbb{R}^{d}$, the local number variance $\sigma^{2}_{N}(R)$
is determined exactly by the pair statistics [34]:

$\displaystyle\sigma_{N}^{2}(R)$ $\displaystyle=$
$\displaystyle\bar{\rho}v_{1}(R)\Big{[}1+\rho\int_{\mathbb{R}^{d}}h(\bm{r})\alpha(r;R)\mathrm{d}\bm{r}\Big{]}$
(7) $\displaystyle=$
$\displaystyle\bar{\rho}v_{1}(R)\Big{[}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\mathcal{S}(\bm{k}){\tilde{\alpha}}(k;R)\mathrm{d}\bm{k}\Big{]},$

where $v_{1}(R)=\pi^{d/2}R^{d}/\Gamma(1+d/2$ is the volume of a
$d$-dimensional sphere of radius $R$, and $\alpha(r;R)$ is the scaled
intersection volume, the ratio of the intersection volume of two spherical
windows of radius $R$ whose centers are separated by a distance $r$ to the
volume of a spherical window, known analytically in any space dimension [1,
37]. Its Fourier transform is the nonnegative function given by

${\tilde{\alpha}}(k;R)=2^{d}\pi^{d/2}\Gamma(1+d/2)\frac{[J_{d/2}(kR)]^{2}}{k^{d}},$
(8)

where $J_{\nu}(x)$ is the Bessel function of order $\nu$.

Consider translationally invariant point configurations that are characterized
by a structure factor with a radial power-law form in the vicinity of the
origin, i.e.,

$\displaystyle\mathcal{S}(\bm{k})\sim|\bm{k}|^{\alpha}\quad\text{for
}|\bm{k}|\to 0.$ (9)

For hyperuniform systems, the exponent $\alpha$ is positive ($\alpha>0$) and
its value determines three hyperuniformity classes corresponding to different
large-$R$ scaling behaviors of the number variance [34, 38, 35]:

$\displaystyle\sigma^{2}(R)$ $\displaystyle\sim\left\\{\begin{array}[]{l
l}R^{d-1},&\alpha>1\text{ (class I)}\\\ R^{d-1}\ln R,&\alpha=1\text{ (class
II)\,.}\\\ R^{d-\alpha},&\alpha<1\text{ (class III)}\end{array}\right.$ (13)

Classes I and III describe the strongest and weakest forms of hyperuniformity,
respectively. States of matter that belong to class I include all perfect
crystals [34, 38], many perfect quasicrystals [38, 39, 40], and ‘randomly’
perturbed crystal structures [41, 42, 43, 44], classical disordered ground
states of matter [34, 45, 46] as well as systems out of equilibrium [47, 48].
Class II hyperuniform systems include some quasicrystals [40], the positions
of the prime numbers [49], and many disordered classical [50, 51, 52, 53, 47]
and quantum [54, 55, 56] states of matter. Examples of class III hyperuniform
systems include classical disordered ground states [57], random organization
models [58] and perfect glasses [47]. Certain disordered hyperuniform systems
are poised at an ‘inverted’ critical point in which the volume integral of the
total correlation function $h(\bm{r})$ is quasi-long-ranged but its volume
integral is bounded [34, 35].

By contrast, for any nonhyperuniform system, the local variance has the
following large-$R$ scaling behaviors [59]:

$\displaystyle\sigma^{2}(R)\sim\begin{cases}R^{d},&\alpha=0\quad\text{(typical
nonhyperuniform)}\\\ R^{d-\alpha},&-d<\alpha<0\quad\text{(anti-
hyperuniform)}.\\\ \end{cases}$ (14)

For a ‘typical’ nonhyperuniform system, the structure factor $\mathcal{S}(0)$
is bounded [35]. In anti-hyperuniform systems, $\mathcal{S}(0)$ is unbounded,
i.e.,

$\lim_{|\bm{k}|\to 0}\mathcal{S}(\bm{k})=+\infty,$ (15)

and hence are diametrically opposite to hyperuniform systems. Anti-
hyperuniform systems include fractals, systems at thermal critical points
(e.g., liquid-vapor and magnetic critical points) [60, 61, 62, 63, 64] as well
as certain substitution tilings [65].

###  II.2 Order Metric

Given the richness of the spectrum of possible microstructures that can arise
in condensed phase systems, an outstanding challenging task has been the
quantification of their degree order/disorder. Scalar order/disorder metrics
have been profitably employed to quantify the degree of order in many-particle
systems, including sphere packings; see [1] & [66] and references therein. Any
scalar order metric $\Psi({\bf R})$ is a well-defined nonnegative scalar
function of a many-particle configuration ${\bf R}$ and if, for any two
configurations ${\bf R}_{A}$ and ${\bf R}_{B}$, $\Psi({\bf R}_{A})>\Psi({\bf
R}_{B})$, we say that configuration ${\bf R}_{A}$ is to be considered more
ordered than configuration ${\bf R}_{B}$. It has been suggested that a good
scalar order metric should have the following additional properties [67]: (1)
sensitivity to any type of ordering without bias toward any reference system;
(2) ability to reflect the hierarchy of ordering between prototypical systems
given by common physical intuition (e.g., perfect crystals with high symmetry
should be highly ordered, followed by quasicrystals, correlated disordered
packings without long-range order, and finally spatially uncorrelated or
Poisson distributed particles); (3) capacity to detect order at any length
scale; and (4) incorporation of both the variety of local coordination
patterns and the spatial distribution of such patterns.

The recently introduced $\tau$ order metric [46] fulfills these requirements
and has been fruitfully employed to characterize the degree of order across
length scales of a diverse set of disordered media [46, 68, 69, 70]. This
metric, which we compute here for the first time for the galaxies, is defined
as

$\displaystyle\tau$ $\displaystyle=$
$\displaystyle\frac{1}{D^{d}}\int_{\mathbb{R}^{d}}\mathrm{d}\bm{r}\,h^{2}(\bm{r})$
(16) $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{d}D^{d}{\bar{\rho}}^{2}}\int_{\mathbb{R}^{d}}\mathrm{d}\bm{k}\,[\mathcal{S}({\bm{k}})-1]^{2},$

where $D$ is a characteristic ‘microscopic’ length scale. This scalar metric
measures deviations of two-particle statistics from that of the Poisson
distribution. Since both positive and negative correlations contribute to the
integral, due to the fact that $h(\bm{r})$ is squared, $\tau$ measures the
degree of translational order across length scales. It clearly vanishes for
the uncorrelated Poisson distribution, diverges for an infinite crystal and is
a positive bounded number for correlated disordered systems without long-range
order (i.e., Bragg peaks). It is interesting to note that the $\tau$ order
metric is closely related to the negative of the excess two-particle entropy
of the system [71].

###  II.3 Nearest-Neighbor Functions

Another well-known set of statistical descriptors that arise in rigorous
bounds on the macroscopic physical properties of disordered heterogeneous
media, such as suspensions of spheres, and employed in the statistical
mechanics of many-particle systems are nearest-neighbor functions [72, 1, 56].
There two types of such functions: ‘void’ and ‘particle’ quantities. The void
and particle nearest-neighbor probability density functions $H_{V}(r)$ and
$H_{P}(r)$, respectively, are defined as follows:

$\displaystyle\begin{array}[]{ccp{3.7in}}H_{V}(r)\,\mathrm{d}r&=&Probability
that a point of the point configuration lies at a distance between $r$ and
$r+\mathrm{d}r$ from an arbitrary point in the space.\end{array}$ (18)
$\displaystyle\begin{array}[]{ccp{3.7in}}H_{P}(r)\,\mathrm{d}r&=&Probability
that a point of the point configuration lies at a distance between $r$ and
$r+\mathrm{d}r$ from another point of the point configuration.\end{array}$
(20)

The associated dimensionless ‘exclusion’ probabilities $E_{V}(r)$ and
$E_{P}(r)$ are defined as follows:

$\displaystyle\begin{array}[]{ccp{3.7in}}E_{V}(r)&=&Probability of finding a
spherical cavity of radius $r$ empty of any points in the point
configuration.\end{array}$ (22)
$\displaystyle\begin{array}[]{ccp{3.7in}}E_{P}(r)&=&Probability of finding a
spherical cavity of radius $r$ centered at an arbitrary point of the point
configuration empty of any other points.\end{array}$ (24)

It follows that the exclusion probabilities are complementary cumulative
distribution functions associated with the density functions and thus are
related to the latter via

$\displaystyle E_{V}(r)=1-\int_{0}^{r}H_{V}(x)\,\mathrm{d}x,\qquad
H_{V}(r)=-\partial_{r}E_{V}(r)$ (25)

and

$\displaystyle E_{P}(r)=1-\int_{0}^{r}H_{P}(x)\,\mathrm{d}x,\qquad
H_{P}(r)=-\partial_{r}E_{P}(r).$ (26)

The moments of $H_{V}(r)$ and $H_{P}(r)$, defined by

$\ell_{V}^{(k)}=\int_{0}^{\infty}r^{k}H_{V}(r)dr,$ (27)
$\ell_{P}^{(k)}=\int_{0}^{\infty}r^{k}H_{P}(r)dr,$ (28)

are particularly useful integral nearest-neighbor measures, with the $k=1$
version of the latter representing the mean nearest-neighbor distance between
particles. The void nearest neighbor function, $H_{V}$, has received some
attention in cosmology, both historically [73, 74, 75], and in recent works,
in particular via the ‘$k$NN’ statistics [76, 77, 78, 79], generalizing the
above to the $k$-th nearest neighbor. This has been shown to yield strong
constraints on cosmological parameters (cf. §VI), and can be modeled semi-
analytically. It is noteworthy that $k$NN statistics and related quantities
have been studied and fruitfully applied in the field of statistical mechanics
[80, 81, 82, 83].

Both the void and particle nearest-neighbor functions generally involve
integrals over all the $n$-particle correlation functions, $\\{g_{n}\\}$
($n=2,3,4,\ldots)$ [72, 76]. While the void and particle nearest neighbor
functions are identical to one another for a Poisson point configuration
(e.g., $E_{V}(r)=E_{P}(r)=\exp(-4\pi r^{3}{\bar{\rho}}/3$) in three
dimensions), they are generally different from one another for correlated
systems, as manifested by their different series representations [72]. Both
the void and particle quantities arise in rigorous bounds on the effective
transport and mechanical properties of heterogeneous media [1]. It is
noteworthy that the void nearest-neighbor functions play a deep role in the
covering problem of discrete geometry [84]. The covering problem asks for the
point configuration that minimizes the radius of overlapping spheres
circumscribed around each of the points required to cover $d$-dimensional
Euclidean space $\mathbb{R}^{d}$ [85]. The above void statistic also bears
some similarities to the ‘void size function’ used in cosmology [e.g., 10,
86]. The latter quantity is usually constructed from smoothed density fields,
with voids identified as spherical regions with a mean density below some
critical threshold, usually $30\%$ of the system’s mean density. This differs
from $H_{V}$ in two key ways: (a) the voids defined by (18) contain no
particles, thus are equivalent to requiring a critical density of zero, (b)
the cosmologists’ void contains no sub-voids: i.e. any empty area of space
within a void cannot be classified as a smaller void, unlike for $H_{V}$.

###  II.4 Clustering and Connectedness Functions

To quantify the geometrical and topological properties of the class of
disordered heterogeneous media consisting of particles distributed throughout
a matrix phase, it is often useful to statistically characterize particle
clusters that are defined according to some connectivity criterion [e.g., 87,
88, 89, 1, 90, 91] For point (zero-dimensional) particles, this can be
achieved by circumscribing each point by spheres of diameter $D$, which
generally may overlap with one another. Such a decoration of the points by
possibly overlapping spheres divides the space into two disjoint regions or
‘phases,’ encompassing points that do and do not lie within a distance $D/2$
of at least one point. Two spheres are deemed to be connected if they overlap.
Defining the reduced density $\eta\equiv\bar{\rho}\pi D^{3}/6$ (in
$\mathbb{R}^{3}$), it is clear that as the diameter $D$ at fixed mean density
$\bar{\rho}$ increases from zero, $\eta$ and the fraction of space occupied by
the spheres will increase and clusters of various sizes will form and grow
[1].

Once clusters have been identifed, one can determine the pair-connectedness
function, $P_{2}(r,\eta)$, where $P_{2}(r,\eta)\times 4\pi r^{2}dr$ is the
conditional probability of finding a particle in a shell of radius $dr$ at
radial distance $r$ from another particle in the same cluster (assuming
statistical homogeneity and isotropy). Equivalently, this quantity gives the
probability that there exists a path from the first to the second point that
never leaves the particle phase, i.e. one that is always within a distance
$D/2$ of at least one particle (cf. Fig. 2).222Formally, this can be defined
as
$P_{2}(r,\eta)\equiv\mathbb{E}\left[\hat{\rho}(\bm{x})\hat{\rho}(\bm{x}+\bm{r})\hat{\Phi}(\bm{x},\bm{x}+\bm{r},\eta)\right]$,
for clustering function
$\displaystyle\Phi(\bm{x},\bm{x}+\bm{r},\eta)=\begin{cases}1&\text{if}\,\,\exists\,\,\bm{\Gamma}:[0,1]\to\mathbb{R}^{3}\,\text{
s.t.
}\,\left\\{\bm{\Gamma}(0)=\bm{x},\quad\bm{\Gamma}(1)=\bm{x}+\bm{r},\quad\int_{0}^{1}\mathrm{d}\gamma\,\varphi_{D}(\bm{\Gamma}(\gamma))\geq\ell(\bm{\Gamma})\right\\}\\\
0&\text{else},\end{cases}$ (29) where we consider all paths $\bm{\Gamma}$
connecting $\bm{x}$ and $\bm{x}+\bm{r}$ for which the integral of
$\varphi(\bm{r})$ (defined as the particle phase, i.e.
$\varphi(\bm{x})=\Theta_{H}\left[\int\mathrm{d}\bm{y}\,\hat{\rho}(\bm{y})\Theta_{H}(D-|\bm{x}-\bm{y}|)\right]$
for Heaviside $\Theta_{H}$) is at least the line length $\ell(\bm{\Gamma})$,
i.e. those passing only through connected regions. This is the connected
contribution to the full pair correlation function, $g_{2}(r)$:

$\displaystyle g_{2}(r)\equiv P_{2}(r,\eta)+B_{2}(r,\eta),$ (30)

where the pair-blocking function $B_{2}(r,\eta)$ gives the correlation between
pairs of particles which do not lie in the same cluster. For $r<D(\eta)$, any
pair of points must be within the same cluster, thus
$P_{2}(r<D,\eta)=g_{2}(r)$.

A related quantity is the direct-connectedness function $C_{2}(r,\eta)$ (also
known as the non-nodal correlation function) [87]. This is the probability
that two points separated by a distance $r$ are connected by a path through
the set of random particles that does not involve nodes (i.e. one that cannot
be broken by a single cut, as in Fig. 2). A general path between two points
contains either zero or at least one node: this permits the Ornstein-Zernike
(OZ) decomposition [87, 88]

$\displaystyle
P_{2}(\bm{r}_{12},\eta)=C_{2}(\bm{r}_{12},\eta)+\bar{\rho}\int\mathrm{d}\bm{r}_{3}\,C_{2}(\bm{r}_{13},\eta)P_{2}(\bm{r}_{32},\eta),$
(31)

labelling $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$. The first quantity on the RHS
contains paths with no nodes, thus involves the first factor of $C_{2}$,
whilst for the second, we integrate over the position of the node closest to
$\bm{r}_{1}$ (assuming a statistically homogeneous field $\bar{\rho}$), noting
that the path from $\bm{r}_{1}$ to $\bm{r}_{3}$ contains no nodes by
definition, yielding another function of $C_{2}$. Finally, the path from
$\bm{r}_{3}$ to $\bm{r}_{2}$ can contain nodes, leading to the final factor of
$P_{2}$. In Fourier-space this gives a simple relation between the pair-
connectedness and direct-connectedness functions:

$\displaystyle\tilde{C}_{2}(k,\eta)=\frac{\tilde{P}_{2}(k,\eta)}{1+\bar{\rho}\,\tilde{P}_{2}(k,\eta)},\qquad\tilde{P}_{2}(k,\eta)=\frac{\tilde{C}_{2}(k,\eta)}{1-\bar{\rho}\,\tilde{C}_{2}(k,\eta)}.$
(32)

To gain intuition for the pair-connectedness function (and related
statistics), it is instructive to consider its form for a Poissonian system
(noting that there is no Gaussian limit, given that we are dealing with
discrete systems). At low densities, it can be computed as a perturbation
series in $\bar{\rho}$ (or, more strictly, in $\eta/\eta_{c}$), first
considering pairs of particles that are directly linked by the covered phase
(i.e. their centers lie within $D(\eta)$), then moving to pairs linked via a
third particle and so on. This leads to the decomposition

$\displaystyle P_{2}^{\rm Poiss}(\bm{r}_{12},\eta)$ $\displaystyle=$
$\displaystyle\Theta_{D}(\bm{r}_{12})+\bar{\rho}\,\left[1-\Theta_{D}(\bm{r}_{12})\right]\int\mathrm{d}\bm{r}_{3}\,\Theta_{D}(\bm{r}_{13})\Theta_{D}(\bm{r}_{32})$
$\displaystyle\,+\,\bar{\rho}^{2}\left[1-\Theta_{D}(\bm{r}_{12})\right]\int\mathrm{d}\bm{r}_{3}\rm\mathrm{d}\bm{r}_{4}\,\Theta_{D}(\bm{r}_{13})\Theta_{D}(\bm{r}_{34})\Theta_{D}(\bm{r}_{24})\left[1-\Theta_{D}(\bm{r}_{14})\right]\left[1-\Theta_{D}(\bm{r}_{23})\right]+\ldots,$

where the Heaviside function $\Theta_{D}(\bm{r})\equiv\Theta_{H}(r-D)$ selects
pairs with separations below $D$. In this expansion, successive terms
integrate over progressively more particle positions with, for example, the
second term averaging over the position of $\bm{r}_{3}$, which must be within
a distance of $D$ from both $\bm{r}_{1}$ and $\bm{r}_{2}$. As such, this
expression is difficult to compute beyond second order (which is
convolutional) and thus rarely used in practice, unless $\eta\ll\eta_{c}$ and
we restrict to small scales. In practice, approximate treatments are usually
adopted, such as via the OZ equation (31) combined with heuristic ‘closure’
relations such as the Percus-Yevick form [92]. These give accurate predictions
for $P_{2}(r,\eta)$ in low density regimes at relatively small $r$. This
stands in contrast to the case familiar from cosmology, when the modeling of
$g_{2}(r)$ becomes progressively more accurate as $r$ increases.

For a general system, a similar decomposition to (II.4) is possible, and takes
the form:

$\displaystyle P_{2}(\bm{r}_{12},\eta)$ $\displaystyle=$ $\displaystyle
g_{2}(\bm{r}_{12})\Theta_{D}(\bm{r}_{12})+\bar{\rho}\,\left[1-\Theta_{D}(\bm{r}_{12})\right]\int\mathrm{d}\bm{r}_{3}\,g_{3}(\bm{r}_{13},\bm{r}_{32})\Theta_{D}(\bm{r}_{13})\Theta_{D}(\bm{r}_{32})$
$\displaystyle\,+\bar{\rho}^{2}\,\left[1-\Theta_{D}(\bm{r}_{12})\right]\int\mathrm{d}\bm{r}_{3}\mathrm{d}\bm{r}_{4}\,g_{4}(\bm{r}_{13},\bm{r}_{34},\bm{r}_{42})\Theta_{D}(\bm{r}_{13})\Theta_{D}(\bm{r}_{34})\Theta_{D}(\bm{r}_{24})\left[\Theta_{D}(\bm{r}_{14})\right]\left[1-\Theta_{D}(\bm{r}_{23})\right]+\ldots.$

In this case, the expansion depends on the correlation functions $g_{n}$,
since there exists background correlations in addition to that induced by the
circumscribed spheres around points. As expected, this implies that
$P_{2}(r<D,\eta)=g_{2}(r)$, with the $n$-th order term involving correlators
of the form $g_{n+2}$. Formally, this expression may be extended to all
orders, via the relation

$\displaystyle\bar{\rho}^{2}P_{2}(r,\eta)=\prod_{k=2}^{\infty}\bar{\rho}^{k}\int\left(\prod_{i=1}^{k}\mathrm{d}\bm{r}_{i}\right)\delta_{\rm
D}(|\bm{r}_{1}-\bm{r}_{k}|-r)g_{k}(\bm{r}_{1},\ldots,\bm{r}_{k})\prod_{i=1}^{k-1}\Theta_{D}(\bm{r}_{i(i+1)})\prod_{m=1}^{k-2}\left[\prod_{j=m+2}^{k}(1-\Theta_{D}(\bm{r}_{mj}))\right].$
(35)

Whilst this form is not particularly useful for analytic treatments, due to
the difficulty inherent in performing the high-dimensional integrals present
for $k>3$ it illustrates how the pair-connectedness function is comprised of
all possible correlation functions, and thereby partly resums information
found at all orders. At very low densities, $P_{2}(r)$ tends to
$\Theta_{D}(r)g_{2}(r)$, thus we do not expect this statistic to add
information; however, as $\eta$ increases, the fraction of information
contributed by the higher-order $g_{n}$ increases, until the system becomes
non-perturbative at $\eta\approx\eta_{c}$ (whence the notion of connectedness
breaks down). One may ask whether an OZ-like equation with some closure
relation can be used to provide an approximate analytic form for
$P_{2}(r,\eta)$ in the general case. Unfortunately, this is far from trivial,
since any scheme only involving $g_{2}$ will miss any contributions to $P_{2}$
from $g_{3}$ and above, which are of particular cosmological interest,
especially when performing parameter inference in conjunction with $g_{2}$ (as
is the case below). We leave further treatment of this problem to future work.

###  II.5 Continuum Percolation

Percolation describes the appearance of a phase transition in the system,
which, in the above case, corresponds to the emergence of long-range
connectivity in the point cloud due to arbitrarily large clusters of points,
again defined by spheres of some diameter. This phenomenon is of relevance in
a wide variety of physical settings such as the transport of fluid in porous
systems, the appearance of fractures in geological formations, spread of
diseases, and the collapse of gas into stars [93, 1, 89, 2, 91]. This
crossover from non-percolating clusters to the appearance of the incipient
sample-spanning cluster (infinite in the thermodynamic limit) is characterized
by a critical reduced density (also known as a percolation density),
$\eta_{c}={\bar{\rho}}\pi D_{c}/6$ (in $\mathbb{R}^{3}$), where $D_{c}$ is the
critical sphere diameter at fixed mean density $\bar{\rho}$ with the system
said to have percolated for $\eta>\eta_{c}$. For a Poissonian system in
$\mathbb{R}^{3}$, numerical simulations find $\eta_{c}\approx 0.34$ [1],
implying that the connected phase fills about $29\%$ of the space. In addition
to the connectedness functions described above, percolation theory utilizes a
number of other statistical descriptors, which we outline below.

The mean cluster size, $S(\eta)$, gives a simple manner in which to
characterize clustering and percolation at some reduced density $\eta$ [e.g.,
87]. This quantity is simply the mean number of in a cluster containing a
randomly chosen particle. It is directly related to the pair-connectedness
function by the following relation:

$\displaystyle
S(\eta)=1+\bar{\rho}\int\mathrm{d}\bm{r}\,P_{2}(\bm{r},\eta)=\left[1-\bar{\rho}\tilde{C}_{2}(0,\eta)\right]^{-1},$
(36)

where the second equality follows from a Fourier-transform and using (32).
When $\eta>\eta_{c}$, clusters of infinite extent appear, thus
$S(\eta)\to\infty$, and the volume integral of the pair-connectedness function
diverges. Given a form for $C_{2}(r,\eta)$, the second relation provides a
useful manner in which to estimate $\eta_{c}$, by solving
$\bar{\rho}\tilde{C}_{2}(0,\eta_{c})=1$. Furthermore, the behavior close to
the phase transition can be expressed in terms of critical exponents of the
field: in particular, $S(\eta)\sim(\eta-\eta_{c})^{-\gamma}$ and
$C_{2}(r,\eta)\sim r^{-\alpha}$ at large $r$, for $\eta\to\eta_{c}^{-}$, where
$\alpha$ and $\gamma$ are found to be universal for a broad class of physical
models [94, 91].

The mean cluster size can also be written in terms of so-called $s$-mer
cluster statistics for densities below the percolation threshold $\eta_{c}$
[95, 1] , namely,

$\displaystyle
S(\eta)=\frac{\displaystyle\sum_{s=1}^{\infty}s^{2}n_{s}}{\displaystyle\sum_{s=1}^{\infty}sn_{s}},\qquad\eta<\eta_{c}.$
(37)

Here $n_{s}$ is the average number of $s$-mers, clusters containing $s$
particles, per unit number of particles. This representation will be employed
to estimate $S(\eta)$ from simulations.

For a finite (aperiodic) system of size $L$, percolation may be studied by
considering the existence of sample spanning clusters, i.e. clusters of
connected points which reach from the top to the bottom of the system (in some
dimension). In the $L\to\infty$ limit, these clusters will appear only for
$\eta>\eta_{c}$: for finite systems, the behavior can be characterized using
the percolation probability, $\Pi(\eta,L)$, which is the probability that a
realization of size $L$ will contain a sample-spanning cluster. This will be
used to compute the percolation threshold, $\eta_{c}$, in §V, via a finite-
scaling analysis. If such a cluster exists, its mass (i.e. the number of
constituent particles, denoted $M(L,\eta)$) can be used to ascertain the
effective mass fractal dimension, $d_{F}$, of the field. In particular:

$\displaystyle M(L,\eta)\sim L^{d_{F}},\qquad\eta\to\eta_{c}^{-}$ (38)

where $d_{F}\approx 2.42$ for Poisson systems, and any other processes in the
same universality class [cf., 1, 96, 97].

##  III Galaxy Surveys as Point Clouds

Although the statistics described in §II have been principally applied to
study the properties of physical materials, their applicability extends far
beyond the terrestrial regime. In this work, we consider their application to
spectroscopic galaxy surveys, such as those of the upcoming Dark Energy
Spectroscopic Instrument (DESI) and Euclid projects [7, 8]. Such projects will
measure the angular positions and redshifts of $\mathcal{O}(10^{7})$ bright
galaxies from ground- and space-based telescopes, providing a three-
dimensional map of the Universe with unprecedented resolution. Fundamentally,
galaxy surveys measure a set of $N$ galaxy positions with some associated
weights, representing experimental effects. In many typical analyses [e.g.,
98, 99, 100], these are assigned to some coarse grid in $\mathbb{R}^{3}$, and
the associated field taken to be an inhomogeneous Poisson sample of an
underlying continuous field. This is itself modelled as a non-linear
transformation of the underlying dark matter distribution, whose correlation
functions (particularly $g_{2}$ and $g_{3}$) encode early Universe physics,
with the field obeying Gaussian statistics on sufficiently large scales,
before a perturbative (and well understood, [e.g., 101]) regime takes hold.
Explicitly, the microscopic density $\hat{\rho}(\bm{x})$ satisfies:

$\displaystyle\hat{\rho}(\bm{x})\sim\mathrm{Poisson}\left(\bar{\rho}\left[1+\delta_{g}(\bm{x})\right]\right),\qquad\tilde{\delta}_{g}(\bm{k})\sim\mathcal{N}(0,\tilde{h}_{2,g}(k))$
(39)

on sufficiently large scales, where $\delta_{g}$ is a normally-distributed
continuous background field with variance $\tilde{h}_{2,g}(k)$. On small, non-
linear, scales, a variety of galaxy formation processes become important and
the above approach is known to be insufficient. This has led to a flurry of
interest in additional statistics beyond the simple correlation functions.

An alternative to the standard approach is to consider the point cloud traced
by the galaxies as the fundamental object, facilitating direct application of
the clustering techniques described in §II. Rather than working with
observational data directly, this will work will make use of simulated data,
drawn from the publicly available Quijote suite [102], which is a collection
of $40\,000$ realizations of the Universe, each contained within a cubic
volume of size $L=1000h^{-1}\mathrm{Mpc}$.333Following cosmologists
convention, we work in $h^{-1}\mathrm{Mpc}$ units, where
$1\,\mathrm{Mpc}\equiv 10^{6}\,\mathrm{parsec}$ and $h^{-1}\approx 1.4$ is
used to remove a leading scaling. In particular, we use dark matter
simulations that have been evolved down to redshift zero (today), and contain
a set of $\sim 10^{5}$ dark matter halos: spheroidal agglomerations of matter
in which galaxies are known to form.444In this work, we use only halos
containing at least 64 dark matter particles to avoid discreteness effects;
these have masses $M\gtrsim 5\times 10^{12}h^{-1}M_{\odot}$ in our baseline
simulations.

Rather than dealing with the complexities of assigning galaxies to dark matter
halos as a function of their mass (for example using a halo occupation
distribution [103]), we use the positions of the dark matter halos as a direct
proxy for the galaxy positions, which is sufficient for this initial study. As
such, we do not require the Poisson-Gaussian assumptions of (39), and will
utilize the galaxy catalog only as a discrete point cloud. In most scenarios,
we will use the ‘galaxy’ catalogs extracted from $1000$ Quijote
‘High–Resolution’ simulations, each run with the same underlying physical
model, but with varying realizations of the (stochastic) initial conditions. A
section of a typical simulation is shown in Fig. 2. We caution that these
simulations do not fully represent observational data, in particular due to
their limited mass resolution and lack of (magneto-)hydrodynamic effects.
However, their simplified nature makes them ideal for the proof-of-concept
study considered herein, since it allows for a large number of simulations
(and thus determination of accurate covariances). Further work will
necessarily require application of the above tools to higher-resolution
simulations, though these are fewer in number.

##  IV Phenomenological Clustering Statistics

(a) Pair Correlation Function

(b) Structure Factor

Figure 3: Measurements of the pair correlation function and structure factor
from 1000 Quijote cosmological simulations (blue) alongside 1000 Poisson
realizations (red). Both sets of simulations have the same number density
($\approx 1\times 10^{-4}h^{3}\mathrm{Mpc}^{-3}$) and a volume of
$1h^{-3}\mathrm{Gpc}^{3}$, with a mean pairwise particle separation of
$\approx 8h^{-1}\mathrm{Mpc}$. The shaded regions show the statistical
variance between realizations (which grows large on large scales), and we note
that the two statistics are related by a Fourier transform. In Fig. 3a, the
inset shows the pair correlation function in the typical cosmologists’
normalization, plotting $r^{2}h_{2}(r)\equiv r^{2}[g_{2}(r)-1]$; this clearly
brings out the structure imprinted by early Universe physics. At large $r$, we
find $h_{2}(r)\sim r^{-(3+n_{s})}$, where $n_{s}\approx 0.96$: this indicates
the presence of quasi-long-range correlations. Similarly, the inset of Fig. 3b
shows $k\tilde{h}(k)$, equal to the cosmologists’ $kP(k)$ (with a slope of
$k^{n_{s}}$ on large scales). The oscillatory features at $k\sim
0.1h\,\mathrm{Mpc}^{-1}$ arise from acoustic waves in the early Universe.

###  IV.1 Pair Correlation Function

We begin by considering the pair correlation function, $g_{2}(r)$, of the
Quijote simulation suite. This is estimated from the array of galaxy positions
using the corrfunc code [104], which computes the statistic in a set of bins
with centers $\\{r_{a}\\}$ via:

$\displaystyle\widehat{g}_{2}(r_{a})=\frac{1}{N\bar{\rho}\,v_{a}}\sum_{i,j=1}^{N}\begin{cases}1&\text{if
}|\bm{r}_{i}-\bm{r}_{j}|\text{ in bin $r_{a}$}\\\ 0&\text{else,}\end{cases}$
(40)

where $v_{a}$ is the volume of bin $a$ and $\bm{r}_{i}$ is the 3D position of
galaxy $i$ (accounting for periodic wrapping). In Fig. 3a, we display the
obtained $g_{2}(r)$ functions, alongside corresponding results from a Poisson
random sample with the same number density ($\approx 1\times
10^{-4}h^{3}\mathrm{Mpc}^{-4}$) and volume. As expected, the latter is simply
unity everywhere, whilst the former shows considerable structure, and is quite
different to that expected from most simple heterogeneous media [e.g., 1]. On
small scales (with $r\lesssim 0.2h^{-1}\mathrm{Mpc}$, considerably less than
the average galaxy separation of $\approx 8h^{-1}\mathrm{Mpc}$), $g_{2}(r)$
decays to zero; this is as expected, since the galaxies are of finite size and
cannot overlap, enforcing some minimum separation. 555Physically, galaxies can
overlap; however, they would be classed as a single object in this paradigm.
The fact the mean galaxy (particle) separation is about an order of magnitude
greater than the minimum pair separation is atypical behavior for most
condensed phase systems in which these two length scales are comparable to one
another (see Refs. [3] and [1]). At large scales, $g_{2}(r)$ decays to zero
slightly slower than $r^{-4}$ (in fact, $r^{-3.96}$), implying that large-
scale correlations are suppressed; since the decay is between $r^{-d}$ and
$r^{-(d+1)}$ for dimension $d=3$, the system is said to have quasi-long-range
correlations. This has a physical origin: the large-scale behavior of
$g_{2}(r)$ arises from correlations in the Universe’ quantum initial
conditions, imprinted before cosmological inflation. Due to the dynamics of
expansion (and slight breaking of time invariance in ‘slow-roll’ inflation),
these are suppressed on the largest scales. From the inset of Fig. 3a, we note
that $g_{2}(r)$ has considerable structure on intermediate scales, with a
prominent peak at $r\approx 100h^{-1}\mathrm{Mpc}$ sourced by acoustic
oscillations fourteen billion years ago [105]. Clearly, the correlation
properties of the galactic point cloud are very different to those for most
media; this arises due to the combination of Poisson-like placements of
galaxies and an underlying background stochastic field from the early
Universe.

###  IV.2 Structure Factor

The structure factor, $\mathcal{S}(k)$, tells a similar story as the pair
correlation function. This is computed by first assigning the galaxies to a
grid, then computing $\tilde{h}(k)$ via fast Fourier transforms, here
implemented using nbodykit [106]. From Fig. 3b, we observe a super-Poissonian
signature on all scales, with a characteristic decline following a peak at
$k\approx 0.05h\,\mathrm{Mpc}^{-1}$. This peak (known as the ‘equality peak’)
corresponds to a change in the Universe’s expansion rate at early times, with
matter starting to drive the expansion rather than radiation pressure. At
larger $k$, we again see the characteristic acoustic features, here shown by
oscillations in $\tilde{h}(k)$. On the largest scales, the power spectrum or
structure factor $\mathcal{S}(k)\sim k^{0.96}$ in the infinite-wavelength
limit $k\to 0$ (just visible in this plot), with a slope set by the physics of
inflation, and hence because the exponent $\alpha$ in (9) is 0.96, the
Universe belongs to class III hyperuniformity, as defined in relation (13). If
the Universe was scale-invariant according to the Peebles-Harrison-Zeldovich
spectrum with $\mathcal{S}(k)\sim k$ [4], then, because $\alpha=1$, it would
be hyperuniform of class II. Of course, either scenario implies that the
structure factor vanishes in the limit $k\to 0$, the system is hyperuniform
[e.g., 107, 108]. In contrast to many terrestrial media, the large-scale
behavior is well understood, and can be predicted using a variety of
cosmological codes; this occurs since it is an imprint of underlying dark
matter physics, rather than a true pairwise interaction.

###  IV.3 $\tau$ Order Metric

To quantify the degree of order/disorder of the the galaxies, we compute the
metric $\tau$, which is defined by (16), as discussed in §II.2. We take $d=3$
and the unclustered interparticle separation $D=\bar{\rho}^{-1/3}\approx
20h^{-1}\mathrm{Mpc}$ to be the characteristic length-scale. Here $\tau=4.85$
for the cosmological sample, which is to be compared to $\tau=8.37\times
10^{-6}$ for the (finite-volume) Poisson realizations, close to the infinite-
volume expectation of $\tau=0$. This result supports the well-known results
that the galaxy distribution is not purely random (uncorrelated), but instead
is a correlated disordered system. To place the magnitude of $\tau$ for the
galaxies in the context of other models of correlated disordered media, we
compute $\tau$ for the random sequential addition (RSA) process, which is a
time-dependent (nonequilibrium) procedure that generates disordered sphere
packings in $\mathbb{R}^{d}$ [109, 110]. Starting with an empty but large
volume in $\mathbb{R}^{d}$, the RSA process is produced by randomly,
irreversibly, and sequentially placing nonoverlapping spheres into the volume.
If a new sphere does not overlap with any existing spheres, it will be added
to the configuration; otherwise, the attempt is discarded. This procedure is
repeated for ever-increasing volumes; then, an appropriate infinite-volume
limit is obtained. One can stop the addition process at any time $t$,
obtaining RSA configurations with a range of packing fractions $\phi(t)$ up to
the maximal ‘saturation’ value $\phi(\infty)$ in the infinite-time limit,
which for three dimensions is about 0.3812 [110]. Using the data for pair
statistics given in Ref. [110], we find $\tau=6.17$ for saturated RSA packings
in $\mathbb{R}^{3}$, which is close in value to that of the galaxies.

###  IV.4 Local Number Variance

As discussed in §II, the local number variance, $\sigma_{N}^{2}(R)$, can also
be computed from the measured correlation function, and provides a useful tool
with which to assess the system’s order. Here, this is computed from the
measured $g_{2}(r)$ values via (7), and plotted in Fig. 5, alongside its
extrapolation to large $R$, using the well-known large-scale limit,
$\mathcal{S}(k)\sim k^{0.96}$. Notably, we find the number variance to
increase faster than the Poisson case at small $R$, roughly up to the scale
corresponding to the second peak in $r^{2}h(r)$ (arising from the imprint of
acoustic oscillations from the early Universe), then fall to sub-Poisson
values by scales corresponding to the peak in $\mathcal{S}(k)$. A variance
that increases much faster than that for Poisson systems at small $R$ is
unusual for typical correlated disordered systems that have been investigated
in condensed matter physics. Of course, that the large-scale variance
approaches zero indicates that the system is hyperuniform; however, these
scales are difficult to measure with most cosmological surveys.

Figure 4: Local number variance for the cosmological (blue) and Poisson (red)
simulations as a function of scale $R$. This is defined as the variance of the
number of particles found in spheres of radius $R$, and computed directly from
the correlation function shown in Fig. 3. The dashed lines show an
extrapolation to large $R$, using the theoretical asymptotic structure factor
form (with $\mathcal{S}(k)\sim k^{0.96}$). The vertical lines show two
characteristic scales: the sound horizon at recombination (dashed), which
sources acoustic wave in the early Universe, giving the bump in $h(r)$, and
the size of the Universe as it transitioned from radiation- to matter-
dominated (dotted).

Figure 5: Comparison of the cluster size distribution in the cosmological and
Poisson simulations. We plot the mean number of $s$-mers per unit particle
(defined as clusters containing $s$ member) as a function of $s$, normalizing
to the Poisson prediction. Results are shown for various values of the reduced
density $\eta$ and we assume an aperiodic simulation volume of size
$L=800h^{-1}\mathrm{Mpc}$, averaging over $100$ realizations. At large $\eta$,
the two sets of simulations have a similar $n_{s}$-distribution (at least for
small $s$), whilst at low $\eta$, the enhanced clustering in the cosmological
simulations leads to significantly more $s$-mers, with the mean number of
particles per cluster increasing with $\eta$ (cf. Fig. 9).

(a) Void Nearest-Neighbor PDF (b) Particle Nearest-Neighbor PDF

Figure 6: Comparison of the void and particle nearest-neighbor probability
density functions for the cosmological and Poissonian data-set, as defined in
(18) & (20). The left panel shows the probability distribution of finding a
void of radius $r$ in the cosmological (blue) and Poisson (red) datasets,
averaged over 100 realizations, whilst the right panel gives the distribution
function of the distance of a given particle from its nearest neighbor. For
the Poisson case, the dashed lines show theory curves (defined in §II.3),
which are in excellent agreement with the simulations. The void distribution
in the Quijote simulations follows the (mean-density-matched) Poisson
distribution at small $r$, but has an excess of large voids (shown by much
broader tails), due to the quasi-long-range correlations. In contrast, the
particle distribution $H_{P}(r)$ differs between the cosmological and
Poissonian simulations on all scales, notably with an absence of small
separations (due to halo exclusion effects) and an enhancement on large
scales.

###  IV.5 Void and Particle Nearest-Neighbor Functions

In Fig. 6, we depict the nearest-neighbor functions of the two sets of
simulations, which provide an alternative description of the system’s
geometrical and topological properties, as discussed in §II.3 (see also [76]
for a previous discussion of the void function of galaxies, yielding similar
results). These are obtained from the simulations by histogramming the minimum
distance between each pair of particles (for $H_{P}$) or a pair of particles
and a Poisson random particle (for $H_{V}$, determining if this particle lies
within a void). For the Poisson system, we find identical results for the void
and particle nearest-neighbor density functions, as expected, but significant
differences for the Quijote simulations. Whilst the cosmological case has a
similar distribution of small ($r\lesssim 5h^{-1}\mathrm{Mpc}$) voids to that
found in the Poisson realizations, it boasts significantly broader tail
towards large $r$, and thus a somewhat larger mean void size. Specifically,
the first moment of $H_{V}(R)$, $\ell_{V}^{(1)}$, defined by (27), is equal to
$14h^{-1}\mathrm{Mpc}$ ($11h^{-1}\mathrm{Mpc}$) for the cosmological (Poisson)
simulations. Interestingly, the maximal void size (averaged over realizations)
for the Quijote simulations is $42h^{-1}\mathrm{Mpc}$, which is is almost
$50\%$ larger than that for the Poisson system with a maximal void size of
$30h^{-1}\mathrm{Mpc}$. Furthermore, the variance of $H_{V}$, defined as
$\ell_{V}^{(2)}-\left(\ell_{V}^{(1)}\right)^{2}$ is much larger for the
cosmological case: $2.7h^{-2}\mathrm{Mpc}^{2}$ instead of
$1.5h^{-2}\mathrm{Mpc}^{2}$. In particular, the above results suggest that the
galaxies will also boast a lower percolation threshold, foreshadowing what we
describe below.

For the particle distribution, we note that (a) the cosmological simulations
have enhanced large-scale clustering, and thus a broad tail to the nearest-
neighbor distance at large $r$, (b) there is a sharp cut at low $r$, with no
galaxies found within a separation of $\sim 1h^{-1}\mathrm{Mpc}$. This is a
consequence of ‘halo exclusion’; a pair of galaxies cannot be arbitrarily
close, else they would be identified as a single object in the simulation
code. Between these two effects, we find a reduced mean particle nearest-
neighbor distribution in Quijote, indicating that galaxies are more likely to
be found in large-scale clusters. This matches theoretical expectations.
Specifically, the mean nearest-neighbor distance between particles,
$\ell_{P}^{(1)}$, defined by (28),) is equal to $8h^{-1}\mathrm{Mpc}$
($11h^{-1}\mathrm{Mpc}$) for the cosmological (Poisson) simulations, with a
minimum distance of $1.1h^{-1}\mathrm{Mpc}$ ($0.29h^{-1}\mathrm{Mpc}$). In
addition, the variance of the cosmological $H_{P}$, is again larger than the
Poisson case, finding $2.5h^{-2}\mathrm{Mpc}^{2}$ instead of
$1.5h^{-2}\mathrm{Mpc}^{2}$.

###  IV.6 Pair-Connectedness and Direct-Connectedness Functions

The ‘astrophysical’ pair-connectedness function $P_{2}(r)$ has not been
previously studied in the literature, and is of particular interest to both
cosmology and condensed matter physics. To construct this, we first take the
set of $N\sim 10^{5}$ galaxy positions in each Quijote (or Poisson)
simulation, and assign clusters via a ‘burning’ algorithm (often known as
‘friends-of-friends’ in cosmology) [e.g., 111, 112], here using the nbodykit
implementation [106]. This finds sets of points for which each member is
connected to each other member via a path through the clustered phase formed
of spheres of radius $D(\eta)=\left[6\eta/\pi\bar{\rho}\right]^{1/3}$ around
each point, where $\eta$ is the reduced density. Given the set of particles
and cluster memberships (visualized in Fig. 2), we compute the pair-
connectedness function in bins with centers $\\{r_{a}\\}$ via

$\displaystyle\widehat{P}_{2}(r_{a},\eta)=\frac{1}{N\bar{\rho}v_{a}}\sum_{i,j=1}^{N}\begin{cases}1&\text{if
}|\bm{r}_{i}-\bm{r}_{j}|\text{ in bin $r_{a}$}\textbf{ and }\text{$i$, $j$ in
same cluster}\\\ 0&\text{else.}\end{cases},$ (41)

analogous to (40). This is achieved using a custom modification of the
corrfunc code [104], which accepts pairs only if they have the same cluster
index.

Figure 7: Pair-connectedness function, $P_{2}(r,\eta)$, for the Poisson (left)
and cosmological (right) simulations. Results are shown for a variety of
reduced densities $\eta$, corresponding to clustering distances $D(\eta)$ in
the range $[10,20]h^{-1}\mathrm{Mpc}$, and shaded regions show the $1\sigma$
deviations expected from statistical fluctuations. The Quijote simulations
show significantly enhanced correlations on large scales, due to the
underlying correlations of matter imprinted in the early Universe. This
additionally suggests that the galaxy sample will percolate at lower $\eta$:
this will be explored in §V.

Figure 7 displays the pair correlation functions from the Quijote simulations
alongside the more familiar Poisson case. The latter match our expectations:
$P_{2}(r)=1$ for $r<D$ (since all particles with this separation must be in
the same cluster), and $P_{2}(r)$ falls sharply with $r$ for $r>D$ (due to an
absence of large-scale clusters), with an enhanced decline at low $\eta$. For
large $\eta$, the volume integral of $P_{2}(r)$ appears to diverges (at least
in the infinite volume limit), indicating percolation.666Note that the
simulations are computed in periodic boxes, which are known to be suboptimal
for computing $P_{2}(r)$ on the largest scales [94]. This will be addressed in
§V in the context of finite scaling analyses. For the cosmological
simulations, we firstly note that $P_{2}(r)=g_{2}(r)$ for $r<D$, as expected.
At larger $r$, we find that $P_{2}^{\rm sim}(r)>P_{2}^{\rm Poiss}(r)$ for all
choices of $\eta$, indicating that our galaxy catalogs contain more long-range
correlations than a Poisson random field of the same density, and suggesting
that the system will also percolate quicker. In the large $\eta$ limit (i.e.
above percolation), $P_{2}(r)\to g_{2}(r)$, since all points belong to the
same cluster. The cosmological utility of $P_{2}(r,\eta)$ will be discussed in
§VI.

Figure 8: Direct-connectedness function, $C_{2}(r,\eta)$, for the Poisson
(left) and cosmological (right) simulations. This follows Fig. 7, but focusses
on smaller scales, and additionally includes analytic predictions from the
Percus-Yevick model. We additionally normalize all quantities by $g_{2}(r)$,
and exclude values of $\eta$ for which $C_{2}(r,\eta)$ is not well-behaved
(beyond the percolation threshold).

The direct-connectedness function also plays an important role in the analysis
of connected systems, in part due to its appearance in the Ornstein-Zernike
equation (31). Given $P_{2}(r)$, this can be computed using (32), performing
the Fourier transforms numerically via the FFTLog prescription [113]. Figure 8
shows $C_{2}(r,\eta)$ for both the cosmological and Poisson simulations,
alongside the analytic ‘Percus-Yevick’ model, which solves the OZ equation by
asserting that $C_{2}(r>D,\eta)=0$ and $P_{2}(r<D,\eta)=1$ [92].777This is
computed for point objects by using the correspondence with the known (cubic)
form for hard spheres via $C_{2}^{\rm PY,Poiss}(r,\eta)=-C_{2}^{\rm
PY,hard\,sphere}(r,-\eta)$ [114]. For the Poisson case, we find good agreement
between theory and simulations for small $\eta$ (far from the percolation
threshold of $\eta_{c}\approx 0.34$), particularly away from the boundary at
$r=D(\eta)$. We observe very little power from the region with $r>D(\eta)$,
since most intra-cluster path-ways with $r>D(\eta)$ contain at least one node,
and thus do not contribute to $C_{2}(r)$. The cosmological simulations show a
very different behavior, with two peaks observed, with one in similar location
to the Poisson system and one at smaller $r$. This statistic represents the
complexities of the clustering on smaller scales than that typically seen in
$P_{2}$ (with $D(\eta)\sim 10h^{-1}\mathrm{Mpc}$), and the differences arise
primarily due to small-scale physics, such as the restriction that galaxies
cannot be arbitrarily close together. We also note that $C(r,\eta)$ was found
to be ill-behaved for $\eta\gtrsim 0.3$ (due to $P_{2}(r,\eta)$ not being
square integrable), indicating that the cosmological simulations have
percolated by around this value of $\eta$ (cf. §V). In practice, we expect the
percolation threshold to depend on the peculiarities of the galaxy sample in
question: this will be discussed further below.

##  V Percolation and Fractal Dimensions

We now turn to the issue of percolation, following the discussion in §II.5. As
noted earlier, determining the mean cluster size $S(\eta)$ in a system as a
function of $\eta$ is a useful way in which to test whether a system has
reached percolation. We utilize the representation of $S(\eta)$ in terms of
$s$-mer cluster statistic, $n_{s}$, as defined by (37). In principle, we
expect $S(\eta)\to\infty$ as $\eta\to\eta_{c}$; in practice, $S(\eta)\leq N$,
where $N$ is the total number of particles in the box. To account for this, it
is useful to analyze a number of different configurations with different
boxsizes, $L$ (and thereby $N\equiv\bar{\rho}L^{3}$). Here, we construct
(aperiodic) subboxes from the Quijote simulations, with $L$ in the range
$[400,800]h^{-1}\mathrm{Mpc}$ (noting that the majority of our analyses are
restricted to $r<200h^{-1}\mathrm{Mpc}$), and construct analogous Poisson
realizations for each. To examine percolation at each choice of boxsize, we
generate clusters for various values of $\eta$ by varying the sphere radius
$D(\eta)$, and utilizing burning (‘friends-of-friends’) algorithms, as
described above.

Figure 9: Mean number of particles per cluster, $S(\eta,L)$, as a function of
the reduced density $\eta$ and (cubic) simulation boxsize $L$. We show results
for both cosmological data and Poisson realizations, plotting the reciprocal
$1/S(\eta,L)$ averaged over 1000 realizations. Note that $S(\eta)$ is bounded
by the number of particles in the data-set ($N$), shown by dotted lines. We
additionally show the Percus-Yevick prediction in black, which is inaccurate
for all but the smallest $\eta$.

Figure 9 shows the mean cluster size for the two data-sets as a function of
$L$ and $\eta$. In both cases, we observe that $S(\eta,L)$ begins to approach
its asymptotic limit as $\eta$ increases, and, moreover, the limit is
approached faster as the boxsize increases. Extrapolating the Poissonian
results to large $L$, the percolation threshold (whence $S(\eta,L)\approx N$)
appears to be around $\eta_{c}=0.35$, matching previous studies [95]. For the
cosmological simulations, we find a generally slower approach to $\eta_{c}$
(corresponding to a different critical exponent), and additionally a lower
percolation threshold, around $\eta_{c}=0.28$ for $L\to\infty$. As above, this
arises since the galaxy sample contains a stochastic background inhomogeneity,
leading to various areas being super- or sub-Poisson populated in a correlated
manner. The difference is evident even low $\eta$: for $\eta=0.1$, clusters in
the full-volume cosmological simulation contain an average of $\approx$ 15
particles, whilst those in the Poisson realizations contain $\approx 2$ only.

It is further instructive to consider the size distribution of clusters, via
the average number of $s$-mers, $n_{s}$, as defined in (37). This is shown in
Fig. 5 for a suite of cosmological and Poissonian boxes at $L=800$ with a
variety of values of the reduced density $\eta$. At low $\eta$, we find that
the ratio of cosmological and Poissonian simulations is a strongly increasing
function of $s$, with the largest slopes seen for small reduced densities. In
this limit, the system is far from percolation, thus large clusters are rare
in both systems. The enhanced correlations in the galaxy distribution seen in
the cosmological case increase the probability of an $s$-mer forming (at fixed
$s$), giving this stark difference in behavior. As $\eta$ approaches the
percolation threshold, the $n_{s}$-ratio becomes roughly constant with $s$;
this indicates that the additional galaxy correlations impacts only the
largest $s$-mers, as we are dominated by the clustering signal imprinted by
the circumscribed spheres, rather than any intrinsic effects.

To measure the percolation threshold of the cosmological simulations in a
robust fashion, we perform a finite scaling analysis, following the approach
of [115], originally formulated in [116]. In essence, this computes the
percolation probability ($\Pi(\eta,L)$, defined as the fraction of
realizations containing a cluster for which the circumscribed spheres overlap
with both the top and bottom of the box) for the simulations at various values
of $L$ and $\eta$ and extrapolate using asymptotic scaling relations to find
the $L\to\infty$ limit. Figure 10 shows the obtained percolation probability
distribution for both sets of simulations as a function of the volume filling
fraction $\phi$. This is computed numerically for each simulation from the
probability that a randomly chosen point within the box is within a distance
$D/2$ from the nearest particle, i.e. whether it is within the sphere phase;
for the Poisson case, this is asymptotically equal to $1-e^{-\eta}$. The
behavior seen in Fig. 10 is qualitatively similar for the Poisson and
cosmological system: the percolation probability is small for low $\phi$
(whence the typical extent of the cluster is far below $L$), and asymptotes to
unity at large $\phi$. As the boxsize increases, the transition becomes
sharper, and asymptotes to a Heaviside function in the $L\to\infty$ limit. It
is also clear that the cosmological simulations percolate at smaller values of
$\eta$ than the Poisson realizations; this is as expected, and indicates their
enhanced clustering due to background inhomogeneities.

Figure 10: Percolation probability, $\Pi(\phi,L)$, obtained from 1000 Quijote
and Poisson simulations, for various non-periodic boxsizes, $L$, and volume
filling fractions, $\phi$. This is defined as the fraction of realizations
containing a sample spanning cluster, and asymptotes to a Heaviside function
centered at the percolation threshold for $L\to\infty$ (shown as a dotted
line). Points represent the values computed from simulations, whilst the solid
lines show a fit using the sigmoid function of (42). Using a finite scaling
analysis, we find the percolation thresholds of $\eta_{c}=0.252$ ($0.343$) for
the Quijote (Poisson) simulations, though we caution that the cosmological
result depends on the sample density.

To extract the percolation thresholds from $\Pi(\phi,L)$, we fit the data to
the phenomenological sigmoid model of [115], as shown in Fig. 10:

$\displaystyle\Pi(\phi,L)\approx\frac{1}{2}\left[1+\mathrm{tanh}\left(\frac{\phi-\phi_{c}(L)}{\Delta(L)}\right)\right]$
(42)

where $\phi_{c}(L)$ and $\Delta(L)$ are the percolation volume fraction and
width at boxsize $L$ and $\phi$ is the volume filling fraction obtained as
described above. Asymptotically, $\Delta(L)\sim L^{-1/\nu}$ and
$\phi_{c}(L)-\phi_{c}\sim L^{-1/\nu}$ for critical exponent [e.g., 1]; by
fitting for $\nu$ from the obtained values of $\Delta(L)$, we can thus obtain
$\phi_{c}\equiv\lim_{L\to\infty}\phi_{c}(L)$. Here, we find a critical
exponent of $\nu=0.85\pm 0.03$ ($0.88\pm 0.03$) for the Quijote (Poisson)
simulations, with a corresponding percolation threshold of $\phi_{c}=0.223$
($0.290$) or $\eta_{c}=0.252$ ($0.343$), each with a statistical error around
$0.001$. The Poissonian case matches standard results [e.g., 1], and, as
foreshadowed in Figs. 7 & 9, the cosmological system percolates at lower
densities, due to the additional clustering signature imprinted by early-
Universe and galaxy formation physics. In addition, the fact that the two sets
of simulations appear to share the same critical exponent $\nu$ suggests that
they belong to the same universality class, as do other correlated disordered
systems [117].

It is important to note that the percolation thresholds found herein are not a
universal property of galaxy distributions; rather, they depend on the galaxy
sample in question. To explore this, we have repeated the analysis using a
galaxy sample with half the density of the fiducial sample, and another
including dark matter halos (i.e. galaxies) down to half of the aforementioned
minimum size. For the former case (with $\bar{\rho}\sim 0.5\times
10^{-4}h^{3}\mathrm{Mpc}^{-3}$, we find that the percolation threshold for the
cosmological sample increases to $\eta_{c}=0.285$, whilst remains the same for
the density-matched Poissonian sample (as expected). This can be rationalized
by noting that the galaxies roughly follow Poisson statistics above a
stochastic background, caused by the matter density; if $\bar{\rho}$ is
reduced, the Poisson part of the stochasticity becomes more dominant, thus
$\eta_{c}$ tends towards its Poisson limit. In the second scenario, we find
$\eta_{c}=0.207$, significantly lower than the fiducial analysis. In this
case, we have both a sample of almost twice greater density, and one that is
more biased with respect to the continuous dark matter density (such that
$\hat{\rho}_{\rm galaxy}/\hat{\rho}_{\rm dark\ matter}$ is larger, smoothed on
sufficiently large scales). In both cases, however, we find a similar critical
exponent, ($0.82\pm 0.03$) to the above.

Finally, we consider the effective fractal dimension of the system, $d_{F}$.
As discussed in §II, this may be computed from the dependence of the sample
spanning cluster mass, $M(\eta,L)$ (i.e. its number of constituent particles)
on the simulation boxsize $L$ at the percolation threshold $\eta_{c}$ (38). To
explore this, we repeat the above analysis for the fiducial sample, computing
the mass of the sample spanning cluster (when it exists) for five boxsizes in
the range $[700,800]h^{-1}\mathrm{Mpc}$ and five reduced densities in the
range $\eta_{c}\pm 0.1$. For the Poisson system, fitting for the relationship
$M(\eta,L)\sim L^{d}$ and interpolating to $\eta_{c}$ gives $d_{F}\equiv
d(\eta_{c})\approx 2.40\pm 0.07$, matching that predicted from theory [e.g.,
1]. For the galaxy sample, we find $d_{F}=2.36\pm 0.08$, which is consistent
with the Poisson realizations, even though the percolation threshold differs.
This is an important result: the cosmological sample lies in the same
universality class as simple Poisson realizations, for the range of scales
considered: $r\sim 500-1000h^{-1}\mathrm{Mpc}$. This is broadly consistent
with previous results on smaller scales; [97, 96, 118] describes a variety of
methods to ascertain the effective fractal index, with galaxy counts yielding
$d_{F}=2.2\pm 0.2$ on $\lesssim 10h^{-1}\mathrm{Mpc}$ scales, and correlation
functions finding the same on $\lesssim 100h^{-1}\mathrm{Mpc}$ scales.

##  VI The Pair-Connectedness Function as a Cosmological Descriptor

###  VI.1 Background

A crucial problem in modern-day cosmology is the extraction of physical
parameters from observed statistics, such as the distribution of galaxies. In
the standard paradigm (dubbed $\nu\Lambda$CDM), six parameters are of
relevance: (1) the Universe’s current expansion rate, $H_{0}$, (2) the density
of baryonic matter, $\omega_{b}$, (3) the combined density of dark matter and
baryonic matter, $\Omega_{m}$, (4) the amplitude of clustering in the
Universe, $\sigma_{8}$, (5) the slope of the primordial power spectrum (i.e.
structure factor), $n_{s}$, (6) the sum of the neutrino masses, $\sum
m_{\nu}$. Whilst $n_{s}$ and $\Omega_{b}$ are well constrained by observations
of the cosmic microwave background [e.g., 119], the remaining parameters are a
key target for upcoming galaxy surveys. Traditionally, they are constrained
through summary statistics such as the two-particle and three-particle
correlation functions. Such an analysis proceeds by the fitting measured
statistics to analytic models depending on the above physical parameters [e.g.
100].

In this section, we consider the utility of alternative statistics (described
above) in this effort. Three metrics by which we judge a statistic to be
useful are: (a) ease of computation, (b) dimensionality, (c) information
content. Here, we will principally concentrate on the pair-connectedness
function, $P_{2}(r)$, since this has not been previously been used in
cosmology, unlike void probability or nearest-neighbor functions [79]. As
shown above, the statistic is simple and fast to measure from the data, and
has a low-dimensional form, satisfying two of the above criteria. We now
proceed to quantify its ability to constrain cosmological parameters.

An alternative approach to the above prescription is to model the entire
galaxy distribution directly (without compressing to statistics such as the
correlation functions), either with perturbative methods [e.g., 120, 121, 122,
123] or machine learning approaches [e.g., 124, 125, 126, 127, 128, 129]. In
principle, this approach enables one to obtain optimal constraints on all
parameters of interest, though is non-trivial to implement in practice, due to
the huge dimensionality of the galaxy distribution and the necessity to run a
large number of expensive simulations.

###  VI.2 Quantifying Information Content

Standard cosmological analyses proceed by measuring a set of statistics from a
dataset, then comparing them to accurate physical models depending on
cosmological parameters, including those discussed above. If the noise
properties of the statistics are known (for example, if we assert that the
distribution of $g_{2}$ is a multivariate Gaussian), this comparison can be
used to place constraints on the underlying parameters via Bayes theorem. A
useful estimate of the constraining power of some statistic $X$ (e.g.,
$g_{2}$) can be obtained using a Fisher matrix [130], defined as

$\displaystyle
F^{X}_{\alpha\beta}=\left(\frac{\mathrm{d}X}{\mathrm{d}\theta_{\alpha}}\right)^{T}\mathsf{C}_{X}^{-1}\left(\frac{\mathrm{d}X}{\mathrm{d}\theta_{\beta}}\right),$
(43)

where $\\{\theta_{\alpha}\\}$ are the set of cosmological parameters of
interest, and $\mathsf{C}$ is the covariance matrix of $X$ (treated as a
vector), i.e. $\mathsf{C}_{X}=\mathbb{E}\left[X\,X^{T}\right]$, averaging over
realizations of the underlying microscopic density at fixed $\theta$.
According to the Cramér-Rao theorem, $\left(F^{-1,X}\right)_{\alpha\alpha}$
gives the best possible constraint on $\theta_{\alpha}$ from a measurement
$X$, i.e.
$\mathrm{var}(\theta_{\alpha})\geq\left(F^{-1,X}\right)_{\alpha\alpha}$.888This
limit is saturated if $X$ obeys Gaussian statistics, i.e.
$\hat{X}\sim\mathcal{N}\left(X(\theta),\mathsf{C}\right)$. To assess the
utility of statistics such as $P_{2}$ and $g_{2}$, we need simply compute the
covariance matrix and the parameter derivatives appearing in (43), both of
which can be done using a set of simulations. Explicitly, given a set of $n$
realizations $\hat{X}^{(i)}$ with varying initial conditions, the two can be
computed via

$\displaystyle\widehat{\mathsf{C}}_{X}$ $\displaystyle=$
$\displaystyle\frac{1}{n-1}\sum_{i=1}^{n}\left(\hat{X}^{(i)}-\overline{X}\right)\left(\hat{X}^{(i)}-\overline{X}\right)^{T},$
(44) $\displaystyle\widehat{\frac{\mathrm{d}X}{\mathrm{d}\theta_{\alpha}}}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{X}^{(i)}(\theta_{\alpha}+\delta\theta_{\alpha})-\hat{X}^{(i)}(\theta_{\alpha}-\delta\theta_{\alpha})}{2\,\delta\theta_{\alpha}},$

using finite-difference for the parameter derivatives, and denoting
$\bar{X}=(1/n)\sum_{i=1}^{n}\hat{X}^{(i)}$ [e.g., 102].

Whilst the Fisher forecast is appealing in its simplicity, it is not without
limitations. Firstly, it gives accurate bounds on cosmological properties only
if the statistics are Gaussian distributed, which can break down in the case
of large correlations between bins, and the parameter posterior is Gaussian,
which fails for non-negative parameters, for example. Secondly, a large number
of simulations may be required to compute the quantities in (44), and, if too
few are used, the constraining power of a given statistic will be
overestimated.999This occurs since noise in the parameter derivatives add a
positive definite contribution to the Fisher matrix, and thus reduce the size
of the inverted matrix, i.e. the output parameter variances. An alternative
approach is to use simulation-based inference (also known as ‘likelihood-free
analysis’) [e.g., 131, 132, 133]. In essence, this draws a set of cosmological
parameters $\theta$ from some input prior, computes a realization
$\hat{X}(\theta)$ for each, and compares a ‘true’ data-set to the empirical
distribution from the simulations. This does not make assumptions on the
statistics’ noise properties, and, in the case of too few simulations, will
only underestimate the cosmological utility.

Here, we examine the constraining power of various summary statistics using
both the Fisher matrix formalism (which has become commonplace in cosmology)
and simulation based inference (which is far less common, though more
accurate). In particular, we consider the cosmological parameters $H_{0}$,
$\Omega_{m}$, $\sigma_{8}$ and $\sum m_{\nu}$ and the following descriptors:
$g_{2}$, $g_{3}$, $P_{2}$, and $H_{V}$, all of which can be defined for
discrete point clouds, such as the galaxy density used in this work. For the
pair-connectedness function, we fix the reduced density to
$\eta\equiv\bar{\rho}\pi D^{3}/6=0.2$, which is a useful balance between the
uninformative case ($\eta=0$) and the percolated limit discussed in §V, though
we note that other choices may yield somewhat different results. Additional
statistics could be straightforwardly added, though we caution that
descriptors such as the number variance are fully described by $g_{2}$, and
will thus not add additional information. Using the fastpm code [134], we run
$n=512$ Quijote-like dark-matter simulations with the following fiducial
parameters:
$\\{H_{0}=68\,\mathrm{km}\,\mathrm{s}^{-1}\mathrm{Mpc}^{-1},\Omega_{m}=0.31,\omega_{b}=0.0227,n_{s}=0.96,\sigma_{8}=0.8178,\sum
m_{\nu}=0\,\mathrm{eV}\\}$, and compute the four statistics for each
realization.

In all cases, we consider only scales above $20h^{-1}\mathrm{Mpc}$ (where the
simulations are accurate, given the mean pairwise separation of $\approx
8h^{-1}\mathrm{Mpc}$), and choose the radial bin sizes to keep the
dimensionality fixed to $\mathcal{O}(100)$ elements. $g_{2}$, $P_{2}$, and
$H_{V}$ are computed as before, with $g_{3}$ computed using the approach of
[135], involving a decomposition into a Legendre multipole basis, with
components $g_{3,\ell}(r_{1},r_{2})$ for $\ell\in\\{0,1,2,3\\}$. In accordance
with §III, we compute this statistic in configuration-space (rather than as a
Fourier-space bispectrum), which obviates the need to grid the particles.
These simulations are used to compute the covariance matrix $\mathsf{C}$ of
the statistics using (44),101010Initial testing demonstrated that this is a
sufficient number of simulations to estimate $\mathsf{C}$ and its inverse
robustly, after including the correction factor of [136]. the structure of
which is visualized in Fig. 11. We find significant correlations both within
and between a number of statistics. In particular, the individual bins of
$g_{3}$ and $P_{2}$ are highly correlated, indicating that their noise
properties may not be Gaussian. In contrast, the void nearest-neighbor
function has an almost diagonal correlation matrix, and is seen to be largely
independent from other statistics. This suggests that it can add significant
information compared to analyses using $g_{2}$-alone.

Figure 11: Correlation matrix of the pair correlation function ($g_{2}$), the
pair-connectedness function ($P_{2}$, with reduced density $\eta=0.2$), the
three-particle correlation function ($g_{3}$) and the void nearest-neighbor
function ($H_{V}$) obtained from $512$ cosmological simulations run with a
fiducial set of cosmological parameters. The correlation matrix is defined as
the covariance matrix normalized by its leading diagonal, i.e.
$C_{ij}/\sqrt{C_{ii}C_{jj}}$, with all values lying in the range $[-1,1]$. The
labels show the statistic of interest, which are demarcated by the dotted
lines. In each statistic, the bins are ordered from small-scale (bottom left)
to large-scale (top-right). For $g_{3}$, we use four Legendre multipoles,
$g_{3,\ell}$ with $\ell\in\\{0,1,2,3\\}$. Notably, the individual bins of
$g_{3,0}$ and $P_{2}$ are highly correlated, and there are a number of non-
trivial correlations between observables, though few with $H_{V}$.

The other ingredient required for Fisher forecasting is the set of parameter
derivatives. These are computed using (44), with $n=512$ simulations (again
computed using fastpm, with a total cost of $\sim 10^{4}$ CPU-hours),
utilizing finite difference in each of the eight sets of parameters. For the
neutrino mass, we have the bound $\sum m_{\nu}>0$, thus we instead utilize
one-sided derivatives, following [102], and using the method of [137] to
emulate the effects of massive neutrinos by modifying the initial conditions.
Following this, we compute the Fisher matrix via (43) for various combinations
of statistics. We caution that this result appears to retain some dependence
on $n$ due to residual noise in the parameter derivatives. This will lead to
the constraints being artificially tightened somewhat; however, it is
computationally impractical to increase $n$ by a significant amount.

For the simulation-based inference (hereafter SBI), we utilize a set of $8192$
galaxy simulations computed using fastpm with the method of [137] at random
locations in parameter space, according to the flat priors:
$H_{0}\in[0,100]\mathrm{km}\,\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$,
$\Omega_{m}\in[0,0.3]$, $\sigma_{8}\in[0.4,1.2]$, $\sum
m_{\nu}\in[0,4]\mathrm{eV}$.111111The neutrino mass limit is significantly
weaker than the bound from the latest probes [119], but is appropriate given
the small volume of the simulations. Summary statistics for each are computed
as before, and fed into the sbi code, which uses neural networks (via the
‘Sequential Neural Posterior Estimation’ method) to compute the parameter
posterior, given a ‘true’ observation, taken from the mean of the fiducial
simulations discussed above.

Figure 12: Forecasted constraints on key cosmological parameters using the
pair correlation function $g_{2}$, and the pair-connectedness function,
$P_{2}$, using a Fisher forecast (left) and simulation-based inference
(right). The dark (light) ellipses represent our $68\%$ ($95\%$) confidence
intervals on the various parameters possible from a single galaxy simulation
(as in Fig. 2), with the grey regions show the prior used in the simulation-
based approach (from 8192 simulations). We show results for four parameters:
$H_{0}$, giving the Universe’s expansion rate, the matter density
$\Omega_{m}$, the clustering amplitude $\sigma_{8}$ and the total neutrino
mass $\sum m_{\nu}$. The diagonal figures give the marginalized constraints on
single parameters, whilst the off-diagonals show the correlations, e.g., the
bottom left panel shows the correlation between $h$ and $\sigma_{8}$. We find
that $P_{2}$ is a poor predictor of $H_{0}$ and $\Omega_{m}$ (shown by the
broad contours), but can tightly constrain $\sigma_{8}$. The combination
$g_{2}+P_{2}$ significantly increases the precision by which this parameter
can be measured, and we find similar results from the Fisher forecasts (which
are, in general, overoptimistic), and the simulation-based inference (which is
usually conservative).

###  VI.3 Results

Figure 12 shows the constraints on cosmological parameters from analyses using
the pair correlation function and pair-connectedness function, both via the
Fisher and SBI forecasts. From the Fisher forecast, we observe that $P_{2}$ is
a poor predictor of the expansion rate and matter density, but gives much
tighter constraints on the clustering amplitude and neutrino mass than
$g_{2}$. This is unsurprising: the slope of $P_{2}(r)$ is strongly sensitive
to the galaxy clustering properties (set by $\sigma_{8}$, and, on small
scales, $\sum m_{\nu}$), but (being a monotonic function) contains little
information on other properties such as early Universe physics. The SBI
forecasts give qualitatively similar results, with the $P_{2}$ constraints on
$H_{0}$ and $\Omega_{m}$ being largely dominated by the prior, with the
combined $g_{2}+P_{2}$ constraints reproducing those of $g_{2}$ alone, For
$\sigma_{8}$, $P_{2}$ is again shown to be of considerable use, with a
significant (factor of a $\approx 2.6\times$, equivalently to observing a
seven times greater volume of space) tightening in the one-dimensional
posterior found by adding $P_{2}$, driven by the differing degeneracy
directions in the $\sigma_{8}-H_{0}$ and $\sigma_{8}-\Omega_{m}$ planes. In
contrary to the Fisher result, the SBI forecast suggests that the pair-
connectedness function does not give significant additional information on the
neutrino mass, however, it is shown to change the $\sigma_{8}-\sum m_{\nu}$
degeneracy direction considerably.

The disagreement between Fisher and SBI forecasts both quantitatively (in
terms of the reduction in width of the $\sigma_{8}$ posterior) and
qualitatively (whether the $\sum m_{\nu}$ posterior is affected) may appear a
little unsettling. We attribute this to a number of reasons: (1) as mentioned
above, the Fisher forecast will give artificially narrow constraints if
insufficient simulations have been run, (2) the Fisher forecast is inaccurate
for parameters whose posterior is non-Gaussian (such as the neutrino mass, due
to the $\sum m_{\nu}>0$ constraint), (3) the SBI forecasts can be artificially
broadened by insufficient simulations being run. However, the results of Fig.
12 are enough to convince us that $P_{2}$ contains significant information
regarding the clustering amplitude $\sigma_{8}$, and its inclusion greatly
aids cosmological analyses, including via degeneracy breaking with $\sum
m_{\nu}$. Although we present results only for reduced density $\eta=0.2$
here, a similar story holds also for $\eta=0.1$; in this case, the
improvements in cosmological parameters are somewhat weaker, due to the
higher-order correlator contributions to $P_{2}(r,\eta)$ being suppressed
(35). We expect that combining measurements of the pair-connectedness function
with multiple values of $\eta$ could further increase the constraining power,
again at little computational cost.

Figure 13: As Fig. 12 but comparing the constraining power from the
combination of the pair correlation function with the pair-connectedness
function (green), the three-particle correlation function (red) and the void
nearest-neighbor function (yellow). The pair-connectedness function and void
nearest-neighbor functions are found to perform almost as well as the three-
particle function in this scenario (in terms of constraining $\sigma_{8}$),
and are much faster to compute and analyze.

It is interesting to compare the cosmological utility of the pair-
connectedness function to that of other higher-order statistics. Before doing
so, let us briefly outline our predictions. In this test, we are limited to
relatively large scales ($r\gtrsim 20h^{-1}\mathrm{Mpc}$, due to simulation
resolution effects), where the galaxy distribution (if treated as a continuous
field) is close to Gaussian. As such, we expect the majority of the
information content on cosmological parameters to be encapsulated by the two-
and three-point functions, $g_{2}$ and $g_{3}$, with only a small amount
leaking into higher-order statistics. In this case, the combination of $g_{2}$
with alternative statistics will likely perform worse than that of $g_{2}$ and
$g_{3}$; our question is whether there are statistics that are able to recoup
most of the information present in $g_{3}$ in a simpler form (for example in
the unidimensional $H_{V}$ and $P_{2}$ statistics). If such a statistic
exists, it is likely that it also contains significant information on small
scales (as probed by future surveys such as that of the Subaru Prime Focus
Spectrograph and the DESI Bright Galaxy Survey [138, 139]), where the
perturbative hierarchy described above breaks down. Indeed, statistical
physics provides examples of small-scale systems where $g_{2}+P_{2}$
outperforms $g_{2}+g_{3}$ [e.g., 90]; it will be interesting to study such
effects further in the future.

In Fig. 13 we show the Fisher and SBI constraints on the same parameter set as
above for $g_{2}$ in combination with $g_{3}$, $P_{2}$, and $H_{V}$.121212We
recall that our $H_{V}$ void statistic follows a somewhat different definition
to the void size function often used in cosmology, and is restricted to
comparatively large scales ($r>20h^{-1}\mathrm{Mpc}$, as for the other
statistics), with only spherical voids. This differs from the approach used in
several cosmological studies [e.g., 11, 140, 141] and explains the reduced
utility found herein, and the different correlation properties seen in Fig.
11. From both the Fisher and SBI forecasts, we find that no additional
statistics lead to significant improvements in the expansion rate constraints,
except for a slight tightening from $g_{3}$. This is not surprising: $H_{0}$
is primarily measured from an oscillatory feature in $g_{2}$ arising from
acoustic waves in the early Universe, which is generally absent in other
statistics. In the SBI forecast, the same is true for the matter density,
$\Omega_{m}$ and the neutrino mass, though the Fisher forecasts disagree on
this aspect, as above, and should therefore be taken with a grain of salt
(especially given the larger dimensionality of $g_{3}$). For the clustering
amplitude, $\sigma_{8}$, we find similar improvement when combining $g_{2}$
with any other statistic, with a slight preference for $g_{3}$ in the SBI
analysis (or a significant one for the Fisher forecast). This matches the
above predictions.

Our conclusion from this exercise is the following: if one wishes to constrain
the Universe’s clustering amplitude (a key target of modern-day cosmology),
the addition of $P_{2}$ or $H_{V}$ into cosmological analyses provides an
excellent route [cf., 76], and contains similar information to $g_{3}$.
Importantly, the alternative statistics are of much lower dimension than
$g_{3}$ and $P_{2}$ is much less computationally expensive to measure
(requiring $\sim 5$ CPU-minutes per simulation, instead of $\sim 1$ CPU-hour
for $g_{3}$ or $H_{v}$). Furthermore, if performs the analysis using a
combination of pair-connectedness functions with different values of $\eta$,
the results may be stronger still. Whilst this analysis is necessarily
simplistic and limited to comparatively large scales (due to the nature of the
simulation suite), it nevertheless suggests that the pair-connectedness
function is a new statistic of significant potential, and could carry
important information also on small scales. In contrast to correlators such as
$g_{2}$ and $g_{3}$ [e.g., 142], this is difficult to model analytically even
at large $r$, due to its inherent dependence on short-scale physics including
the connection between galaxies and dark matter. For this reason it will
likely prove useful to adopt a simulation-based methodology to analyze
$P_{2}$, such as the SBI techniques discussed above, and marginalize over
parameters controlling galaxy formation. On small scales, a similar approach
is required for any statistic, due to the breakdown of perturbative modeling.

##  VII Summary

In this work, we have considered the application of the theory of disordered
heterogeneous media and statistical mechanics to cosmology, and of cosmology
to the former. By treating the distribution of galaxies in the present-day
Universe as a point process, we can analyze the data using techniques
developed to characterize heterogeneous media, such as the correlation
functions and nearest-neighbor distributions. Furthermore, by augmenting the
data-set with some concept of ‘connectedness’ (here defined by circumscribing
the galaxies with spheres), we may utilize various clustering diagnostics and
pair-connectedness functions, which encode a different subset of the
information present within the distribution and additionally allows the
percolation properties to be determined. Such a framework (a) provides a novel
method for understanding the galaxy distribution, whose importance will only
grow in the next decade with the plurality of upcoming telescopes, and (b)
demonstrates the applicability of heterogeneous media and statistical
mechanical techniques in a very different regime to that usually explored.

Our main conclusions are the following:

  * •

The galaxy distribution exhibits very different physical properties to those
of conventional materials, leading to distinct signatures in a wide variety of
clustering and correlation descriptors. On the largest scales, the system
approaches hyperuniformity, whilst on the the smallest, it becomes almost
antihyperuniform and strongly inhomogeneous; this dichotomy arises from the
fact that the minimum separation between galaxies is much smaller than the
mean interparticle distance, with localized groups of galaxies separated by
vast cosmic distance.

  * •

Physically, the cosmological system has two peculiarities: (a) although we
treat the galaxies as point objects, they have some physical scale in
practice, and cannot overlap, (b) the distribution carries the signatures of a
large-scale stochastic background that modulates the quasi-Poissonian
distribution; this is sourced by early Universe physics and gravitational
evolution,

  * •

The galaxy pair correlation function shows this behavior clearly, with the
expected hyperuniform tail appearing only at gargantuan scales ($r\gtrsim
200h^{-1}\mathrm{Mpc}$), and with a sharp peak at the mean pairwise particle
separation of $r\approx 8h^{-1}\mathrm{Mpc}$. This scale separation induces a
large number variance, which is unusually super-Poissonian on small scales,
yet sub-Poissonian on the largest. These results are consistent with the order
metric, $\tau$, which we determine for the galaxy sample for the first time:
its value ($\tau=4.85$) implies that the system is strongly correlated and
disordered. The nearest-neighbor functions are again consistent with this
picture, with considerably extended tails, a lack of particle pairs below some
critical galaxy size, and much enhanced variance relative to the Poissonian
case (and most other common scenarios).

  * •

Analysis of pair-connectedness functions, mean particle numbers, and sample
spanning clusters indicate that the galaxy sample percolates at significantly
lower reduced densities than corresponding Poisson realizations. For the
fiducial galaxy simulations, finite-scaling analysis gives $\eta_{c}=0.25$ in
the former case compared to $\eta_{c}=0.34$ in the latter, a difference which
is amplified by increasing the sample density. This is again supported by the
above evidence: the scale separation is a consequence of the extra small-scale
clustering in the galaxy distribution which leads to faster percolation. Both
scenarios appear to have the same critical exponents and fractal dimensions,
implying that they live in the same universality class, despite very different
physics operating.

  * •

The pair-connectedness function is a conceptually straightforward and easy-to-
measure statistic that carries useful and accessible large-scale information
about the underlying physical parameters of the Universe, and can be trivially
extended to small scales. This could enhance the cosmological utility of
future galaxy surveys, in combination with conventional techniques. Using
simulation-based analysis techniques, we forecast that constraints on
amplitude of clustering improve by a factor of $\approx 5$ (or $\approx 25$ in
terms of survey volume) when performing inference using the large-scale pair-
connectedness and pair correlation functions as opposed to the pair
correlation function alone (which is standard in cosmology). This provides a
useful alternative to the three-particle correlation function $g_{3}$, which
is of significantly lower dimension and much faster to model, and is shown to
be a resummation of correlation functions of all order. Unlike the large-scale
three-particle function, it seems unlikely that $P_{2}(r)$ can be modelled
analytically; simulation-based treatments will likely be required in this
case.

The galaxy samples used in this work are purposefully simplified, in order to
provide a proof-of-concept study capturing the essential physical attributes
of the cosmological set-up. More work is needed before the statistics can be
applied to real data, and will require the following: (a) higher resolution
simulations containing more particles, allowing smaller scales to be probed,
(b) inclusion of real galaxies in the simulations, rather than dark matter
halos, and the associated physical uncertainties with their formation [e.g.,
143, 144], (c) anisotropic distortions in the Universe created by transforming
from redshifts to physical coordinates [145], (d) inhomogeneities in the field
induced by observational effects, such as the limited field-of-view of the
telescope. However, all of these complexities have been overcome a number of
times before for other statistics (such as the correlation functions [e.g.,
100]) and we expect can be similarly surmounted in this case. Furthermore, it
is important to characterize how the statistics depend on the galaxy sample:
for the correlation functions, this is well understood (and encapsulated by
‘bias parameters’, which depend on galaxy mass and luminosity), but should be
explored further for nearest-neighbor and pair-connectedness functions, as
well as the percolation threshold.

Finally, we consider the broader extensions of this work. Although we have
restricted our gaze to galaxy distributions, this is far from being the only
stochastic distribution in the Universe. One additional application could
include a more principled treatment of cosmic voids [e.g., 11, 10]: these are
low-density regions in the galaxy distribution that form a partition of the
space, and could be described by the same mathematics as that invoked for
percolation. Even more relevant is the distribution of ‘bubbles’ of ionized
gas around the first galaxies [e.g., 146, 147]. The growth of such bubbles
likely led to the Universe’s reionization approximately one billion years
after the Big Bang, the time of which is set by percolation itself. Finally,
we note that there are a wealth of techniques from the theory of disordered
heterogeneous media that have not been considered in this work. It would be
interesting to consider the utility of the various descriptors in the context
of ‘simulated annealing’ [e.g., 148, 90], to understand the extent to which
any statistic can capture the full complexities of the field, though we
caution that conventional approaches will likely need to be modified to
account for the peculiarities of the galaxy distribution, in particular its
significant scale separation. Further still, we may consider how annealing
techniques allow us to recover ‘effective pair interactions’ between
individual galaxies [149], and thus learn more about the Universe’s average
dynamics.

###### Acknowledgements.

We thank Jim Peebles, Paul Steinhardt, Robert Scherrer, Arka Banerjee, Tom
Abel, and Haina Wang for insightful comments on this manuscript. OHEP is
additionally grateful to Alice Pisani and Will Coulton for useful discussions
regarding cosmic voids and Fisher forecasts respectively. We additionally
thank the anonymous referees for an insightful report. OHEP is a Junior Fellow
of the Simons Society of Fellows and thanks the Institute for Advanced Study
for their hospitality and abundance of baked goods. ST thanks the Institute
for Advanced Study for their hospitality during his sabbatical leave there.
The authors are pleased to acknowledge that the work reported in this paper
was substantially performed using the Princeton Research Computing resources
at Princeton University, which is a consortium of groups led by the Princeton
Institute for Computational Science and Engineering (PICSciE) and the Office
of Information Technology’s Research Computing Division. Additional
computations were performed on the Helios cluster at the Institute for
Advanced Study, Princeton.

## References

  * [1] S. Torquato, _Random Heterogeneous Materials: Microstructure and Macroscopic Properties_. Springer-Verlag, New York, 2002.
  * [2] M. Sahimi, _Heterogeneous Materials I: Linear Transport and Optical Properties_. Springer-Verlag, New York, 2003.
  * [3] J. P. Hansen and I. R. McDonald, _Theory of Simple Liquids_. Academic Press, New York, 4th ed., 2013.
  * [4] P. J. E. Peebles, _The large-scale structure of the universe_. Princeton University Press, 1980.
  * [5] W. C. Saslaw, _The Distribution of the Galaxies_. Cambridge University Press, Cambridge, England, 2000.
  * [6] A. Gabrielli, F. S. Labini, M. Joyce and L. Pietronero, _Statistical Physics for Cosmic Structures_. Springer-Verlag, New York, 2005.
  * [7] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _arXiv e-prints_ (2016) arXiv:1611.00036 [1611.00036].
  * [8] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _arXiv e-prints_ (2011) arXiv:1110.3193 [1110.3193].
  * [9] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, _Large scale structure of the universe and cosmological perturbation theory_ , _Phys. Rept._ 367 (2002) 1 [astro-ph/0112551].
  * [10] R. K. Sheth and R. van de Weygaert, _A Hierarchy of voids: Much ado about nothing_ , _Mon. Not. Roy. Astron. Soc._ 350 (2004) 517 [astro-ph/0311260].
  * [11] A. Pisani et al., _Cosmic voids: a novel probe to shed light on our Universe_ , 1903.05161.
  * [12] R. K. Sheth, _The halo-model description of marked statistics_ , _MNRAS_ 364 (2005) 796 [astro-ph/0511772].
  * [13] M. White, _A marked correlation function for constraining modified gravity models_ , _JCAP_ 2016 (2016) 057 [1609.08632].
  * [14] E. Massara, F. Villaescusa-Navarro, S. Ho, N. Dalal and D. N. Spergel, _Using the Marked Power Spectrum to Detect the Signature of Neutrinos in Large-Scale Structure_ , _Phys. Rev. Lett._ 126 (2021) 011301 [2001.11024].
  * [15] O. H. E. Philcox, E. Massara and D. N. Spergel, _What does the marked power spectrum measure? Insights from perturbation theory_ , _Phys. Rev. D_ 102 (2020) 043516 [2006.10055].
  * [16] D. H. Weinberg, _Reconstructing primordial density fluctuations. I - Method_ , _MNRAS_ 254 (1992) 315.
  * [17] M. C. Neyrinck, I. Szapudi and A. S. Szalay, _Rejuvenating the Matter Power Spectrum: Restoring Information with a Logarithmic Density Mapping_ , _ApJ_ 698 (2009) L90 [0903.4693].
  * [18] X. Wang, M. Neyrinck, I. Szapudi, A. Szalay, X. Chen, J. Lesgourgues et al., _Perturbation Theory of the Cosmological Log-density Field_ , _ApJ_ 735 (2011) 32 [1103.2166].
  * [19] H. Rubira and R. Voivodic, _The effective field theory and perturbative analysis for log-density fields_ , _JCAP_ 2021 (2021) 070 [2011.12280].
  * [20] D. J. Eisenstein, H.-j. Seo, E. Sirko and D. Spergel, _Improving Cosmological Distance Measurements by Reconstruction of the Baryon Acoustic Peak_ , _Astrophys. J._ 664 (2007) 675 [astro-ph/0604362].
  * [21] F. Elsner, F. Schmidt, J. Jasche, G. Lavaux and N.-M. Nguyen, _Cosmology inference from a biased density field using the EFT-based likelihood_ , _JCAP_ 2020 (2020) 029 [1906.07143].
  * [22] G. Cabass and F. Schmidt, _The EFT likelihood for large-scale structure_ , _JCAP_ 2020 (2020) 042 [1909.04022].
  * [23] F. Schmidt, _Sigma-eight at the percent level: the EFT likelihood in real space_ , _JCAP_ 2021 (2021) 032 [2009.14176].
  * [24] I. Gott, J. Richard, D. H. Weinberg and A. L. Melott, _A Quantitative Approach to the Topology of Large-Scale Structure_ , _ApJ_ 319 (1987) 1.
  * [25] T. Matsubara, _Analytic expression of the genus in weakly non-gaussian field induced by gravity_ , _Astrophys. J. Lett._ 434 (1994) L43 [astro-ph/9405037].
  * [26] T. Matsubara, _Nonlinear Evolution of Topology of Large Scale Structure_ , in _Cosmological Constant and the Evolution of the Universe_ , K. Sato, T. Suginohara and N. Sugiyama, eds., p. 45, Jan., 1996\.
  * [27] J. Schmalzing, _Minkowski functionals in cosmology: An Overview_ , in _12th Potsdam Cosmology Workshop: Large Scale Structure: Tracks and Traces: International Workshop_ , 9, 1997, astro-ph/9710302.
  * [28] SDSS collaboration, _Minkowski Functionals of SDSS galaxies. 1\. Analysis of excursion sets_ , _Publ. Astron. Soc. Jap._ 55 (2003) 911 [astro-ph/0304455].
  * [29] P. Parihar, M. S. Vogeley, I. Gott, J. Richard, Y.-Y. Choi, J. Kim, S. S. Kim et al., _A Topological Analysis of Large-Scale Structure, Studied Using the CMASS Sample of SDSS-III_ , _ApJ_ 796 (2014) 86.
  * [30] M. Biagetti, J. Calles, L. Castiblanco, A. Cole and J. Noreña, _Fisher Forecasts for Primordial non-Gaussianity from Persistent Homology_ , _arXiv e-prints_ (2022) arXiv:2203.08262 [2203.08262].
  * [31] R. Meester, R. Roy and A. Sarkar, _Nonuniversality and continuity of the critical covered volume fraction in continuum percolation_ , _J. Stat. Phys._ 75 (1994) 123.
  * [32] K. R. Mecke, T. Buchert and H. Wagner, _Robust morphological measures for large-scale structure in the universe_ , _Astronomy Astrophys._ 288 (1994) 697.
  * [33] T. Sousbie, C. Pichon, S. Colombi, D. Novikov and D. Pogosyan, _The three dimensional skeleton: Tracing the filamentary structure of the Universe_ , _Mon. Not. Roy. Astron. Soc._ 383 (2008) 1655 [0707.3123].
  * [34] S. Torquato and F. H. Stillinger, _Local density fluctuations, hyperuniform systems, and order metrics_ , _Phys. Rev. E_ 68 (2003) 041113.
  * [35] S. Torquato, _Hyperuniform states of matter_ , _Physics Reports_ 745 (2018) 1.
  * [36] A. Gabrielli, M. Joyce and F. S. Labini, _Glass-like universe: Real-space correlation properties of standard cosmological models_ , _Phys. Rev. D_ 65 (2002) 083523.
  * [37] S. Torquato and F. H. Stillinger, _New conjectural lower bounds on the optimal density of sphere packings_ , _Experimental Math._ 15 (2006) 307.
  * [38] C. E. Zachary and S. Torquato, _Hyperuniformity in point patterns and two-phase heterogeneous media_ , _J. Stat. Mech.: Theory & Exp._ 2009 (2009) P12015.
  * [39] C. Lin, P. J. Steinhardt and S. Torquato, _Hyperuniformity variation with quasicrystal local isomorphism class_ , _J. Phys.: Cond. Matter_ 29 (2017) 204003.
  * [40] E. C. Oğuz, J. E. S. Socolar, P. J. Steinhardt and S. Torquato, _Hyperuniformity of quasicrystals_ , _Phys. Rev. B_ 95 (2017) 054119.
  * [41] A. Gabrielli, _Point processes and stochastic displacement fields_ , _Phys. Rev. E_ 70 (2004) 066131.
  * [42] A. Gabrielli and S. Torquato, _Voronoi and void statistics for superhomogeneous point processes_ , _Phys. Rev. E_ 70 (2004) 041105\.
  * [43] A. Gabrielli, M. Joyce and S. Torquato, _Tilings of space and superhomogeneous point processes_ , _Phys. Rev. E_ 77 (2008) 031125\.
  * [44] J. Kim and S. Torquato, _Effect of imperfections on the hyperuniformity of many-body systems_ , _Phys. Rev. B_ 97 (2018) 054105.
  * [45] O. U. Uche, F. H. Stillinger and S. Torquato, _Constraints on collective density variables: Two dimensions_ , _Phys. Rev. E_ 70 (2004) 046122.
  * [46] S. Torquato, G. Zhang and F. H. Stillinger, _Ensemble theory for stealthy hyperuniform disordered ground states_ , _Phys. Rev. X_ 5 (2015) 021020.
  * [47] G. Zhang, F. H. Stillinger and S. Torquato, _The perfect glass paradigm: Disordered hyperuniform glasses down to absolute zero_ , _Sci. Rep._ 6 (2016) 36963.
  * [48] Q.-L. Lei and R. Ni, _Hydrodynamics of random-organizing hyperuniform fluids_ , _Proc. Nat. Acad. Sci._ 116 (2019) 22983.
  * [49] S. Torquato, G. Zhang and M. de Courcy-Ireland, _Hidden multiscale order in the primes_ , _J. Phys. A: Math. & Theoretical_ 52 (2019) 135002.
  * [50] A. Donev, F. H. Stillinger and S. Torquato, _Unexpected density fluctuations in disordered jammed hard-sphere packings_ , _Phys. Rev. Lett._ 95 (2005) 090604.
  * [51] C. E. Zachary, Y. Jiao and S. Torquato, _Hyperuniform long-range correlations are a signature of disordered jammed hard-particle packings_ , _Phys. Rev. Lett._ 106 (2011) 178001.
  * [52] Y. Jiao and S. Torquato, _Maximally random jammed packings of Platonic solids: Hyperuniform long-range correlations and isostaticity_ , _Phys. Rev. E_ 84 (2011) 041309.
  * [53] S. Atkinson, G. Zhang, A. B. Hopkins and S. Torquato, _Critical slowing down and hyperuniformity on approach to jamming_ , _Phys. Rev. E_ 94 (2016) 012902.
  * [54] R. P. Feynman and M. Cohen, _Energy spectrum of the excitations in liquid helium_ , _Phys. Rev._ 102 (1956) 1189.
  * [55] L. Reatto and G. V. Chester, _Phonons and the properties of a Bose system_ , _Phys. Rev._ 155 (1967) 88.
  * [56] S. Torquato, A. Scardicchio and C. E. Zachary, _Point processes in arbitrary dimension from Fermionic gases, random matrix theory, and number theory_ , _J. Stat. Mech.: Theory Exp._ 2008 (2008) P11019.
  * [57] C. E. Zachary and S. Torquato, _Anomalous local coordination, density fluctuations, and void statistics in disordered hyperuniform many-particle ground states_ , _Phys. Rev. E_ 83 (2011) 051133.
  * [58] D. Hexner and D. Levine, _Hyperuniformity of critical absorbing states_ , _Phys. Rev. Lett._ 114 (2015) 110602.
  * [59] S. Torquato and J. Kim, _Nonlocal effective electromagnetic wave characteristics of composite media: Beyond the quasistatic regime_ , _Phys. Rev. X_ 11 (2021) 021002.
  * [60] B. Widom, _Equation of state in the neighborhood of the critical point_ , _J. Chem. Phys._ 43 (1965) 3898.
  * [61] L. P. Kadanoff, _Scaling laws for Ising models near $T_{c}$_, _Physics_ 2 (1966) 263.
  * [62] M. E. Fisher, _The theory of equilibrium critical phenomena_ , _Rep. Prog. Phys._ 30 (1967) 615.
  * [63] K. G. Wilson and J. Kogut, _The renormalization group and the $\epsilon$ expansion_, _Phys. Rep._ 12 (1974) 75.
  * [64] J. J. Binney, N. J. Dowrick, A. J. Fisher and M. E. J. Newman, _The Theory of Critical Phenomena: An Introduction to the Renormalization Group_. Oxford University Press, Oxford, England, 1992.
  * [65] E. C. Oğuz, J. E. S. Socolar, P. J. Steinhardt and S. Torquato, _Hyperuniformity and anti-hyperuniformity in one-dimensional substitution tilings_ , _Acta Cryst. Section A: Foundations & Advances_ A75 (2019) 3.
  * [66] S. Torquato, _Perspective: Basic understanding of condensed phases of matter via packing models_ , _J. Chem. Phys._ 149 (2018) 020901\.
  * [67] A. R. Kansal, S. Torquato and F. H. Stillinger, _Diversity of order and densities in jammed hard-particle packings_ , _Phys. Rev. E_ 66 (2002) 041109.
  * [68] S. Atkinson, F. H. Stillinger and S. Torquato, _Static structural signatures of nearly jammed disordered and ordered hard-sphere packings: Direct correlation function_ , _Phys. Rev. E_ 94 (2016) 032902\.
  * [69] M. A. Klatt, J. Lovrić, D. Chen, S. C. Kapfer, F. M. Schaller, P. W. A. Schönhöfer et al., _Universal hidden order in amorphous cellular geometries_ , _Nature Comm._ 10 (2019) 811.
  * [70] F. Martelli, S. Torquato, N. Giovambattista and R. Car, _Large-scale structure and hyperuniformity of amorphous ices_ , _Phys. Rev. Lett._ 119 (2017) 136002.
  * [71] E. Lomba, J.-J. Weis and S. Torquato, _Disordered hyperuniformity in two-component non-additive hard disk plasmas_ , _Phys. Rev. E_ 96 (2017) 062126.
  * [72] S. Torquato, B. Lu and J. Rubinstein, _Nearest-neighbor distribution functions in many-body systems_ , _Phys. Rev. A_ 41 (1990) 2059.
  * [73] B. S. Ryden and E. L. Turner, _A statistical comparison of voids in the galaxy distribution and n-body simulations_ , _ApJ_ 287 (1984) L59.
  * [74] S. D. M. White, _The hierarchy of correlation functions and its relation to other measures of galaxy clustering_ , _Mon. Not. Roy. Astron. Soc._ 186 (1979) 145.
  * [75] R. Balian and R. Schaeffer, _Scale-invariant matter distribution in the universe._ , _A &A_ 220 (1989) 1.
  * [76] A. Banerjee and T. Abel, _Nearest neighbour distributions: New statistical measures for cosmological clustering_ , _MNRAS_ 500 (2021) 5479 [2007.13342].
  * [77] A. Banerjee and T. Abel, _Cosmological cross-correlations and nearest neighbour distributions_ , _MNRAS_ 504 (2021) 2911 [2102.01184].
  * [78] A. Banerjee, N. Kokron and T. Abel, _Modelling nearest neighbour distributions of biased tracers using hybrid effective field theory_ , _MNRAS_ 511 (2022) 2765 [2107.10287].
  * [79] Y. Wang, A. Banerjee and T. Abel, _Detection of spatial clustering in the 1000 richest SDSS DR8 redMaPPer clusters with nearest neighbor distributions_ , _MNRAS_ 514 (2022) 3828 [2112.04502].
  * [80] D. J. Vezzetti, _A new derivation of some fluctuation theorems in statistical mechanics_ , _J. Math. Phys._ 16 (1975) 31.
  * [81] R. M. Ziff, _On the bulk distribution functions and fluctuation theorems_ , _J. Math. Phys._ 18 (1977) 1825.
  * [82] T. M. Truskett, S. Torquato and P. G. Debenedetti, _Density fluctuations in many-body systems_ , _Phys. Rev. E_ 58 (1998) 7639.
  * [83] S. Torquato, J. Kim and M. A. Klatt, _Local number fluctuations in hyperuniform and nonhyperuniform systems: Higher-order moments and distribution functions_ , _Phys. Rev. X_ 11 (2021) 021028.
  * [84] S. Torquato, _Reformulation of the covering and quantizer problems as ground states of interacting particles_ , _Phys. Rev. E_ 82 (2010) 056109.
  * [85] J. H. Conway and N. J. A. Sloane, _Sphere Packings, Lattices and Groups_. Springer-Verlag, New York, 1998.
  * [86] A. Pisani, E. Massara, D. N. Spergel, D. Alonso, T. Baker, Y.-C. Cai et al., _Cosmic voids: a novel probe to shed light on our Universe_ , _BAAS_ 51 (2019) 40 [1903.05161].
  * [87] A. Coniglio, U. D. Angelis and A. Forlani, _Pair connectedness and cluster size_ , _Journal of Physics A: Mathematical and General_ 10 (1977) 1123.
  * [88] G. Stell, _Exact equation for the pair-connectedness function_ , _J. Phys. A: Math. Gen._ 17 (1984) L855.
  * [89] S. Torquato, J. Beasley and Y. Chiew, _Two-point cluster function for continuum percolation_ , _The Journal of chemical physics_ 88 (1988) 6540.
  * [90] Y. Jiao, F. H. Stillinger and S. Torquato, _A superior descriptor of random textures and its predictive capacity_ , _Proceedings of the National Academy of Sciences_ 106 (2009) 17634 [https://www.pnas.org/doi/pdf/10.1073/pnas.0905919106].
  * [91] S. Torquato, _Effect of dimensionality on the continuum percolation of overlapping hyperspheres and hypercubes_ , _The Journal of Chemical Physics_ 136 (2012) 054106 [https://doi.org/10.1063/1.3679861].
  * [92] J. K. Percus and G. J. Yevick, _Analysis of classical statistical mechanics by means of collective coordinates_ , _Phys. Rev._ 110 (1958) 1.
  * [93] H. Hadwiger, _Über Treffanzahlen bei translationsgleichen Eikörpern_ , _Arch. Math._ 8 (1957) 212.
  * [94] S. B. Lee and S. Torquato, _Monte carlo study of correlated continuum percolation: Universality and percolation thresholds_ , _Phys. Rev. A_ 41 (1990) 5338.
  * [95] S. B. Lee and S. Torquato, _Pair connectedness and mean cluster size for continuum‐percolation models: Computer‐simulation results_ , _The Journal of Chemical Physics_ 89 (1988) 6427 [https://doi.org/10.1063/1.455411].
  * [96] A. Gabrielli, F. Sylos Labini, M. Joyce and L. Pietronero, _Statistical Physics for Cosmic Structures_. Springer-Verlag, 2005, 10.1007/b138455.
  * [97] Y. Baryshev and P. Teerikorpi, _Fractal approach to large-scale galaxy distribution_ , .
  * [98] M. Tegmark, A. J. S. Hamilton, M. A. Strauss, M. S. Vogeley and A. S. Szalay, _Measuring the Galaxy Power Spectrum with Future Redshift Surveys_ , _ApJ_ 499 (1998) 555 [astro-ph/9708020].
  * [99] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek et al., _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample_ , _MNRAS_ 470 (2017) 2617 [1607.03155].
  * [100] O. H. E. Philcox and M. M. Ivanov, _BOSS DR12 full-shape cosmology: $\Lambda$ CDM constraints from the large-scale galaxy power spectrum and bispectrum monopole_, _Phys. Rev. D_ 105 (2022) 043517 [2112.04515].
  * [101] D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga, _Cosmological non-linearities as an effective fluid_ , _JCAP_ 2012 (2012) 051 [1004.2488].
  * [102] F. Villaescusa-Navarro et al., _The Quijote simulations_ , _Astrophys. J. Suppl._ 250 (2020) 2 [1909.05273].
  * [103] Z. Zheng, A. A. Berlind, D. H. Weinberg, A. J. Benson, C. M. Baugh, S. Cole et al., _Theoretical Models of the Halo Occupation Distribution: Separating Central and Satellite Galaxies_ , _ApJ_ 633 (2005) 791 [astro-ph/0408564].
  * [104] M. Sinha and L. H. Garrison, _corrfunc – a suite of blazing fast correlation functions on the CPU_ , _Mon. Not. Roy. Astron. Soc._ 491 (2020) 3022 [1911.03545].
  * [105] D. J. Eisenstein and W. Hu, _Baryonic features in the matter transfer function_ , _Astrophys. J._ 496 (1998) 605 [astro-ph/9709112].
  * [106] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak et al., _nbodykit: an open-source, massively parallel toolkit for large-scale structure_ , _Astron. J._ 156 (2018) 160 [1712.05834].
  * [107] S. Torquato and F. H. Stillinger, _Local density fluctuations, hyperuniformity, and order metrics_ , _Phys. Rev. E_ 68 (2003) 041113.
  * [108] S. Torquato, _Hyperuniform states of matter_ , _Physics Reports_ 745 (2018) 1.
  * [109] B. Widom, _Random sequential addition of hard spheres to a volume_ , _J. Chem. Phys._ 44 (1966) 3888.
  * [110] S. Torquato, O. U. Uche and F. H. Stillinger, _Random sequential addition of hard spheres in high Euclidean dimensions_ , _Phys. Rev. E_ 74 (2006) 061308.