# DarkMix: Mixture Models for the Detection and Characterization of Dark
Matter Halos

Lluís Hurtado-Gil eDreams ODIGEO
C/ Bailèn 67-69,
08009 Barcelona, Spain.  Observatori Astronòmic
Universitat de València
C/ Catedràtic José Beltrán 2
E-46980, Paterna, Spain  Michael A. Kuhn California Institute of Technology
Pasadena, CA 91125, USA  Pablo Arnalte-Mur Observatori Astronòmic
Universitat de València
C/ Catedràtic José Beltrán 2
E-46980, Paterna, Spain  Departament d’Astronomia i Astrofísica
Universitat de València
E-46100, Burjassot, Spain  Eric D. Feigelson Department of Astronomy &
Astrophysics,
Penn State University
University Park, PA 16802, USA  Center for Astrostatistics,
Penn State University, University Park, PA 16802, USA  Vicent Martínez
Observatori Astronòmic
Universitat de València
C/ Catedràtic José Beltrán 2
E-46980, Paterna, Spain  Departament d’Astronomia i Astrofísica
Universitat de València
E-46100, Burjassot, Spain  Unidad Asociada Observatorio Astronómico (IFCA-UV)
E-46980, Paterna, Spain

(Accepted July 26, 2022)

###### Abstract

Dark matter simulations require statistical techniques to properly identify
and classify their halos and structures. Nonparametric solutions provide
catalogs of these structures but lack the additional learning of a model-based
algorithm and might misclassify particles in merging situations. With mixture
models, we can simultaneously fit multiple density profiles to the halos that
are found in a dark matter simulation. In this work, we use the Einasto
profile (Einasto, 1965, 1968, 1969) to model the halos found in a sample of
the Bolshoi simulation (Klypin et al., 2011), and we obtain their location,
size, shape and mass. Our code is implemented in the R statistical software
environment and can be accessed on https://github.com/LluisHGil/darkmix.

Dark matter distribution (356), Galaxy dark matter halos (1880), Spatial point
processes (1915), Mixture model (1932)

††journal: ApJ

##  1 Introduction

###  1.1 Dark Matter Halos

In 1952, Neyman & Scott proposed the first statistical model of large-scale
galaxy distribution (Neyman & Scott, 1952, 1953, 1954). This model interpreted
luminous matter to be distributed in the universe according to a stochastic
process, where a discrete distribution of galaxies is aggregated into
overdense clusters that themselves are distributed in space as a Poisson
process. This model needs three main components to be built: the distribution
of galaxies within each cluster, the size distribution of the clusters, and a
description of the clustering of clusters. While their original formulation
was too simple, this framework is still a valid approach for the baryonic
distribution in the universe and it can be extended to dark matter
distribution (Cooray & Sheth, 2002). Here, dark matter is a continuous field
of particles that collapses from overdensities into strongly clustered
structures, which are called halos. The identification and morphology of these
halos is the main topic of the present work.

Analytic models and numeric simulations show how the initial dark matter
field, which is made of particles, evolves from an initially smooth state to a
highly clustered final condition, which results in a complex cosmic web of
knots (the halos), filaments, sheets and voids (Bertschinger, 1985).
Simulations show that the halo mass (Moore et al., 1999; Navarro et al.,
1996), abundance and spatial distribution (Church et al., 1997; Jenkins et
al., 2001) are highly dependent on the initial conditions. The final structure
within a halo can be reasonably assumed to be in virial equilibrium. The
nature and evolution of the galaxies that have formed inside halos is strongly
dependent on the parent halo’s properties (Kauffmann et al., 1999; Somerville
& Primack, 1999; Benson et al., 2001; Colberg et al., 2000), with more
galaxies formed in more massive and clustered halos. However, the galaxy
distribution is biased toward stronger clustering conditions (Madgwick et al.,
2003; Zehavi et al., 2011; Hurtado-Gil et al., 2016).

The spherical collapse model is a classic approximation for initial conditions
leading to dark matter halos (Gunn & Gott, 1972; Fillmore & Goldreich, 1984;
Bertschinger, 1985). Here, dark matter overdensity collapses from a tophat
density perturbation into a halo that, depending on the overdensity mass and
density, virializes when a certain size is reached. The final density of the
halo is much higher that the prediction of a linear model because the
evolution of the halo clustering is governed by nonlinear processes, where the
densest regions are populated by the most massive halos (Bardeen, 1980;
Kaiser, 1984) and the dark matter follows a log-normal distribution (Coles &
Jones, 1991; Hurtado-Gil et al., 2017; Clerkin et al., 2017).

The shape of the collapsed halos, Bertschinger (1985) and Fillmore & Goldreich
(1984) suggests that their density profile around the center depends on the
initial density distribution of the parent overdense region. Although more
massive halos arise from denser peaks in the initial fluctuation field
(Kaiser, 1984; Hoffman & Shaham, 1985), these dense peaks are also less
centrally concentrated (Bardeen, 1980). Massive virialized halos are thereby
less centrally concentrated than less massive halos (Navarro et al., 1996).

Several halo density profiles have been proposed, including the Hernquist
(1990), Navarro-Frenk-White Navarro et al. (1996, 1997), and Einasto (1968)
profile. In sections 2.1 and 2.2, we introduce the two later profiles and we
justify our final selection of the Einasto profile.

###  1.2 Finding Structure with Mixture Models and Other Methods

The growth in volume and detail of astronomical observations and simulations
requires automated tools to characterize the properties and evolution of halo
structures. These tools should be robust and objective, and should not depend
on heuristic choices or subjective judgment.

Most of the widely-used clustering algorithms for galaxy clustering analysis
are nonparametric and are based on dissimilarity distances, which is a metric
that is used to decide if two particles are sufficiently close to belong to
the same cluster. Astronomers commonly use the Friends-of-Friends algorithm,
which is known in statistics as single-linkage hierarchical agglomerative
clustering (Gower & Ross, 1969). This method is highly performant and is
scalable, which is crucial for large volume data sets, such as dark matter
numerical simulations (Willis et al., 2020). However, this method is prone to
‘chaining’ of unrelated groupings (Everitt et al., 2011). The resulting
clusters depend strongly on the choice of density or size threshold value.
Other nonparametric clustering procedures—such as Ward’s hierarchical
clustering, $k$-means algorithm, DBSCAN, kernel density estimation bump
hunting, and their many extensions—also have heuristic thresholds or stopping
rules with a corresponding loss of objectivity and robustness. These
procedures also have the limitation of giving ‘hard’ classifications where
data points are strictly classified on one cluster or another without regard
to the reliability of this decision. Therefore, the characterization of
merging or blended clusters is limited. In contrast, parametric methods assume
a particular shape to the structures in the population. They also have the
advantage that modeling can be based on maximum likelihood estimation or (if
prior information is available) Bayesian inference, without requiring the
addition of heuristic thresholds or stopping rules. These methods are
typically based on probability density functions. This allows us to perform
statistical tests (e.g., significance tests on the parameters) and goodness of
fit validations (e.g., the coefficient of determination) (Rao et al., 1973).
Furthermore, they give ‘soft’ probabilities for each point belonging to each
cluster. Hard membership classifications can be decided afterwards using
heuristic decision rules.

These benefits motive us to use finite mixture models to detect and
characterize dark matter structures (Peel & MacLahlan, 2000; McLachlan & Peel,
2000; Everitt, 2005; Frühwirth-Schnatter, 2006; McLachlan & Krishnan, 2007;
Everitt et al., 2011). This technique has been widely used in astronomy and
astrophysics with considerable success (Kuhn et al., 2014; Fruhwirth-Schnatter
et al., 2019; Kuhn et al., 2019).

We understand the dark matter distribution to be structured in halos, which
can be described by a parametric surface-density distribution. The halos and
other structures are described by the mixture model components, which are
summed to create the surface density function. Mixture models are widely used
in parametric modeling of point processes, usually with Gaussian function
components (Fraley & Raftery, 2002). As previously explained, we will instead
use an astrophysically motivated function, such as the Einasto profile.

The mixture model is then estimated following a three step method: first, the
number of components is chosen by the user; second, the properties of these
components are obtained through maximum likelihood parameter estimation (MLE);
and third, the component assignments for dark matter particles are determined
using posterior probabilities from the fitted models (Everitt et al., 2011).

Mixture models are dependent on three main choices: the chosen number of
components, the convergence criteria towards the surface density distribution,
and the selected profile. The first is assessed using model selection, such as
Bayesian and Akaike Information criteria (Schwarz et al., 1978; Akaike, 1998).
The model goodness-of-fit will be measured using well-known statistics, such
as the coefficient of determination, or residual maps between the input data
and model predictions (Kuhn et al., 2014; de Souza et al., 2017). We will
justify the election of the Einasto profile in section 2.2. Once the model is
estimated, we may want to give a final particle classification in the
clusters. Mixture models allow for a probability based classification, but as
we explain in section 7.2 we recommend the application of a heuristic
threshold for background particles. This threshold should be based on the
merging conditions of our data.

Aside from the mixture model alternatives, several algorithms are widely used
by the community to detect structures in dark matter N-body simulations and
classify their particles. We mention the Bound Density Maximum algorithm
(Klypin & Holtzman, 1997; Riebe et al., 2013) in section 7.4, which is based
on density maximums and spherical halos. Another alternative is the ROCKSTAR
algorithm (Behroozi et al., 2013), ‘based on adaptive hierarchical refinement
of friends-of-friends groups in six phase-space dimensions and one time
dimension’. With this method, the algorithm can provide an effective particle
classification and can even detect small subhalos in merging conditions. The
halo finder VELOCIraptor (Elahi et al., 2019) has also been proven to be able
to ‘identify sub-halos deep within the host that have negligible density
contrasts to their parent halo’. Both methods are meant to be used on large
N-body simulations, providing a full catalog of halos, sub-halos and even
tidal features.

In contrast, our method is a parametric based approach. Although this limits
its use to small-sized samples of particles, it provides a parametric fitting
of the halos density profile. In this work, we will focus on the advantages of
obtaining such a profile-based description.

This paper is organized as follows. Section 2 presents the finite mixture
modelling technique, the Einasto profile and different applications for the
estimated results. Section 3.1 illustrates the MLE calculation and section 3.2
describes the tools that we used to validate the best-fit model. In section 4
we summarize the steps of our code darkmix and we apply it in section 5 to a
set of generated realizations of a simulated dark matter distribution with
Einasto profile and fit them with our software to validate it. Section 6
presents a data set from the Bolshoi simulation, and section 7 shows the
analysis and results of our mixture modeling and validation. In section 8, we
summarize our main conclusions and outline future work.

##  2 Finite Mixture Models for Dark Matter Halos

Finite mixture densities are a family of probability density distributions
that combine multiple components into a single probability function. Each
component has a probability density function (e.g., the Einasto profile) and
the final mixture model is the weighted sum of $c$ components. All components
are continuous, and can be evaluated at any location of the data space,
occupied or not by a particle. Mixture models admit as many different
components as desired, as long as they can be defined as probability
distributions. In this work we will consider two kinds of components: $k$
halos, and a single background component containing all of the particles that
are not associated with a halo. The total number of the mixture model is
$c=k+1$.

