#  Unravelling the role of cosmic velocity field in dark matter halo mass
function using deep learning

Saba Etezad-Razavi1 Erfan Abbasgholinejad1 Mohammad-Hadi Sotoudeh1,2,3 Farbod
Hassani4 Sadegh Raeisi1 Shant Baghram1<EMAIL_ADDRESS>1Department of
Physics, Sharif University of Technology, Tehran 11155-9161, Iran
2Department of Physics, Univesité de Montréal, Montréal, Canada
3Mila - Quebec Artificial Intelligence Institute, Montréal, Canada
4Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway

###### Abstract

We discuss an implementation of a deep learning framework to gain insight into
dark matter (DM) structure formation. We investigate the impact of initial
velocity and density field information on the construction of halo mass
function (HMF) through cosmological $N$-body simulations. We train a
Convolutional Neural Network (CNN) on the initial snapshot of a DM-only
simulation to predict the HMF that individual particles fall into at $z=0$, in
the halo mass range of $10.5<\log(M/M_{\odot})<14$. Our results show a
negligible improvement from including the velocity in addition to the density
information when considering simulations based on ($\Lambda$CDM) with the
amplitude of initial scalar perturbations $A_{s}=2\times 10^{-9}$. To
investigate the effect of the initial velocity field in constructing the halo
mass function, we increase the initial power spectrum such that we see the
effect of velocities in larger halos visible to the resolution of our
simulations. The CNN model trained on the simulation snapshots with large
$A_{s}$ shows a considerable improvement in the HMF prediction when adding the
velocity field information. Eventually, for the simulation with $A_{s}=8\times
10^{-8}$, the model trained with only density information shows at least
$80\%$ increase in the mean squared error relative to the model with both
velocity and density information, which indicates the failure of the density-
only model to predict the halo mass function in this case. Our work shows the
interpretability and ability of CNNs to read higher-order information from
simple images, making them an excellent tool for cosmological studies.

##  I Introduction

