# Differentiating small-scale subhalo distributions in CDM and WDM models
using persistent homology

Jessi Cisewski-Kehe<EMAIL_ADDRESS>Department of Statistics, University of
Wisconsin-Madison, 1300 University Ave, Madison, WI 53706, USA  Brittany
Terese Fasy  Gianforte School of Computing, Montana State University, 357
Barnard Hall, P.O. Box 173880, Bozeman, MT, 59717-3880, USA  Wojciech Hellwing
Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotników
32/46, 02-668 Warsaw, Poland  Mark R. Lovell  Science Institute, University of
Iceland, Dunhaga 5, 107 Reykjavík  Paweł Drozda  Center for Theoretical
Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland
Mike Wu  Computer Science Department, Stanford University, 353 Jane Stanford
Way, Stanford, CA, 94305, USA

###### Abstract

The spatial distribution of galaxies at sufficiently small scales will encode
information about the identity of the dark matter. We develop a novel
description of the halo distribution using persistent homology summaries, in
which collections of points are decomposed into clusters, loops and voids. We
apply these methods, together with a set of hypothesis tests, to dark matter
haloes in MW-analog environment regions of the cold dark matter (CDM) and warm
dark matter (WDM) Copernicus Complexio $N$-body cosmological simulations. The
results of the hypothesis tests find statistically significant differences
(p-values $\leq$ 0.001) between the CDM and WDM structures, and the functional
summaries of persistence diagrams detect differences at scales that are
distinct from the comparison spatial point process functional summaries
considered (including the two-point correlation function). The differences
between the models are driven most strongly at filtration scales $\sim 100$
kpc, where CDM generates larger numbers of unconnected halo clusters while WDM
instead generates loops. This study was conducted on dark matter haloes
generally; future work will involve applying the same methods to realistic
galaxy catalogues.

##  I Introduction

The large scale structure (LSS)—as defined by the spatial distribution of
galaxies—encodes information on many vital aspects of the standard model of
cosmology that remain open questions in physics [1, 2, 3, 4]. For example, the
LSS is sensitive to the characteristics of dark energy, the unexplained
phenomenon that drives the accelerated expansion of the Universe [5, 6] and
also holds clues as to the nature of dark matter (DM). Typical LSS observables
that are relevant for DM studies include the abundance of low mass galaxies
[7, 8], the paucity of galaxies in voids [9] and the spatial distribution of
MW satellite galaxies [10]. An additional, as yet largely untapped, method for
analysing LSS models is the application of topological methods to the
distribution of galaxies and haloes. These methods describe the spatial
distribution of points as different dimensional holes with clusters, filaments
loops, and voids in dimensions 0, 1, and 2, respectively, and it is possible
to envisage that the imprint of DM physics on the primordial density field may
be detectable in their topological statistics [5, 11]. In this paper we will
apply topological methods in order to identify differences between two
competing DM models. The simplest viable model of DM is the cold DM matter
model (CDM), in which the DM particle has a negligible velocity dispersion at
early times and thus DM halos are able to start collapsing early and in large
quantities. The combination of CDM with the cosmological constant model of
dark energy is known as $\Lambda$CDM. This model has enjoyed success in
predicting the properties of the cosmic microwave background (CMB) radiation
[12] and the distribution of galaxies at large scales ($>$2 Mpc) [13].
However, at smaller scales ($<$1 Mpc) there are tensions among others with the
densities of dwarf galaxies that may hint at problems for the CDM model [3].
Given the simultaneous failure to detect the particle physics candidates that
correspond to CDM in direct detection experiments [14, 15] or in indirect
detection observations [16], it is important to consider alternatives.

One compelling alternative to CDM is the warm DM (WDM) model, in which the DM
particles have a significant velocity dispersion in the early Universe [17].
The effects of this velocity dispersion include a drastic reduction in the
number of low mass DM halos. In this study we compare simulations of these two
models to determine whether persistent homology can detect differences in the
DM halo spatial distribution. In Fig. 1, we present images of two realizations
of the Copernicus Complexio (COCO) cosmological volume [18], one simulated
with the CDM model, and the second with the WDM model [19]. The large scale
distribution of matter is nearly identical in the two images—thus WDM
preserves the large scale successes in explaining the distribution of massive
galaxies of CDM—but at smaller scales the abundance of WDM halos is strongly
suppressed relative to CDM, and the distribution of the remaining subhalos is
much less homogeneous.

In this work, we investigate differences in the spatial distribution of DM
haloes as described in CDM and WDM. The primary goal is to ascertain whether
topological methods are sensitive to differences between the models and the
second goal is to interpret the differences to determine whether topological
methods have the potential to discern which DM model most accurately describes
the properties of our own Universe. We restrict our analysis to the
distribution of haloes, which will work as a proof of concept. A comprehensive
comparison with observations will require a mock galaxy catalogue and we defer
this step to future work.

(a) COCO-CDM

(b) COCO-WDM

Figure 1: Illustration of the DM distributions in the COCO-CDM (a) and COCO-
WDM (b) simulations. Each image is a slice through the simulation of 23 Mpc on
a side with an image depth of 10 Mpc. The image intensity encodes the DM
column density and the image color indicates velocity dispersion. Eight of the
77 volumes used in this study (see §II.2) are included in this slice, and
their locations are indicated as follows. The MW-analog halo on which each
volume is centered is enclosed by a green circle, and the full extent of the
analysis volume (a radius 3 Mpc) is shown with a white circle. Note that the
apparent overlap of the white circles in this projection does not imply that
the volumes overlap: there can still be considerable separation between the
volumes in the depth direction. See §II for details on the COCO data.

The persistent homology formulation of topology offers a novel way to
represent, visualize, and interpret complex data by extracting homological
features, which can be used to infer properties of the underlying structures.
Homological features include the decomposition of halo distributions into
clusters, filaments loops, and voids at different scales controlled by a
parameter that is analogous to halo linking lengths—which in statistics is
known as a filtration parameter—and persistent homology in particular tracks
how the number of such features changes as the filtration parameter is
increased. It has been successfully applied to problems in astronomy (e.g.,
Refs. [20, 21, 5, 22, 23, 24, 25, 26, 27]), along with other areas of science
(e.g., Refs. Ref. [28, 29, 30, 31]). There have been proposals for hypothesis
testing using persistent homology (e.g., Refs. [32, 33, 34, 31, 35]), which we
build on as we construct tests that can detect differences between DM model
predictions in the LSS.

We investigate several test statistics to discriminate between the CDM and WDM
halo spatial distributions that are based on persistent homology functional
summaries. Each functional summary is a different transformation of
information to a function that approximates a property of the topological
features, and is a function of the filtration parameter. We also consider
different visualizations in order to investigate detected differences.

This paper is organised as follows. We begin with background on the
cosmological simulation data we use in the analysis (§II), then we introduce
persistent homology and functional summaries of persistence diagrams that are
used in the proposed test statistics (§III). Then the hypothesis testing
framework is presented (§IV), followed by the investigation of the
cosmological simulation data (§V). We end with concluding remarks (§VI).

##  II Cosmological Simulation Data

This section begins with a description of the COCO simulations, and then
continues with our procedure for selecting MW halo-analog sample regions.

###  II.1 The Copernicus Complexio (COCO) cosmological simulations

The COCO simulation volume constitutes a high resolution spherical region of
space with a comoving radius of approximately 25 Mpc; the full (low-
resolution) simulation volume is a periodic box 100 Mpc on a side.111All
distances are in comoving Mpc. The numerical integration of the gravitational
forces begins at redshift 127. The cosmological parameters are consistent with
the 7-year results from the WMAP satellites: matter density
$\Omega_{0}=0.272$, dark energy density $\Omega_{\Lambda}=0.728$,
$\Omega_{b}=0.04455$, Hubble parameter $h_{0}=0.704$, spectral index
$n_{s}=0.967$, and power spectrum normalization $\sigma_{8}=0.81$. The mass of
the simulation particle is $1.135\times 10^{5}$ $M_{\odot}$. DM halos and
subhalos were identified using the SUBFIND algorithm [36], and the smallest
permitted number of particles to identfy a subhalo is 20 particles. Our
definition of halo mass is the total mass bound gravitationally to each halo
as determined by the halo finder.

Two copies of this volume were run, the first applying CDM [18] and the second
WDM [19]. Both simulations use the same initial phases, and differ in that the
WDM simulation had wave amplitudes rescaled using the transfer function of a
3.3 keV thermal relic DM particle, with the relic mass chosen to be in
agreement with the Lyman-$\alpha$ forest constraints of Ref. [37]. This
results in the suppression of structure on the scale of dwarf galaxies. These
large-scale structure similarities between the WDM and CDM data due to the
same initial phases are shown in Fig. 1. One issue peculiar to WDM simulations
is the spurious numerical fragmentation of filaments into halos; these so-
called spurious subhalos are identified and removed from the halo catalog
using the algorithm described in Ref. [38].

###  II.2 Milky Way-analog DM halos and their associated halo samples

Given that we intend to use future work to compare the models with
observations of galaxies around our own MW, we identify MW-analog halos and
their surrounding regions in the two simulations. The criteria for our MW-
analog halos were that they must be located within 21 Mpc of the center of the
simulation222The central high-resolution sphere of COCO extends out to about
25 Mpc. and have a mass in the range $[0.5,2]\times 10^{12}$ $M_{\odot}$. We
also required that there be no other halo with a mass greater than $0.5\times
10^{12}$ $M_{\odot}$ within 0.7 Mpc. This procedure resulted in 77 MW-analog
DM halos in each of the COCO CDM and WDM simulations; for each MW-analog DM
halo in the WDM data, there is a matching MW-analog DM halo in the CDM.

We now discuss our selection of halos in the vicinity of the 77 MW halo-
analogs. In both CDM and WDM realizations we identify halos that are within 3
Mpc of the MW-analog333Objects with $\sim$0.2 Mpc of each analog are typically
referred to as ‘subhalos’ that orbit within the analog ‘host halo.’ In this
study we refer to all bound DM objects simply as ‘halos’ and include all of
them in our analysis, not drawing any distinction between ‘subhalos’ and other
classes of object.. For the WDM case we include all halos in the 3 Mpc region.
However, CDM forms hundreds of times more halos than WDM in our resolved mass
range. If we were to include all CDM halos, the abundance difference would
dominate our statistical results. Therefore, the CDM samples were downsampled
to match the number of DM halos in the corresponding WDM sample. The
downsampling was accomplished by selecting the most massive DM halos from each
of the CDM samples. An example of one of the MW-analog DM halo neighborhoods
from COCO-CDM and COCO-WDM is displayed in Fig. 2.