Given a data sample of $N$ dark matter particles from a simulation, we define
$\mathbf{X}$ as a three column matrix containing the coordinates of our
particles. Each component will be modeled by a profile function
$\rho_{j}(\mathbf{r}_{i},\vec{\theta}_{j})$, $j=1,\dotsc,c$, where
$\vec{\theta}_{j}$ is the parameters vector for component $j$ and $\Theta$ is
the matrix of the $c$ vectors $\vec{\theta}_{j}$. The sum of these components
is weighted by the mixing proportions $\vec{w}=\\{w_{j}\\}$, $j=1,\dotsc,c$,
which are non-negative. In point process statistics, the model is defined over
a window $W$ containing the sample of $N$ points
$\mathbf{X}=\\{\mathbf{r}_{i}\\}$, where $\mathbf{r}_{i}=(x_{i},y_{i},z_{i})$.
Together, the finite mixture model $\Sigma$ can be written as

$\Sigma(\mathbf{X}|\vec{w},\Theta)=\sum_{j=1}^{c}w_{j}\cdot\rho_{j}(\mathbf{X}|\vec{\theta}_{j})=\sum_{j=1}^{c}\sum_{i=1}^{N}w_{j}\cdot\rho_{j}(\mathbf{r}_{i}|\vec{\theta}_{j})$
(1)

Our implementation of the mixture model for the dark matter distribution will
follow the strategy of (Kuhn et al., 2014), and therefore we chose our
notation after this paper.

Functions $\rho_{j}$ represent the profile functions of our dark matter
structures, which are summed into the surface density function $\Sigma$ and
the probability density function that we use as a model and will fit to our
data. These functions typically create a multimodal distribution that matches
the clusters that are present in our data. Therefore, individual dark matter
particles are used to obtain the overall density distribution of the structure
that they belong to. No learning can be obtained from the internal
distribution (i.e., the relative position of the particles inside their
structure). The total finite mixture model $\Sigma$ is the weighted sum of
these components and models the dark matter density distribution as generated
by the N-body simulation.

The profile of a mixture model component can be as irregular as we are able to
model it. However, since this method is generally used as a cluster
classification method, we will model two different components—the halos and
the background.

We will introduce the NFW and the Einasto profiles in the next section.
However, as we justify, we will only make use of the latter. We will also
introduce the final version of our mixture model before deducing several
quantities and functions of interest.

###  2.1 The NFW Profile

The NFW profile (Navarro et al., 1996, 1997) has the following shape

$\rho(r|M)=\frac{\rho_{s}(M)}{\left[c(M)\frac{r}{r_{\rm
vir}(M)}\right]\left[1+c(M)\frac{r}{r_{\rm vir}(M)}\right]^{2}}\;,$ (2)

truncated at the virial radius $r=r_{\rm vir}(M)$. This profile has a
logarithmic slope of $-1$ at small scales ($r\ll r_{\rm vir}/c$) and of $-3$
at large scales ($r\gg r_{\rm vir}/c$). Here, $\rho_{s}(M)$ is a normalization
factor that is fixed from the condition that the total mass integrated to
$r_{\rm vir}(M)$ must be equal to $M$ and $c(M)$ is the concentration
parameter. We note that the NFW profile or power-law profiles cannot be used
in this context because they do not have a finite integral. Instead, we will
use the Einasto profile for our study.

###  2.2 The Einasto Profile

The Einasto profile (Einasto, 1965, 1968, 1969) for three-dimensional density
profiles is similar to the Sérsic profile used for two-dimensional galaxy
brightness profiles (Sérsic, 1963; Sersic, 1968). Assuming spherical symmetry,
the Einasto profile describes the density $\rho$ of halo $j$ as a function of
distance from the halo center $\mathbf{r_{0,j}}$. The Sérsic index $n$ is a
free parameter that is used to measure the shape of the halo density profile:
larger values of $n$ create centrally concentrated profiles. As we will see in
section 7, large values of this parameter (above $n\sim 8$) degenerate the
profile into a power law-like function. The other free parameter of the
profile is Einasto’s radius $r_{e}$, which defines a volume containing half of
the total mass that can be used to understand the size of the halo.

Given these parameters, we can define $d_{n}$, which is a function of $n$ such
that $\rho_{e}$ is the density at the radius $r_{e}$, and we have defined
$\rho_{e}$. The factor $d_{n}$ can be obtained by solving

$\Gamma(3n)=2\gamma(3n,d_{n})$ (3)

where $\Gamma$ is the complete gamma function and $\gamma$ is the lower
incomplete gamma function111The gamma functions are extensions of the
factorial function to complex numbers..

We are now ready to define the Einasto profile as

$\rho(r)=\rho_{e}\exp{\Big{(}-d_{n}\Big{[}(r/r_{e})^{1/n}-1\Big{]}\Big{)}}$
(4)

where $r=||\mathbf{r-r_{0}}||$ is the distance between a chosen location
$\mathbf{r}$ and the center of the halo $\mathbf{r_{0}}$. As we see in eq. 1,
this profile will be multiplied by parameter $p_{j}$, which is the mixture
coefficient. This makes it impossible for us to estimate $p_{j}$ and $p_{e}$
separately. The results that we provide for $p_{j}$ in Tables 3 and 5 should
be read as $p_{j}\cdot p_{e}$.

Before moving on to the next section, we define the profile of the background
component as a constant density, having value 1 at all locations of the window
$W$. Since the background profile is multiplied by its our mixture coefficient
$p_{b}$, we expect it to be close to the mean density of the data set. In a
large data set, where conditions are comparable to the universe, we would
expect it to be around the mean density of the universe. This is equivalent to
defining the background as a homogeneous Poisson distributed point process,
with profile function

$\rho_{b}(\mathbf{r})=1$ (5)

###  2.3 Dark Matter Halo Mixture Model

The mixture model for our problem is a weighted sum of the background
component ($\rho_{b}$) plus the $k$ halo components ($\rho$), hence $c=k+1$
components (see eq. 6). The parameters for each halo profile are the three
vectors representing the halo center $\mathbf{r_{0}}$, the size parameter
$r_{e}$, and the shape parameter $n$. These parameters are collected into a
five-dimensional parameter vector $\vec{\theta}=(\mathbf{r_{0}},r_{e},n)$.

The contribution of each component to the final density distribution is
uneven. The mixture proportions $\vec{w}=\\{w_{j}\\}$, $j=1,\dotsc,c$ are used
as weights to normalize each halo’s contribution to the mixture model. Based
on equation (1), the resulting model is

$\Sigma(\mathbf{r}|\vec{w},\Theta)=\frac{N}{M}\Big{(}w_{b}\cdot\rho_{b}(\mathbf{r})+\sum_{j=1}^{k}w_{j}\cdot\rho(\mathbf{r}-\mathbf{r_{0,j}}|r_{e,j},n_{j})\Big{)}$
(6)

where the weights $w_{j}$ give the mixture proportions $p_{j}=N\cdot w_{j}/M$
in equation (1) and $M$ is the total mass given by

$M=\int_{W}\Big{(}w_{b}+\sum_{j=1}^{k}w_{j}\cdot\rho(\mathbf{r}-\mathbf{r_{0,j}}|r_{e,j},n_{j})\Big{)}d\mathbf{r}$
(7)

The term $N/M$ works as a normalizing constant so that the integral of the
model $\Sigma$ is always the number of particles $N$. Note that we ignore
function $\rho_{b}(\mathbf{r})$ because it is constant 1. Function $\Sigma$ is
our probability density function for a mixture model problem as in equation 1.
Since the density of the universe is constrained to not be infinite, we can
model the density of the components relative to each other. In equation 6, the
mixture proportions are defined so that $\sum_{i=1}^{c}p_{i}=N$, and by
definition we can set $w_{1}=1$.

This normalization of the model by the total number of objects $N$ has other
advantages. The integral of $\Sigma(\mathbf{r})$ in a region $A\subset W$
gives the number of model objects in $A$. The integration of a chosen model
component over the entire volume $W$ gives the number of model objects
belonging to this component. The statistical model can thus be understood as
an inhomogeneous Poisson distribution with intensity $\Sigma(\mathbf{r})$.

Notice how the quantity $\rho_{e}$ in the Einasto profile will be masked into
the mixture proportion as a consequence of the mixture model because both are
factors to the profile.

###  2.4 Halo Identification and Characterization

Once the mixture model is fitted (section 3.1), analysis of each Einasto-
shaped component representing a halo can follow. Due to the additive nature of
a mixture model, it is easy to determine the dominant component at every
location. Here, we apply a simple decision rule to assign individual data
points to a halo: given a particle, we calculate its probability to belong to
each component. We then randomise the membership of the particle using the
probabilities.

Each model halo component can be compared to the empirical halo profile in the
context of the total mixture model of overlapping halos. Component $j$ has a
three-vector $\mathbf{r_{0,j}}$, giving its location, and three Einasto
parameters $w_{j}$, $r_{e,j}$ and $n_{j}$. We define the empirical halo
profile as the number of real particles per volume, the particle density, in
concentric shells centred at $\mathbf{r_{0,j}}$ with volumes
$S_{j}(r)=V_{j}(r+\text{d}r)-V_{j}(r)$. The empirical halo profile is

$\hat{\delta}_{j}(r)=\frac{1}{S_{j}(r)}(n_{j}(r+\text{d}r)-n_{j}(r))$ (8)

where $n_{j}(r)$ is the number of particles in a sphere of radius $r$ and
center $\mathbf{r_{0,j}}$. Note this has the same units, number of particles
per unit volume, as the mixture model in equation 6 by $N/M$. The estimated
density profile of the mixture model centered at $\mathbf{r_{0,j}}$ is

$\hat{P}(r|\mathbf{r_{0,j}})=\int_{S_{j}(r)}\Sigma(\mathbf{r}|\vec{p},\Theta)d\mathbf{r}.$
(9)

The integral of $\Sigma$ at any volume $V$ gives us the estimated number of
particles. Since profile $\hat{P}(r)$ is evaluated at the concentric shells
$V_{j}(r)$ with thickness $\text{d}r$, it should be understood as a density
profile over radius $r$. Component $j$ in the previous integral can be
isolated from the total model profile as

$\hat{\rho}(r|\mathbf{r_{0,j}})=\int_{S_{j}(r)}p_{j}\rho(\mathbf{r-r_{0,j}}|r_{e,j},n_{j})d\mathbf{r}.$
(10)

We can now compare the observed profile of $\hat{\delta}_{j}$ with the full
model estimator $\hat{P}(r|\mathbf{r_{0,j}})$ and with the isolated estimated
profile $\hat{\rho}(r|\mathbf{r_{0,j}})$, which may show important departures
from the total profile. For short distances, $\hat{\rho}$ and $\hat{P}$ should
have similar values. However, for distances at which our component $j$ starts
to increasingly overlap with other components, these functions will start to
diverge as other structures contribute more to $\hat{P}$.

Parametric models offer the possibility of easily generating new samples
following the model distribution. For a mixture model, this is done with an
inverse transform sampling for the estimated number of objects per component
$N_{j}$.

$N_{j}=\int_{W}p_{j}\rho(\mathbf{r-r_{0,j}}|r_{e,j},n_{j})d\mathbf{r}.$ (11)

For the possible interest of the user, software implementations of these
functions are included in the $R$ language and are included in our repository.

##  3 Model Fitting

###  3.1 Estimation of Model Parameters and Model Selection

The optimal mixture model for a data set is calculated by maximum likelihood
estimation (MLE) for the log-likelihood

$\log{L(\vec{p},\Theta|\mathbf{X})}=\sum_{i=1}^{N}\log{\Sigma(\mathbf{r}_{i}|\vec{p},\Theta)}-\int_{W}\Sigma(\mathbf{r^{\prime}}|\vec{p},\Theta)d\mathbf{r}^{\prime}$
(12)