The standard model of cosmology known as $\Lambda$CDM (cosmological constant
plus cold dark matter) is established in the last two decades through
observation of the cosmic microwave background radiation (CMB) [1] and large-
scale structure (LSS) surveys [2, 3, 4]. Despite many successes of this model
in describing the evolution of the Universe, and the statistical distribution
of structures, the nature of dark energy, dark matter (DM) and the physics of
the early Universe is still unknown. The future LSS surveys, for example, the
Vera Rubin observatory, Euclid satellite, and Nancy Groce Roman mission are
planned to map the Universe in more detail over a wider range of redshifts
[5]. Accordingly, the study of LSS in the Universe employing high precision
$N$-body simulations becomes crucial, especially to compare different models
against precise observational data. In the standard picture of structure
formation, the galaxies, group of galaxies and galaxy clusters reside in DM
halos [6]. These DM halos are the nodes of the cosmic web, which are
structured by hierarchical formation. In this sense, the study of DM halo
formation is essential for understanding galaxy formation and evolution and
testing our cosmological models. The evolution of DM halos from an initial
condition to the late time structures is a complicated non-linear process
based on gravitational instability, mainly through mass accretion and mergers
of DM structures [7]. Also, the knowledge of complicated baryon physics
becomes essential to populate DM halos. In this direction cosmological
$N$-body simulations are used to model these non-linear processes [8, 9, 10].
On the other hand, the analytical models can give insights into the process of
structure formation. One of the old but enlightened ideas of these models of
structure formation is related to the nominal work of Press and Schechter (PS)
[11]. They propose that by looking at the statistics of the density field in
higher redshifts, where the perturbations are almost Gaussian and linear, we
can predict the resulting halo mass function by linear evolution of the
density contrast field to the late time. This raises the question that what
are the main features of the early Universe which are important to predict the
DM halo mass function. The two ideas of peak theory [12] and excursion set
theory (EST) [13] are theoretical developments to use the idea of using
initial conditions to predict the late time distribution of matter. The two
main features of these analytical models are the collapse model (spherical
collapse, ellipsoidal collapse [14, 15] and the smoothing functions of the
initial density field [16]. Regardless of the massive success of analytical
frameworks and simulation in shedding light on our understanding of structure
formation, they still are not complete. Analytical models fail to accurately
predict the formation of structures for the smallest halo ranges, and the
simulations lack the physical insight that we could get from analytical models
and leaves unanswered questions like, Is the initial information enough to
predict halo mass function (HMF) accurately? What are the main features of the
initial cosmic fields affecting the halo formation?

It is worth mentioning, that due to the enormous amount of data going to be
delivered in near future from the cosmological surveys and simulations, we
have already entered an era of big data in cosmology. Thus, the Machine
Learning (ML) techniques become a complementary way to study the LSS of the
Universe [17]. ML can be used to have a more precise prediction for the LSS
observables and to shed light on the physics and process of structure
formation [18, 19].
In this work, we use ML methods to investigate the proposition of the
analytical models of structure formation in predicting the late-time
distribution of matter. Mainly, we aim to investigate the ideas of spherical
collapse and ellipsoidal collapse. The first one is related to the initial
smoothed density and the second one incorporates the physics of velocity field
and tidal shear of the collapse region. We use Convolutional Neural Networks
(CNN) to study the effect of density contrast and velocity field in LSS
formation. The organisation of the paper is as follows: In section II we
review the theoretical background of this work, and in section III we discuss
the $N$-body simulation setups and the data being used in our work, namely we
describe our simulations. Section IV is devoted to the technical description
of our method and the deep learning framework. We present our results and
interpretations in section V, and finally, we have our conclusion and future
remarks in section VI.

##  II Theoretical Background

In the context of the standard model, DM halos are the host of the luminous
matter [11, 13]. Accordingly, the statistics and distribution of galaxies are
promising tracers of DM distribution in the Universe [20]. There are two
different approaches to studying the formation and evolution of DM halos,
first through $N$-body simulations and second, utilising theoretical ideas to
study the mass profile of DM halos. The analytical models are mainly based on
the idea of PS [11]. In this framework, we study the Universe in high redshift
where the perturbations are almost Gaussian and linear. Then, we predict the
number density of DM halos by integrating the probability distribution (PDF)
of density contrast from the critical density of collapse in spherical model
$\delta_{\text{sc}}\simeq 1.69$ to infinity [11]. The PDF is constructed using
a specific smoothing scale. The mass enclosed in the smoothing region is equal
to the mass of the DM halo formed later during the evolution of the overdense
region. Then the two ideas of peak theory [12] and excursion set theory [13]
were developed accordingly. An important ingredient of these frameworks is the
collapse model. The spherical collapse model [21] only depends on the value of
the over-density in the smoothed scale. However, a more realistic model of
ellipsoidal collapse [15] indicates the importance of a strong shearing field
of velocity which makes the collapse of smaller mass halos more difficult than
the larger halo mass. This effect manifests itself in a mass-dependent
collapse barrier. To quantify the process we show how the non-linear
distribution of DM structures is related to the initial condition. The initial
power spectrum of curvature perturbation is defined as

${\cal{P}}_{{\cal{R}}}=A_{s}(\frac{k}{k_{p}})^{n_{s}-1},$ (1)

where $A_{s}$ is the amplitude of the curvature perturbation, $n_{s}$ is the
spectral index and $k_{p}$ is the pivot wavenumber. The amplitude and the
spectral index are fixed by the CMB observations to the values $A_{s}\simeq
2\times 10^{-9}$ and $n_{s}\simeq 0.96$ [1]. We can find the matter power
spectrum via the evolution of the scalar Bardeen potential through cosmic time
for different modes as [22]

$P_{m}(k,z)=A_{l}k^{n_{s}}T^{2}(k)D^{2}(z),$ (2)

where $A_{l}$ is the late time matter power spectrum amplitude, $T(k)$ and
$D(z)$ are transfer function and growth function, respectively. The late time
amplitude of matter power spectrum is related to the initial curvature
amplitude as
$A_{l}=A_{s}\times(8\pi^{2}k_{p}^{1-n_{s}})/({25\Omega_{m}^{2}H_{0}^{4}})$. We
can obtain the scale of non-linearity which is a function of redshift and mass
using the definition of variance $S$

$S(M,z)\equiv\sigma^{2}(M,z)=\int\frac{dk}{2\pi^{2}}k^{2}P_{m}(k,z)\tilde{W}^{2}(kR),$
(3)

where $\tilde{W}(kR)$ is the Fourier transfer of the top-hat window function
with smoothing scale $R$ related to the mass scale $M$. For the standard cold
dark matter model, the variance is a monotonic decreasing function of mass.
Larger mass corresponds to lower variance scales. This means that by
increasing the variance, we will probe the smaller mass scales. In this work,
we increase the primordial curvature perturbation, which results in a higher
variance mimicking the small mass range. We will discuss this idea in the
simulation section as well. The number density of dark matter halos $n(M,z)$
can be expressed in terms of variance and the collapse threshold as [15]

$n(M,z)=\frac{\bar{\rho}}{M^{2}}f(\nu)\frac{d\ln\nu}{d\ln M},$ (4)

where $\bar{\rho}$ is the mean density of matter in the present time,
$\nu=(\delta_{\text{sc}}/\sigma(M))^{2}$ and $f(\nu)$ is a universal function,
related to the first up-crossing statistics in the EST [23]. For the
ellipsoidal collapse, the universality function is given

$\nu
f(\nu)=\bar{A}[1+(\bar{a}\nu)^{-p}](\frac{{\bar{a}}\nu}{2})^{1/2}\frac{e^{-\bar{a}\nu/2}}{\sqrt{\pi}},$
(5)

where $p=0.3$, $\bar{a}=0.7$ and $\bar{A}$ is the normalisation factor. The DM
halo number density obtained from the universality function defined in
equation 5 is known as the Sheth-Tormen (ST) number density. In the case of
$\bar{a}=1$ and $p=0$, we will find the extended Press-Schechter mass function
for spherical collapse [13]. Note that the ratio of critical density in
ellipsoidal collapse $\delta_{\text{ell}}$ to spherical collapse
$\delta_{\text{sc}}$ model is defined in ST model as below [15]:

$\frac{\delta_{\text{ell}}}{\delta_{\text{sc}}}=\sqrt{\bar{a}}[1+\beta(\bar{a}\nu)^{-\alpha}],$
(6)

where $\bar{a}\approx 0.7$, $\alpha\approx 0.615$ and $\beta\approx 0.485$.
The halo mass function we used in this work for comparison with $N$-body
simulation data is from [24], which calibrates the mass function at $z=0$ for
virial masses in the range $10^{11}h^{-1}M_{\odot}\leq M\leq
10^{15}h^{-1}M_{\odot}$ to $5\%$:

$\frac{dn}{dM}=f(\sigma)\frac{\bar{\rho}_{m}}{M}\frac{d\ln\sigma^{-1}}{dM},$
(7)

where

$f(\sigma)=A\left[(\frac{\sigma}{b})^{-a}+1\right]e^{-c/\sigma^{2}},$ (8)

with the fitting functions for constants

$\displaystyle A$ $\displaystyle=$
$\displaystyle\begin{cases}0.1(\log\Delta)-0.05~{}~{}~{}&\Delta<1600,\\\
0.26~{}~{}~{}&\Delta\geq 1600,\end{cases}$ $\displaystyle a$ $\displaystyle=$
$\displaystyle 1.43+(\log\Delta-2.3)^{1.5},$ (9) $\displaystyle b$
$\displaystyle=$ $\displaystyle 1.0+(\log\Delta-1.6)^{-1.5},$ $\displaystyle
c$ $\displaystyle=$ $\displaystyle 1.2+(\log\Delta-2.35)^{1.6},$

where all logarithms are base 10 and $\Delta$ is the over-density, where in
this work, we set it equal to $200$. An important question to ask: what is the
crucial information in the initial condition which helps to predict the
distribution of dark matter? In linear perturbation theory through the
continuity and Euler equations, there is a relation between the density
contrast field and velocity. However, this assumption can be checked by
studying the normalized cross-power spectrum of velocity and density contrast
field. We define the normalized cross-power spectrum NCP as

$\rm{NCP}=\frac{P_{\delta\theta}(k)}{\sqrt{P_{\delta\delta}(k)P_{\theta\theta}(k)}},$
(10)

where $P_{\delta\delta}(k)$, $P_{\theta\theta}(k)$ are the density and
momentum $(p=\rho V)$ divergence power-spectra and $\theta$ is defined as
$\nabla.(1+\delta)V$. Using the momentum field instead of the velocity field
is to avoid complexities rises in defining a velocity field in simulation
voxels without any particle. $P_{\delta\theta}(k)$ is the cross power spectrum
of density and momentum divergence. In a linear regime by definition, we have
$NCP=1$, which means that all the information is in the density-contrast field
and the velocity field information is entirely constructible from the density
field using the Poisson equation. However, 1 shows a deviation from unity in
small scales, as well as when we’re increasing $A_{s}$ values. It means that
in N-body simulations, because of non-linear effects, we have access to
additional information in the velocity field. This is a point we use
throughout this work to justify our results.

Figure 1: The normalized cross power spectrum (NCP) is plotted versus
wavenumber for different values of $A_{s}$. The unity shows that all
information in the velocity field is present in the density field as well. The
deviation from unity shows that the non-linear effect becomes important.

In 1, we plot the NCP for different values of initial conditions. It shows
that by increasing the amplitude the scale of non-linear effects goes to
larger scales and we have more deviation than unity.

Figure 2: The variance of linear matter perturbations is plotted versus mass
(left panel) and versus the ratio of the critical density of ellipsoidal to
spherical (right panel). In the left panel, the variance is plotted for two
different initial power spectrum amplitudes. The green dotted vertical lines
show the scope of the mass range under study.

To conclude this section, we refer to figure 2, in which we plot the variance
of perturbation in terms of DM halo mass in the left panel, for two amplitudes
of the initial power spectrum. The vertical dashed green lines show the mass
range in our simulation. In the right panel, we plot the variance (same y-axes
for both panels) in terms of the ratio of the critical density of ellipsoidal
to the spherical model introduced in equation 6. Note that the variance is a
monotonic decreasing function of mass, and the critical density ratio is a
monotonic increasing function of mass. This is a crucial point to assert that
by increasing the amplitude of perturbations for a specific mass range (orange
dash-dotted line) we obtain higher variances. It means that by fixing the
variance, a mass range in the cosmology with higher initial amplitude
corresponds to a smaller mass scale in the $\Lambda\text{CDM}$ paradigm, where
the velocity shear field is crucial in the collapse model. [15, 25].
Accordingly, due to the right panel, we are probing the scales in which the
ratio of critical densitis are higher than unity. The unity ratio of critical
densities means that we have almost spherically collapsed halos, however large
ratio indicates that we should consider ellipsoidal collapse, in which the
velocities play an important role. This is a key idea in this work, which we
will discuss further in upcoming sections.

##  III $N$-body Simulations

Figure 3: Halo mass function for simulation boxes with different sizes and
theoretical halo mass function from the Tinker profile is plotted. For box
size $50{\text{ Mpc}}/h$ we have the validity range for halo masses up to
$2\%$ in comparison to the Tinker profile. This is for the mass range of
$10^{10.5}M_{\odot}$ to $10^{14}M_{\odot}$. For all the box sizes, the number
of simulating grids is 600, so the number of dark matter particles in the box
will be $600^{3}$.

To study the structure formation using deep learning methods, we use data
obtained from gevolution, a relativistic cosmological $N$-body simulation [26,
27]. In gevolution, the full set of Einstein’s equations is solved to update
the positions and velocities of particles. One of the main advantages of
gevolution is its natural extendability to consider relativistic particles
(e.g., massive neutrinos [28]) and non-standard cosmologies (e.g., different
dark energy models and modified gravity theories) [29, 30, 31]. The presence
of a horizon in relativistic equations compared with the Newtonian ones
provides a more accurate analysis of the structure formation at large scales.
The gevolution’s potential for covering the relativistic particles and seeing
their effects in the matter power spectrum on small scales as well as large
scales shows us a wide range of exciting possibilities to investigate the
effects of non-standard cosmological models (like assuming massive neutrinos
or considering different dark energy models) on the perturbation theory and
formation of structures in our future works.
Assuming the standard model of cosmology, $\Lambda$CDM, we run our simulations
with a standard cosmological model configuration. Our main training data
consists of 10 gevolution boxes with 50 ($\text{Mpc }h^{-1}$) box length and
$600^{3}$ DM particles within the box. We use the cosmological parameters
$\Omega_{\text{cdm}}h^{2}=0.122$, $\Omega_{\text{b}}h^{2}=0.021$,
$T_{\text{CMB}}=2.725$, $h=0.67$ , $A_{s}=2\times 10^{-9}$, $n_{s}=0.96$ and
$N_{\text{eff}}=3.046$. We made use of the public halo-finder code Rockstar
[32] to identify halos at $z=0$ in the box.
Due to the fixed mesh and finite resolution in our simulation data sets, there
is an inaccuracy in finding small-mass halos. Thus, to find the best mapping
between the initial condition of dark matter particles and their final halo
mass, we should recognise the resolution effect on the halo mass function. The
resolution limitation is not physical and is only due to the simulation or
halo-finder imprecisions. So we have to choose our training data-set from the
particles that end up in the halos with proper mass in our specific training
box. To find the resolution effect, we run four simulations with different
resolutions, where the number of particles is fixed to $600^{3}$ but different
box sizes (25, 50, 100, and 300 $\text{Mpc }h^{-1}$) are considered. The
overlapped region of the halo mass function for the different resolutions
shows the validity range of each simulation. In figure 3 we show the validity
range for each simulation. Based on the figure the validity range for the bbox
size50 $\text{ Mpc }h^{-1}$, which is the main simulation being used in this
study, within 2% accuracy is $10^{10.5}M_{\odot}$ to $10^{14}M_{\odot}$.

##  IV Method: Convolutional Neural network

In this section, we discuss the ML methods we use in this work. In the first
subsection IV.1, the convolutional neural network is discussed. In the
upcoming subsections, the image construction IV.2, architecture IV.3, loss
function IV.4, image size IV.5 and image resolution IV.6 are described
respectively.

Figure 4: We train our CNN on a set of cosmological simulations with different
phase space information. Cosmological simulations alone cannot provide us with
insight into the physical processes in structure formation and its relation to
the main parameters in the initial condition, which are important in halo mass
function. We use our neural net, which is trained on different sets of initial
condition information of the simulation box to provide an understanding of the
structure formation. We demonstrated our two kinds of training in this plot.
We once train our model on the images with the velocity field information
($15^{3}\times 4$ input shape) and once without($15^{3}\times 1$ input shape)
to investigate the effect of the velocity field in our machinery’s learning.

###  IV.1 Convolutional Neural Networks

Convolutional neural networks (CNNs) [33] are one of the most powerful types
of neural network models and are widely used for image processing
applications. Inspired by Neocognitrons [34], they benefit from locality-based
operations to hierarchically detect translation-invariant patterns and make
predictions based on them. Also, they use parameter sharing, which reduces
their complexity compared to traditional neural networks. This allows for
building extremely deep neural networks that can learn complex patterns
efficiently. From the early applications of CNNs in [35] and [36] to the
modern applications in life sciences [37], medical imaging [38, 39] and
autonomous driving [40]. CNN models have always been one of the driving forces
in the field of deep learning.

In recent years, applications of CNNs in cosmology have attracted a lot of
attention. This is partly motivated by the big data from cosmological
observations and large-scale simulations [41]. As a more specific example of
data-driven cosmology, we follow the line of the research in the LSS, using
$N$-body simulations as an arena to study the distribution of matter.
For instance, [42] used a CNN to estimate cosmological parameters directly
from zero-redshift density field images and outperformed conventional methods.
In [43], U-net architecture has been used to infer baryonic matter
distribution from less costly dark matter only $N$-body simulations. [44] used
a CNN model to learn moments of margthe inalized joint posterior of
cosmological and astrophysical parameters, and thereby performed likelihood-
free inference to work out $\Omega_{m}$ and $\sigma_{8}$ from 2D field maps.

Like other deep learning models, CNN-learned feature maps are not easy to
interpret. Thus, deriving theoretical conclusions from such models is
challenging. Some works have been dedicated to developing methods to interpret
neural networks [45, 46, 47] or to implement easily interpretable networks
[48]. It is also possible to shed some light on deep neural networks by
training them using dissimilar or feature-reducdata setsets and comparing the
results. For example, [19] tried to interpret the learned feature maps of a
CNN to assess the importance of anisotropic features in the prediction of dark
matter halo masses. We employed the latter method to conclude our CNN’s
trained parameters.

We develop a deep CNN to learn about the physical processes of structure
formation directly from $N$-body simulations. We train our CNN to learn the
halo mass that each particle resides in at the redshift zero. It is worth
mentioning that the halo mass is defined based on the number of bound
particles, thus there is no assumption on the shape of halos in our analysis.
Our CNN’s input is the density and velocity field information at $z=100$. We
use our trained network to investigate the effect of the velocity field in
halo collapse and the final halo mass function. In the upcoming subsections,
we first describe our data preparation and image construction and then give an
overview of the architecture of our CNN model.

###  IV.2 Image Construction

We use one of the most basic features of the simulation box, namely the phase
space information of the DM particles. We make a 3D cubic image of the local
neighbourhood of each particle at the initial snapshot of simulation boxes at
$z=100$. Each image consists of the particle itself at the central voxel. We
use $0.5\text{ Mpc }h^{-1}$ comoving size for resolution and an image size of
$15$ voxels corresponding to $7.5\text{ Mpc }h^{-1}$ comoving size for our
model. We prepare the density and velocity field from the phase space by
averaging the information of the particles for each voxel. Images consist of
four channels, one for the normalized number of particles in each image pixel
and three for the net velocity of particles in the three Cartesian
coordinates. We want our model to predict the mass of each particle’s final
halo destination at the last snapshot $z=0$. For this task, first, we should
fix the image resolution and the image size. The combination of $7.5\text{ Mpc
}h^{-1}$ image size and $0.5\text{ Mpc }h^{-1}$ resolution is the conservative
image properties that give us a reasonable run time while preserving the
essential information. We will study in more detail the choice of these two
parameters in the next section after describing our deep learning model and
method. The input of our ML model is the image constructed for each particle,
and the output is the particle’s final destination halo mass.

ITodecrease our machine’s sensitivity to different simulation realisations and
to reduce the effect of Poisson error, we use 10 independent simulations.
Since the number of large halos in the simulation dramatically decreases due
to small boxsize and mass resolution consideration, our model may face a
sample variance problem as it can not see enough halos at a large mass range.
Considering the independent simulations decreases the Poisson error due to the
small number of large halos by a factor of one. To reduce it by a larger
factor of two, we need 100 independent simulation boxes, which is challenging
based on our computational power. We pick 300,000 particles randomly from all
these boxes as our training sample. Increasing this number from 300,000 to
400,000 showed minor improvement in the result while adding to the
computational cost extensively.

We make the images for each particle in the way described above. We exclude
all the particles ending into a halo with masses less than the allowed halo
mass range by simulations in both training and testing sets. This choice of
minimum halo mass is compatible with the main goal of this work, namely,
studying the effect of the velocity field information. Note that the particles
ending in cosmic web structures other than halos such as voids and filaments
are not considered in this study; however, they are present in the image of
the particles. One will need a simulation with a larger resolution to broaden
the allowed halo mass range for a comprehensive study of structure formation.
Based on the physical properties of structure formation, we expect the
inferred halo masses to be invariant under reflection and rotation of the
input subregions. In order Forl to meet such symmetries, we use 48 members of
the cube symmetry group to replicate our input images. This has been shown to
effectively reduce the variance of predictions within the replicates of each
subregion in [42]. However, in order to speed up the training process, we do
not use the whole replicates to augment our data-set. Instead,data setthe
generation of each mini-batch, each input image is transformed by a randomly
selected member of the symmetry group. While preserving the possibility of
feeding our model with different transformations of each subregion, this
approach considerably reduces our training time.

###  IV.3 Architecture

Input: | Subregion around a particle | shape: | (15, 15, 15, 4)
---|---|---|---
layer | type | kernel size | padding | activation | param # | output shape
1 | Conv3D | (5, 5, 5) | same | sigmoid | 16032 | (15, 15, 15, 32)
2 | Conv3D | (3, 3, 1) | valid | relu | 9248 | (13, 13, 15, 32)
3 | Maxpooling3D | (2, 2, 1) | same | - | 0 | (7, 7, 15, 32)
4 | Conv3D | (3, 3, 1) | same | relu | 18496 | (7, 7, 15, 64)
5 | Conv3D | (4, 4, 1) | valid | relu | 65600 | (4, 4, 15, 64)
6 | Maxpooling3D | (2, 2, 1) | valid | - | 0 | (2, 2, 15, 64)
7 | Dropout (0.25) | - | - | - | 0 | (2, 2, 15, 64)
8 | Flatten | - | - | - | 0 | (3840)
9 | Dense | - | - | relu | 460920 | (120)
10 | Dense | - | - | relu | 7260 | (60)
11 | Dense | - | - | linear | 61 | (1)
Output: | Destination halo mass | shape: | (1)

Table 1: The architecture of the sequential Convolutional Neural Network
model.

Our model is a deep convolutional neural network which its architecture is
described in table 1. The model has two parts: a contracting path and a fully
connected component. The contracting path consists of two scales. At each
scale, the input is passed through two convolutional layers, followed by a max
pooling layer. In each convolutional layer, a set of filters are convolved
with the input and the result is passed through an activation function. These
operations are used to extract features from the input. Max pooling layers
down-sample the feature maps on each scale while preserving their important
structural information. This is done by sliding a kernel across each image
channel and putting the maximum value of each resulting region in the
corresponding location of the respective output channel. Fully connected
layers use the feature maps extracted by the contracting path to infer the
scalar value of the halo mass to which the central particle of the input image
will finally belong (see figure 4).

###  IV.4 Loss Function

The objective loss function used for the model is mean squared errora (MSE)
which is minimised using the Adam optimisation algorithm [49]. As we will
demonstrate later, our model with this symmetric loss function introduces a
systematic error, which causes us to over-estimate thoverestimatet smaller
mass ranges.
In figure 5, we see an over-prediction for the mass of smaller halos. This
feature is present in all samples, indicating an intrinsic error in the ML
method. Accordingly, as the over-prediction increased with a higher minimum
mass range, we will not restrict our sample with this consideration. Even
more, we use the smaller mass range which is not in the range of our model
prediction. This is done in order to transfer the over-prediction error to a
smaller mass range which is out of the scope of the study. As demonstrated in
figure 5, we can see that changing the minimum halo mass given to the CNN
model does not affect the prediction of the model more than the systematic
error due to the symmetric loss function that we discussed above. This means
that predictions of the halo mass are not affected by slightly changing the
minimum used in our training and thus it is independent of our choice of halo
mass range in the training set.

Figure 5: The mass of the predicted dark matter halos to the real mass sample
is plotted for three minimum mass ranges of the machine entry data. The error
bars show the mean and the $\pm 25\%$ of the distribution. Different colours
show the minimum mass which is chosen for ML. The presence of the non-
symmetric outliers shows the systematic overprediction of our model at the
smallest halo ranges. We resolve this issue by keeping the complete halo mass
range of the box while only trusting the results on an interval. This
technique pushes the over-prediction to the mass ranges, which are out of the
range of interest.

###  IV.5 Image size

Figure 6: The mass of predicted dark matter halo concerning real one is
plotted versus mass for different image sizes. The different colours
correspond to different image sizes ina units of $\text{Mpc }h^{-1}$. The
error bars show the $\pm 25\%$ of the prediction distribution in each mass
bin.

One of the critical hyperparameters of our model is the physical size of the
input images. The initial density contrast field at $z=100$ is almost
homogeneous and less than unity. Furthermore, for a desired final halo, we can
approximate the initial proto-halo (ph) size as
$V_{\text{ph}}=M_{\text{halo}}/\bar{\rho}_{i}$, where $\bar{\rho}_{i}$ is the
average density of the simulation box at $z=100$. As a result, the image size
should be at least the volume needed to contain all the most massive halo
particles. There is a trade-off in choosing our image size. Our mapping
between these images to the halo mass only aims to find the future host halo’s
mass of the central part of the image. Besides their substantial computational
cost, large-volume images may mislead the machine when we are not on the
density field peaks, and many proto-halo regions exist within a single image.
Also, we can not choose our images too small, as we miss some information
necessary to find the future halo mass of the specified region in the initial
box. We do some scale tests for choosing the best possible image size. We
change the image size slightly around the value we obtained as described and
compare our model’s performance with different image sizes. Figure 6 indicates
that the error due to image size is almost the same for all sizes that we
investigate. In general, we can categorise the systematic errors involved in
our machinery as a) error due to the resolution of the simulation (see figure
3), b) the systematic over-prediction error (see figure 5), and c) error due
to the characteristic of the image (figure 6, figure 7). We set an upper limit
of $20\%$ for our uncertainty, and for the rest of the work, we choose a
specific set of hyperparameters for our images concerning the three sources of
error mentioned earlier. This means the minimum mass is set to the minimum
halo mass in simulations as described in section IV.4, and the image size is
set to $7.5\text{ Mpc }h^{-1}$ ($15$ voxels).