(a) CDM

(b) WDM

Figure 2: A MW-analog DM halo neighborhood sample for (a) the COCO-CDM data
and (b) the COCO-WDM data. The red triangles indicate the MW-analog DM halo
that was selected along with the other DM halos that are within 3 Mpc from it.

The two sets of samples for the CDM and WDM data are defined as

$\mathbb{Y}_{c}=\\{\mathbf{Y}_{c,1},\ldots,\mathbf{Y}_{c,77}\\},\mathbb{Y}_{w}=\\{\mathbf{Y}_{w,1},\ldots,\mathbf{Y}_{w,77}\\}$
(1)

where $\mathbb{Y}_{c}$ and $\mathbb{Y}_{w}$ represent the set of 77 CDM and 77
WDM samples, respectively. Each $\mathbf{Y}_{k,i}\in\mathbb{R}^{n_{i}\times
3}$ for $k=c,w$ and $i=1,\ldots,77$ where $n_{i}$ indicates the number of DM
halos in sample $i$; an individual DM halo in simulation $k$ of sample $i$ is
indicated by $Y_{k,i,j}$ for $j=1,\ldots,n_{i}$.

##  III Topological Data Analysis Methods for Quantifying LSS

Homology is one way to study the features of topological spaces (e.g.,
manifolds), specifically the dimensional “holes” in the space (e.g., connected
components, loops, voids). Persistent homology studies the spatial structure
of a parameterized family of topological spaces that keeps track of the so-
called _births_ and _deaths_ of homological features as a topological space
changes with a filtration parameter. In particular, we focus on point cloud
data, where each point can represent some unit of mass or an object (e.g., a
point may represent a center of a DM halo). In this section, we provide a
brief overview of the necessary concepts; however, see, e.g., Refs. [39, 40,
41] for a more thorough introduction to algebraic and computational topology.
The homological features that are tracked in the filtration have cosmological
interpretations in dimensions zero, one, and two. Before providing more
details about persistent homology, we explain the interpretation of different
dimensional holes with respect to the distribution of DM halos.

#### Clusters

A _connected component_ , or zeroth-dimensional homology feature ($H_{0}$), is
a maximal subspace of a topological space that cannot be covered by two
disjoint open sets; that is, a connected component is a whole piece of the
space. For example, under some assumptions on the topological space, a
connected component is a cluster of data points.444There is no clear or
established relationship between the definition of homological clusters and
galaxy clusters. In this paper, the term ‘cluster’ is only used for the
homology definition. In cosmology, the connected components represent clusters
of halos or galaxies. Persistent homology tracks the appearance of new
connected components and the merging of distinct components.

#### Filaments and Loops

A _loop_ , or one-dimensional homology feature ($H_{1}$), provides information
about the connectivity of data. As the filtration parameter increases, nearby
connected components can merge together in such a way that a loop is formed.
For DM halos, this would appear as filaments of halos joined together in a
loop.

#### Cosmological Voids

A _void_ , or two-dimensional homology feature ($H_{2}$), represents the
boundary of three-dimensional empty regions within the topological space
(e.g., the boundary of a football). In cosmology, these are the thin walls
surrounding the low-density regions that are typically referred to as
cosmological voids.

###  III.1 Persistent homology

Persistent homology is a framework for computing the homology of a data set at
different scales. Given a data set, one defines a filtration (that is, a
sequence of nested topological spaces) of intermediate structures, on which
the homology is computed at different values of the filtration parameter.
Homology group generators (e.g., $H_{0}$, $H_{1}$, $H_{2}$ features) are
tracked as they form and die as the filtration parameter changes. Various
methods can be used in order to transform a discrete point set into a
connected topological space. For example, simplicial complexes (see below)
such as the Vietoris–Rips complex (VR complex) can be used, or a function can
be defined over the domain of the data using an empirical distance function or
a kernel density estimate (KDE) of the point cloud. In this work, we use a VR
complex to construct the filtration (discussed below). Next, we introduce some
key components of persistent homology.

#### Simplicial complexes

A geometric $k$-simplex is the convex hull of $k+1$ affinely independent
points. For our 3D halo data, the simplices we use are zero-simplices
(vertices), one-simplices (edges), two-simplices (triangles), and three-
simplices (tetrahedrons). A face of a simplex is another simplex obtained by
removing zero or more points (e.g., a triangle has seven faces: itself, three
edges, and the three vertices). A simplicial complex, $\mathcal{K}$, is
defined as a finite set of simplices such that (i) if $\sigma\in\mathcal{K}$,
then every face of $\sigma$ is also in $\mathcal{K}$, and (ii) if
$\sigma_{1},\sigma_{2}\in\mathcal{K}$, then either
$\sigma_{1}\cap\sigma_{2}\in\mathcal{K}$ or
$\sigma_{1}\cap\sigma_{2}=\emptyset$. In this work, we use simplicial
complexes to represent topological spaces, as they are the standard input to
code to compute homology.

#### Filtrations

A filtration is a sequence of nested topological spaces. Given dataset
$y_{1},y_{2},\ldots,y_{n}\in\mathcal{Y}\subseteq\mathbb{R}^{3}$, one common
way to create a simplicial complex is to choose some $t\in\mathbb{R}$ such
that $t\geq 0$ and replace each $y_{i}\in\mathcal{Y}$ with a ball of diameter
$t$. The Vietoris–Rips complex at scale $t$ (the $t$-VR complex) is created by
representing each of these balls as a vertex, and creating a $k$-simplex
anytime there are $k+1$ balls that pairwise intersect. Specifically:

$VR_{t}(S)=\\{\sigma\subseteq S\mid d(x,z)\leq t,\forall x,z\in\sigma,\\}$ (2)

where $d(\cdot,\cdot)$ is the Euclidean distance [41, 42]. That is,
$VR_{t}(S)$ is a simplicial complex containing the vertex set $\mathcal{S}$,
edges between all the vertices that are separated by at most $t$, and
triangles for sets of three vertices that have pairwise distances of at most
$t$.

We obtain the VR filtration by increasing $t$ from $0$ to $\infty$ (here, $t$
is referred to as the filtration parameter).555In practice, the maximum
filtration value $t$ we consider correspond to the largest scales encompassed
by a given galaxy/halo catalog. Note that $VR_{t_{1}}(S)$ is a subset of
$VR_{t_{2}}(S)$ (i.e., $VR_{t_{1}}(S)\subseteq VR_{t_{2}}(S)$) for $t_{1}\leq
t_{2}$. Sometimes, for the right selection of $t$ and a dense enough sample,
we can recover the homology of $\mathcal{Y}$ (see, e.g., Ref. [43]); however,
using the whole sequence of complexes, we can recover information about
$\mathcal{Y}$ with more relaxed sampling conditions.

To derive the persistent homology for a VR filtration, the homology of
$VR_{t}(S)$ is computed as $t$ changes. If $t$ is initialized at $0$, then
only the data points contribute to the homology. The evolving topological
space is characterized by its homology as $t$ increases toward $\infty$. For a
$\mathcal{Y}\subseteq\mathbb{R}^{3}$, the persistent homology would then track
the connected components ($H_{0}$), loops ($H_{1}$), and voids ($H_{2}$) that
appear and disappear in the VR filtration. An example of a VR filtration with
a 2-dimensional domain is presented in Fig. 3: Figures 3(a) and 3(b) display
the data points with balls of diameter t = 0.5 and 1, respectively, along with
the one- and two-simplices of the corresponding VR complex. Fig. 3(c) shows
the persistence diagram for the data points using the VR filtration, which is
discussed next.

#### Tracking Homology Generators

The birth and death times of the homology group generators are displayed in a
persistence diagram. These times correspond to values of the filtration
parameter, which is the diameter of the balls $t$ when considering a VR
filtration. Suppose a filtration is defined over some data points
$y_{1},y_{2},\ldots,y_{n}\in\mathcal{Y}\subseteq\mathbb{R}^{3}$, then a
persistence diagram, $\mathbf{D}$, can be written as a multiset of points:

$\mathbf{D}=\\{(r_{j},b_{j},d_{j}):j=1,\ldots,|\mathbf{D}|\\}\cup\Delta$ (3)

where $(r_{j},b_{j},d_{j})$ are the homology group dimension, the birth time,
and the death time, respectively, of feature $j$, $|\mathbf{D}|$ indicates the
number of homology group generators with $d_{j}>b_{j}$, and $\Delta$
represents a set of points on the diagonal (birth time = death time) with
infinite multiplicity. The persistence diagram is a nice summary because small
changes in the input data $\mathcal{Y}$ will result in only small changes in
the diagram [44, 45].

Figures 3(a) and 3(b) show an example where the filtration parameter, $t$,
increases from 0.5 to 1. In that interval, the homology changed from having 15
connected components ($H_{0}$’s) and zero loops ($H_{1}$’s) to having 5
connected component and 1 loop. The time in the filtration when homology
features appear, the _birth_ of the feature, and the time when a feature joins
other features, the _death_ of the feature, are captured in a persistence
diagram. Fig. 3(c) displays the persistence diagram, where the location of
each point represents the birth height (x-axis) and death height (y-axis) of a
homological feature for a VR filtration, and the shape and color represent the
homology group dimension. A point $(\cdot,x,x)$ on the diagonal represents a
feature with a zero-length lifespan. The _persistence_ of a point
$(\cdot,b,d)$ is the length of the interval of the persistence parameter that
supports that feature: $|d-b|$. In the persistence diagram, the distance from
$(\cdot,b,d)$ to the diagonal is proportional to this value; in fact, the
(Euclidean) distance to the diagonal is $\frac{|d-b|}{\sqrt{2}}$. When working
with empirical experimental or observational data it is necessary to be
concerned with the associated intrinsic noise of such measurements. This is
especially important in the context of a spatial distribution of objects
derived either from N-body simulations or galaxy catalogs. N-body simulations
are limited by their spatial resolution, where their Monte Carlo sampling
nature starts to breaks down and is overrun by the shot-noise. The
astronomical observations are limited by imperfections including those related
to involved instruments which contribute to measurements errors. For these
reasons it is important in a persistent homology analysis to be able to
distinguish between real features really present in the target and those that
are noise-induced transients. A notion of a topological significance can be
derived in this context by considering features with longer lifetimes as more
significant, and those with short lifetimes (i.e., closer to the diagonal) as
topological noise [46]. Distinguishing between topological signal and noise is
a problem of a great interest in real applications of TDA (e.g., Ref. [26]).