where $\mathbf{X}=\\{\mathbf{r}_{i}\\}$ contains the point process distributed
in the window $W$ with surface density distribution $\Sigma$. Note that the
right-hand side term is the mass $M$ for the parameters $\Theta$. The MLE and
Bayesian best-fit model parameters $\Theta$ are calculated for a chosen number
of $c$ components. Model selection among models of different complexities is
based on minimising two commonly-used penalized likelihood measures, the
Bayesian Information Criterion (BIC) and the Akaike Information Criterion
(AIC) (Schwarz et al., 1978; Akaike, 1998),

$\displaystyle\text{BIC}(k)$ $\displaystyle=-2\log L+6(c-1)\log N$ (13)
$\displaystyle\text{AIC}(k)$ $\displaystyle=-2\log L+12(c-1)$ (14)

where $6(c-1)$ is the number of parameters in $c-1$ halo components plus the
background component and $\text{log}L$ is the log-likelihood for the best fit
parameters.

The relative strengths of AIC and BIC is widely debated, although both are
founded on powerful theorems (Lahiri, 2001; Konishi & Kitagawa, 2008; Burnham
& Anderson, 2002; Kass & Wasserman, 1995; Everitt et al., 2011). The BIC has a
well-accepted valuation for relative model merit: one model is strongly (very
strongly) favored over another when $\Delta(BIC)>6$ ($>10$) (Kass & Raftery,
1995). For the mixture model problem, Kuhn et al. (2014) found that the AIC
was more sensitive to the presence of sparse clusters in the presence of rich
clusters because its penalty for complexity is weaker when $N$ is large.

###  3.2 Goodness-of-Fit

A best-fit model with optimum complexity selected with a likelihood-based
criterion is not guaranteed to be a good fit to the data. A complex clustered
spatial distributions cannot be effectively fitted with a mixture model of a
few halos. It is therefore necessary to assess the overall quality of the fit
for the entire pattern by a study of the residuals to identify departures of
the model from the real data density distribution. Several such tests are
outlined here. Studies involving residual analysis of astronomical spatial
point processes for goodness-of-fit evaluation of maximum likelihood mixture
models include Kuhn et al. (2014) and de Souza et al. (2017).

The residual analysis that is used here is described in Baddeley et al. (2005,
2015) as a ‘raw residuals’. Raw residuals are defined as the absolute
difference between the real number of points in a region $A$ and our
estimation for the same region. For our mixture model,

$R(A)=n(\mathbf{X}\cap A)-\int_{A}\Sigma(u|\vec{p},\Theta)\text{d}u$ (15)

where $n(\mathbf{X}\cap A)$ is the number of data points in the region $A$.
For well-fitted models, the sum of the raw residuals should approach zero when
integrated over $W$, and the residual map should approach a random spatial
distribution with no correlations between the values of the residuals and the
locations of the data points (spatial white noise). Residuals should have low
amplitude when compared with the surface density function based on the data.

The calculation of the residuals is made using a quadrature or grid of dummy
points $\mathbf{Q}=\\{u_{i}\\}\subset W$, $i=1,\dotsc,T$ (Baddeley et al.,
2005). Each point defines a small cell where the residuals will be calculated.
These residuals create a sparse distribution that is 1 when the cell contains
a data point $\mathbf{r}$ and $-\Sigma(u|\mathbf{X},\vec{p},\Theta)$ at empty
locations in $W$, which is negative and typically close to zero.

To effectively visualise the spatial distribution of the residuals, they have
to be smoothed. With grid $\mathbf{Q}$ dense enough to approximate an
integral, we can obtain the smoothed residual map $s(u)$ with

$s(u)=\int_{W}\kappa_{\omega}(u-v)\text{d}R(v)=\sum_{i=1}^{N}\kappa_{\omega}(u-\mathbf{r}_{i})-\int_{W}\kappa_{\omega}(u-v)\cdot\sum(u|\mathbf{X},\vec{p},\Theta)\text{d}v$
(16)

where $\kappa_{\omega}$ is a kernel function and $R_{\Theta}(v)$ is the raw
residual in a cell $v$.

It is also useful to define the relative residuals $e(u)$, which are defined
as the residual map normalized by the model intensity. The model intensity
$\Sigma^{{\dagger}}_{\omega}$ is the right-hand term of eq. 16. Hence the
definition of these two functions is

$\Sigma^{{\dagger}}_{\omega}(u|\vec{p},\Theta)=\int_{W}\kappa_{\omega}(u-v)\Sigma(u|\vec{p},\Theta)\text{d}v$
(17)

$e(u)=s(u)/\Sigma^{{\dagger}}_{\omega}(u|\vec{p},\Theta)$ (18)

As we will see in the following sections, this function can be used to detect
data structures that have not been modeled by any model component. For any
structure that is properly mapped by the model, $s(u)$ will be a small
quantity and the values of $e(u)$ in its region will also be small. However,
if a structure in the data is not included in the model, then the value of
function $\Sigma^{{\dagger}}_{\omega}(u|\vec{p},\Theta)$ will be close to zero
for any location close to that structure and the error value in $s(u)$ will be
amplified in $e(u)$ and unfitted structures can be easily detected.

The kernel function that is used in this work is the Gaussian filter, which is
commonly used in cosmology (Martínez & Saar, 2002), where $\omega$ is the
smoothing radius or bandwidth. The appearance of the smoothing strongly
depends on the choice of this quantity. The bandwidth can be selected using
cross-validation or other techniques and we select it heuristically to give
informative residual maps.

Finally, we use the coefficient of determination (Rao et al., 1973) to assess
the global goodness of fit. Using the expected proportionality between the
model density ($\Sigma_{\omega}^{{\dagger}}$) and the data density
($\Sigma_{\omega}^{*}$) functions, we estimate a simple linear regression
between them and use the resulting $R^{2}$ coefficient as a measure of
goodness of fit (de Souza et al., 2017).

The data $X$ can be similarly smoothed with the same kernel. This function is
not used in the fitting or model validation process, which would be incorrect
given the loss of data structure detail. However, we include it in our code
repository for completeness and visualization purposes (as in Fig 6 top left-
hand panel):

$\Sigma^{*}_{\omega}(u)=\sum_{i=1}^{N}\kappa_{\omega}(u-\mathbf{r}_{i})\qquad\mathbf{r}_{i}\in
W\\\ $ (19)

###  3.3 Available Software Packages

It is not our aim to compare the performance of our code with other packages
available for users. In contrast, our solution is adapted to the particular
scenario of dark matter halos with Einasto profile, while most of the public
solutions use GMM or only accept two-dimensional data. However, we find it
interesting to provide a small summary of the available solutions.

Substantial code devoted for mixture model analysis are available in the R
public domain statistical software environment such as CRAN packages mixtools,
mclust, and EMCluster (Benaglia et al., 2009; Scrucca et al., 2016; Chen &
Maitra, 2015). The Python Machine Learning library scikit-learn (Pedregosa et
al., 2011) includes Gaussian mixture models (GMM) algorithms. Another option
is $EMMIX$, which is written in Fortran (McLachlan & Peel, 1999), where the
Expectation-Maximization (EM) algorithm is used to find the MLE (Krishnan &
McLachlan, 1997; McLachlan Geoffrey & Krishnan, 2008). However, estimation of
best fits in complex data sets can be difficult to achieve thanks to the
multiple peaks. A more robust procedure makes use of the stochastic EM
algorithm (Celeux & Govaert, 1992), where randomization of the steps seeks to
avoid trapping in the first found local maximum. Even with this technique, it
is advisable to repeatedly run estimation algorithms using different starting
values to avoid convergence to a non-optimal local maximum. If the same best
fit solution is achieved every time, then we can be more confident of reaching
an absolute maximum.

###  3.4 Software Repository

Our solution necessarily departs from the previous packages. As explained,
while most mixture model software assumes Gaussian shapes and use the EM
algorithm for parameter estimation, we will adopt different algorithms.
Mixture model solutions with non-Gaussian functions such as the Einasto
profile are not so easily achieved with this method, due to the reduced
curvature in the derivative.

Following Kuhn et al. (2014), who calculate two-dimensional mixture models
with the isothermal ellipsoid shapes, we will maximize the log-likelihood
function using the Nelder-Mead simplex algorithm (Nelder & Mead, 1965) as
implemented in function optim within CRAN package stats (Team, 2015). Testing
with different initial values is recommended to avoid trapping in local
maxima. Since the models have high dimensionality, a strategy of systematic
freezing and thawing parameters during estimation can be useful. Once the
Nelder-Mead algorithm has finished a first attempt at estimation, the results
can often be improved by freezing some of the parameters and repeating the
calculation for the remaining free parameters. This technique is valuable in
mixture models because the parameters of distant halos tend to be
uncorrelated. It can also be helpful to obtain a good estimates of halo
centers $\mathbf{r_{0}}$ before attempting the fitting of the rest of the
parameters. Confidence intervals on MLE parameters can be estimated from the
Fisher Information Matrix.

MLE can present additional problems that depend on the profile functions
$\rho_{i}(\mathbf{X},\Theta)$. It is possible for singularities to exist for
certain values of $\mathbf{X}$ or $\Theta$ where the likelihood becomes
infinite. This might happen when the number of parameters to be estimated is
high when compared with the sample size. This problem can be addressed with a
Bayesian approach where the posterior distribution of the parameters is mapped
instead of maximizing the log-likelihood function. As priors, we use Gaussian
functions with large variances for the centers of the halos, and log-normal
distributions for the size and shape parameters $r_{e}$ and $n$. These
calculations are performed using Markov chain Monte Carlo (MCMC). However,
even if aided by MLE results that were previously obtained by the Nelder-Mead
algorithm, sufficient MCMC evaluations of high-dimensional parametric models
can be computationally expensive.

The software that we have used for our Bayesian calculations is the CRAN
package LaplacesDemon (Statisticat & LLC., 2013a, b, c, d), which provides
more than 40 different MCMC algorithms. These algorithms make different
decisions regarding the next combination of parameters $\Theta$ to be tested
to efficiently map the normalized log-likelihood function 12. Adaptive MCMC
algorithms, which use the previous evaluations to choose the next $\Theta$,
are often more efficient in finding the overall distribution of the function
but we must always finish with a long run of a non adaptive algorithm to
ensure convergence. In this work, we make use of the Adaptive Metropolis-
Within-Gibbs (AMWG, Roberts & Rosenthal (2009)) and the twalk (Christen et
al., 2010) algorithms. These calculations start with a preliminary MLE
solution that is obtained with the Nelder-Mead algorithm.

Our software implementation of this mixture model decomposition of the dark
matter distribution is named DarkMix and is written in the R statistical
software environment, which has many tools relating to spatial point
processes, mixture models, likelihood calculations, model assessment, and
graphics. Our work extends that developed by Kuhn et al. (2014) to identify
star clusters in two dimensions. We present a new library to model three-
dimensional data sets (Hurtado-Gil et al., 2022), which is publicly available
on Github222https://github.com/LluisHGil/darkmix . Its use is explained in
sections 4 and 7. Additional documentation for the code can be found
online333https://darkmix.readthedocs.io/.

##  4 Darkmix Overview

With the equations from section 3 the darkmix code has the following
capabilities:

  1. 1.

Create data and model objects that are compatible with the R library spatstat.
This is a powerful library to work with point processes and its capabilities
can be combined with our code.

  2. 2.

Parametric estimation with two available methods: the Nelder-Mead simplex
algorithm (Nelder & Mead, 1965) and MCMC.

  3. 3.