###  IV.6 Image resolution

Figure 7:  The relative error and R2-score of the predicted final halo mass
with respect to the real one is plotted versus mass for different image
resolution of $100^{3}$ and $125^{3}$ voxels for a $50\text{ Mpc }h^{-1}$ box.

After fixing the appropriate image size, we specify the image resolution for
the CNN input. The volume of the minimum halo mass in the initial condition is
the lower bound for the resolution. Thus we pick a value using this minimum
bound, and then we change our model’s resolution around this specific value to
reach the best resolution. As we show in figure 7, increasing the image’s
resolution by a factor of $1.25$ does not change the overall accuracy of the
model in the desired mass ranges noticeably. Note that this is done while
keeping the physical length of the image constant. However, this increase in
the resolution adds to the computational cost by more than a factor of two.
Furthermore, resolution $100$ (physical voxel size of $0.5\text{ Mpc }h^{-1}$)
is a convenient choice.

##  V Results

We use the basic phase space information of the initial simulation box with a
finite resolution smoothed with a cubic window function, as the input of our
CNN model. We make images of the particle’s local neighbourhood and train our
machine to predict the mass of the final halo to which this particle ends up
from its local neighbourhood picture. We use a training set of $300,000$
images randomly chosen from 10 different realisations of the same cosmological
parameters. Using the results described in section IV to fix the
hyperparameters of the CNN model, we set the image size to be $7.5\text{ Mpc
}h^{-1}$ or $15$ pixels and with a resolution of $0.5\text{ Mpc }h^{-1}$. Our
results can be separated into two main parts. First, we use our model to
recover the halo mass function of a completely new box but with the same
cosmological parameters as the training boxes in subsection V.1. Then we
present a physical interpretation of the structure formation using our
framework based on the effect of the initial velocity field in the prediction
of HMF in subsection V.2.