(a) Data with $t=0.5$ balls

(b) Data with $t=1$ balls

(c) Persistence diagram

Figure 3: Persistence diagram example where observations were sampled around
three circles with noise. The data are displayed in (a) and (b) as black
points (zero-simplices) with cyan balls with diameters of 0.5 and 1,
respectively, along with the one- and two-simplices of the corresponding VR
complexes. The persistence diagram for the VR filtration of the points is
displayed in (c) with the the three circles indicated by the cyan triangles
($H_{1}$). The $H_{0}$ features represent the connected components, which all
have birth times at 0.

###  III.2 Persistence diagram summaries

While persistence diagrams and their individual features provide useful
information about the topology of a data set, persistence diagrams are not
easy objects to work with directly for statistical analyses. For example, the
distance between two persistence diagrams can be calculated using metrics such
as the bottleneck distance or the $p$-Wasserstein distance, but both are
computationally expensive because they require finding a certain optimal
matching between the features on each diagram; see Equation (21) in the
Appendix for the definition of the bottleneck distance. Fréchet means and
medians have been defined for spaces of persistence diagrams [47], but are
also computationally expensive and not necessarily unique (although ways
around this exist as addressed in Ref. [48]). Instead, we consider
transformations and summaries of persistence diagrams that make computations
more tractable [31]. Below are several approaches that transform a persistence
diagram into a functional summary, which are used in §IV to formulate test
statistics for hypothesis tests.

#### Landscape Functions

Landscape functions [33] are popular functional summaries of persistence
diagrams [49, 31, 50], which are defined as follows. Let
$\mathbf{D}_{r}=\\{(b_{j},d_{j})\\}_{j=1}^{n_{r}}$ be the finite set of off-
diagonal points of a homology dimension $r$ persistence diagram. Next, rotate
the persistence diagram such that each point $(b_{j},d_{j})\in\mathbf{D}_{r}$
is mapped to
$p_{r,j}=\left(\frac{b_{j}+d_{j}}{2},\frac{d_{j}-b_{j}}{2}\right)\in\widetilde{\mathbf{D}}_{r}$.
Isosceles right triangles are formed from each $p_{r,j}$ to the base as

$\Lambda_{p_{r,j}}(t)=\begin{cases}t-b_{j}&\quad
t\in[b_{j},\frac{d_{j}+b_{j}}{2}]\\\ d_{j}-t&\quad
t\in[\frac{d_{j}+b_{j}}{2},d_{j}]\\\ 0&\quad\text{ otherwise},\\\ \end{cases}$
(4)

where $t\in[t_{\min},t_{\max}]$. The persistence landscape is then defined as
the following collection of functions

$\lambda_{\mathbf{D}_{r}}(k,t)=\underset{p_{r,j}\in\widetilde{\mathbf{D}}_{r},}{\text{kmax
}}\Lambda_{p_{r,j}}(t),t\in[t_{\min},t_{\max}],k=1,\ldots,n_{r},$ (5)

where kmax is the $k$-th largest value. An example of a persistence landscape
function is displayed in Fig. 4. Rather than working with each $k$ of
$\lambda_{\mathbf{D}_{r}}(k,t)$ individually, a subset of the landscape layers
can be concatenated to a long vector as

$\mathcal{F}_{\text{land}}(\mathcal{I},r,t)=\oplus_{k\in\mathcal{I}}\lambda_{\mathbf{D}_{r}}(k,t),$
(6)

where $\mathcal{I}$ is the index set of the included landscape layers.

(a) Persistence Diagram

(b) Summaries

Figure 4: A persistence diagram (a) along with its landscape functions and a
weighted silhouette (b) for an arbitrary homology dimension $r$. The dotted
pink curve is the weighted silhouette function with tuning parameter $p=1$;
the other four curves correspond to landscape functions
$\lambda_{\mathbf{D}_{r}}(k,t)$ for $k=1,\ldots,4$.

#### Weighted Silhouette Functions

Rather than working with each $k$ of $\lambda_{\mathbf{D}_{r}}(k,t)$ from
Equation (5) individually, weighted silhouette functions provide a way of
combining the information in the collection of landscape functions.
Silhouettes are weighted averages of the individual functions for homology
dimension $r$ defined as

$\mathcal{F}_{\text{sil}}(r,t\mid
p)=\frac{\sum_{j=1}^{n_{r}}|d_{r,j}-b_{r,j}|^{p}\Lambda_{p_{r,j}}(t)}{\sum_{j=1}^{n_{r}}|d_{r,j}-b_{r,j}|^{p}},$
(7)

where the $|d_{r,j}-b_{r,j}|^{p}$ act as weights that can give more emphasis
or less emphasis to features with longer lifetimes depending on the user-
specified parameter $p$. The form of these weights are suggested in Ref. [51].
An example of a weighted silhouette function is provided in Fig. 4(b). More
details and theoretical properties of landscapes and silhouettes can be found
in Ref. [51].

#### Euler Characteristic and Betti Functions

The Euler characteristic (EC) is a topological invariant and can be defined as
the alternating sum of the rank of the homology groups, where the rank of the
$r$th homology group is the $r$th Betti number. As the persistent homology
filtration parameter $t$ changes and new features are born or old ones die,
the Betti numbers and EC changes, allowing for the definition of Betti
functions and an EC function. The Betti functions can be defined as

$\mathcal{F}_{\text{betti}}(r,t)=|\\{(r,b_{j},d_{j}):b_{j}\leq t,d_{j}>t\\}|,$
(8)

which indicates the number of dimension $r$ homology group generators that
persist in the filtration at time $t$. The only non-trivial homology groups
for data in $\mathbb{R}^{3}$ are in dimensions 0, 1, and 2; thus, the Euler
characteristic equation we use is

$\mathcal{F}_{\text{ec}}(t)=\sum_{r=0}^{2}(-1)^{r}\mathcal{F}_{\text{betti}}(r,t).$
(9)

Betti and EC functions have been used in applications [e.g., 52, 53, 54, 23,
25, 55, 56] and some of their theoretical properties have been explored [e.g.,
57, 58, 59, 60, 35].

There are a number of other summary functions of persistence diagrams that
have been defined [e.g., 61, 62, 34]. For a general discussion of summary
functions of persistence diagrams, including some theoretical properties, see
Ref. [31].

##  IV Methods: topological hypothesis tests for LSS

A primary goal in this work is to develop a framework that can inferentially
discriminate between different realizations of web-like geometric data
structures such as the Cosmic Web. Our TDA-based framework allows for
extracting information, encoded in large-scale galaxy/halo distribution, that
goes beyond methods commonly used in cosmology N-point clustering statistics.
The motivation is to detect differences between the DM halo spatial
distributions (i.e., 3-manifolds) evolved in cosmological simulations where
initial conditions were set to be either that of CDM or WDM-type. In this
section we present a hypothesis testing framework using test statistics
derived from the summaries of persistence diagrams presented in §III.2. These
topological hypothesis tests build on the work outlined in Ref. [31],
including their notation.

The proposed hypothesis tests rely on permutation methods to compute the
p-values. There have been some central limit theorem results for summary
functions of persistence diagrams (e.g., Betti functions) and persistent
homology-based hypothesis tests statistics using asymptotic theory (e.g.,
Refs. [63, 58, 35]), but they generally assume the data were drawn from a
homogeneous Poisson point process. Both the large and small-scale halo/galaxy
distribution in the Universe cannot be described by a homogeneous Poisson
process. Owing to the nature of initial conditions (i.e., an adiabatic
Gaussian random field) and the gravitational instability (a mechanism
responsible for the growth and evolution of the cosmic structures) halos
spatial distribution is clustered with non-Gaussian features on small (non-
linear) scales. Naturally, also the COCO DM simulations provide halo samples
which are not close to resembling homogeneous Poisson point processes, which
is discussed in Appendix §A. Furthermore, as we present below, the WDM and CDM
samples are not independent of one another due to the cosmological simulation
design: it is therefore necessary to use matched-pairs hypothesis tests.

###  IV.1 Test statistic and p-value computations

For the proposed hypothesis tests, we consider two samples of observations,

$\mathbf{Y}_{1}=\\{Y_{1,1},\ldots,Y_{1,n_{1}}\\},\text{ and
}\mathbf{Y}_{2}=\\{Y_{2,1},\ldots,Y_{2,n_{2}}\\}$ (10)

where each $Y_{j,i}$, $i=1,\ldots,n_{j}$ and $j=1,2$ is a data set of which a
persistence diagram can be computed. For our cosmological simulation data,
each $Y_{j,i}$ will have a set of points in $\mathbb{R}^{3}$, but, in general,
the $Y_{j,i}$’s could take different forms; for example, each $Y_{j,i}$ could
be an image of a fibrin network [31] or a brain artery tree [29].

Each observation from Equation (10) will have a corresponding persistence
diagram

$\mathbf{D}_{1}=\\{\mathbf{D}_{1,1},\ldots,\mathbf{D}_{1,n_{1}}\\},\text{ and
}\mathbf{D}_{2}=\\{\mathbf{D}_{2,1},\ldots,\mathbf{D}_{2,n_{2}}\\}.$ (11)

These samples of diagrams can be used to test the hypotheses,

$H_{0}:\mathcal{P}_{1}=\mathcal{P}_{2}\text{ vs.
}H_{1}:\mathcal{P}_{1}\neq\mathcal{P}_{2},$ (12)