Functions are incorporated for model selection (through AIC and BIC
calculations) and goodness-of-fit (through $R^{2}$ and residuals plots).

  4. 4.

Generation of model outputs, which comprise the soft classification of the
particles, the extraction of individual profiles for components, the
generation of model realizations and the visualization of the components’
sizes.

##  5 Model Validation

Before we introduce our real case data set in section 6 and estimate a mixture
model, we must first validate the performance of our algorithm. To do so, we
will estimate the mixture model in three different configurations: 3, 9 and 16
dark matter halos with background particles in a cube of 25 length units side.
For each configuration, we generate 20 realizations with three different
particle densities: $0.25$, $0.5$ and $0.75$ particles per volume unit.

The halos are modeled using the Einasto profile and the background components
is uniformly distributed, as assumed by our model. The true parameters
defining the location, size and shape of the halos, together with the number
of particles, are generated randomly (see Tables 1 and 2). We then obtained
the mixture coefficients using eq. 7 knowing that $w_{0}=1$. When using three
halos, the configuration is a simple set of components around the center of
the volume, and its fitting presented no problem. We will use the more complex
cases of 9 and 16 halos to study merging halos or locate halos near the volume
window boundary.

Table 1: True parameters for model validation: 9 halos

k | $x_{0}$ | $y_{0}$ | $z_{0}$ | $r_{e}$ | $n$ | $\log w$ | $N$
---|---|---|---|---|---|---|---
$1$ | $2.9$ | $21.0$ | $21.7$ | $1.1$ | $2.4$ | $0.0$ | $256$
$2$ | $8.2$ | $6.5$ | $18.6$ | $0.9$ | $1.4$ | $0.68$ | $544$
$3$ | $8.7$ | $14.9$ | $16.0$ | $1.4$ | $1.7$ | $-0.86$ | $66$
$4$ | $10.1$ | $16.2$ | $4.4$ | $1.3$ | $1.9$ | $-0.64$ | $92$
$5$ | $16.1$ | $7.7$ | $5.9$ | $2.0$ | $1.7$ | $-0.42$ | $518$
$6$ | $16.4$ | $22.8$ | $19.5$ | $2.2$ | $2.9$ | $-0.65$ | $454$
$7$ | $18.4$ | $16.3$ | $22.3$ | $1.4$ | $1.9$ | $-0.08$ | $403$
$8$ | $20.3$ | $6.1$ | $13.9$ | $0.7$ | $1.7$ | $1.09$ | $717$
$9$ | $21.7$ | $14.9$ | $7.9$ | $1.1$ | $2.3$ | $0.20$ | $415$

Note. — These parameters have been randomly generated in a cube of side 25 to
describe a configuration of 9 Einasto halos (here they are sorted by $x_{0}$)
with density 0.25 particles per volume unit. The background contains 442
particles, which gives $\log w_{b}=-2.4$. The populations $N$ for densities
0.5 and 0.75 can obtained multiplying the values in the table by 2 and 3.

Table 2: True parameters for model validation: 16 halos

k | $x_{0}$ | $y_{0}$ | $z_{0}$ | $r_{e}$ | $n$ | $\log w$ | $N$
---|---|---|---|---|---|---|---
$1$ | $1.3$ | $6.2$ | $5.2$ | $1.3$ | $1.9$ | $0.0$ | $204$
$2$ | $2.9$ | $21.0$ | $21.7$ | $1.1$ | $2.4$ | $-0.21$ | $83$
$3$ | $5.4$ | $15.7$ | $5.2$ | $1.1$ | $1.8$ | $0.76$ | $729$
$4$ | $5.7$ | $20.0$ | $6.9$ | $2.0$ | $2.1$ | $-0.57$ | $212$
$5$ | $7.9$ | $15.4$ | $12.8$ | $1.0$ | $1.4$ | $0.33$ | $179$
$6$ | $8.2$ | $6.5$ | $18.6$ | $0.9$ | $1.4$ | $0.25$ | $108$
$7$ | $8.7$ | $14.9$ | $16.0$ | $1.4$ | $1.7$ | $-0.56$ | $69$
$8$ | $9.2$ | $21.3$ | $19.4$ | $0.9$ | $1.6$ | $0.29$ | $130$
$9$ | $10.1$ | $16.2$ | $4.4$ | $1.3$ | $1.9$ | $-0.58$ | $55$
$10$ | $15.5$ | $10.6$ | $5.7$ | $1.0$ | $1.8$ | $0.71$ | $481$
$11$ | $16.1$ | $7.7$ | $5.9$ | $2.0$ | $1.7$ | $-0.4$ | $293$
$12$ | $16.4$ | $22.8$ | $19.5$ | $2.2$ | $2.9$ | $-0.84$ | $154$
$13$ | $18.4$ | $16.3$ | $22.3$ | $1.4$ | $1.9$ | $-0.5$ | $83$
$14$ | $20.3$ | $6.1$ | $13.9$ | $0.7$ | $1.7$ | $0.91$ | $255$
$15$ | $21.7$ | $14.9$ | $7.9$ | $1.1$ | $2.3$ | $0.22$ | $230$
$16$ | $23.8$ | $15.5$ | $16.8$ | $1.2$ | $2.3$ | $0.55$ | $571$

Note. — These parameters have been randomly generated in a cube of side 25 to
describe a configuration of 16 Einasto halos (here they are sorted by
$x_{0}$). Components 2, 6, 7, 9, 11, 12, 13, 15, and 16 are copied from the 9
halos configuration (see Table 1). The background contains 68 particles, which
gives $\log w_{b}=-2.85$. The populations $N$ for densities 0.5 and 0.75 can
obtained multiplying the values in the table by 2 and 3.

Once these realizations are generated, we estimate the best fit parameters
using our darkmix algorithm. As explained in section 4, this algorithm starts
with an approximation of initial values, which is later optimized using the
Nelder-Mead algorithm. The three-halos case is easy to model and the true
parameters lie inside the confidence intervals of our best fit. The samples
with nine-halos are more challenging, and we will present several results and
plots to show the performance of our model, which we consider acceptable. With
16 halos, the model is no longer able to find and correctly model some of the
extra halos because they are too small and faint to be distinguished from the
background.

Figure 1: Left: distribution of best fit parameters for the 20 realizations of
9 halos with low density (black) and true parameters from Table 1 (red
crosses). The validation of the center of the halos is showed as the distance
(first panel) to the true center (red line). Right: same for the 16 halos
configuration (see Table 2).

As we can see in Fig. 1, the parameters for a population of nine components is
generally correctly estimated. The center of the halos is found. This allows
for correct estimations of the radius $r_{e}$ and the Sérsic index $n$, only
in two cases do we see one of these parameters to be outside the 1-sigma error
bars. The mixture coefficients tend to be slightly overestimated, but a
correction of these values does not greatly improves the maximum likelihood.
Consequently, the estimation of the components population is rather accurate,
which allows for a correct classification of the particles. Even in the case
of 16 halos, when several components are incorrectly estimated, the halos’
population is close to the true values.

An additional validation analysis was performed with the nine-halo
configuration and low density: the same realization of particles (one of the
20 samples generated for the previous analysis) was estimated using mixture
models of different components, from 6 to 12 halos (see Fig. 2). The results
were shown to be satisfactory: for $k\leq 8$, the model finds the real halos
and tries to estimate the Einasto profile. However, it is not until we input
the right value of $k=9$ halos that the halos are not only found but correctly
estimated, as seen in Fig. 1. When $k\geq 10$, the model tries to estimate
spurious over-densities in the background population, which creates halos with
large radius and Sérsic index. With such values, the mass of the halo is
concentrated around a few close particles, which adds a minimal contribution
to the maximum likelihood of the model. This contribution is negligible, and
the AIC and BIC clearly show a regression for such models.

Figure 2: Relative change in the log-likelihood (gray band), and in the
information criteria BIC and AIC (symbols), as a function of the number of
components included in the model, $c$. The model is fit over the nine-halo
configuration. Changes are shown with respect to the best model in each case.
The inset shows a zoom of the range $10\leq c\leq 12$. The horizontal dotted
(dashed) line shows the threshold for strong (very strong) evidence,
$\Delta(BIC)>6$ ($>10$) according to Kass & Raftery (1995).

We conclude that our model is fully validated for configurations of at least
nine halos with background component and densities equal or greater than 0.25
particles per volume unit. In Fig. 3, we provide a four-panel plot with the
data and model densities, plus the raw and relative residuals. The plot shows
the close agreement between data and model. Only two halos (top left-hand a
bottom left-hand panels) seem to be overestimated (blue raw residuals) and
underestimated (red in raw residuals). The relative residuals show an
uncorrelated pattern with the data, which means that no component is excluded
from the model.

The case of 16 halos is still interesting and valid conclusions can be
obtained for the richer halos, although the fainted halos might add biased
results. All of the results and related data files can be found in our Git Hub
repository.

Figure 3: Comparison of one realization of a nine-halo sample and density 0.25
with the estimated model. Top left-hand panel: data kernel particle density
field $\Sigma_{\omega}^{*}(\mathbf{r}|\mathbf{X})$. Top right-hand panel:
smoothed raw residuals $s(\mathbf{r})$. Bottom left-hand panel: model particle
density field
$\Sigma_{\omega}^{{\dagger}}(\mathbf{r}|\mathbf{X},\vec{p},\Theta)$. Bottom
right-hand panel: relative residuals $e(\mathbf{r})$. In the data and model
density fields, gray intensity is in logarithmic scale. Both images have been
normalized over the number of particles. Dark areas indicate a denser region
or higher probability of being occupied by a point, while lighter areas are
unlikely to be occupied by a particle. Smoothing bandwidth for all maps is
$\omega=1$, and the values summed over the dimension Z. Data points from
Figure 4 are shown as small black dots.

##  6 Data

We perform this analysis for a case study data set from the Bolshoi simulation
(Klypin et al., 2011). This cosmological N-body simulation offers the
necessary conditions and a high particle resolution that allows us to apply
our methodologies over a sample of scientific interest. The MultiDark Database
(Riebe et al., 2013) hosts two 8.6 billion particle cosmological N-body
simulations: the already mentioned Bolshoi simulation (Klypin et al., 2011)
and MultiDark Run1 simulation (MDR1, or BigBolshoi) (Prada et al., 2012). The
Bolshoi simulation can be used to study both the large scale structure of the
universe and the properties of dark matter halos. In this work we focus on the
latter, which agreed with the assumptions made for our model.

The MultiDark Database allows us to use a SQL (Structured Query Language)
query interface to extract the desired sample. In this paper, data has been
extracted from the Bolshoi simulation. The simulation has been performed in a
volume of $250\,h^{-3}\text{Mpc}^{3}$, having a mass resolution of $1.35\times
10^{8}\,h^{-1}$ M⊙, and a force resolution of physical (proper) scale of
$1\,h^{-1}\text{kpc}$ (Klypin et al., 2011). The cosmological parameters that
we have used for this simulation are $\Omega_{m}=0.27$, $\Omega_{b}=0.0469$,
$\Omega_{\Lambda}=0.73$, $\sigma_{8}=0.82$, spectral index $n_{s}=0.95$, and
$H_{0}=100\,h$ km $\text{s}^{-1}\text{Mpc}^{-1}$ with $h=0.70$. The snapshot
that we have used is at redshift $z=0$.