###  V.1 Halo Mass Function Prediction

The aim of constructing an accurate predictor for halo masses is to eventually
predict the halo mass function, given an initial distribution of particles.
Furthermore, we use the distance of our predicted halo mass function from the
real halo mass function as a metric to measure the performance of our model.
To predict the halo mass function, we first need to apply a binning on our
regression problem and then find the predicted number of halos corresponding
to each halo mass bin. As the prediction the of halo mass function is strongly
dependent on the bin number, we should be careful about our binning
definition. We define 40 bins between $10.5\leq\log(M/M_{\odot})\leq 14$
linear in the logarithm of mass. By applying the trained model on a set of
$280,000$ randomly chosen particles in our test box and predicting the halo
mass for each of them, we already acquired the number of particles in each
halo mass bin. Next, we use the average density of each halo mass bin to
convert the number of particles to the number of halos in that bin. One also
must scale this result by the relative ratio of the total number of randomly
chosen particles used in this prediction over the total number of particles in
the box. The main two contributing ftotors in the error of this distribution
are first, the systematic error in the prediction of halo masses arising from
our CNN model’s imprecision, and the second error comes from the lack of
enough data in chosen random point. Each error is shown in figure 8, where we
demonstrated the predicted halo mass function of the box. At our ending mass
bins, where the sample variance due to the lack of most massive halos
increases, our confidence in prediction drops. The colour bar in this plot
corresponds to the statistical error with an n-sigma deviation of our
prediction from the real HMF of the box which is our metric. Our prediction
for the HMF is in agreement with the real one from the simulations in the most
mass ranges. However, our prediction ability drops considerably for the mass
ranges of $\log(M/M_{\odot})<11$, which is because of the effect we described
in section IV. This result shows that our model prediction is trustable the
most in the range of $11<\log(M/M_{\odot})<13$.