where $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are the true underlying
distributions of persistence diagrams from group 1 and 2,
respectively.666Probability measures can be theoretically defined on a space
of persistence diagrams (with a Wasserstein metric) as presented in Ref. [64].

Given two samples of persistence diagrams, there are a number of possible ways
to derive test statistics; we consider the functional versions of persistence
diagrams presented in §III.2 as test statistics. These functional summaries
can be understood as a map between the space of persistence diagrams,
$\mathcal{P}$, to the space of functions, $\mathcal{F}$, defined as
$\mathbb{F}:\mathcal{P}\longrightarrow\mathcal{F}$. Therefore, the diagrams
from above can be used to define the collection of functional summaries with
$j=1,2$ as

$\mathbf{F}_{j}=\\{F_{j,1}=\mathbb{F}(\mathbf{D}_{j,1}),\ldots,F_{j,n_{j}}=\mathbb{F}(\mathbf{D}_{j,n_{j}})\\}.$
(13)

A test statistic for the two-sample hypothesis test of Equation (12) can be
derived using estimates of functional summaries. Letting
$F_{j,i}=\mathbb{F}(\mathbf{D}_{j,i}),i=1,\ldots,n_{j},\text{ and }j=1,2$ (see
Equation (13)), the mean functional summaries are defined as

$\bar{F}_{j}(t)=n_{j}^{-1}\sum_{i=1}^{n_{j}}F_{j,i}(t).$ (14)

Then our test statistic for the different functional summaries is based on the
following distance between mean functional summaries,

$d(\bar{F}_{1},\bar{F}_{2})=\int_{\mathbb{T}}|\bar{F}_{1}(t)-\bar{F}_{2}(t)|dt,$
(15)

where $\mathbb{T}$ defines the domain of the functions. Note that the
$\mathbb{T}$ is related to the range of values of the filtration parameter,
which depends on the functional summary. For example, for the Euler
characteristic function, the $\mathbb{T}$ covers the range of the filtration
parameter, but for the landscape and silhouette functions it represents the
range of a transformed filtration parameter since the persistence diagram is
rotated.

#### Matched-Pairs Permutation Test.

Since the distributions of the test statistics of Equation (15) for the
different functional summaries we consider are unknown, the p-values for the
two-sample hypothesis tests can be computed using the usual permutation
testing framework. The general procedure is to randomly assign the
$n_{1}+n_{2}$ functional summaries into two groups, because this random
assignment is consistent with the null hypothesis ($H_{0}$) from Equation (12)
where the two groups follow the same distribution. Using the random group
assignments, the new mean functional summaries are estimated using Equation
(14), $\widetilde{F}_{1}^{(l)}\text{ and }\widetilde{F}_{2}^{(l)}$, which are
used to compute the distance
$d(\widetilde{F}_{1}^{(l)},\widetilde{F}_{2}^{(l)})$ from Equation (15), for
$l=1,\ldots,n_{l}$ random permutations [31]. The resulting (approximate)
permutation p-value can then be computed as

$p_{\text{perm}}=n_{l}^{-1}\sum_{l=1}^{n_{l}}I(d(\widetilde{F}_{1}^{(l)},\widetilde{F}_{2}^{(l)})\geq
d(\bar{F}_{1},\bar{F}_{2})),$ (16)

where $I(A)$ is an indicator function that takes the value 1 if $A$ is true
and 0 if $A$ is false.

When the two sets of samples are independent, the above permutation p-value is
reasonable. However, as explained in §II, the CDM and WDM COCO data,
$\mathbb{Y}_{c}$ and $\mathbb{Y}_{w}$, are not independent due to the initial
conditions of the simulations. Instead, the samples $\mathbf{Y}_{c,i}$ and
$\mathbf{Y}_{w,i}$ for $i=1,\ldots,77$ have a similar spatial structure which
should be accounted for in the computation of the permutation p-values.
Therefore, we consider a matched-pairs version of the permutation p-values.
The difference between this matched-pairs version and the permutation test
outlined above is in how the two groups are randomly assigned for each
permutation. The matched-pairs permutation involves randomly selecting one of
the two matched samples to go into each of the two groups (e.g., one of
$\mathbf{Y}_{c,i}$ or $\mathbf{Y}_{w,i}$ will be randomly assigned to group 1,
and the other will be assigned to group 2). The matched-pairs permutation
p-value is then defined in the same manner as above, as

$p_{\text{matched}}=\sum_{l=1}^{n_{l}}I(d(\widetilde{F}_{1,\text{matched}}^{(l)},\widetilde{F}_{2,\text{matched}}^{(l)})\geq
d(\bar{F}_{1},\bar{F}_{2})),$ (17)

where $\widetilde{F}_{j,\text{matched}}^{(l)},j=1,2,$ are the mean functional
summary for permutation $l$ using the matched-pairs random assignment. This
matched-pairs permutations p-value computation accounts for correlations
between the COCO CDM and WDM samples by including one of the two matched
samples within each group for each permutation, but randomizing which label
(CDM or WDM) is assigned.

##  V Investigation of COCO Simulation Data

In order to investigate differences between the CDM and WDM COCO samples of
MW-analog halo neighborhoods described in §II.2, $\mathbb{Y}_{c}$ and
$\mathbb{Y}_{w}$, respectively, we carry out the two-sample hypothesis tests
defined in Equation (12) and described in the previous section. The test
statistics are based on the functional summaries of persistence diagrams
outlined in §III.2, along with several other methods discussed below. The
comparison methods include a test statistic that uses persistence diagrams
directly (rather than a functional summary of them) and non-TDA functional
summaries that capture second-order properties of spatial point processes. The
collective goals of the test statistics considered are (i) to detect
differences between CDM and WDM MW-analog halo neighborhoods, and (ii) to
understand and interpret any detected differences (e.g., the distance scale at
which differences occur).

In addition to the test statistics using the functional summaries presented in
§IV, we also consider other approaches. One method is the persistence diagram-
based test (PDT) of Ref. [32] which has a test statistic defined using
distances between persistence diagrams. We also consider two functional
summaries of spatial point processes which do not use persistence diagrams,
namely the G-function and the two-point correlation function (2PCF). The
G-function gives the distribution function of the nearest-neighbor distances,
and 2PCF uses the Landy-Szalay estimator [65] and is one of the most basic and
fundamental objects used to study clustering in cosmology [66]. These methods
are described in more detail in Appendix §B. The p-values for these additional
tests are also carried out using permutations, and we have also adapted them
to work for our matched-pairs design.

###  V.1 Hypothesis testing results

All the hypothesis tests that use statistics derived from persistence diagrams
(including the PDT) use the same persistence diagrams, which were computed
using a VR filtration. These computations were carried out with _Ripser_ [67].
P-values were computed for the test statistics discussed previously based on
20,000 permutations using the traditional and matched-pairs permutation
methods of §IV.1. The results are displayed in Table 1. Below we discuss the
resulting p-values, and in the next section we investigate and interpret where
the differences are most pronounced.

Overall, statistically significant differences with p-values $<0.001$ are
apparent between the CDM and WDM MW-analog DM halo neighborhoods samples, and
$p_{\text{perm}}\geq p_{\text{matched}}$ for all test statistics
considered.777If our test statistics were Gaussian distributed, a p-values
$<0.001$ would correspond to $>3\sigma$ significance. Since the two sets of
samples are from different populations (CDM vs. WDM COCO data), it is a
positive result that our proposed tests are able to detect differences. For
$H_{0}$ and $H_{1}$, all the function-based tests had
$p_{\text{matched}}<0.001$, and this was also the case for the G-function and
2PCF test statistics. PDT has $p_{\text{perm}}$ and $p_{\text{matched}}\leq
0.003$ for $H_{0}$, but higher p-values for $H_{1}$ with
$p_{\text{perm}}=0.255$ and $p_{\text{matched}}=0.036$. Because the PDT uses
the bottleneck distance, only one $H_{1}$ feature on each of the persistence
diagrams contribute to the test statistic for each MW-analog halo neighborhood
sample, while the functional summary-based test statistics considers all the
features on the persistence diagrams (except for the landscape functions which
only includes features that contribute to the first 10 layers).

Aside from the silhouette function tests, the $H_{2}$ p-values are $<0.01$ for
both $p_{\text{perm}}$ and $p_{\text{matched}}$. The $p_{\text{perm}}$ for the
silhouette function tests are $>0.10$, but then drop below $0.01$ for
$p_{\text{matched}}$. The EC function test statistic, similar to the related
Betti function test statistics, has both $p_{\text{perm}}$ and
$p_{\text{matched}}\leq 0.001$. For the TDA-based test statistics, the EC,
Betti, and Landscape function test statistics appear to be best able to detect
differences between the CDM and WDM MW-analog halo neighborhood samples for
both the traditional and matched-pairs permutation tests across the three
homology dimensions ($H_{0}$, $H_{1}$, $H_{2}$). Tuning could be carried out
for the landscape function tests to find which landscape function layers are
most informative at detecting differences. Since the results with the first 10
layers performed well, we did not consider tuning for this analysis.