Since this work is a case study, we select a small region of interest from the
table Bolshoi.Particles416 at https://www.cosmosim.org/. We are interested in
a volume $W$ containing an interesting structure with halos of different sizes
and merging cases. With this purpose, we select a flat cuboid with an squared
face to facilitate the two-dimensional examination. This sample contains three
halos that are among the 100 most massive, plus several other halos of smaller
size. The final sample contains 2081 particles in a volume of $4375\,h^{-3}$
Mpc3 (defined by a box of $25\,h^{-1}$ Mpc$\times 25\,h^{-1}$ Mpc$\times
7\,h^{-1}$ Mpc). We wish to remark here that we tested the method on other
more sparse galaxy samples. However, success is not reached when the data
sparsity of the sample is high. Evaluating which is the minimum needed
structure density to be properly described by our model is outside the scope
of this work, but this can serve as a reference. The Bolshoi simulation also
provides a catalog of halos that have been categorized with the BDM algorithm
(Klypin & Holtzman, 1997; Riebe et al., 2013). In section 7.4, we make a
comparison between these halos and our findings. It is worth noting here that
neither the BDM catalog of halos, nor any other catalog, has been used in this
work apart from at the appendix.

An image of the selected sample can be seen in Figure 4, notice the abundant
structure and variations in the shape and size of its clusters.

Figure 4: Sample of dark matter particles from the Bolshoi simulation.

##  7 DarkMix Application to the Bolshoi Simulation

###  7.1 The Fitting Procedure

This section serves as an example of how our code can be used to estimate a
model for the data set presented in section 6. All of the calculations and
results presented in this section have been obtained with darkmix code and a
full walk-through can be found in the code documentation
444https://darkmix.readthedocs.io/en/latest/darkmix_steps.html.

The data can be easily loaded into a spatsat point process object, while the
model and the parameters are defined as arrays of functions and values. The
code contains functions for Einasto profiles and the background component (a
constant function) but the user can easily define additional functions to
similarly model a structure of interest. Regarding the integration grid, this
object is used to calculate the model mass and in eq. 7. After testing
different grid sizes, the authors recommend using a relation of 2:1 (i.e., we
divide the sides of length 25 of our window into 50 parts). Thinner grids do
not produce a different result. We only recommend using a thinner grid for
plotting purposes (in this work, we use 128 grid points per side of 25 units).

Once these objects are created, we can proceed to estimate the parameters. It
is advisable to start with an initial guess of the number of halos and their
centers. We recommend using function centers, which outputs a list with the
$k$ densest locations in the data set and will be used as the initial guess
for the halo components.

The remaining parameters are defaulted to $1$ for the radius $r_{e}$ and $3$
for the Sérsic index $n$. The mixing coefficients are in logarithmic form: $1$
for the halos weights and $-2$ for the background component.

We will start with a $c=11$ components model and the Nelder-Mead simplex
algorithm via the R function optim, which maximizes the likelihood function of
our mixture model. However, the model needs several iterations to achieve our
best fit. To improve the estimation, additional functions have been created to
freeze some parameters while estimating the rest. While algorithms spend
little time finding the center of the halos, most of the computation is
devoted to the estimation of the three Einasto model parameters. The radius of
a component is generally independent of the center of a distant location and
reducing the total number of parameters per optimization can greatly improve
the fit. Each function fixed one of the different parameter types: centers,
radii, Sérsic index and mixture coefficients. It is advisable to repeat this
procedure until convergence and ending by calling the optim function with no
frozen parameters for a final fit. For the data set that we have used in this
work, the full procedure might take around 30 minutes to complete on a
commonly-available laptop.

The results for our data set in Figure 4 for an assumed $c=11$ number of
components are shown in Table 3. The expected number of particles per
component can be obtained directly from the mixture model: the component’s
mass is integrated independently and normalized to the total number of
particles. The last column of Table 3 shows the expected number of particles
per component.

Function optim and the other R routines that are designed to estimate the
maximum likelihood provide the hessian matrix of the best fit set of
parameters. This matrix can be used to obtain confidence intervals for the
parameters. However, the numerical approximation of the hessian matrix
obtained in our problem was ill-defined and produced negative variances. Given
the impossibility of calculating the confidence intervals, we decided to run a
MCMC routine with a simpler model and estimate them here. We explain this with
more detail in §7.5 .

The mixture model should be fitted for several values of $c$, and the AIC and
BIC (eq. 14 and 13) functions should be calculated for model selection. Figure
5 shows the log-likelihood and information criteria values from fits of the
Bolshoi simulation data set for models with $c=4$ to $15$ components. As
expected, the log-likelihood increases monotonically with the number of
components, while both the BIC and the AIC reach a minimum value for $c=11$.
Consequently, this is the preferred model according to both criteria. We can
interpret the evidence in favour of this model according to both criteria
using the scale introduced by Kass & Raftery (1995). In the case of BIC, the
$c=11$ model is very strongly favoured in all cases, as $\Delta(BIC)>10$ for
all other models. However, when considering the AIC, the models with $c=12,13$
cannot be discarded. This is a consequence of the AIC’s lower penalty for
increased model complexity (see eqs. 13, 14). Following the BIC result, and
choosing the most parsimonious among the accepted models by the AIC, we stick
to the model with $c=11$.

As we discuss in §7.5, it is possible for scientific considerations to prefer
a non-optimal model, such as $k=6$ rather than $k=10$. Generally, the most
populous halos will appear in all models and increasing $k$ will identify
small, sparser halos. However, it is also possible that small changes in $k$
may lead to major changes in the structure of the best-fit model. In
particular, spatial mixture models can identify a large diffuse halo that
encompasses or overlaps smaller, denser halos. This will occur when the
clustering of points has a complicated hierarchical structure, rather than
exhibiting distinct halo structures. Since the model for each value of $k$ is
a maximum likelihood fit, they are all statistically valid. Whether to use the
BIC and AIC for model selection or another choice of $k$ that gives a more
parsimonious or more complete model of the particle distribution is a
scientific decision.

Figure 5: Relative change in the log-likelihood (gray band), and in the
information criteria BIC and AIC (symbols), as a function of the number of
components included in the model, $c$. The model is fit over the Bolshoi
dataset. Changes are shown with respect to the best model in each case. The
inset shows a zoom of the range $c\geq 10$. The horizontal dotted (dashed)
line shows the threshold for strong (very strong) evidence, $\Delta(BIC)>6$
($>10$) according to Kass & Raftery (1995).  Table 3: Maximum Likelihood Fit
to the 10 Halo Model

| Initial estimates |  | Best fit parameters |
---|---|---|---|---
k | x | y | z | $\rho$ |  | $x_{0}$ | $y_{0}$ | $z_{0}$ | $r_{e}$ | $n$ | $\log w$ | N
1 | $33.5$ | $178.5$ | $99.5$ | $0.078$ |  | $33.5$ | $178.9$ | $99.7$ | $2.1$ | $14.1$ | $0.00$ | $688$
2 | $36.5$ | $192.5$ | $98.5$ | $0.060$ |  | $36.7$ | $192.4$ | $98.8$ | $0.9$ | $2.9$ | $0.95$ | $354$
3 | $39.5$ | $174.5$ | $97.5$ | $0.023$ |  | $39.5$ | $174.3$ | $97.5$ | $1.7$ | $28.5$ | $-0.37$ | $138$
4 | $25.5$ | $189.5$ | $98.5$ | $0.018$ |  | $26.0$ | $189.4$ | $98.9$ | $2.0$ | $26.2$ | $-0.14$ | $449$
5 | $20.5$ | $192.5$ | $100.5$ | $0.015$ |  | $20.4$ | $192.8$ | $100.9$ | $1.6$ | $28.3$ | $-0.53$ | $94$
6 | $38.5$ | $193.5$ | $96.5$ | $0.011$ |  | $38.7$ | $193.9$ | $96.3$ | $1.4$ | $4.0$ | $-0.06$ | $119$
7 | $37.5$ | $175.5$ | $99.5$ | $0.005$ |  | $37.4$ | $175.4$ | $99.5$ | $8.7$ | $9.2$ | $-2.75$ | $47$
8 | $33.5$ | $172.5$ | $100.5$ | $0.005$ |  | $33.1$ | $190.4$ | $99.4$ | $2.0$ | $26.9$ | $-1.43$ | $23$
9 | $32.5$ | $185.5$ | $99.5$ | $0.005$ |  | $32.8$ | $185.5$ | $100.0$ | $3.8$ | $30.0$ | $-1.92$ | $48$
10 | $21.5$ | $190.5$ | $98.5$ | $0.004$ |  | $25.1$ | $181.9$ | $95.3$ | $4.1$ | $25.3$ | $-1.91$ | $48$
Bk |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | $0.13$ | $73$

Note. — Input point process shown in Figure 4. $k$ identifies the halo
component. Initial $\mathbf{r_{0}}$ values from kernel density estimator in
order of decreasing maximum density $\rho$. Best fit parameters give the halo
center, Einasto parameters $r_{e}$ and $n$, mixing coefficient, and number of
dark matter particles.

###  7.2 Particle Membership in Halos

Once we have our estimated mixture model, we may want to classify the data set
particles into its different components. This comes naturally with a soft
classifier method such as a mixture model: given a particle, we evaluate each
model component at the particle location and we then normalize the obtained
quantities to one. The resulting values can be understood as membership
probabilities and used in a multinomial distribution to assign one particle to
$c$ components with the $c$ probabilities. If the model is correct, then the
number of particles assigned per component should match, on average, to that
shown in Table 3.

We recommend adapting this criteria to our scientific interests. Halos with
heavy tails might populate areas with particles far from the real halo
boundaries, hence overestimating the halo population. In contrast, under
merging circumstances, several halos might be competing for the same particle
and the model will assign low probabilities to each of them. Assuming same
size halos, the border between two merging halos will be populated by
particles with 0.5 probability of belonging to each halos, 0.33 for three
halos, etc. The multinomial criteria might end up assigning these particles to
the background instead, especially when the number of merging halos is high.
To compensate for this, we provide two modifications to the multinomial
criteria. First, the background component will not compete with the halos and
the particle will be directly assigned to the background whenever the
background probability is higher than that of any halo. Second, the user can
input a threshold value such that the particle is assigned directly to the
background component if all probabilities are below this value.

In our case study data set, we detect cases of the merging halos (see Figure
8). Following the rule given above, we recommend a threshold of 0.3, which is
sufficient to assign to the background all particles in an undecided
situation.

In Table 4, we provide an example of this classification procedure for the
first five particles in our data set, showing the probability of each particle
(rows) to belonging to each component (column). A final column gives the
identifying number of the halo component for each particle according to the
chosen decision rule. This table is used to plot Figure 8, with each component
given in a different color.

Table 4: Halo Membership Probabilities for Bolshoi Dark Matter Particles

Part | Location |  | Halo Component | Bkgd | Memb
---|---|---|---|---|---
| x | y | z |  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |  |
1 | 39.125 | 192.997 | 95.635 |  | 0.001 | 0.020 | 0.000 | 0.001 | 0.000 | 0.969 | 0.000 | 0.000 | 0.000 | 0.000 | 0.007 | 6
2 | 35.139 | 191.504 | 98.906 |  | 0 004 | 0.942 | 0.000 | 0.007 | 0.000 | 0.010 | 0.000 | 0.022 | 0.003 | 0.000 | 0.010 | 2
3 | 34.817 | 191.624 | 98.213 |  | 0.006 | 0.896 | 0.000 | 0.013 | 0.001 | 0.021 | 0.000 | 0.038 | 0.004 | 0.002 | 0.018 | 2
4 | 32.199 | 189.673 | 99.982 |  | 0.028 | 0.029 | 0.001 | 0.112 | 0.003 | 0.004 | 0.002 | 0.730 | 0.044 | 0.002 | 0.042 | 8
5 | 31.930 | 189.982 | 99.307d |  | 0.024 | 0.029 | 0.001 | 0.127 | 0.003 | 0.005 | 0.001 | 0.732 | 0.032 | 0.001 | 0.021 | 8