Figure 8: The number density of predicted and real dark matter halos plotted
versus mass in the upper panel. The colour bar shows the $\sigma$-deviation of
our predicted model compared to the real data. In the lower panel, we show the
ratio of the number density, and the shaded region shows the $3\%$ confidence
level. In both panels, we show the model which considers the velocity as well
as density field as an input.

###  V.2 The role of the velocity field in structure formation

To investigate the effectthe of velocity field in the process of dark matter
halo formation, we train our model on different sets of input features. One
only consists of the density field information, while the other also consists
of each particle’s velocities in the three Cartesian coordinates. We change
our CNN model in each case to fit the input image shape ( 15,15,15,1) for the
case without velocity information and (15,15,15,4) for when the velocity
information is also provided. Except for the first layer, which is changed to
fit the input’s shape, all the other layer’s architecture remains the same.

Figure 9: The predicted mass of dark matter halos versus the real mass is
plotted. The violins show the distribution of our prediction for two different
sets of training data with velocity (blue/left) information and the other
without velocity (orange/right) information.

In figure 9, we show the prediction of halo mass with and without considering
the velocity information. Violins represent two different sources of error in
each halo mass bin. One is the distance of the solid line to the distribution
average for the predicted halo mass in each mmass bin(the accuracy of the
result), and the other is the variance of the distribution for the predicted
halo mass (precision of our prediction). To consider these two error sources,
we use the R2-score, which would be a reasonable metric to see whether there
is an improvement in the prediction with and without velocity in each halo
mass bin.
The advantage of this metric is that by dividing the model’s standard
deviation, it measures the score relative to the model itself and thus enables
us to compare the score of the two different models. We use R2-score
measurement to quantify the performance of our model, where the R2-score is
measured in each mass scale as:

$\displaystyle\text{R2-score}(y_{\text{real}},y_{\text{pred}})=1-\frac{\sum_{i}(y_{\text{real}}^{(i)}-y_{\text{pred}}^{(i)})^{2}}{\sum_{i}(y_{\text{real}}^{(i)}-\bar{y})^{2}}\
,$

where $y_{\text{real}}^{(i)}$ is the real mass of the particle’s destination
halo and $y_{\text{pred}}^{(i)}$ is the prediction of our model for it, and
the $\bar{y}$ is the average $y_{\text{pred}}$ on that mass bin. We measure
the R2-score on different mass scales for the two configurations. One with and
the other without velocity information in the case of standard initial power
amplitude $A_{s}=2\times 10^{-9}$. Note that the overall R2-score for the
models would be a value between 0 and 1, as expected from its definition,
while when we measure the R2-score on individual bins, it can get other values
as well. Being negative on a bin would only mean we use small bin sizes
compared to the complete mass range.

We use the cross-validation technique to find the error bars on our model. In
each case, we train our model five separate times on randomly chosen training
sets of the same size. We measure the R2-score for each of these different
trained models at each mass scale and use the mean and standard deviation of
the bootstrapped sample as our model’s mean and uncertainty in terms of the
R2-score.

Figure 10:  R2-score for the two with and without velocity models. The result
shows that there is a negligible difference between the performance of these
two models in the standard model of cosmology with initial curvature
perturbations of $A_{s}=2\times 10^{-9}$.

Based on figure 10 the mean R2-score value does not change considerably for
the model without velocity in comparison with the one with velocity field
information. However, the uncertainty of prediction (standard deviation of the
bootstrapped sample) reduces for the model with velocity field information.
Our result in figure 10, demonstrates that in the mass range
$10.5\lesssim\log(M/M_{\odot})\lesssim 14$, there is no compelling evidence on
the importance of velocities in our predictions. However, due to our
discussion in section II, we anticipate that the velocities should be
important in small mass scales where the non-linear effects are no longer
negligible. Nonetheless, due to the mass resolution limitations of our
simulations, this regime is not accessible to our current computational power
and configuration. To overcome this obstacle, we use an idea to check the
velocity dependence of the structure formation. We increase the initial power
spectrum of the perturbations in the $N$-body simulation runs. This causes
dark matter halos in our simulation mimic the velocity field in smaller mass
range halos due to their higher over-density, as discussed in section II.

Figure 11:  The R2-score statistics iareplotted versus mass for three
different initial power spectrums in the left panel. In the right panel, we
show the distribution ofthe dark matter particle’s velocThe rightRight panel
indicates that the peak of the velocity distribution of particles moves to
larger values when we increase the initial curvature perturbations.

Again, we train our model five times for the new set of simulations with a
higher initial power spectrum, and we compare the R2-score of bootstrapped
cross-validated samples. The result is demonstrated in figure 11. As before,
the solid lines are the extrapolation of the mean of the R2-score in each mass
bin and the shaded regions are the 1$\sigma$ uncertainty obtained from
bootstrap. The right panel in figure 11 indicates that by increasing the
initial power spectrum, the distribution of the velocities shifts to higher
values. This implies that by increasing the initial power spectrum we can
bring the ethe ffect of velocity field to larger mass scales which effectively
means that by looking at a fixed halo with a specific mass in this non-
standard case we are looking at smaller mass dark matter halos compared to the
standard model case. The left panel of figure 11 shows that the precision of
the model without velocity field drops significantly when the initial power
spectrum increases. Moreover, the difference between the model trained with
velocity field information and the one without it increases noticeably when
increasing the initial power. To show this more explicitly, we measure the
difference between the two models for each case. As both models (with and
without velocity channels) are trained on the same data-set, the standard
deviation in each mass bin is the same. Thus using the mean squared error
(MSE) of each model should be sufficient to compare their performances. We
measure the relative difference of the two models in each mass bin,
$d(M_{\text{bin}})$, using:

$\displaystyle
d(M_{\text{bin}})=\frac{|\text{MSE}_{\text{without}}(M_{\text{bin}})-\text{MSE}_{\text{with}}(M_{\text{bin}})|}{\text{MSE}_{\text{with}}(M_{\text{bin}})},$

where,

$\displaystyle\text{MSE}(M_{\text{bin}})=\frac{1}{N}\sum_{i=1}^{N}(M_{\text{pred}}^{(i)}-M_{\text{real}}^{(i)})^{2},$

and $N$ is the total number in particles on that specific mass bin. Figure 12
shows the result of this the distance for different amplitude of initial
curvature perturbations, where the errors are obtained from cross-validation.
The figure implies that the prediction of the model with velocity information
is much better than the one without velocity for the large amplitude of
initial curvature perturbations. This is consistent with our previous
discussion that higher initial amplitudes bring the effect of velocities into
larger mass halos compared to the $\Lambda\text{CDM}$ case. This is done by
increasing the ratio $\delta_{\text{ell}}/\delta_{\text{sc}}$ at larger mass
halos to mock the effect of velocity shear field which is visible in the small
mass halos in the $\Lambda\text{CDM}$ paradigm. As a result, we assert that
the velocity information is important on small scales. Our result obtained
from a deep learning framework which is independent of the collapse mdel,
illustrates that the ellipsoidal collapse model will diverge from the
spherical collapse model on small scale in the standard model.