Table 1: COCO data results. Permutation p-values ($p_{\text{perm}}$) and matched permutation p-values ($p_{\text{matched}}$) for tests comparing the CDM and WDM MW-analog halo neighborhood samples. The p-values are rounded to three decimal places and are based on 20,000 permutations as described in §IV.1.  Test statistic | Notation | Homology dimension | $p_{\text{perm}}$ | $p_{\text{matched}}$
---|---|---|---|---
Landscape | $\mathcal{F}_{\text{land}}(1:10,0,t)$ | 0 | 0 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(0,t\mid p=0.5)$ | 0 | 0.009 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(0,t\mid p=1)$ | 0 | 0.001 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(0,t\mid p=2)$ | 0 | 0 | 0
Betti | $\mathcal{F}_{\text{betti}}(0,t)$ | 0 | 0.001 | 0
PDT | $\mathcal{T}_{\text{PDT}}(D_{1,\cdot\mid 0},D_{2,\cdot\mid 0}\mid\infty,1)$ | 0 | 0.003 | 0
Landscape | $\mathcal{F}_{\text{land}}(1:10,1,t)$ | 1 | 0 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(1,t\mid p=0.5)$ | 1 | 0.009 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(1,t\mid p=1)$ | 1 | 0.007 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(1,t\mid p=2)$ | 1 | 0.014 | 0
Betti | $\mathcal{F}_{\text{betti}}(1,t)$ | 1 | 0 | 0
PDT | $\mathcal{T}_{\text{PDT}}(D_{1,\cdot\mid 1},D_{2,\cdot\mid 1}\mid\infty,1)$ | 1 | 0.255 | 0.036
Landscape | $\mathcal{F}_{\text{land}}(1:10,2,t)$ | 2 | 0 | 0
Silhouette | $\mathcal{F}_{\text{sil}}(2,t\mid p=0.5)$ | 2 | 0.123 | 0.003
Silhouette | $\mathcal{F}_{\text{sil}}(2,t\mid p=1)$ | 2 | 0.158 | 0.009
Silhouette | $\mathcal{F}_{\text{sil}}(2,t\mid p=2)$ | 2 | 0.135 | 0.008
Betti | $\mathcal{F}_{\text{betti}}(2,t)$ | 2 | 0.001 | 0
PDT | $\mathcal{T}_{\text{PDT}}(D_{1,\cdot\mid 2},D_{2,\cdot\mid 2}\mid\infty,1)$ | 2 | 0.009 | 0
Euler characteristic | $\mathcal{F}_{\text{ec}}(t)$ | 0-2 | 0.001 | 0
G-function | $\mathcal{F}_{\text{G}}(t)$ | N/A | 0 | 0
2PCF | $\mathcal{F}_{\text{2PCF}}(t)$ | N/A | 0 | 0

Given that we only have one COCO-CDM and one COCO-WDM realization, and that we
seek to evaluate the performance of the proposed test statistics when the null
hypothesis is true, we consider bootstrap samples of the data from the CDM
data and from the WDM data. The distribution of the p-values when the null
hypothesis is true should follow a uniform distribution. Details of this
simulation study and the results are presented in Appendix §C. Overall, we
find the p-values resulting from proposed test statistics based on the
functional summaries of persistence diagrams under the null hypothesis are
generally consistent with uniform distributions.

###  V.2 Interpretation of results

In this section, we explore the Betti functions in more detail and develop
other visualizations to aid in the interpretation of the results in order to
investigate the scales at which the differences between the CDM and WDM MW-
analog halo neighborhood samples occur and are significant. Since our interest
is in where the test statistics diverge, the mean difference function is
displayed where the signal is based on the matched data in the CDM and WDM
COCO MW-analog halo neighborhood samples using

$\bar{F}_{\text{diff}}(t)=n_{s}^{-1}\sum_{i=1}^{n_{s}}\left(F_{c,i}(t)-F_{w,i}(t)\right)$
(18)

where $F_{c,i}(t)$ and $F_{w,i}(t)$ are functional summaries for CDM and WDM
sample $i$, respectively, and $n_{s}=77$. Additionally, 95% global confidence
bands are computed using the bootstrap approach outlined in Section 3.2 of
Ref. [31], with 1000 bootstrap samples.888Note that the hypothesis tests use
$L_{1}$ distances between functions (see Equation (15)) while the confidence
bands are investigating differences across the functions.

The CDM and WDM MW-analog halo neighborhood samples’ persistence diagrams were
generated using using a VR filtration. For example, Fig. 5 displays the
persistence diagrams for the COCO CDM and WDM samples of Fig. 2(a) and 2(b),
respectively. The CDM and WDM persistence diagrams in this example share a
similar pattern where generally the $H_{0}$ features are all connected by
around a filtration parameter value of 1, $H_{1}$ features persist longer than
the $H_{2}$ features across the range of birth times. The $H_{0}$ feature
plotted on both diagrams at (0, 2.57) represents an $H_{0}$ feature that in
fact persists indefinitely and should, technically, be plotted at a death time
of infinity.

(a) CDM

(b) WDM

Figure 5: Persistence diagrams for the MW-analog halo neighborhood sample for
(a) the CDM data of Fig. 2(a), and (b) the WDM data of Fig. 2(b).

The mean differences (WDM - CDM) of the Betti functional summaries are
displayed in Fig. 6 along with the corresponding 95% confidence bands.
Overall, these summaries suggest that the CDM and WDM samples differ on
shorter distance scales, but then start to resemble each other at longer
distance scales in keeping with Fig. 1. Recall that the Betti functions count
the number of features that are persistent at the filtration parameter values
(i.e., the x-axis) so by considering the average difference of the Betti
functions we observe at which scales the number of features differ between the
CDM and WDM. For $H_{0}$, the number of features, on average, for the CDM data
is larger than the number for the WDM for distances until scales of around 0.4
Mpc, and then the number of WDM features is slightly higher than the number of
CDM features until distances of $\sim$0.75 Mpc. The number of $H_{1}$ features
is greater, on average, for the WDM data over the CDM data when $t\leq 0.13$
Mpc, and then the CDM has more $H_{1}$ features until around 0.9 Mpc. A
similar pattern is observed with the $H_{2}$, but the average differences
between the CDM and WDM are within only two $H_{0}$ features.

(a) $\mathcal{F}_{\text{betti}}(0,t)$

(b) $\mathcal{F}_{\text{betti}}(1,t)$

(c) $\mathcal{F}_{\text{betti}}(2,t)$

Figure 6: Mean differences (CDM-WDM) of the Betti functional summaries along
with 95% confidence bands (shaded regions) for the noted functional summaries.
The x-axis limits were set to highlight the non-zero mean differences regions.

While Betti functions capture the number of features that persist across the
filtration parameter values, we defined analogous functions that instead
capture the maximum persistence (MaxPers) and average persistence (AvePers),
which are displayed in Fig. 7(a) and Fig. 7(b), respectively. Similar to the
plots in Fig. 6, the mean difference (CDM-WDM) of these MaxPers and AvePers
functions for the matched samples were computed. However, for Fig. 7, in order
to visualize the variability in the mean differences, pointwise error bars
($\pm$ one standard error) are included. The filtration parameter grid ranges
from 0 to 2.5 Mpc with a spacing of 0.05. This is a lower resolution than the
Betti function figures, which we adopt here in order to be able to improve the
visibility of the individual error bars. There are larger differences between
CDM and WDM MaxPers in $H_{0}$, $H_{1}$, and $H_{2}$ for
$t\stackrel{{\scriptstyle<}}{{{}_{\sim}}}1.85$ Mpc: generally the $H_{0}$
MaxPers are greater for WDM than CDM, the $H_{1}$ MaxPers is greater for CDM
than WDM at scales $\stackrel{{\scriptstyle<}}{{{}_{\sim}}}$1.1 Mpc when this
tendency switches and WDM has greater MaxPers, and the $H_{2}$ MaxPers are
higher for CDM than WDM. A similar pattern is apparent with the AvePers
functions except the $H_{1}$ AvePers are similar for CDM and WDM until scales
around 1 Mpc, after which WDM generally has greater AvePers until around 2
Mpc.

(a) Max persistence

(b) Average persistence

Figure 7: Mean differences (CDM-WDM) of the maximum (a) and average (b)
persistences $\pm$ one standard error. Means and standard errors were computed
every 0.05 Mpc between 0 and 2.5. Gray dotted vertical lines are plotted every
0.10 Mpc. Note that the $H_{0}$ feature that persists indefinitely has been
removed from this analysis.

Basic spatial point process summary functions, such as the 2PCF, are commonly
employed tools in cosmological large-scale structure study. To quantify the
degree to which our persistence diagrams provide new information over these
standard statistics, we calculate and show mean difference functions for the
G-functions and 2PCFs in Fig. 8(a) and Fig. 8(b), respectively. The plotted
data indicate that the WDM functions take, on average, greater values than the
CDM variants for $t\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.5$ Mpc. This
result points toward a similar direction as what we observed with the
differences in the $H_{0}$ Betti function mean differences displayed in Fig.
6(a). This is not surprising since the $H_{0}$ Betti functions, the
G-functions, and the 2PCFs have different ways of assessing the closeness of
the halos within the samples. However, the $H_{1}$ and $H_{2}$ Betti
functions, together with the MaxPers and AvePers methods, appear to detect
differences between the CDM and WDM at different scales, suggesting that they
provide distinct information from the spatial point process functions. In
particular, the $H_{1}$ and $H_{2}$ functions suggest that as the halos become
connected (i.e., the death of $H_{0}$ features), the CDM and WDM models are
forming loops and voids (i.e., $H_{1}$ and $H_{2}$ features, respectively) in
different ways. Also, the MaxPers and AvePers of the $H_{0}$ features differ
between the CDM and WDM data on different scales than those of the $H_{0}$
Betti functions.

(a) $\mathcal{F}_{\text{G}}(t)$

(b) $\mathcal{F}_{\text{pcf}}(t)$

Figure 8: Mean differences (CDM-WDM) of the point process functions along with
95% confidence bands (shaded regions).

##  VI Conclusion

The LSS contains valuable information about the composition, and evolution,
and the physical nature of the Universe. TDA tools such as persistent homology
provide a novel opportunity to extract this information from cosmological
data. While TDA-based approaches have been applied in various fields of
statistical studies, its application to cosmological data and analysis is
still in its infancy. In this paper, we introduced a hypothesis testing
framework built on persistent homology that extends the work of Ref. [31] in
order to compare topological summaries of MW-analog halo neighborhoods (3 Mpc
spheres) evolved under two different DM models: CDM and WDM (WDM thermal relic
mass: 3.3 keV). Next, we have assessed the sensitivity and robustness of this
framework in the context of differentiating between CDM and WDM variants. The
proposed collection of test statistics based on persistence diagrams uses
summaries that were recently proposed in the literature [51, 33, 49, 23, 31,
50], and are easier to work with than the original persistence diagrams. The
results of the persistence diagram-based functional summaries were compared to
two spatial point process functional summaries (G-functions and 2PCF) and test
statistics that use persistence diagrams directly (PDT) [32].