Note. — This only shows a portion of the membership table. The full table is
available in
https://github.com/LluisHGil/darkmix/blob/master/Output/membership.txt.

###  7.3 Goodness-of-Fit and Residual Analysis

Once our best fit parameters are estimated, they can be displayed and analyzed
for astrophysical properties. Some halo components show degenerate values for
$r_{e}$ and $n$. As can be seen in Figure 8, each component is numbered
according to Table 3 and the circles have radius $r_{e}$. Halos 7 to 10 have
extremely large radius. From Table 3, we also note that these components have
mixing coefficients $>20$ times lower than the other components—they have low
central densities with large radii and few particles. This is a sign that
these components do not follow an Einasto profile and cannot be correctly
modeled by our model. In addition, the Sérsic index $n$ is also
unreliable—while values around 3 are expected in astrophysics (Merritt et al.,
2006), much greater values are found instead. In §7.5, we provide confidence
intervals for the parameters and a covariance matrix can be used to evaluate
the reliability of the coefficients.

An Einasto profile with a high $n$ mimics a power law with a strong
concentration of points around the center and weak tails. Optimally, we would
use a different function to model this kind of structure, but we can choose
between including degenerated halos or neglecting them and incorporate their
particles to the background.

For our best fit model of 11 components, we have a coefficient of
determination $R^{2}=0.92$. In this calculation, the kernel density estimator
of the data (eq. 19) is compared to the model density field (eq. 17). Even if
these fields are expressed in different units, they must follow a linear
relation if the fitting is perfect. A linear regression between the densities
is fitted in each spatial tile and $R^{2}$ is reported as a measure of
goodness-of-fit.

These two fields, $\Sigma_{\omega}^{*}$ and $\Sigma_{\omega}^{{\dagger}}$, can
be seen in the left-hand panel of Fig. 6. We can see a clear agreement in both
the distribution, size and range of densities of the halos. The main
differences arise from asymmetries in the real structure, which are difficult
to capture with spherical halos.

Figure 6: Comparison of smoothed Bolshoi sample data and the 10-halo model.
Top left-hand panel: data kernel particle

density field $\Sigma_{\omega}^{*}(\mathbf{r}|\mathbf{X})$. Top right-hand
panel: smoothed raw residuals $s(\mathbf{r})$. Bottom left-hand panel: model
particle density field
$\Sigma_{\omega}^{{\dagger}}(\mathbf{r}|\mathbf{X},\vec{p},\Theta)$. Bottom
right-hand panel: relative residuals $e(\mathbf{r})$. In the data and model
density fields, gray intensity is in logarithmic scale. Both images have been
normalized over the number of particles. Dark areas indicate a denser region
or higher probability of being occupied by a point, while lighter areas are
unlikely of being occupied by a particle. The smoothing bandwidth for all maps
is $\omega=1$, and the values summed over the dimension Z. Data points from
Figure 4 are shown as small black dots.

Figure 7: Probability density distribution of the smoothed raw residuals.

We can see the probability density distribution of the smoothed raw residuals
in Fig. 7. The residuals are skewed ($g_{1}=-5.4$) at low values away from the
halo peaks and exhibit heavy tails at high values near halo peaks. The reason
for the heavy tails is evident in the lower left-hand panel, which shows the
absolute ($s(\mathbf{r})$) from equations 16. While errors have mean zero by
construction, a strong red-blue pattern around one of the biggest halos is
responsible for the largest residuals. This suggests that the shape is poorly
fitted by an Einasto profile. It is difficult to assess if these errors are
within the expected shot-noise errors of the model or if we should consider a
change in the model, such as a parameter fit refining or a different profile.
However, the bottom right-hand panel of Fig. 6, which maps the relative raw
residuals ($e(\mathbf{r})$), shows that these model errors are small compared
to the background, which is covered by red areas when data is present. This
implies that, in relative terms, the background has stronger fitting problems
that the halo components.

This panel is particularly effective at revealing structures that are not
included in the model. The most prominent structure missing in the model is at
the bottom where, due to truncation by the window edge, an Einasto profile of
a sparse cluster did not fit well. Finally, on Fig. 8 we present a randomized
classification of the particles based on our model. The plot shows how a soft
classifier mixes the particles of different components, specially under
merging conditions.

Figure 8: Particle classification based on Table 4. Clusters are labeled as in
Table 3 and the circle radii are the fitted Einasto radii $r_{e}$.

Another interesting visualization analysis is the profile extraction. With eq.
9 and 10 we can plot the one-dimensional profile of the whole mixture model
and also that of any individual halo component. In Fig. 9, we see the case for
components 2 and 4 in our classification. We can see that the Einasto profiles
(red curves) match the empirical profile (black dots) out to $\sim 2$ Mpc,
beyond which additional halos appear. The full model (green curves) follows
the empirical profile remarkably well, with peaks associated with other
clusters. With this method, mixture models can be used to disentangle the real
profile of merging halos, recovering the distribution of these structures
where the empirical profile does not allow it.

Figure 9: Radial density profiles of clusters 2 (left-hand panel) and 6
(right-hand panel). Black dots represent the particle number densities. The
red curve is the density profile of an isolated component
($\hat{\rho}(r|\mathbf{r_{0,j}})$) and the green curve is the density profile
of the full mixture model ($\hat{P}(r|\mathbf{r_{0,j}})$). Note the close
fitting of the full mixture model profile (in green).

###  7.4 Comparison with BDM Calculations

The tables provided by the MultiDark Database (Riebe et al., 2013), including
the Bolshoi simulation used in this work, include a catalog of halos that have
been detected with the Bound Density Maximum (BDM) algorithm (Klypin &
Holtzman, 1997; Riebe et al., 2013). This algorithm detects local density
maxima and defines a spherical halo that removes unbound particles. This
algorithm allows for the detection of subhalos, which are smaller structures
inside parent halos.

Our work makes use of a low density data sample, which does not permit the
detection of small structures, such as the less massive halos. Therefore, we
compare our detected halos with the most massive BDM structures. In table
Bolshoi.BDMV from the Multidark Database, and inside the region of our data
set, 26 halos can be found with more than 30,000 particles. Since we have 10
halo components in our best fit model, we select the 10 halos laying closer to
our centers. The resulting list can be found in the Github repository, file
bdm_halo.txt. In Figure 10 we compare both sets of halos. As can be seen, the
selected BDM halos lie very close to our halo centers and have been calculated
to contain a number of particles that highly correlates with our estimation
(see Table 3, column $N$). Since we are using just a sample of particles to
illustrate the methods, our numbers are only a fraction of the real number of
particles.

Figure 10: Left-hand panel: layout of dark matter particles with BDM halo
centers (squares) and our estimations (triangles). Right-hand panel: number of
particles found in the BDM catalog vs estimated number of particles in our
mixture model. Line indicates the best linear regression between both
variables.

###  7.5 A Simpler, Bayesian Model

In Fig. 5 we see how the AIC and BIC showed a plateau after $c=6$, Fig. 8
shows how the 7th component has a radius $r_{e}$ that extends around other
halos and does not follow an Einasto profile, and components $8-10$ are
sparse. It may thus be scientifically useful to examine a simpler model with 7
rather than 11 components.

This simplified model still has 36 parameters and was estimated using the MCMC
sampling approach with maximum a posterior best fit parameters, where
3,600,000 iterations are needed to satisfactorily map the distributions.
Although it is not included in darkmix.R and shown here, we recommend the user
to examine graphical and scalar (e.g., the Gelman-Rubin statistic) diagnostics
of convergence of the MCMC chains. CRAN package CODA is widely used for this
purpose. The long calculation time can be a considerable handicap and the MLE
procedure that was given earlier may be operationally more feasible for large
data sets. The advantage of a Bayesian approach is to map non-Gaussianities in
the posterior distributions and curved relationships in bivariate parameter
confidence intervals. However, for many science applications these
capabilities are not needed and the MLE approach would be preferred.

Table 5: Bayesian Fit to a Simpler c=7 Model

$k$ | $x_{0}$ | $y_{0}$ | $z_{0}$ | $r_{e}$ | $n$ | $\log w$ | $N$
---|---|---|---|---|---|---|---
1 | $33.46\pm 0.01$ | $178.80\pm 0.02$ | $99.70\pm 0.02$ | $1.05\pm 0.04$ | $2.4\pm 0.2$ | $1.5\pm 0.9$ | $599\pm 22$
2 | $36.67\pm 0.02$ | $192.33\pm 0.01$ | $98.75\pm 0.02$ | $0.77\pm 0.05$ | $3.0\pm 0.3$ | $2.0\pm 1.0$ | $413\pm 19$
3 | $26.00\pm 0.03$ | $189.35\pm 0.03$ | $98.87\pm 0.02$ | $1.11\pm 0.07$ | $2.5\pm 0.3$ | $0.8\pm 0.4$ | $376\pm 19$
4 | $38.60\pm 0.04$ | $193.88\pm 0.03$ | $96.18\pm 0.03$ | $0.75\pm 0.08$ | $2.2\pm 0.4$ | $0.7\pm 0.5$ | $108\pm 11$
5 | $39.47\pm 0.02$ | $174.22\pm 0.02$ | $97.44\pm 0.03$ | $0.90\pm 0.10$ | $3.2\pm 0.6$ | $0.4\pm 0.3$ | $127\pm 12$
6 | $20.40\pm 0.02$ | $192.72\pm 0.02$ | $100.81\pm 0.02$ | $0.56\pm 0.08$ | $3.2\pm 0.9$ | $1.0$ | $74\pm 8$
Bk |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | $0.14\pm 0.004$ | $384\pm 26$

Note. — Values are the mean and standard variation of the posterior
distribution. Component 6 has been used as the first component in the model,
and therefore $w_{6}=1$.

Figure 11: Residuals from the Bayesian fitting of the $k=6$ halo mixture
model. Absolute raw residuals $s(\mathbf{r})$ (left-hand panel) and relative
residuals $e(\mathbf{r})$ (right-hand panel) of the Bolshoi sample. Smoothing
bandwidth is $\omega=1$ and values summed over dimension Z.

In Table 5 we summarize the estimated maximum a posterior estimated parameters
with error bars.

Compared to the $c=11$ mixture model, the removal of components 7 to 10 has an
impact on the estimation of the other parameters. The $R^{2}$ coefficient is
0.936, which is similar to the more elaborate model. We show in Figure 11 the
absolute and the relative residuals ($s(\mathbf{r})$ and $e(\mathbf{r})$). As
can be seen, the structures that were previously modeled by components 7 and 8
appear now in red colors. The absence of model components leave these
structures underestimated. However, its intensity is similar to that of some
other regions in the data set and a visual inspection alone might not be
enough to detect them. We use the relative residuals in Figure 11 bottom panel
to detect their presence.

##  8 Conclusions

We have demonstrated the use of finite mixture models to detect and
characterize a sample of simulated dark matter halos. The maximum-likelihood
solution to a parametric model of the particle distribution produces a maximum
likelihood probability density function. Two results emerge from the model: a
list of halos with their properties (Table 3), and the probability that each
particle is a member of each cluster (Table 4). This is a ‘soft
classification’ measure and additional decision rules are needed to assign
membership probabilities to each particle. A variety of graphical and
statistical measures of goodness-of-fit are provided. This mixture model
approach has several important differences from other conventional clustering
techniques.

  1. 1.