Figure 12: The difference between the two (with and without velocity) models
for different initial curvature perturbations are plotted versus mass. By
increasing the power to $A_{s}=8\times 10^{-8}$ the difference betthe ween two
models increases to more than $80\%$. This implies that for higher $A_{s}$
values the density field information is not enough to predict the halo mass
function accurately.

##  VI Conclusions and Future Remarks

We developed an interpretable deep-learning framework to shed light on the
process of structure formation. We train a CNN model on an initial simulation
box at $z=100$ to find the final destination halo mass for each DM particle.
Our method, in comparison with its previous studies, has the advantage of
being simple in terms of the initial features. We use our prediction for the
halo mass function as the indicator of the performance of our machinery. Our
results in terms of the halo mass function are in agreement with the one from
the simulation box up to $1\sigma$ in the halo mass range of
$11<\log(M/M_{\odot})<13$. We use our framework to find the importance of the
velocity field in the structure formation process. We train our CNN on two
different sets of samples: one on the images with velocity channels and the
other without velocity channels. Our results imply that the effect of the
velocity field in the standard model of cosmology given the initial curvature
perturbations $A_{s}=2\times 10^{-9}$ is negligible in the mass range of
$10.5<\log(M/M_{\odot})<14$. To study the effect of the velocity field on
smaller scales we use a novel idea. We increase the initial curvature
perturbations to increase the variance of perturbations in the scales that are
in the scope of our simulation. This causes massive halos to mimic the effect
of velocity fields of a small mass range. Accordingly, the effect of velocity
becomes visible in larger mass halos of our simulations with the given
resolution. We observe that increasing the initial curvature perturbations
increases the effect of the velocity field in the structure formation process
in the given halo mass range, correspondingly. We also discuss, why we
anticipate that in small scales and non-linear regimes, we have more
information in the velocity field added to the density contrast field.
Although, at first glance, they are related by the continuity equation. In
this direction, we show that the normalized cross power-spectrum deviates from
unity on small scales. By increasing the amplitude of the initial conditions
the deviation from unity occurs in larger scales in comparison. This result
implies that the effect of the velocity field in structure formation is scale-
dependent and increases in the smaller mass scales. Our work shows the vast
potential of Convolutional Neural Networks to extract higher-order information
from basic phase space information and their capability to help us investigate
fundamental physical questions due to their simplicity and interpretability.
It is noteworthy that the halo mass, in our analysis, is defined based on the
number of bound particles and we do not have any pre-assumptions on the shape
of halos. So, there is no constraint on our learning process to select any
specific collapse model. Accordingly, our results indicate the full
consideration of an ellipsoidal collapse may be essential to reach a more
viable halo mass function on small scales, and the simple spherical collapse
model would not be sufficient. Accordingly, it will be interesting to study
the effect of different smoothing window functions on the halo mass prediction
in a separate work.
Future work might study the correlation between the accuracy of halo mass
prediction and the particle’s distance from its host halo centre. This is to
integrate the idea of peak theory and excursion set theory of peaks in our ML
methods. We can use our deep learning method to probe other cosmologies; for
example, the models that velocity may play an important role in large mass
halo formation. Also, this method can be used to investigate the effect of the
velocity field on the structure formation in extensions of the standard model
(e.g., massive neutrinos) and alternative models (e.g., $k$-essence models).

Acknowledgments: We would like to thank Francisco Villaescusa-Navarro for his
helpful comment and discussion on the normalized cross-power spectrum. We also
thank Arya Farahi for many useful discussions. This work is supported by Oslo
University computational facilities and a grant from the Swiss National
Supercomputing Centre (CSCS) under project ID s1051. SB is partially supported
by the Abdus Salam International Center of Theoretical Physics (ICTP) under
the junior associateship scheme. This research is supported by Sharif
University of Technology Office of Vice President for Research under Grant No.
G960202.