We showed empirically that such a framework is able to infer differences
between CDM and WDM, and investigated the scales at which differences occur.
While most of the test statistics were able to detect statistically
significant differences with $p_{\text{matched}}\leq 0.009$ for all tests
considered except the PDT for $H_{1}$ (§V.1, especially Table 1), the
persistent homology-based functional summaries appear to detect differences
between the CDM and WDM data on different scales from the spatial point
process functional summaries (§V.2).

Our results imply that the homology properties of clustered CDM and WDM haloes
distributions are very different on small scales. CDM haloes are distributed
across a larger number of clusters (homology dimension 0) than WDM haloes,
especially at the $\sim 80$ kpc filtration scale (Fig. 6(a)), although the
clusters that form in WDM are more persistent on average (Fig. 7(b)). The $80$
kpc scale is also where the two process functions return the biggest
difference between the models—in both cases an excess of clustering in CDM
relative to WDM—plus the filtration scale at which WDM features more loops
(homology dimension 1) than CDM. We thus build a picture in which CDM rapidly
builds up a large number of small clusters, whereas WDM builds a smaller
number of clusters, many of which will be rapidly converted into loops.

This picture is consistent with the formation of haloes in and around
cosmological filaments. In CDM, the distribution of filaments extends to near
arbitrarily small scales and fills much of configuration space, whereas the
WDM cutoff restricts WDM haloes to lie along large filaments and so their
spatial distribution is much more constrained. Therefore, the dispersed CDM
haloes form large numbers of small, isolated clusters, whereas WDM haloes are
quickly joined up along cosmological filaments into loops.

The question remains as to whether this difference between the models can be
detected in the spatial distribution of observed Local Group galaxies. One
will have to select haloes that are likely to form a galaxy, where most of the
haloes that we included in this study will be below the HI cooling limit and
thus dark [68, 69]. Reducing the number of haloes available in this manner
will likely lead to a reduction in the statistical significance of differences
between the models’ persistent homology properties: it is therefore imperative
to make halo selections based on, for example, peak halo mass or a semi-
analytic models [70] to confirm the potential persistent homology has for
understanding which DM models best describe the Local Group.

##  Appendix A Distributional Assumptions of DM Samples

In this section, we carry out tests to show that the spatial distributions of
the halo samples do not follow a homogeneous Poisson point process (i.e.,
_complete spatial randomness_ , CSR), which then precludes the use of many
theoretical results that rely on that assumption; see §IV for a brief
discussion about some asymptotic results in persistent homology. Ref. [71]
describe a straightforward Monte Carlo test for checking CSR. Using the same
number of observations (i.e., the number of halos in the MW-analog halo
neighborhoods) and the same window volume (i.e., a sphere with radius 3 Mpc),
$N_{\text{MC}}$ Monte Carlo realizations are generated assuming CSR and then
their G-functions are estimated; see Appendix §B and Equation (23) for
background on G-functions. Then a global envelope is defined using the
$N_{\text{MC}}$ summary functions based on the maximum absolute deviation of
the simulated summary functions from the (known) theoretical summary function
(assuming CSR). If the summary function for the observations are outside the
band, then that is evidence against CSR for those data.

The global envelope was computed by generating $N_{\text{MC}}$ realizations of
a homogeneous Poisson process within a sphere of radius 3 Mpc using rejection
sampling. The number of points was set to match the number of halos in each of
the MW-analog halo neighborhood samples. Then a G-function was estimated for
each sample using the _G3est_ function in the _spatstat_ R package. For each
simulated G-function, the maximum absolute deviation was computed using the
true G-function of the corresponding Poisson process, defined as

$\mathcal{F}_{G}(t)=1-e^{-\frac{4}{3}\pi\hat{\lambda}t^{3}}$ (19)

where $\hat{\lambda}$ is the intensity estimated as the number of points
divided by the volume of the sphere. The interpretation is that if the
observed G-function is outside the envelope for any value $t$ Mpc, then we can
reject CSR at a significance level of $1/(1+N_{\text{MC}})$ (Ch. 10, Ref.
[71]). Using $N_{\text{MC}}=19$, these global envelopes were computed for all
77 samples of the CDM and WDM data, and all observed G-functions have regions
outside the envelopes. For illustration purposes, the resulting global
envelope and observed G-functions for the CDM and WDM MW-analog halo
neighborhood samples from Fig. 2 are displayed below in Fig. 9.

Figure 9: The estimated spatial point process summary functions for the CDM
(dashed blue line) and WDM (dotted red line) MW-analog halo neighborhood
samples displayed in Fig. 2, along with the theoretical summary function
assuming CSR (solid black line) and its global envelope (gray region) using
$N_{\text{MC}}=19$ samples. Because this is a global envelope, we can reject
the hypothesis that the CDM and WDM samples were generated from homogenous
Poisson point process at the $1/(1+N_{\text{MC}})=0.05$ level of significance.

##  Appendix B Comparison Methods in COCO Analysis

In addition to the test statistics proposed based on functional summaries of
persistence diagrams, we include three comparison test statistics in our
investigation of the COCO simulation data presented in §V. The comparison
methods are the Persistence Diagram Test (PDT), and test statistics derived
using the G-function and two-point correlation function (2PCF) which are
popular functional summaries of spatial point processes. The comparison
methods are described below.

#### Persistence Diagram Test (PDT)

Ref. [32] developed a two-sample test that compares persistence diagrams
rather than functional summaries of persistence diagrams. The PDT test
statistic takes the following form,

$\displaystyle\mathcal{T}_{\text{PDT}}$ $\displaystyle(D_{1,1:n_{1}\mid
r},D_{2,1:n_{2}\mid r}\mid p,q)=$ (20)
$\displaystyle\sum_{l=1}^{2}\frac{1}{2n_{l}(n_{l}-1)}\sum_{i=1}^{n_{l}}\sum_{j=1}^{n_{l}}W_{p}(D_{l,i\mid
r},D_{l,j\mid r})^{q}$

where $D_{l,1:n_{l}\mid r}$ is a set of $n_{l}$ persistence diagrams for
homology dimension $r$ from population $l=1,2$, $q$ satisfies $1\leq
q<\infty$, and $W_{p}(\cdot,\cdot)$ is the $p$-Wasserstein distance, with
$1\leq p\leq\infty$. In this work, we set $q=1$ and $p=\infty$. The
$W_{\infty}$ distance is also known as the _bottleneck distance_ , and is
defined as

$W_{\infty}(D_{1},D_{2})=\inf_{\eta:D_{1}\rightarrow D_{2}}\sup_{x\in
D_{1}}\|x-\eta(x)\|_{\infty}$ (21)

where $D_{1}$ and $D_{2}$ are persistence diagrams, $\eta$ defines a bijection
between the two persistence diagrams that allows for matches to the diagonal
$\Delta$, and $\|\cdot\|_{\infty}$ is the $L_{\infty}$ norm in
$\mathbb{R}^{2}$ computed between the birth and death coordinates of $x$ and
$\eta(x)$ for a fixed homology dimension $r$.

#### Spatial Point Process Functions

In order to investigate properties of the spatial distributions of the data,
which in our setting is the location of the DM halos, we consider a popular
functional summary of spatial point processes, the G-function999Also referred
to as the “nearest-neighbor distance distribution function”, along with the
2PCF which is commonly used in cosmology research. The G-function is defined
below and estimated using the implementation for three-dimensional point
patterns in the R package spatstat [72, 73];101010The R function from the
spatstat package is G3est.  see Ref. [71] for more details. The 2PCF uses the
Landy-Szalay estimator [65]

The G-function and 2PCF are estimated for each sample of the DM halos,
$\mathbf{Y}_{k,i}\in\mathbb{R}^{n_{i}\times 3}$ for $k=w,c$ and
$i=1,\ldots,77$. Given a sample $\mathbf{Y}\in\mathbb{R}^{n\times 3}$, let
each point be denoted by $Y_{i}=(Y_{i,1},Y_{i,2},Y_{i,3})$ for $i=1,\ldots,n$.
Define a distance function, $\rho$, as

$\rho(x,\mathbf{A})=\inf\\{\|x-a\|:a\in\mathbf{A}\\}$ (22)

which represents the shortest distance between some point $x\in\mathbb{R}^{3}$
and a closed set $A\subset\mathbb{R}^{3}$. The G-function gives the
distribution function of the nearest neighbor distances, and can be defined as

$\mathcal{F}_{G}(t)=\mathbb{P}(\rho(Y_{i},\mathbf{Y}_{-Y_{i}})\leq t\mid
Y_{i}\in\mathbf{Y})$ (23)

where $\mathbf{Y}_{-Y_{i}}$ is the set of points $\mathbf{Y}$ excluding the
point $Y_{i}$. The Kaplan-Meier estimator of Ref. [74] is used to address the
edge effects (i.e., boundary issues).

##  Appendix C Distribution of p-Values Under the Null Hypothesis

The results of the proposed hypothesis tests are presented in §V.1. Many of
the test statistics find statistically significant differences between the CDM
and WDM models with p-values $\leq 0.001$. In order to verify that the test
statistics do not inappropriately reject the null hypothesis when the null
hypothesis is true (i.e., when both groups come from either CDM or WDM), we
carry out the following experiment. We repeatedly generate two sets of
boostrap realizations from either the CDM or WDM samples, and then compute
permutation p-values for the test statistics presented in the main text. The
distribution of the p-values in this setting where the null hypothesis is true
should follow a uniform distribution. To compute one p-value, two bootstrap
samples (with replacement) of 77 MW-analog halo neighborhoods are selected
from the CDM (WDM) data. Then the hypothesis testing framework presented in
§IV.1 is used to compute a traditional permutation p-value (since the matched
pairs design is not present in this setting) using 20,000 permutations. This
computation is repeated for 100 independent iterations for the CDM (WDM) data
with the same sampled indexes used for the CDM and WDM bootstrap samples. Fig.
10 and Fig. 11 display the results for the CDM and WDM samples, respectively,
as uniform quantile-quantile plots with 99% pointwise bands based on the
distribution of order statistics of uniform random variables (i.e., Beta($k$,
$n+1-k$) where $k$ is the order and $n=100$). The resulting p-values for each
test statistic are generally consistent with uniform distributions. The CDM
$H_{1}$ landscape function p-values (Fig. 10(b)) have some values that are not
within the 99% confidence band, but this does not occur with the WDM $H_{1}$
landscape function p-values (Fig. 11(b)) nor with the other landscape function
p-values so it appears to not be a reason for concern about the landscape
function-based test statistics.