The user specifies a radial profile function based on previous astronomical
experience or astrophysical insight. Here, we chose an Einasto profile, which
is the three-dimensional analog of the Sérsic profile and a generalization of
the de Vaucouleurs profile for elliptical galaxies. This permits an estimation
of parameters with a physical interpretation, such as the size, shape and
particle abundance of a halo. In addition, this estimation is done directly on
the three-dimensional data set, without the need to integrate the one-
dimensional profile and the loss of information that it involves.

  2. 2.

No threshold densities, maximum distance scales, minimum membership
population, or other arbitrary parameters are needed, as in nonparametric
clustering procedures such as the ‘friends-of-friends’ algorithm (single-
linkage hierarchical clustering) or DBSCAN. A unique MLE for each value of $k$
is calculated without any free parameters. Model complexity (i.e., the value
of $k$) is chosen based on penalized likelihood information criteria combined
with scientific judgment. The location, size and population of each halo is
not subject to arbitrary parameter choices but are outcomes of a likelihood
calculation based on an astrophysically reasonable profile function. Only the
membership of each particle can be enhanced using a threshold determined by
the number of merging halos. When two or three halos overlap, a value of 0.3
is advised.

  3. 3.

As a soft classification method, we have a measure of the reliability that a
given dark matter particle is a member of a given halo. Researchers who are
interested in merging scenarios might use these results in various ways to
obtain more informative descriptions of the interactions. For example, sparse
satellite halos can be identified and their merging history into large
protogalaxies can be traced. The ability of the mixture model to identify
clustered ensembles of halos (see Figure 8) may be useful. Mixture modeling of
time sequences of the dark matter distribution can reveal, in a statistically
rigorous fashion, when and where merging occurs. Alternatively, decision rules
can be constructed to identify the most isolated dark matter particles to
understand the evolution of particles that have not merged into equilibrated
halos.

  4. 4.

The parametric mixture model would be challenged in situations where halos
with a wide range of populations $N$ are present. If, for example, the data
set here was dominated by a large halo with $N\sim 10,000$ particles, then the
likelihood may not be sufficiently improved by the addition of a small halo
with $N\sim 100$. However, most nonparametric clustering techniques have
similar difficulties. Methods such as hierarchical density-based clustering
seek to address this difficulty (Campello et al., 2013).

Mixture models have been discouraged for use to understand the galaxy large-
scale structure distribution because it is dominated by interconnected curved
filamentary structures rather than centrally condensed clusters (Kuhn et al.,
2019). However, the dark matter collapses into distinct equilibrated halos
early in the evolution of cosmic structures and is more amenable to mixture
analysis. The biggest handicap of these models is the calculation costs for
large data sets and large $k$. This method is best applied to the
understanding of dark matter evolution within small boxes, rather than
analysis of a full multi-billion particle simulation. There is also the
possibility that convergence will be difficult, either for the optimization of
the MLE or for the MCMC chains of the Bayesian calculation. Other algorithms,
such as MultiNest (Skilling, 2004; Feroz & Skilling, 2013) or PolyChord
(Handley et al., 2015a, b), that are specially meant to estimate high-
dimensional multimodal likelihood functions might be tried. Finally, we
reiterate that all of the R functions used in this work can be downloaded from
our GitHub repository https://github.com/LluisHGil/darkmix.

## Acknowledgements

This work has been funded by the project
PID2019-109592GB-I00/AEI/10.13039/501100011033 from the Spanish Ministerio de
Ciencia e Innovación - Agencia Estatal de Investigación, by the Project of
excellence Prometeo/2020/085 from the Conselleria d’Innovació, Universitats,
Ciència i Societat Digital de la Generalitat Valenciana, and by the Acción
Especial UV-INV-AE19-1199364 from the Vicerrectorado de Investigación de la
Universitat de València.

The CosmoSim database used in this paper is a service by the Leibniz-Institute
for Astrophysics Potsdam (AIP). The MultiDark database was developed in
cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. The
Bolshoi and MultiDark simulations have been performed within the Bolshoi
project of the University of California High-Performance AstroComputing Center
(UC-HiPACC) and were run at the NASA Ames Research Center. The MultiDark-
Planck (MDPL) and the BigMD simulation suite have been performed in the
Supermuc supercomputer at LRZ using time granted by PRACE. E.D.F. thanks Penn
State’s Center for Astrostatistics for an environment where cross-disciplinary
research can be effectively pursued.

## References

  * Akaike (1998) Akaike, H. 1998, in Selected papers of hirotugu akaike (Springer), 199–213
  * Baddeley et al. (2015) Baddeley, A., Rubak, E., & Turner, R. 2015, Spatial point patterns: methodology and applications with R (CRC press)
  * Baddeley et al. (2005) Baddeley, A., Turner, R., Møller, J., & Hazelton, M. 2005, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 617, doi: 10.1111/j.1467-9868.2005.00519.x
  * Bardeen (1980) Bardeen, J. M. 1980, Phys. Rev. D, 22, 1882, doi: 10.1103/PhysRevD.22.1882
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Benaglia et al. (2009) Benaglia, T., Chauveau, D., Hunter, D. R., & Young, D. 2009, Journal of Statistical Software, 32, 1.  http://www.jstatsoft.org/v32/i06/
  * Benson et al. (2001) Benson, A. J., Frenk, C. S., Baugh, C. M., Cole, S., & Lacey, C. G. 2001, MNRAS, 327, 1041, doi: 10.1046/j.1365-8711.2001.04824.x
  * Bertschinger (1985) Bertschinger, E. 1985, ApJS, 58, 1, doi: 10.1086/191027
  * Burnham & Anderson (2002) Burnham, K. P., & Anderson, D. R. 2002, Model selection and multimodel inference: a practical information-theoretic approach (Springer Science & Business Media)
  * Campello et al. (2013) Campello, R. J., Moulavi, D., & Sander, J. 2013, in Pacific-Asia conference on knowledge discovery and data mining, Springer, 160–172
  * Celeux & Govaert (1992) Celeux, G., & Govaert, G. 1992, Computational statistics & Data analysis, 14, 315
  * Chen & Maitra (2015) Chen, W.-C., & Maitra, R. 2015, EMCluster: EM Algorithm for Model-Based Clustering of Finite Mixture Gaussian Distribution
  * Christen et al. (2010) Christen, J. A., Fox, C., et al. 2010, Bayesian Analysis, 5, 263
  * Church et al. (1997) Church, S. E., Ganga, K. M., Ade, P. A. R., et al. 1997, ApJ, 484, 523, doi: 10.1086/304362
  * Clerkin et al. (2017) Clerkin, L., Kirk, D., Manera, M., et al. 2017, MNRAS, 466, 1444, doi: 10.1093/mnras/stw2106
  * Colberg et al. (2000) Colberg, J. M., White, S. D. M., Yoshida, N., et al. 2000, MNRAS, 319, 209, doi: 10.1046/j.1365-8711.2000.03832.x
  * Coles & Jones (1991) Coles, P., & Jones, B. 1991, MNRAS, 248, 1, doi: 10.1093/mnras/248.1.1
  * Cooray & Sheth (2002) Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1, doi: 10.1016/S0370-1573(02)00276-4
  * de Souza et al. (2017) de Souza, R. S., Dantas, M. L. L., Costa-Duarte, M. V., et al. 2017, MNRAS, 472, 2808, doi: 10.1093/mnras/stx2156
  * Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Einasto (1968) —. 1968, Publications of the Tartu Astrofizica Observatory, 36, 414
  * Einasto (1969) —. 1969, Astrofizika, 5, 137
  * Elahi et al. (2019) Elahi, P. J., Cañas, R., Poulton, R. J. J., et al. 2019, PASA, 36, e021, doi: 10.1017/pasa.2019.12
  * Everitt et al. (2011) Everitt, B., Landau, S., Leese, M., & Stahl, D. 2011, Cluster Analysis (¿ Wiley Series in Probability and Statistics), Chichester, UK: Wiley
  * Everitt (2005) Everitt, B. S. 2005, Encyclopedia of statistics in behavioral science
  * Feroz & Skilling (2013) Feroz, F., & Skilling, J. 2013, in American Institute of Physics Conference Series, Vol. 1553, American Institute of Physics Conference Series, ed. U. von Toussaint, 106–113, doi: 10.1063/1.4819989
  * Fillmore & Goldreich (1984) Fillmore, J. A., & Goldreich, P. 1984, ApJ, 281, 9, doi: 10.1086/162071
  * Fraley & Raftery (2002) Fraley, C., & Raftery, A. E. 2002, Journal of the American statistical Association, 97, 611
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter, S. 2006, Finite mixture and Markov switching models (Springer Science & Business Media)
  * Fruhwirth-Schnatter et al. (2019) Fruhwirth-Schnatter, S., Celeux, G., & Robert, C. P. 2019, Handbook of mixture analysis (Chapman and Hall/CRC)
  * Gower & Ross (1969) Gower, J. C., & Ross, G. J. 1969, Journal of the Royal Statistical Society: Series C (Applied Statistics), 18, 54
  * Gunn & Gott (1972) Gunn, J. E., & Gott, J. Richard, I. 1972, ApJ, 176, 1, doi: 10.1086/151605
  * Handley et al. (2015a) Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015a, MNRAS, 453, 4384, doi: 10.1093/mnras/stv1911
  * Handley et al. (2015b) —. 2015b, MNRAS, 450, L61, doi: 10.1093/mnrasl/slv047
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hoffman & Shaham (1985) Hoffman, Y., & Shaham, J. 1985, ApJ, 297, 16, doi: 10.1086/163498
  * Hurtado-Gil et al. (2022) Hurtado-Gil, L., Kuhn, M. A., & Arnalte-Mur, P. 2022, LluisHGil/darkmix: v1.0.3, v1.0.3, Zenodo, doi: 10.5281/zenodo.6855703
  * Hurtado-Gil et al. (2017) Hurtado-Gil, L., Martínez, V. J., Arnalte-Mur, P., et al. 2017, A&A, 601, A40, doi: 10.1051/0004-6361/201629097
  * Hurtado-Gil et al. (2016) Hurtado-Gil, L., Arnalte-Mur, P., Martínez, V. J., et al. 2016, ApJ, 818, 174, doi: 10.3847/0004-637X/818/2/174
  * Jenkins et al. (2001) Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372, doi: 10.1046/j.1365-8711.2001.04029.x
  * Kaiser (1984) Kaiser, N. 1984, ApJ, 284, L9, doi: 10.1086/184341
  * Kass & Raftery (1995) Kass, R. E., & Raftery, A. E. 1995, Journal of the american statistical association, 90, 773
  * Kass & Wasserman (1995) Kass, R. E., & Wasserman, L. 1995, Journal of the american statistical association, 90, 928
  * Kauffmann et al. (1999) Kauffmann, G., Colberg, J. M., Diaferio, A., & White, S. D. M. 1999, MNRAS, 303, 188, doi: 10.1046/j.1365-8711.1999.02202.x
  * Klypin & Holtzman (1997) Klypin, A., & Holtzman, J. 1997, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/9712217
  * Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102, doi: 10.1088/0004-637X/740/2/102
  * Konishi & Kitagawa (2008) Konishi, S., & Kitagawa, G. 2008, Information criteria and statistical modeling (Springer Science & Business Media)
  * Krishnan & McLachlan (1997) Krishnan, T., & McLachlan, G. 1997, Wiley, 1, 58
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32, doi: 10.3847/1538-4357/aaef8c
  * Kuhn et al. (2014) Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2014, ApJ, 787, 107, doi: 10.1088/0004-637X/787/2/107
  * Lahiri (2001) Lahiri, P. 2001, in Model selection, IMS
  * Madgwick et al. (2003) Madgwick, D. S., Hawkins, E., Lahav, O., et al. 2003, MNRAS, 344, 847, doi: 10.1046/j.1365-8711.2003.06861.x
  * Martínez & Saar (2002) Martínez, V. J., & Saar, E. 2002, Statistics of the Galaxy Distribution (Chapman & Hall/CRC Press)
  * McLachlan & Krishnan (2007) McLachlan, G., & Krishnan, T. 2007, The EM algorithm and extensions, Vol. 382 (John Wiley & Sons)
  * McLachlan & Peel (1999) McLachlan, G., & Peel, D. 1999, Journal of Statistical Software, 004.  https://EconPapers.repec.org/RePEc:jss:jstsof:v:004:i02
  * McLachlan & Peel (2000) —. 2000, Finite Mixture Models, 238
  * McLachlan Geoffrey & Krishnan (2008) McLachlan Geoffrey, J., & Krishnan, T. 2008, The EM algorithm and extensions, Wiley Series in Probability and Statistics. Wiley-Interscience
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685, doi: 10.1086/508988
  * Moore et al. (1999) Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999, MNRAS, 310, 1147, doi: 10.1046/j.1365-8711.1999.03039.x
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * Navarro et al. (1997) —. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Nelder & Mead (1965) Nelder, J. A., & Mead, R. 1965, The computer journal, 7, 308
  * Neyman & Scott (1952) Neyman, J., & Scott, E. L. 1952, ApJ, 116, 144, doi: 10.1086/145599
  * Neyman & Scott (1953) —. 1953, Proceedings of the National Academy of Science, 39, 737, doi: 10.1073/pnas.39.8.737
  * Neyman & Scott (1954) —. 1954, Proceedings of the National Academy of Science, 40, 873, doi: 10.1073/pnas.40.10.873
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Peel & MacLahlan (2000) Peel, D., & MacLahlan, G. 2000, John & Sons
  * Prada et al. (2012) Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018, doi: 10.1111/j.1365-2966.2012.21007.x
  * Rao et al. (1973) Rao, C. R., Rao, C. R., Statistiker, M., Rao, C. R., & Rao, C. R. 1973, Linear statistical inference and its applications, Vol. 2 (Wiley New York)
  * Riebe et al. (2013) Riebe, K., Partl, A. M., Enke, H., et al. 2013, Astronomische Nachrichten, 334, 691, doi: 10.1002/asna.201211900
  * Roberts & Rosenthal (2009) Roberts, G. O., & Rosenthal, J. S. 2009, Journal of Computational and Graphical Statistics, 18, 349, doi: 10.1198/jcgs.2009.06134
  * Schwarz et al. (1978) Schwarz, G., et al. 1978, The annals of statistics, 6, 461
  * Scrucca et al. (2016) Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. 2016, The R Journal, 8, 205.  https://journal.r-project.org/archive/2016-1/scrucca-fop-murphy-etal.pdf
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Sersic (1968) Sersic, J. L. 1968, Atlas de Galaxias Australes, Vol. 1 (Observatorio Astronomico, Universidad Nacional de Cordoba)
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405, doi: 10.1063/1.1835238
  * Somerville & Primack (1999) Somerville, R. S., & Primack, J. R. 1999, MNRAS, 310, 1087, doi: 10.1046/j.1365-8711.1999.03032.x
  * Statisticat & LLC. (2013a) Statisticat, & LLC. 2013a, LaplacesDemon: Complete Environment for Bayesian Inference.  http://cran.r-project.org/web/packages/LaplacesDemon/index.html
  * Statisticat & LLC. (2013b) —. 2013b, Bayesian Inference.  http://cran.r-project.org/web/packages/LaplacesDemon/index.html
  * Statisticat & LLC. (2013c) —. 2013c, LaplacesDemon Examples.  http://cran.r-project.org/web/packages/LaplacesDemon/index.html
  * Statisticat & LLC. (2013d) —. 2013d, LaplacesDemon Tutorial.  http://cran.r-project.org/web/packages/LaplacesDemon/index.html
  * Team (2015) Team, R. C. 2015, URL http://www. R-project. org
  * Willis et al. (2020) Willis, J. S., Schaller, M., Gonnet, P., & Helly, J. C. 2020, arXiv e-prints, arXiv:2003.11468.  https://arxiv.org/abs/2003.11468
  * Zehavi et al. (2011) Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59, doi: 10.1088/0004-637X/736/1/59