## References

  * [1] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [2] W. J. Percival et al. [SDSS], Mon. Not. Roy. Astron. Soc. 401 (2010), 2148-2168 doi:10.1111/j.1365-2966.2009.15812.x [arXiv:0907.1660 [astro-ph.CO]].
  * [3] S. Alam et al. [BOSS], Mon. Not. Roy. Astron. Soc. 470 (2017) no.3, 2617-2652 doi:10.1093/mnras/stx721 [arXiv:1607.03155 [astro-ph.CO]].
  * [4] H. Camacho et al. [DES], Mon. Not. Roy. Astron. Soc. 487 (2019) no.3, 3870-3883 doi:10.1093/mnras/stz1514 [arXiv:1807.10163 [astro-ph.CO]].
  * [5] J. Rhodes, D. Alonso, B. Ansarinejad, R. Armstrong, J. Asorey, A. Avelino, J. Blazek, F. J. Castander, R. R. Chary and X. Chen, et al. Bull. Am. Astron. Soc. 51 (2019) no.3, 201 FERMILAB-PUB-19-741-PPD.
  * [6] S. D. M. White and M. J. Rees, Mon. Not. Roy. Astron. Soc. 183 (1978), 341-358 doi:10.1093/mnras/183.3.341
  * [7] P. J. E. Peebles, Large-Scale Structure of the Universe, Princeton University Press, 1980.
  * [8] V. Springel, N. Yoshida and S. D. M. White, New Astron. 6 (2001), 79 doi:10.1016/S1384-1076(01)00042-2 [arXiv:astro-ph/0003162 [astro-ph]].
  * [9] V. Springel, S. D. M. White, A. Jenkins, C. S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton and J. Helly, et al. Nature 435 (2005), 629-636 doi:10.1038/nature03597 [arXiv:astro-ph/0504097 [astro-ph]].
  * [10] M. Kuhlen, M. Vogelsberger and R. Angulo, Phys. Dark Univ. 1 (2012), 50-93 doi:10.1016/j.dark.2012.10.002 [arXiv:1209.5745 [astro-ph.CO]].
  * [11] W. H. Press and P. Schechter, Astrophys. J. 187 (1974), 425-438 doi:10.1086/152650
  * [12] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304 (1986), 15-61 doi:10.1086/164143
  * [13] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, Astrophys. J. 379 (1991), 440 doi:10.1086/170520
  * [14] R. K. Sheth, H. J. Mo and G. Tormen, Mon. Not. Roy. Astron. Soc. 323 (2001), 1 doi:10.1046/j.1365-8711.2001.04006.x [arXiv:astro-ph/9907024 [astro-ph]].
  * [15] R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 329 (2002), 61 doi:10.1046/j.1365-8711.2002.04950.x [arXiv:astro-ph/0105113 [astro-ph]].
  * [16] F. Nikakhtar, M. Ayromlou, S. Baghram, S. Rahvar, M. R. Rahimi Tabar and R. K. Sheth, Mon. Not. Roy. Astron. Soc. 478 (2018) no.4, 5296-5300 doi:10.1093/mnras/sty1415 [arXiv:1802.04207 [astro-ph.CO]].
  * [17] S. He, Y. Li, Y. Feng, S. Ho, S. Ravanbakhsh, W. Chen and B. Póczos, Proc. Nat. Acad. Sci. 116 (2019) no.28, 13825-13832 doi:10.1073/pnas.1821458116 [arXiv:1811.06533 [astro-ph.CO]].
  * [18] L. Lucie-Smith, H. V. Peiris and A. Pontzen, Mon. Not. Roy. Astron. Soc. 490 (2019) no.1, 331-342 doi:10.1093/mnras/stz2599 [arXiv:1906.06339 [astro-ph.CO]].
  * [19] L. Lucie-Smith, H. V. Peiris, A. Pontzen, B. Nord and J. Thiyagalingam, [arXiv:2011.10577 [astro-ph.CO]].
  * [20] R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 308 (1999), 119 doi:10.1046/j.1365-8711.1999.02692.x [arXiv:astro-ph/9901122 [astro-ph]].
  * [21] J. E. Gunn and J. R. Gott, III, Astrophys. J. 176 (1972), 1-19 doi:10.1086/151605
  * [22] S. Dodelson, Academic Press, 2003, ISBN 978-0-12-219141-1
  * [23] F. Nikakhtar and S. Baghram, Phys. Rev. D 96 (2017) no.4, 043524 doi:10.1103/PhysRevD.96.043524 [arXiv:1606.06666 [astro-ph.CO]].
  * [24] J. L. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. S. Warren, G. Yepes, S. Gottlober and D. E. Holz, Astrophys. J. 688 (2008), 709-728 doi:10.1086/591439 [arXiv:0803.2706 [astro-ph]].
  * [25] H. Kameli and S. Baghram, Mon. Not. Roy. Astron. Soc. 494 (2020) no.4, 4907-4913 doi:10.1093/mnras/staa1058 [arXiv:1912.12278 [astro-ph.CO]].
  * [26] J. Adamek, D. Daverio, R. Durrer and M. Kunz, Nature Phys. 12 (2016), 346-349 doi:10.1038/nphys3673 [arXiv:1509.01699 [astro-ph.CO]].
  * [27] J. Adamek, D. Daverio, R. Durrer and M. Kunz, JCAP 07 (2016), 053 doi:10.1088/1475-7516/2016/07/053 [arXiv:1604.06065 [astro-ph.CO]].
  * [28] J. Adamek, R. Durrer and M. Kunz, JCAP 11 (2017), 004 doi:10.1088/1475-7516/2017/11/004 [arXiv:1707.06938 [astro-ph.CO]].
  * [29] F. Hassani, J. Adamek, M. Kunz and F. Vernizzi, JCAP 12 (2019), 011 doi:10.1088/1475-7516/2019/12/011 [arXiv:1910.01104 [astro-ph.CO]].
  * [30] F. Hassani, B. L’Huillier, A. Shafieloo, M. Kunz and J. Adamek, JCAP 04 (2020), 039 doi:10.1088/1475-7516/2020/04/039 [arXiv:1910.01105 [astro-ph.CO]].
  * [31] F. Hassani and L. Lombriser, Mon. Not. Roy. Astron. Soc. 497 (2020) no.2, 1885-1894 doi:10.1093/mnras/staa2083 [arXiv:2003.05927 [astro-ph.CO]].
  * [32] P. S. Behroozi, R. H. Wechsler and H. Y. Wu, Astrophys. J. 762 (2013), 109 doi:10.1088/0004-637X/762/2/109 [arXiv:1110.4372 [astro-ph.CO]].
  * [33] T. Homma, L. E. Atlas, R. J. Marks, An Artificial Neural Network for Spatiotemporal: Application to Phoneme Classification, MIT Press, 1987.
  * [34] K. Fukushima, Biological Cybernetics 36 (1980), no. 4, 193.
  * [35] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, L. D. Jackel, Neural Computation 1 (1989), no. 4, 541.
  * [36] Y. Lecun, L. Bottou, Y. Bengio, P. Haffner, Proceedings of the IEEE 86 (1998), no. 11, 2278.
  * [37] B. Ramsundar, P. Eastman, P. Walters, V. Pande, Deep learning for the life sciences: Applying deep learning to genomics, microscopy, drug discovery, and more, O’Reilly, 2019.
  * [38] A. S. Lundervold, A. Lundervold, Zeitschrift für Medizinische Physik 29 (2019), no. 2, 102.
  * [39] R. Yamashita, M. Nishio, R. K. G. Do, k. Togashi, Insights into Imaging 9 (2018), no. 4, 611.
  * [40] S. Grigorescu, B. Trasnea, T. Cocias, G. Macesanu, Journal of Field Robotics 37 (2020), no. 3, 362.
  * [41] A. Siemiginowska et. al., [arXiv:1903.06796 [astro-ph.CO]], (2019).
  * [42] S. Ravanbakhsh, J. Oliva, S. Fromenteau, L. C. Price, S. Ho, J. Schneider and B. Poczos, [arXiv:1711.02033 [astro-ph.CO]].
  * [43] X. Zhang, Y. Wang, W. Zhang, Y. Sun, S. He, G. Contardo, F. Villaescusa-Navarro and S. Ho, [arXiv:1902.05965 [astro-ph.CO]].
  * [44] F. Villaescusa-Navarro, D. Anglés-Alcázar, S. Genel, D. N. Spergel, Y. Li, B. Wandelt, A. Nicola, L. Thiele, S. Hassan and J. M. Z. Matilla, et al. [arXiv:2109.09747 [astro-ph.CO]].
  * [45] M. T. Ribeiro, S. Singh, C. Guestrin, [arXiv:1602.04938].
  * [46] A. Mahendran, A. Vedaldi, International Journal of Computer Vision 120 (2016), no. 3, 233.
  * [47] G. Montavon, W. Samek, K. Müller, Digital Signal Processing 73 (2018), 1.
  * [48] Q. Zhang, Y. N. Wu, S. Zhu, [arXiv:1710.00935].
  * [49] D. P. Kingma, J. Ba, [arXiv:1412.6980].