(a) Betti functions

(b) Landscape functions

(c) Silhouette functions

(d) Spatial functions

Figure 10: Uniform quantile-quantile plots of the 100 permutations p-values
calculated for each test statistic using bootstrap realizations of the CDM MW-
analog halo samples. Each bootstrap sample includes 77 MW-analog halo
neighborhoods, and 20,000 permutations were used to compute each p-value.

(a) Betti functions

(b) Landscape functions

(c) Silhouette functions

(d) Spatial functions

Figure 11: Uniform quantile-quantile plots of the 100 permutations p-values
calculated for each test statistic using bootstrap realizations of the WDM MW-
analog halo samples. Each bootstrap sample includes 77 MW-analog halo
neighborhoods, and 20,000 permutations were used to compute each p-value.

###### Acknowledgements.

This research was performed using the compute resources and assistance of the
UW-Madison Center For High Throughput Computing (CHTC) in the Department of
Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing
Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes
for Discovery, and the National Science Foundation, and is an active member of
the Open Science Grid, which is supported by the National Science Foundation
and the U.S. Department of Energy’s Office of Science. This work used the
DiRAC Data Centric system at Durham University, operated by the Institute for
Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure
capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1,
STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part
of the National E-Infrastructure. This project has also benefited from
numerical computations performed at the Interdisciplinary Center for
Mathematical and Computational Modeling (ICM) University of Warsaw under
grants #no GB79-7, GA67-17 and G63-3. JCK and BTF acknowledge support from NSF
under Grant Numbers DMS 2038556 and 1854336. WAH and PD acknowledge the
support from the Polish National Science Center within research projects no.
2018/31/G/ST9/03388, 2020/39/B/ST9/03494. MRL acknowledges support by a Grant
of Excellence from the Icelandic Research Fund (grant number 206930).