$\sigma_{i}=(1-t)\sigma_{i}^{\text{start}}+t(\sigma_{i}^{\text{end}})$ for
$i=1,2$. For filaments,
$\mu_{1}^{\text{start}}=0.3,\mu_{1}^{\text{end}}=0.38,\mu_{2}^{\text{start}}=0.7,\mu_{2}^{\text{end}}=0.62,\sigma_{1}^{\text{start}}=\sigma_{2}^{\text{start}}=0.05,\sigma_{1}^{\text{end}}=\sigma_{2}^{\text{end}}=0.1$.
For foam cells,
$\mu_{1}^{\text{start}}=0.7,\mu_{1}^{\text{end}}=0.55,\mu_{2}^{\text{start}}=0.3,\mu_{2}^{\text{end}}=0.45,\sigma_{1}^{\text{start}}=\sigma_{2}^{\text{start}}=0.05,\sigma_{1}^{\text{end}}=\sigma_{2}^{\text{end}}=0.1$.
These parameters are chosen such that voxels labeled 1 in the binary image
will, in expectation, be brighter than voxels labeled 0 in the binary image.
Moreover, the parameters are chosen such that, for $t=0$, all algorithms we
consider output near perfect reconstructions of the binary images, and, for
$t=1$, all algorithms we consider output binary images with severe mistakes.

###  A.2 Detailed numerical results for renderings

Table 1 and Table 2 report the numerical results corresponding to the examples
that are depicted in Figure 8 and Figure 9, respectively.

$t$ | Alg. | $\theta_{\text{start}}$ | $\theta_{\text{end}}$ | $b$ | $\operatorname{VI-WS}$ | $\operatorname{FC}$ | $\operatorname{FJ}$ | $\operatorname{VI-NS}$ | $\operatorname{FC-NS}$ | $\operatorname{FJ-NS}$
---|---|---|---|---|---|---|---|---|---|---
0.00 | WS | 0.580 | 0.580 |  | 0.056 | 0.008 | 0.048 | 0.000 | 0.000 | 0.000
| MS |  |  | 0.06 | 0.054 | 0.007 | 0.046 | 0.000 | 0.000 | 0.000
0.25 | WS | 0.580 | 0.585 |  | 0.104 | 0.022 | 0.082 | 0.000 | 0.000 | 0.000
| MS |  |  | 0.09 | 0.079 | 0.013 | 0.065 | 0.000 | 0.000 | 0.000
0.50 | WS | 0.580 | 0.580 |  | 0.250 | 0.075 | 0.174 | 0.000 | 0.000 | 0.000
| MS |  |  | 0.08 | 0.144 | 0.059 | 0.084 | 0.000 | 0.000 | 0.000
0.75 | WS | 0.570 | 0.570 |  | 0.750 | 0.421 | 0.330 | 0.435 | 0.435 | 0.000
| MS |  |  | 0.01 | 0.406 | 0.317 | 0.089 | 0.000 | 0.000 | 0.000
1.00 | WS | 0.555 | 0.565 |  | 1.447 | 0.781 | 0.666 | 1.045 | 0.907 | 0.138
| MS |  |  | -0.11 | 1.345 | 1.274 | 0.071 | 0.749 | 0.749 | 0.000
Table 1: This table summarizes the numerical results corresponding to the renderings in Figure 8 for filaments. For the watershed algorithm (WS) the thresholds $\theta_{\text{start}}$ and $\theta_{\text{end}}$, and for the multi-separator (MS) algorithm the biases are chosen so as to minimize the average $\operatorname{VI-WS}$ across those images of the data set with the amount of noise $t$. $t$ | Alg. | $\theta_{\text{start}}$ | $\theta_{\text{end}}$ | $b$ | $\operatorname{VI-WS}$ | $\operatorname{FC}$ | $\operatorname{FJ}$ | $\operatorname{VI-NS}$ | $\operatorname{FC-NS}$ | $\operatorname{FJ-NS}$
---|---|---|---|---|---|---|---|---|---|---
0.00 | WS | 0.380 | 0.470 |  | 0.566 | 0.421 | 0.144 | 0.023 | 0.023 | 0.000
| MS |  |  | -0.10 | 0.480 | 0.354 | 0.126 | 0.000 | 0.000 | 0.000
0.25 | WS | 0.390 | 0.470 |  | 1.083 | 0.787 | 0.296 | 0.209 | 0.101 | 0.108
| MS |  |  | -0.08 | 0.769 | 0.535 | 0.234 | 0.005 | 0.005 | 0.000
0.50 | WS | 0.120 | 0.480 |  | 2.110 | 1.281 | 0.828 | 0.671 | 0.283 | 0.388
| MS |  |  | -0.03 | 1.400 | 0.970 | 0.430 | 0.023 | 0.022 | 0.001
0.75 | WS | 0.120 | 0.470 |  | 3.273 | 2.540 | 0.734 | 0.657 | 0.472 | 0.185
| MS |  |  | 0.06 | 2.552 | 1.565 | 0.987 | 0.085 | 0.080 | 0.005
1.00 | WS | 0.430 | 0.600 |  | 4.526 | 2.810 | 1.716 | 3.275 | 1.941 | 1.334
| MS |  |  | 0.17 | 4.769 | 2.504 | 2.265 | 1.586 | 0.674 | 0.912
Table 2: This table summarizes the numerical results corresponding to the
renderings in Figure 9 for foam cells. For the watershed algorithm (WS) the
thresholds $\theta_{\text{start}}$ and $\theta_{\text{end}}$, and for the
multi-separator (MS) algorithm the biases are chosen so as to minimize the
average $\operatorname{VI-WS}$ across those images of the data set with the
amount of noise $t$.

###  A.3 Analysis of parameters of watershed algorithm

Figures 15 and 16 depict the inaccuracy of multi-separators computed by the
watershed algorithm for volume images of filaments and foam cells,
respectively.

Figure 15: Depicted above is the inaccuracy (distance from truth) of multi-
separators computed by the watershed algorithm (Section 5.3) for volume images
of filaments with five amounts of noise $t$ (columns) for all values of
$\theta_{\text{start}}$ (x-axis) and $\theta_{\text{end}}$ (y-axis). The three
rows depict the medians across those images of the data set with the amount of
noise $t$ of the $\operatorname{VI-WS}$ (top) and the conditional entropies
due to false cuts ($\operatorname{FC}$, middle) and false joins
($\operatorname{FJ}$, bottom). White indicates a value of $0$ and saturated
blue indicates a value of $2$ or greater. Infeasible parameters (i.e.
$\theta_{\text{start}}>\theta_{\text{end}}$) are depicted in gray. The +
symbol indicates the parameters that archive the best median
$\operatorname{VI-WS}$.  Figure 16: Depicted above is the inaccuracy (distance
from truth) of multi-separators computed by the watershed algorithm (Section
5.3) for volume images of foam cells completely analogous to Figure 15.