## References

  * Davis  _et al._ [1985] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. White, The evolution of large-scale structure in a universe dominated by cold dark matter, The Astrophysical Journal 292, 371 (1985).
  * Bull  _et al._ [2016] P. Bull, Y. Akrami, J. Adamek, T. Baker, E. Bellini, J. B. Jimenez, E. Bentivegna, S. Camera, S. Clesse, J. H. Davis, _et al._ , Beyond $\Lambda$CDM: Problems, solutions, and the road ahead, Physics of the Dark Universe 12, 56 (2016).
  * Bullock and Boylan-Kolchin [2017] J. S. Bullock and M. Boylan-Kolchin, Small-scale challenges to the $\Lambda$CDM paradigm, Annual Review of Astronomy and Astrophysics 55, 343 (2017).
  * Perivolaropoulos and Skara [2021] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update, arXiv preprint arXiv:2105.05208 (2021).
  * van de Weygaert  _et al._ [2011] R. van de Weygaert, G. Vegter, H. Edelsbrunner, B. J. T. Jones, P. Pranav, C. Park, W. A. Hellwing, B. Eldering, N. Kruithof, E. G. P. P. Bos, J. Hidding, J. Feldbrugge, E. ten Have, M. van Engelen, M. Caroli, and M. Teillaud, Alpha, Betti and the Megaparsec Universe: On the Topology of the Cosmic Web, in  _Lecture Notes in Computer Science_, Vol. 6970 (Springer, 2011) pp. 60–101.
  * Sánchez  _et al._ [2012] A. G. Sánchez, C. Scóccola, A. Ross, W. Percival, M. Manera, F. Montesano, X. Mazzalay, A. Cuesta, D. Eisenstein, E. Kazin, _et al._ , The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological implications of the large-scale two-point correlation function, Monthly Notices of the Royal Astronomical Society 425, 415 (2012).
  * Papastergis  _et al._ [2011] E. Papastergis, A. M. Martin, R. Giovanelli, and M. P. Haynes, The Velocity Width Function of Galaxies from the 40% ALFALFA Survey: Shedding Light on the Cold Dark Matter Overabundance Problem, ApJ 739, 38 (2011), arXiv:1106.0710 [astro-ph.CO] .
  * Kennedy  _et al._ [2014] R. Kennedy, C. Frenk, S. Cole, and A. Benson, Constraining the warm dark matter particle mass with Milky Way satellites, MNRAS 442, 2487 (2014), arXiv:1310.7739 [astro-ph.CO] .
  * Tikhonov  _et al._ [2009] A. V. Tikhonov, S. Gottlöber, G. Yepes, and Y. Hoffman, The sizes of minivoids in the local Universe: an argument in favour of a warm dark matter model?, Monthly Notices of the Royal Astronomical Society 399, 1611 (2009).
  * Lovell  _et al._ [2021] M. R. Lovell, M. Cautun, C. S. Frenk, W. A. Hellwing, and O. Newton, The spatial distribution of Milky Way satellites, gaps in streams, and the nature of dark matter, MNRAS 507, 4826 (2021), arXiv:2104.03322 [astro-ph.GA] .
  * Watts  _et al._ [2017] A. L. Watts, P. J. Elahi, G. F. Lewis, and C. Power, Large-scale structure topology in non-standard cosmologies: impact of dark sector physics, MNRAS 468, 59 (2017), arXiv:1702.03066 [astro-ph.CO] .
  * Planck Collaboration [2014] Planck Collaboration, Planck 2013 results. XVI. Cosmological parameters, A&A 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Eisenstein  _et al._ [2005] D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, R. C. Nichol, R. Scranton, H.-J. Seo, M. Tegmark, Z. Zheng, S. F. Anderson, J. Annis, N. Bahcall, J. Brinkmann, S. Burles, F. J. Castander, A. Connolly, I. Csabai, M. Doi, M. Fukugita, J. A. Frieman, K. Glazebrook, J. E. Gunn, J. S. Hendry, G. Hennessy, Z. Ivezić, S. Kent, G. R. Knapp, H. Lin, Y.-S. Loh, R. H. Lupton, B. Margon, T. A. McKay, A. Meiksin, J. A. Munn, A. Pope, M. W. Richmond, D. Schlegel, D. P. Schneider, K. Shimasaku, C. Stoughton, M. A. Strauss, M. SubbaRao, A. S. Szalay, I. Szapudi, D. L. Tucker, B. Yanny, and D. G. York, Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies, ApJ 633, 560 (2005), arXiv:astro-ph/0501171 [astro-ph] .
  * LUX Collaboration [2017] LUX Collaboration, Results from a Search for Dark Matter in the Complete LUX Exposure, Phys. Rev. Lett. 118, 021303 (2017), arXiv:1608.07648 [astro-ph.CO] .
  * Xenon Collaboration [2018] Xenon Collaboration, Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Fermi-LAT Collaboration and DES Collaboration [2017] Fermi-LAT Collaboration and DES Collaboration, Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-Lat, ApJ 834, 110 (2017), arXiv:1611.03184 [astro-ph.HE] .
  * Colombi  _et al._ [1996] S. Colombi, S. Dodelson, and L. Widrow, Large-scale structure tests of warm dark matter, Astrophysical Journal 458 (1996).
  * Hellwing  _et al._ [2016] W. A. Hellwing, C. S. Frenk, M. Cautun, S. Bose, J. Helly, A. Jenkins, T. Sawala, and M. Cytowski, The Copernicus Complexio: a high-resolution view of the small-scale Universe, Monthly Notices of the Royal Astronomical Society 457, 3492 (2016).
  * Bose  _et al._ [2016] S. Bose, W. A. Hellwing, C. S. Frenk, A. Jenkins, M. R. Lovell, J. C. Helly, and B. Li, The Copernicus Complexio: statistical properties of warm dark matter haloes, Monthly Notices of the Royal Astronomical Society 455, 318 (2016).
  * Sousbie [2011] T. Sousbie, The persistent cosmic web and its filamentary structure - I. Theory and implementation, Monthly Notices of the Royal Astronomical Society 414, 350 (2011).
  * Sousbie  _et al._ [2011] T. Sousbie, C. Pichon, and H. Kawahara, The persistent cosmic web and its filamentary structure – II. Illustrations, Monthly Notices of the Royal Astronomical Society 414, 384 (2011).
  * Cisewski  _et al._ [2014] J. Cisewski, R. A. Croft, P. E. Freeman, C. R. Genovese, N. Khandai, M. Ozbek, and L. Wasserman, Non-parametric 3D map of the intergalactic medium using the Lyman-alpha forest, Monthly Notices of the Royal Astronomical Society 440, 2599 (2014).
  * Pranav  _et al._ [2017] P. Pranav, H. Edelsbrunner, R. Van de Weygaert, G. Vegter, M. Kerber, B. J. Jones, and M. Wintraecken, The topology of the cosmic web in terms of persistent Betti numbers, Monthly Notices of the Royal Astronomical Society 465, 4281 (2017).
  * Green  _et al._ [2019] S. B. Green, A. Mintz, X. Xu, and J. Cisewski-Kehe, Topology of our cosmology with persistent homology, CHANCE 32, 6 (2019).
  * Pranav  _et al._ [2019] P. Pranav, R. Van de Weygaert, G. Vegter, B. J. Jones, R. J. Adler, J. Feldbrugge, C. Park, T. Buchert, and M. Kerber, Topology and geometry of Gaussian random fields I: on Betti numbers, Euler characteristic, and Minkowski functionals, Monthly Notices of the Royal Astronomical Society 485, 4167 (2019).
  * Xu  _et al._ [2019] X. Xu, J. Cisewski-Kehe, S. B. Green, and D. Nagai, Finding cosmic voids and filament loops using topological data analysis, Astronomy and Computing 27, 34 (2019).
  * Cole  _et al._ [2020] A. Cole, M. Biagetti, and G. Shiu, Topological Echoes of Primordial Physics in the Universe at Large Scales, arXiv preprint arXiv:2012.03616 (2020).
  * Duong  _et al._ [2012] T. Duong, B. Goud, and K. Schauer, Closed-form density-based framework for automatic detection of cellular morphology changes, Proceedings of the National Academy of Sciences 109, 8382 (2012).
  * Bendich  _et al._ [2016] P. Bendich, J. Marron, E. Miller, A. Pieloch, and S. Skwerer, Persistent homology analysis of brain artery trees, The Annals of Applied Statistics 10, 198 (2016).
  * Lawson  _et al._ [2019] P. Lawson, A. B. Sholl, J. Q. Brown, B. T. Fasy, and C. Wenk, Persistent homology for the quantitative evaluation of architectural features in prostate cancer histology, Scientific reports 9, 1 (2019).
  * Berry  _et al._ [2020] E. Berry, Y.-C. Chen, J. Cisewski-Kehe, and B. T. Fasy, Functional summaries of persistence diagrams, Journal of Applied and Computational Topology 4, 211 (2020).
  * Robinson and Turner [2017] A. Robinson and K. Turner, Hypothesis testing for topological data analysis, Journal of Applied and Computational Topology 1, 241 (2017).
  * Bubenik [2015] P. Bubenik, Statistical topological data analysis using persistence landscapes, Journal of Machine Learning Research 16, 77 (2015).
  * Biscio and Møller [2019] C. A. Biscio and J. Møller, The accumulated persistence function, a new useful functional summary statistic for topological data analysis, with a view to brain artery trees and spatial point process applications, Journal of Computational and Graphical Statistics 28, 671 (2019).
  * Krebs and Hirsch [2021] J. Krebs and C. Hirsch, Functional central limit theorems for persistent Betti numbers on cylindrical networks, Scandinavian Journal of Statistics (2021).
  * Springel  _et al._ [2001] V. Springel, S. D. White, G. Tormen, and G. Kauffmann, Populating a cluster of galaxies–I. Results at z= 0, Monthly Notices of the Royal Astronomical Society 328, 726 (2001).
  * Viel  _et al._ [2013] M. Viel, G. D. Becker, J. S. Bolton, and M. G. Haehnelt, Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman-$\alpha$ forest data, Physical Review D 88, 043502 (2013).
  * Lovell  _et al._ [2014] M. R. Lovell, C. S. Frenk, V. R. Eke, A. Jenkins, L. Gao, and T. Theuns, The properties of warm dark matter haloes, Monthly Notices of the Royal Astronomical Society 439, 300 (2014).
  * Munkres [1984] J. R. Munkres,  _Elements of algebraic topology_ , Vol. 2 (Addison-Wesley Menlo Park, 1984).
  * Hatcher [2002] A. Hatcher,  _Algebraic topology_ (Cambridge University Press, 2002).
  * Edelsbrunner and Harer [2010] H. Edelsbrunner and J. Harer,  _Computational topology: an introduction_ (American Mathematical Soc., 2010).
  * Zomorodian [2010] A. Zomorodian, Fast construction of the Vietoris-Rips complex, Computers & Graphics 34, 263 (2010).
  * Niyogi  _et al._ [2008] P. Niyogi, S. Smale, and S. Weinberger, Finding the homology of submanifolds with high confidence from random samples, Discrete & Computational Geometry 39, 419 (2008).
  * Cohen-Steiner  _et al._ [2007] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer, Stability of persistence diagrams, Discrete & Computational Geometry 37, 103 (2007).
  * Chazal  _et al._ [2016] F. Chazal, V. De Silva, M. Glisse, and S. Oudot,  _The Structure and Stability of Persistence Modules_ (Springer, 2016).
  * Fasy  _et al._ [2014] B. T. Fasy, F. Lecci, A. Rinaldo, L. Wasserman, S. Balakrishnan, A. Singh, _et al._ , Confidence sets for persistence diagrams, The Annals of Statistics 42, 2301 (2014).
  * Turner  _et al._ [2014] K. Turner, Y. Mileyko, S. Mukherjee, and J. Harer, Fréchet means for distributions of persistence diagrams, Discrete & Computational Geometry 52, 44 (2014).
  * Munch  _et al._ [2015] E. Munch, K. Turner, P. Bendich, S. Mukherjee, J. Mattingly, J. Harer, _et al._ , Probabilistic Fréchet means for time varying persistence diagrams, Electronic Journal of Statistics 9, 1173 (2015).
  * Bubenik and Dłotko [2017] P. Bubenik and P. Dłotko, A persistence landscapes toolbox for topological statistics, Journal of Symbolic Computation 78, 91 (2017).
  * Bubenik [2020] P. Bubenik, The persistence landscape and some of its properties, in  _Topological Data Analysis_ (Springer, 2020) pp. 97–117.
  * Chazal  _et al._ [2014] F. Chazal, B. T. Fasy, F. Lecci, A. Rinaldo, and L. Wasserman, Stochastic convergence of persistence landscapes and silhouettes, in  _Proceedings of the thirtieth annual symposium on Computational geometry_ (ACM, 2014) p. 474.
  * Mecke  _et al._ [1993] K. R. Mecke, T. Buchert, and H. Wagner, Robust morphological measures for large-scale structure in the universe, arXiv preprint astro-ph/9312028 (1993).
  * Park  _et al._ [2013] C. Park, P. Pranav, P. Chingangbam, R. Van De Weygaert, B. Jones, G. Vegter, I. Kim, J. Hidding, and W. A. Hellwing, Betti numbers of gaussian fields, Journal of The Korean Astronomical Society 46, 125 (2013).
  * Kimura and Imai [2017] Y. Kimura and K. Imai, Quantification of lss using the persistent homology in the sdss fields, Advances in Space Research 60, 722 (2017).
  * Giri and Mellema [2021] S. K. Giri and G. Mellema, Measuring the topology of reionization with Betti numbers, Monthly Notices of the Royal Astronomical Society 505, 1863 (2021).
  * Wilding  _et al._ [2021] G. Wilding, K. Nevenzeel, R. van de Weygaert, G. Vegter, P. Pranav, B. J. Jones, K. Efstathiou, and J. Feldbrugge, Persistent homology of the cosmic web–I. hierarchical topology in $\Lambda$CDM cosmologies, Monthly Notices of the Royal Astronomical Society 507, 2968 (2021).
  * Bauer and Pausinger [2018] U. Bauer and F. Pausinger, Persistent Betti numbers of random Cech complexes, arXiv preprint arXiv:1801.08376 (2018).
  * Hiraoka  _et al._ [2018] Y. Hiraoka, T. Shirai, K. D. Trinh, _et al._ , Limit theorems for persistence diagrams, Annals of Applied Probability 28, 2740 (2018).
  * Krebs and Polonik [2019] J. T. Krebs and W. Polonik, On the asymptotic normality of persistent Betti numbers, arXiv preprint arXiv:1903.03280 (2019).
  * Biscio  _et al._ [2020] C. A. Biscio, N. Chenavier, C. Hirsch, A. M. Svane, _et al._ , Testing goodness of fit for point processes via topological data analysis, Electronic Journal of Statistics 14, 1024 (2020).
  * Chen  _et al._ [2015] Y.-C. Chen, D. Wang, A. Rinaldo, and L. Wasserman, Statistical analysis of persistence intensity functions, arXiv preprint arXiv:1510.02502 (2015).
  * Adams  _et al._ [2017] H. Adams, T. Emerson, M. Kirby, R. Neville, C. Peterson, P. Shipman, S. Chepushtanova, E. Hanson, F. Motta, and L. Ziegelmeier, Persistence images: A stable vector representation of persistent homology, Journal of Machine Learning Research 18 (2017).
  * Yogeshwaran  _et al._ [2017] D. Yogeshwaran, E. Subag, and R. J. Adler, Random geometric complexes in the thermodynamic regime, Probability Theory and Related Fields 167, 107 (2017).
  * Mileyko  _et al._ [2011] Y. Mileyko, S. Mukherjee, and J. Harer, Probability measures on the space of persistence diagrams, Inverse Problems 27, 124007 (2011).
  * Landy and Szalay [1993] S. D. Landy and A. S. Szalay, Bias and variance of angular correlation functions, The Astrophysical Journal 412, 64 (1993).
  * Peebles [1980] P. Peebles, The large-scale structure of the universe, Large-Scale Structure of the Universe by Phillip James Edwin Peebles. Princeton University Press (1980).
  * Bauer [2021] U. Bauer, Ripser: efficient computation of Vietoris-Rips persistence barcodes, Journal of Applied and Computational Topology 10.1007/s41468-021-00071-5 (2021).
  * Benson  _et al._ [2002] A. J. Benson, C. G. Lacey, C. M. Baugh, S. Cole, and C. S. Frenk, The effects of photoionization on galaxy formation - I. Model and results at z=0, MNRAS 333, 156 (2002), arXiv:astro-ph/0108217 [astro-ph] .
  * Sawala  _et al._ [2016] T. Sawala, C. S. Frenk, A. Fattahi, J. F. Navarro, T. Theuns, R. G. Bower, R. A. Crain, M. Furlong, A. Jenkins, M. Schaller, and J. Schaye, The chosen few: the low-mass haloes that host faint galaxies, MNRAS 456, 85 (2016), arXiv:1406.6362 [astro-ph.CO] .
  * Bose  _et al._ [2017] S. Bose, W. A. Hellwing, C. S. Frenk, A. Jenkins, M. R. Lovell, J. C. Helly, B. Li, V. Gonzalez-Perez, and L. Gao, Substructure and galaxy formation in the Copernicus Complexio warm dark matter simulations, MNRAS 464, 4520 (2017), arXiv:1604.07409 [astro-ph.CO] .
  * Baddeley  _et al._ [2015] A. Baddeley, E. Rubak, and R. Turner,  _Spatial point patterns: methodology and applications with R_ (CRC press, 2015).
  * Baddeley  _et al._ [2004] A. J. Baddeley, R. Turner, _et al._ , Spatstat: An R package for analyzing spatial point patterns (2004).
  * Baddeley  _et al._ [2014] A. Baddeley, R. Turner, E. Rubak, and K. K. Berthelsen, Package ‘spatstat’, The Comprehensive R Archive Network (2014).
  * Baddeley and Gill [1997] A. Baddeley and R. D. Gill, Kaplan-Meier estimators of distance distributions for spatial point processes, The Annals of Statistics , 263 (1997).