# Performance Evaluation of Big Data Processing Strategies for Neuroimaging

Valérie Hayot-Sasson1, Shawn T Brown2 and Tristan Glatard1
1Department of Computer Science and Software Engineering, Concordia
University, Montreal, Canada
2Montreal Neurological Institute, McGill University, Montreal, Canada

###### Abstract

Neuroimaging datasets are rapidly growing in size as a result of advancements
in image acquisition methods, open-science and data sharing. However, the
adoption of Big Data processing strategies by neuroimaging processing engines
remains limited. Here, we evaluate three Big Data processing strategies (in-
memory computing, data locality and lazy evaluation) on typical neuroimaging
use cases, represented by the BigBrain dataset. We contrast these various
strategies using Apache Spark and Nipype as our representative Big Data and
neuroimaging processing engines, on Dell EMC’s Top-500 cluster. Big Data
thresholds were modeled by comparing the data-write rate of the application to
the filesystem bandwidth and number of concurrent processes. This model
acknowledges the fact that page caching provided by the Linux kernel is
critical to the performance of Big Data applications. Results show that in-
memory computing alone speeds-up executions by a factor of up to 1.6, whereas
when combined with data locality, this factor reaches 5.3. Lazy evaluation
strategies were found to increase the likelihood of cache hits, further
improving processing time. Such important speed-up values are likely to be
observed on typical image processing operations performed on images of size
larger than 75GB. A ballpark speculation from our model showed that in-memory
computing alone will not speed-up current functional MRI analyses unless
coupled with data locality and processing around 280 subjects concurrently.
Furthermore, we observe that emulating in-memory computing using in-memory
file systems (tmpfs) does not reach the performance of an in-memory engine,
presumably due to swapping to disk and the lack of data cleanup. We conclude
that Big Data processing strategies are worth developing for neuroimaging
applications.

##  I Introduction

Big Data processing engines have significantly improved the performance of Big
Data applications by diminishing the amount of data movement that occurs
during the execution of an application. Locality-aware scheduling, introduced
by the MapReduce [1] framework, reduced the overall costs associated to
network transfer of data by scheduling tasks to the nodes located nearest to
the data. Reduction of data movement was further improved upon through in-
memory computing [2], which ensures that data is maintained in memory between
tasks whenever possible. To further reduce the cost of data movement, lazy
evaluation, the process of performing computations only when invoked, was
leveraged by Big Data frameworks to enable further optimizations such as
regrouping of tasks and computing only what is necessary.

Frameworks such as MapReduce and Spark [2] have become mainstream tools for
data analytics, although many others, such as Dask [3], are emerging.
Meanwhile, several scientific domains including bioinformatics, physics or
astronomy, have entered the Big Data era due to increasing data volumes and
variety. Nevertheless, the adoption of Big Data engines for scientific data
analysis remains limited, perhaps due to the widespread availability of
scientific processing engines such as Pegasus [4] or Taverna [5], and the
adaptations required in Big Data processing engines for scientific computing.

Scientific applications differ from typical Big Data use cases, which might
explain the remaining gap between Big Data and scientific engines. While Big
Data applications mostly target text processing (e.g. Web search, frequent
pattern mining, recommender systems [6]) implemented in consistent software
libraries, scientific applications often involve binary data such as images
and signals, processed by a sequence of command-line/containerized tools using
a mix of programming languages (C, Fortran, Python, shell scripts), referred
to as workflows or pipelines. With respect to infrastructure, Big Data
applications commonly run on clouds or dedicated commodity clusters with
locality-aware file systems such as the Hadoop Distributed File System (HDFS
[7]), whereas scientific applications are usually deployed on large, shared
clusters where data is transferred between data and compute nodes through
shared file systems such as Lustre [8]. Such differences in applications and
infrastructure have important consequences. To mention only one, in-memory
computing requires instrumentation to be applied to command-line tools.

Technological advances of the past decade, in particular page caching in the
Linux kernel [9], in-memory file systems (tmpfs) and memory-mapped files might
also explain the lack of adoption of Big Data engines for scientific
applications. In such configurations, in-memory computing would be a feature
provided by the operating system rather than by the engine itself. The
frontier between these two components is blurred and needs to be clarified.

Our primary field of interest, Neuroimaging, is no exception to the
generalized rise of data volumes in science due to the joint increase of image
resolution and subject cohort sizes [10]. Processing engines have been
developed with neuroinformatics applications in mind, for instance Nipype [11]
or the Pipeline System for Octave and Matlab (PSOM [12]). Big Data engines
have also been used for neuroimaging applications, including the Thunder
project [13] and in more specific works such as [14]. However, no quantitative
performance evaluation has been conducted on neuroimaging applications to
assess the added-value of Big Data engines compared to traditional processing
engines.

This paper addresses the following questions:

  1. 1.

What is the effect of in-memory computing, lazy evaluation and data locality
on current neuroimaging applications?

  2. 2.

Can in-memory computing be effectively enabled by the operating system rather
than the data processing engine?

Answers to these questions have important implications. In [15], a comparative
study of Dask, Spark, TensorFlow, MyriaDB, and SciDB on neuroinformatics use-
cases is presented. It concludes that these systems need to be extended to
better address scientific code integration, data partitioning, data formats
and system tuning. We argue that such efforts should only be conducted if
substantial performance improvements are expected from in-memory computing,
lazy evaluation or data locality. On the other hand, neuroimaging data
processing engines are still being developed, and the question remains whether
these projects should just be migrated to Spark, Dask, or other Big Data
engines.

Our study focuses on performance. We intentionally do not compare Big Data and
scientific data processing engines on the grounds of workflow language
expressivity, fault-tolerance, provenance capture and representation,
portability or reproducibility, which are otherwise critical concerns,
addressed for instance in [16]. Besides, our study of performance focuses on
the impact of data writes and transfers. It purposely leaves out task
scheduling to computing resources, to focus on the understanding of data
writes and movement. Task scheduling will be part of our discussion, however.

In terms of infrastructure, we focus on the case of High-Performance Computing
(HPC) clusters that are typically available through University facilities or
national computing infrastructures such as XSEDE, Compute Canada or PRACE, as
neuroscientists typically use such platforms. We assume that HPC systems are
multi-tenant, that compute nodes are accessible through a batch scheduler, and
that a file system shared among the compute nodes is available. We
intentionally did not consider distributed, HDFS-like file systems, as
initiatives to deploy them in HPC centers, for instance Hadoop on-demand [17],
have not become mainstream yet.

Our methods, including performance models, processing engines, applications
and infrastructure used, are described in Section II. Section III presents our
results which we discuss in Section IV along with the two research questions
mentioned previously. Section V concludes on the relevance of Big Data
processing strategies for neuroimaging applications.

##  II Materials and Methods

The application pipelines, benchmarks, performance data, and analysis scripts
used to implement the methods described hereafter are all available at
https://github.com/big-data-lab-team/paper-in-mem-locality for further
inspection and reproducibility. Links to the processing engines and processed
data are provided in the text.

###  II-A Engines

####  II-A1 Apache Spark

Apache Spark is a well-established Scala-based processing framework for Big
Data, with APIs in Java, Python and R. Its generalized nature allows Spark to
not only be applied to batch workflows, but also SQL queries, iterative
machine learning applications, data streaming and graph processing. Spark’s
main features include data locality, in-memory processing and lazy evaluation,
which it achieves through its principle abstraction, the Resilient Distributed
Dataset (RDD).

An RDD is an immutable parallel data structure that achieves fault-tolerance
through the concept of lineage [18]. Rather than permitting fine-grained
transformations, only coarse-grained transformations (applying to many
elements in the RDD), can be applied, thereby making it simple to maintain a
log of data modifications. This log, known as the lineage, is used to
reproduce any lost data modifications.

Two types of operations can be performed on RDDs: transformations and actions.
Applying a transformation to an RDD produces a new child RDD through a narrow
or wide dependency. A narrow dependency signifies that the child is only
dependent on a single parent partition, whereas a child RDD is dependent on
all parent partitions in a wide dependency. Examples of transformations
include map, filter and join. To materialize an RDD, an action must be
performed, such as a reduce or a collect. Lazy evaluation is represented in
Spark through the use of transformations and actions. A series of
transformations may be defined without the data ever being materialized. Using
this strategy, Spark can optimize data processing throughout the application.

All actions and wide dependencies require a shuffle – Spark’s most costly
operation. Every shuffle begins with each map task saving its data to local
files for fault tolerance. The shuffle operation then redistributes the data
across the partitions as requested. A shuffle marks a stage boundary in Spark,
where reduce-like operations will not begin until all dependent map tasks have
completed.

Although Spark uses in-memory computing, it is not necessary for all the data
to fit in memory. Spark will spill any RDD elements that cannot be maintained
in memory to disk. Moreover, as Spark transformations generate new RDDs and
numerous transformations may occur within a single application, Spark
implements a Least-Recently Used (LRU) eviction policy. If an evicted RDD
needs to be reused, Spark will recompute it using the lineage data collected.
As demonstrated in [13], caching significantly improves processing times of
iterative algorithms where RDDs are reused. It can be of even greater
importance if the RDD is costly to recompute.

Data locality in Spark is achieved through the scheduling of tasks to
partitions which have the data loaded in memory. If the data is instead stored
on HDFS, the scheduler will assign it to one of the preferred locations
specified by HDFS. Spark’s scheduler utilizes delay scheduling to optimize
fairness and locality for all tasks.

Three different types of schedulers are compatible with Spark: 1) Spark
Standalone, 2) YARN [19] and 3) Mesos [20]. The Spark Standalone scheduler is
the default scheduler. The YARN scheduler designed for Hadoop [21] clusters
and is prepackaged with Hadoop installations, whereas in contrast, Mesos was
designed to be used in multi-tenant cluster environments. In our experiments,
we focus on the Spark Standalone cluster.

Executing Spark applications on HPC systems with Spark-unaware schedulers may
be inefficient. The amount of resources requested by Spark may impede Spark-
cluster scheduling time. Using pilot-scheduling strategies to add nodes to the
Spark cluster as they are allocated by the underlying schedulers may speedup
allocation and overall processing time [22]. This, however, is not studied in
the current paper.

As Spark is frequently used by the scientific community, we designed our
experiments using the PySpark API. This is at a cost to performance as the
PySpark code must undergo Python to Java serialization. We used Spark version
2.3.2 installed from https://spark.apache.org.

####  II-A2 Nipype

Nipype is a popular Python neuroimaging processing engine. It aims at being a
solution for easily creating reproducible neuroimaging workflows. Although
Nipype does not employ any Big Data processing strategies, it provides plugins
to numerous schedulers found in most clusters readily available to
researchers, such as the Sun/Oracle Grid Engine (SGE/OGE), TORQUE, Slurm and
HTCondor. It also includes its own scheduler, MultiProc, for parallel
processing on single nodes. Furthermore, Nipype provides many built-in
interfaces to commonly used neuroimaging tools that can be incorporated within
the workflows. Nipype’s ability to easily parallelize workflows in researcher-
available cluster setups, capture detailed provenance information necessary
for reproducibility, and allow users to easily integrate existing neuroimaging
tools, make it preferable over existing Big Data solutions, which would
necessitate modifications to achieve this.

Jobs, or Interfaces in Nipype, are encapsulated by a Node. A Node dictates
that the job will execute on a single input. However, the MapNode, a child
variant of the Node, can execute on multiple inputs. All tasks in Nipype
execute in their own uniquely named subdirectory which facilitates provenance
tracking of inputs and outputs and also enables checkpointing of the workflow.
In the case of Node failure or application modification, only the nodes which
have been modified (verified by hash) or have not successfully completed, are
re-executed.

In order for Nipype to operate as intended, a filesystem shared by all the
nodes is required. However, it is still possible to save to a non-shared local
filesystem, but it may come at the expense of fault-tolerance, as data located
on failed nodes will be permanently lost. Moreover, the user will need to
ensure that the files are appropriately directed to the nodes that require
them as there is no guarantee of data locality in Nipype.

We used Nipype version 1.1.4 installed through the pip package manager.

###  II-B Data Storage Locations

Data storage location is critical to the performance of Big Data applications
on HPC clusters. Data may reside in the engine memory, on a file system whose
contents reside in virtual memory (for instance tmpfs), on disks local to the
processing node, or on a shared file system. Table I summarizes the Big Data
strategies that can be used depending on the data location. In addition, lazy
evaluation is available in Spark regardless of data location. The remainder of
this Section explains this Table and provides related performance models.

Data Location | In-Memory | Data Locality
---|---|---
| Computing |
In-memory | Yes | Yes
tmpfs | Yes | Yes
Local Disk | Page Caching | Yes
Shared File System | Page Caching | No
TABLE I: Big Data strategies on a shared HPC cluster.

####  II-B1 In-Engine-Memory and In-Memory File System

The main difference between storing the data in the engine memory, as in
Spark, and simply writing to an in-memory file system, such as tmpfs, is what
happens when the processed data fills up the available memory. When engine
memory is used, the engine must cleanup unused data to avoid crashes. When an
in-memory file system is used, the user is responsible ensuring that the
filesystem does not reach capacity. Should it be necessary to free up
additional memory, the kernel will swap filesystem memory to disk and the
performance will become that of local disk writes. In our experiments, we will
explore the configuration where the data consumed by the application
approaches the threshold of available memory.

####  II-B2 Local Disk

Storing data on local disks inevitably enables data locality, since data
transfers are not necessary when tasks are executed on the nodes where the
data resides. However, in absence of a more specific filesystem such as HDFS
to handle file replication across computing nodes, data locality comes at the
price of stringent scheduling restrictions, as tasks can only be scheduled to
the single node that contains their input data.

The performance of local disk accesses is strongly dependent on the page
caching mechanism provided by the Linux kernel, described in details in [9].
To summarize, data read from disk remains cached in memory until evicted by an
LRU (Least Recently Used) strategy. When a process invokes the read() system
call, the kernel will return the data directly from memory if the requested
data lies in the page cache, realizing a cache _hit_. Cache hits drastically
speed-up data reads, by masking the disk latency and bandwidth behind a memory
buffer. In effect, page caching provides in-memory computing transparently to
the processing engine. However, page cache eviction strategies currently
cannot be controlled by the application, which prevents processing engines
from anticipating reads by preloading the cache. Scheduling strategies might
be designed to maximize cache hits, however. For instance, lazy evaluation
could result in more cache hits by scheduling data-dependent tasks on the same
node.

Page caching has a more dramatic effect on disk writes, reducing their
duration by several orders of magnitude. When a process calls the write()
system call, data is copied to a memory cache that is asynchronously written
to disk by flusher threads, when memory shrinks, when “dirty” (unwritten) data
grows, or when a process invokes the sync() system call. This asynchronous
flushing of the page cache is called _writeback_.

Page caching is essentially a way to emulate in-memory computing at the kernel
level, without requiring a dedicated engine. The size of the page cache,
however, becomes a limitation when processes write faster than the disk
bandwidth. When this happens, the page cache rapidly fills up and writes are
limited by the disk write bandwidth as if no page cache was involved.

We introduce the following basic model to describe the filling and flushing of
the page cache by an application:

$d(t)=\left(\frac{D}{C}-\frac{\delta}{\gamma}\right)t+d_{0},$

where:

  * •

$d(t)$ is the amount of data in the page cache at time t

  * •

$D$ is the total amount of data written by the application

  * •

$C$ is the total CPU time of the application

  * •

$\delta$ is the disk bandwidth

  * •

$\gamma$ is the max number of concurrent processes on a node

  * •

$d_{0}$ is the amount of data in the page cache at time $t_{0}$

This model applies to parallel applications assuming that (1) concurrent
processes all write the same amount of data, (2) concurrent processes all
consume the same CPU time, (3) data is written uniformly along task execution.
With these assumptions, all the processes will write at the same rate, which
explains why the model does not depend on the total number of concurrent
processes in the application, but only on the max number of concurrent
processes executing on the same node ($\gamma$). While these assumptions would
usually be violated in practice, this simple model already provides
interesting insights on the performance of disk writes, as shown later.
Naturally, the model also ignores other processes that might be writing to
disk concurrently to the application, which we assume negligible here.

In general, an application should ensure that $\dot{d}$ remains negative or
null, leading to the following inequality:

$\frac{D}{C}\leq\frac{\delta}{\gamma}$ (1)

This defines a D/C (data-write) rate beyond which the page cache becomes
asymptotically useless. It should be noted that the transient phase during
which the page cache fills up might last a significant amount of time, in
particular when $\dot{d}$ is positive and small. We intentionally do not model
the transient phase as it requires detailed knowledge of difficult to estimate
parameters such as the page cache size and the initial amount of data in it
($d_{0}$).

We will use Equation 1 to define our benchmarks and interpret the results. It
should be noted that leveraging the page cache, and therefore ensuring that
Equation 1 holds, has important performance implications: with page caching,
the write throughput will be that of memory, while without page caching it
will be that of the disk.

####  II-B3 Shared File System

We model a shared file system using its global realized bandwidth $\Delta$,
shared by all concurrent processes in the cluster. We are aware that such a
simplistic model does not describe at all the intricacies of systems such as
Lustre. In particular, metadata management, RPC protocol optimizations and
storage optimizations are all covered under the realized bandwidth. We do,
however, consider the effect of page caching in shared file systems, since in
Linux writes to network-mounted volumes benefit from this feature too.

As in the local disk model, we note that page caching will only be useful when
the flush bandwidth is greater than the write throughput of the application,
that is:

$\frac{D}{C}\leq\frac{\Delta}{\Gamma},$ (2)

where $\Gamma$ is the max number of concurrent processes _in the cluster_.
Note that $\frac{\Delta}{\Gamma}$ will usually be much lower than
$\frac{\delta}{\gamma}$.

###  II-C Infrastructure

All experiments were executed on Dell EMC’s Zenith cluster, a Top-500 machine
in the Dell EMC HPC and AI Innovation Lab, running Slurm. For the Spark
experiments, a Spark cluster was started on a Slurm allocation comprised of 16
dedicated nodes. Each Compute node has Red Hat Enterprise Linux Server release
7.4 (Maipo) as the base operating system with kernel version
3.10.0-693.17.1.el7.x86_64 (patched for Spectre/Meltdown). Dell EMC PowerEdge
C6420 with dual Intel Xeon Gold 6148/F processors (40 cores per node) and
192GB ($12\times 16$ GB), 2666 MHz memory, serve as the compute nodes. Each
compute has a 120GB M.2 SATA SSD as local disk. A Dell HPC Lustre Solution
with a raw storage of 960TB is accesible on each compute node through a 100
Gb/s Intel OmniPath network. All the nodes connect to a director switch in a
1:1 non-blocking topology. The realized write bandwidth of the local disk,
Lustre file system and tmpfs were measured by sequentially writing various
numbers of image blocks containing random intensities, to avoid caching
effects (see: measure_bandwidth.py). They are reported in Table II.

Data location | Measured write bandwidths (MB/s)
---|---
tmpfs | 1377.18
Local disk ($\delta$) | 193.64
Lustre ($\Delta$) | 504.03
TABLE II: Measured bandwidths

###  II-D Datasets

We used BigBrain [23], a 75GB 40$\mu$m isotropic histological image of a
65-year-old human brain. The BigBrain was selected due to its uniqueness, as
there does not yet exist a higher-resolution image of a human brain. Moreover,
there currently exists a lack of standardized tools for processing the
BigBrain as a consequence of its size. To examine the effects processing the
BigBrain has on Big Data strategies, we partitioned the full $3845\times
3015\times 3470$ voxel image into 30 ($5\times 3\times 2$) chunks, 125
($5\times 5\times 5$) chunks and 750 ($5\times 15\times 10$) chunks.
Additionally, the full image was also split in half ($769\times 603\times 347$
voxels) and the half image was partitioned into 125 chunks.

Processing large images is only considered to be part of the Big Data problem
in neuroscience. The other problem being the processing large MRI datasets,
that is, datasets consisting of many small brain images belonging to various
different subjects. This situation is commonly observed in functional MRI
(fMRI), where it is becoming increasingly common to process data from hundreds
of subjects. Although we have not explored explicitly the processing of large
MRI datasets, the 75GB BigBrain is within the size ballpark [10] of MRI
datasets commonly processed in today’s studies.

Since both small and large datasets may need to be processed using the same
analysis pipeline, we examined the effects of the data management strategies
on small data as well. For this, we selected a 12MB T1W image belonging to
subject 1 of OpenNeuro’s ds000001 dataset version 6. In order to split the
image into 125 equal-sized chunks, it was necessary to zero-pad the image to
the dimensions $165\times 200\times 200$ voxels, which subsequently increased
the total image size to 13MB.

###  II-E Applications

Algorithm 1 Incrementation

1:Input

2: $x$a sleep delay in seconds

3: $n$a number of iterations

4: $C$a set of image chunks

5: $fs$filesystem to write to (mem, tmpfs, local disk, Lustre).

6:for each $chunk\in C$ do

7: read $chunk$ from Lustre

8: for $i\in[1,n]$ do

9: $chunk\leftarrow chunk+1$

10: sleep $x$

11: if $i<n$ then

12: save $chunk$ to $fs$

13: end if

14: end for

15: save $chunk$ to Lustre

16:end for

To effectively investigate how the different strategies impact processing, we
selected a simple incrementation pipeline (Algorithm 1) that consists
exclusively of map stages. A series of map-only stages would enable us to
evaluate the effects of in-memory computing when data locality is preserved.
Incrementation was selected over other applications, such as binarization, as
it ensured that a new image was created at each step (i.e. no caching effects
within the executing application). Each partitioned chunk was incremented by
1, in parallel, by each task. As incrementing images is not a time consuming
process, we added a sleep delay to the tasks to study the effects of tasks
duration. The incremented chunks would be either maintained in-memory (Spark
only) or saved to either tmpfs, local disk or Lustre (Spark and Nipype).
Should more than a single iteration be requested, the incremented chunks would
be incremented again and saved to the same file system. This would repeat
until the number of requested iterations had elapsed. In all conditions, the
first input chunks and final output chunks would reside on Lustre. We chose to
perform our initial input/final output on Lustre as local storage is typically
only accessible to a user for the duration of the execution in HPC
environments.

###  II-F Experiments

We conducted four experiments in which we varied (1) the number of iterations
in the application ($n$ in Algorithm 1), (2) the task duration ($x$), (3) the
chunk size given a constant image size, (4) the total image size. To evaluate
the page-cache model, experiment conditions fell in different regions of
Equations 1 and 2, as summarized in Table III. Among the 16 nodes available, 1
was dedicated to the Spark master and driver, and the remaining 15 were used
as compute nodes. Since data locality is not normally preserved in Nipype (a
new Slurm allocation is requested for each processed task), we instrumented
Nipype to ensure data locality (see: run_benchmarks.py). That is, the chunks
were split into partitions, and for each partition, we requested a Slurm
allocation to process the entire pipeline in parallel using Nipype’s MultiProc
scheduler on a given node. This was possible as no communication was required
between the processed chunks.

Experiment 1: Number of Iterations
---
n | D (GB) | C (s) | D/C (MB/s) | $\gamma$ |  $\delta/\gamma$ (MB/s) | $\Gamma$ |  $\Delta/\Gamma$ (MB/s) | (D/C)/($\delta/\gamma$) | (D/C)/($\Delta/\Gamma$)
1 | 75 | 430 | 178.6 | 9 | 21.5 | 125 | 4.0 | 8.3 | 44.7
10 | 750 | 4,300 | 178.6 | 9 | 21.5 | 125 | 4.0 | 8.3 | 44.7
100 | 7,500 | 43,000 | 178.6 | 9 | 21.5 | 125 | 4.0 | 8.3 | 44.7
Experiment 2: Task Duration
x (s) | D (GB) | C (s) | D/C (MB/s) | $\gamma$ |  $\delta/\gamma$ (MB/s) | $\Gamma$ |  $\Delta/\Gamma$ (MB/s) | (D/C)/($\delta/\gamma$) | (D/C)/($\Delta/\Gamma$)
2.4 | 750 | 3,000 | 256 | 9 | 21.5 | 125 | 4.0 | 11.9 | 64
3.44 | 750 | 4,300 | 178.6 | 9 | 21.5 | 125 | 4.0 | 8.3 | 44.7
7.68 | 750 | 9,600 | 80 | 9 | 21.5 | 125 | 4.0 | 3.7 | 20
320 | 750 | 400,000 | 1.9 | 9 | 21.5 | 125 | 4.0 | 0.09 | 0.48
Experiment 3: Number of Chunks
chunks | D (GB) | C (s) | D/C (MB/s) | $\gamma$ |  $\delta/\gamma$ (MB/s) | $\Gamma$ |  $\Delta/\Gamma$ (MB/s) | (D/C)/($\delta/\gamma$) | (D/C)/($\Delta/\Gamma$)
30 | 750 | 4,400 | 174.6 | 2 | 96.8 | 30 | 16.8 | 1.8 | 10.4
125 | 750 | 4,400 | 174.6 | 9 | 21.5 | 125 | 4.0 | 8.1 | 43.7
750 | 750 | 4,400 | 174.6 | 25 | 7.7 | 375 | 1.3 | 22.7 | 134.3
Experiment 4: Image Size
image | D (GB) | C (s) | D/C (MB/s) | $\gamma$ |  $\delta/\gamma$ (MB/s) | $\Gamma$ |  $\Delta/\Gamma$ (MB/s) | (D/C)/($\delta/\gamma$) | (D/C)/($\Delta/\Gamma$)
BigBrain | 750 | 2,200 | 349.1 | 9 | 21.5 | 125 | 4.0 | 16.2 | 87.3
Half BigBrain | 375 | 2,200 | 174.6 | 9 | 21.5 | 125 | 4.0 | 8.1 | 43.7
MRI | 0.127 | 2,200 | 0.06 | 9 | 21.5 | 125 | 4.0 | 0.003 | 0.015
TABLE III: Experiment conditions. Red cells denote the conditions where the
inequalities in Equations 1 and 2 do not hold, i.e., the page cache is
asymptotically useless. Green cells show the conditions where the page cache
covers all data writes.

For our first incrementation experiment, we investigated the effects of Big
Data strategies on varying total data size. To achieve this, we increased the
number of incrementation iterations from 1, 10 and 100 times. The total data
size would then increase from 75GB, at 1 iteration, to 7,500GB, at 100
iterations. The total number of chunks was 125. Chunks all ran concurrently
($\Gamma$=125) and were equally balanced among 15 nodes, leading to 8 or 9
concurrent jobs per node ($\gamma$=9). Task duration was fixed at 3.44
seconds.

In the second experiment, we evaluated the effects of Big Data strategies on
varying task duration. If the page cache has sufficient time to flush, it
would be expected that in-memory computing and local disk perform
equivalently. We varied the task duration between 2.4 and 320 seconds such
that the D/C falls into different regions of Equations 1 and 2. The number of
chunks was maintained at 125, leading to $\Gamma$=125 and $\gamma$=9. The
number of iterations was fixed to 10.

As a third incrementation experiment, we were interested in the effects of
chunk size on Big Data strategies. Naturally, a greater chunk size signifies a
decrease in parallelism. However, it also signifies an increase in sequential
I/O (increased $\Delta/\Gamma$ and $\delta/\gamma$). For this experiment we
partitioned the complete BigBrain image into 30, 125 and 750 chunks,
corresponding to $\gamma$ values of 2, 9 and 25 respectively. While Spark
attempted to load-balance the data, it used up only 25 of the 40 cores for 750
chunks. In contrast, Nipype tried to use up as many cores as possible. Unlike
the previous experiment, the D/C rate was kept static at 178.6MB/s, however,
this ratio ensured that different regions of the inequality were reached
depending on amount of parallelism. The number of iterations was fixed to 10,
and the task duration was adjusted so that C remained constant at 4,400s.

For our fourth and final incrementation experiment, we investigated the
effects of the strategies on different image sizes. We selected the 75GB
BigBrain, the 38GB half BigBrain and the 13M T1W MRI image for this
experiment. The number of chunks was fixed to 125. Similarly to the previous
experiment, the total sequential compute time was fixed (10 iterations, 1.76
seconds per task), however, due to varying size in total data processed (D),
the D/C rate varied. Once again, we ensured that the D/C rate fell in multiple
different regions of the inequality. The D/C rate ranged from 349.1 MB/s for
BigBrain and 174.6MB/s for half BigBrain, to 0.06MB/s for the MRI image. Only
the 0.06MB/s MRI satisfied the inequality for both Lustre and local disk.

##  III Results

(a) Experiment 1: complete BigBrain, 125 chunks, 3.44-second tasks.

(b) Experiment 2: complete BigBrain, 125 chunks, 10 iterations.

(c) Experiment 3: complete BigBrain, 10 iterations, C= 4,400s.

(d) Experiment 4: 125 chunks, 10 iterations, 1.76-second tasks.

Figure 1: Experiment results: Makespans of Spark and Nipype writing to memory,
tmpfs, local disk and Lustre.

###  III-A Experiment 1: Number of Iterations

Fig. 1(a) shows the difference between the different filesystem choices given
varying number of iterations. At 1 iteration, all filesystems behave the same,
although the application was writing faster than the disk bandwidth. This is
because application data was not saturating the page cache (transient phase).
The page cache, on Zenith, occupies 40% of total memory. With 192 GB of RAM on
each node, 76.8 GB of dirty data could be held in a node’s page cache at any
given time. As the total amount of data written by the application increases
to 750 GB, there is a greater disparity between Lustre and in-memory (2.67 x
slower, on average). Local disk performance, however, is still comparable to
memory (1.38 x slower, on average). Despite local disk and Lustre both being
in transient state, local disk encounters less contention than what would be
found on Lustre.

At 100 iterations, or 7,500 GB, Lustre can be found to be, on average, 3.82 x
slower than Spark in-memory. The slowdown experienced can be explained by the
smaller percentage of total data residing in the page cache at a given time,
compared to 10 iterations. Therefore, the effects of Lustre bandwidth are more
significant in this application. At 100 iterations, the application was
writing 500GB per node (7,500 GB / 15 nodes) and hence could not run on tmpfs
or local disk.

While there is some variability that can be seen in Fig. 1(a) between the two
engines, this believed to be insignificant, and potentially due to SLURM node
allocation delays in our launching of Nipype.

###  III-B Experiment 2: Task Duration

Increasing task duration ensured that all file systems had a comparable
performance (Fig. 1(b)). Lustre, for instance, is approximately 1.01x slower
than Spark in-memory at a task duration of 320 seconds, whereas it is
approximately 3.25x slower that Spark in-memory with 2.4 second tasks. This
pattern corroborated our page-cache model which postulates that data movement
costs will have little impact on compute-intensive tasks. The reasoning behind
this is that longer tasks give the page cache more time to flush between disk
writes.

###  III-C Experiment 3: Image Block Size

As can be seen in Fig. 1(c), makespan decreases with increasing number of
chunks. This is due to the fact that parallelism increases with an increase in
number of chunks. At 30 chunks, only 2 CPUs per node are actively working. At
125 chunks, this changes to a maximum of 9 CPUs per node, and at 750 chunks,
up to 40 CPUs can be active.

Due to a size limitation of 2 GB imposed on Spark partitions, Spark with in-
memory computing processing 30 chunks was not performed.

Local disk and tmpfs perform comparably for all conditions, with Lustre being
significantly slower. As with varying the number of iterations, Lustre is
slower due to increased filesystem contention, which is, at minimum, 15 x
greater than contention on local disk, due to the number of nodes used. With
an increase in number of chunks, local disk and tmpfs makespans begin to
converge. A potential explanation for this may be that tmpfs is utilizing swap
space. As concurrency increases, the memory footprint of the application also
increases. It is possible that at 750 chunks, swapping to disk is required by
tmpfs, thus resulting in similar processing times as local disk.

Swapping may also be an explanation for the variance between Spark in-memory
and tmpfs performance. While Spark may also spill to disk, it only does so
when data does not fit in memory. As none of the RDDs generated throughout the
pipeline were cached and all data concurrently accessed could be mantained in-
memory, spilling to disk was not necessary.

###  III-D Experiment 4: Image Size

Increasing overall data size decreases performance, as can be seen in Fig.
1(d). When the data size is very small (e.g. MRI image) all file system
makespans are comparable. This is due to the fact that page cache can be
leveraged fully regardless of file system. However, this time, Spark in-memory
performed significantly worse than all other filesystems, with a makespan of
2,211 seconds. Upon further inspection, it appeared that Spark in-memory
executed in a sequential order, on a single worker node. Lack of parallelism
for the MRI image may be a result of Spark’s maximum partition size, which is
by default 128 MB – significantly larger than the 13 MB MRI image.

At half BigBrain, the makespan differences become apparent in both local disk
and Lustre, with Lustre becoming 2.4x slower than in-memory. This can be
attributed to page cache saturation, as predicted by the model for both half
the BigBrain image and the complete BigBrain. Only the MRI image was predicted
to fall within the model constraints.

When the complete BigBrain is processed, the disparity between the different
filesystems becomes even greater. Lustre becomes 3.68x slower, whereas local
disk becomes 1.68x slower. An explanation for this is that the page cache
fills up faster due to data size.

###  III-E Page Cache Model Evaluations

(a) Local Disk

(b) Lustre

Figure 2: Page cache model evaluation. Grey regions denote areas that violate
model predictions.

In order to evaluate the page cache model, we compared the observed speedup
ratio provided by in-memory computing to the (D/C) / ($\delta$/$\gamma$) and
(D/C) / ($\Delta$/$\Gamma$) ratios (Fig. 2). Speed-up ratios were computed as
the ratio between the makespan obtained with Spark on local disk or Lustre,
and the makespan obtained with Spark for in-memory computing. Experiments for
which there were no in-memory equivalent (i.e. BigBrain split into 30 chunks)
were not considered.

Results show that, overall, the model correctly predicted the effect of page
cache on processing times for local disk and Lustre. That is, the speed-up
provided by in-memory computing was larger than 1 for D/C rates larger than
$\delta/\gamma$ (local disk) or $\Delta/\Gamma$ (Lustre). Conversely, the
speed-up provided by in-memory computing remained close to 1 for D/C rates
smaller than $\delta/\gamma$ (local disk) or $\Delta/\Gamma$ (Lustre). The two
points close to the origin correspond to the sequential processing of the MRI
image by Spark mentioned previously.

Points which violated model predictions were found at 1 iteration, where page
cache would not have been saturated in spite of a high D/C (transient state).
However, in all cases, the “1” boundary was never trespassed by more than a
factor of 0.19, and is therefore likely a result of system variability.

##  IV Discussion

###  IV-A Effect of In-Memory Computing

We measured the effect of in-memory computing by comparing the runs of Spark
in-memory (yellow bars in Fig. 1) to the ones of Spark on local disk (non-
hatched green bars). The speed-up provided by in-memory computing is also
reported in Fig. 2(a). The speed-up provided by in-memory computing increases
with (D/C) / ($\delta$/$\gamma$), as expected from the model. In our
experiments, it peaked at 1.6, for a ratio of 16.2. This correspond to the
processing of the BigBrain with 125 chunks and 1.76-second tasks in experiment
4 (total computing time C=2,200s), which is typically encountered in common
image processing tasks such as denoising, intensity normalization, etc. The
speed up of 1.6 is also reached with a ratio of 22.7 in experiment 3, obtained
by processing the BigBrain with 750 chunks.

The results also allow us to speculate on the effect of in-memory computing on
the pre-processing of functional MRI, another typical use case in
neuroimaging. Assuming an average processing time of 20 minutes per subject,
which is a ballpark value commonly observed with the popular SPM or FSL
packages, an input data size of 100MB per subject, and an output data size of
2GB (20-fold increase compared to input size), the D/C rate would be 1.8MB/s,
which would reach the $\delta/\gamma$ threshold measured on this cluster for
$\gamma=108$, that is, if 108 subjects were processed on the same node. This
is very unlikely as the number of CPUs per node was 40. We therefore conclude
that in-memory computing is likely to be useless for fMRI analysis. Naturally,
this estimate is strongly dependent on the characteristics of the cluster.

###  IV-B Effect of Data Locality

We measure the effect of data locality by comparing the runs of Spark on local
disk (non-hatched green bars in Fig. 1) to the ones of Spark on Lustre (non-
hatched blue bars). The speed-up provided by local execution peaked at 3.2,
for 750 chunks in experiment 3. Overall, writing locally was usually
preferable over writing to Lustre, as a result of the lower contention on
local disk. Although it may be true that network bandwidths exceed that of
disks [24], locality remains important as contention on a shared filesystem
tends to be much higher than on local disk. The only time writing locally did
not have significant impact over Lustre was in experiments 1 and 4, at 1
iteration and when processing the MRI image, respectively. In both these
scenarios, the Lustre writes did not impact performance as the data was able
to be written to page cache and flushed to Lustre asynchronously.

###  IV-C Combined Effect of In-Memory and Data Locality

We measure the combined effect of data locality and in-memory computing by
comparing the runs of Spark in-memory (yellow bars in Fig. 1) to the ones of
Spark on Lustre (non-hatched blue bars). The speed-up provided by the combined
use of data locality and in-memory computing is also reported in Fig. 2(b).
The provided speed-up increases with (D/C) / ($\Delta$/$\Gamma$), as expected
from the model. In our experiments, it peaked around 5, for ratios of 120.4
and 64. Again, this configuration is likely to happen in typical image
processing tasks performed on the BigBrain.

As for the fMRI speculation, the D/C rate of 1.8MB/s would reach the
$\Delta/\Gamma$ threshold for $\Gamma=280$, which is a realistic number of
subjects to process on a complete cluster. Naturally, this estimate is highly
dependent on the observed bandwidth of the shared file system ($\Delta$).

###  IV-D Effect of Lazy Evaluation

The effects of lazy evaluation can be seen throughout the experiments. Nipype
was found to be slower than Spark in most experiments. While the Nipype
execution graph is generated prior to workflow execution, there are no
optimizations to ensure that the least amount of work is performed to produce
the required results.

During the processing of Experiment 3, 750 chunks were processed in two
batches for both Spark and Nipype due to CPU limitations. Rather than running
each iteration on the full dataset, as with Nipype, Spark opted to perform all
the iterations on the first batch (load, increment, save), and then proceeded
to process the second batch. Such an optimization is important, even when
processing data on disk, as it would presumably increase the occurrence of
cache hits. This may partially explain the speedup seen at 750 chunks in
Figure 1(c).

###  IV-E Can tmpfs and Page Caches Emulate In-Memory Computing?

Although tmpfs and page cache do improve performance, as seen in Figure 1,
they do not always perform equivalently to in-memory. Tmpfs’s main limitation
is that data residing on it may be swapped to disk if the system’s memory
usage is high. When it reaches this point, its performance slows down to swap
disk bandwidth, as observed in Figure 1(d). Page cache suffers a similar
dilemma. I/O blocking writes to disk occur when a given percentage (e.g. 40 %)
of total memory is occupied by dirty data. When the threshold is exceeded,
processes performing writes must wait for dirty data to be flushed to disk.

Furthermore, like memory, tmpfs and page cache are shared resources on a node.
If users on the node are heavily using memory to incite tmpfs to writes to
swap space, or are performing data-intensive operations that fill up the page
cache, tmpfs/page cache performance on will be limited for other users.
However, it is possible to request through the HPC scheduler a certain amount
of available memory. Ultimately, in-memory data will also need to be spilled
to disk if memory usage exceeds amount of available memory, although disk
writes are likely to occur in tmpfs and page cache before requested available
memory is filled.

###  IV-F Scheduling Remarks

A common recommendation in Spark is to limit the number of cores per executors
to 5, to preserve a good I/O throughput (see Cloudera blog), but the rationale
for this recommendation is hardly explained. We believe that throughput
degradation observed with more than 5 cores per executor might be coming from
full page caches.

Spark does not currently include any active management of disk I/Os or page
caches. We believe that it would be beneficial to extend it toward this
direction, to increase the performance of operations where local storage has
to be used, such as disk spills or shuffles. For instance, workflow-aware
cache eviction policies that maximizes page cache usage for the workflows
could be investigated.

An alternative Nipype plugin designed for running on Slurm was not used in the
experiments. The Slurm plugin requests a Slurm allocation for each processed
data chunk. Such a scheduling strategy was not ideal in our environment where
oversubscription of nodes was not enabled.

Unlike Spark, Nipype by default opts to use all available CPUs rather than to
load balance data across the cluster. That is, given 50 chunks and 40 cores,
Spark will only use up 25 cores and process in two batches. Nipype will also
have no choice but to split up the processing into two batches, but will first
process 40 chunks, immediately followed by the remaining 10. While both are
reasonable strategies for data processing, Spark may end up benefiting more
from the page cache, as less data is written in parallel (25 vs 40), giving
more time for the page cache to flush.

Nipype’s MapNodes, which apply a given function to each element in a list,
were found to be slower than the Node, which apply a function to a single
element, due to a blocking mechanism. For this reason, we selected to iterate
through a series of Nodes in our code despite MapNodes being easier to use.

###  IV-G Other Comments

Writing to node-local storage in a cluster environment comes at a cost, for
both Nipype and Spark without HDFS. When a node is lost, the node-local data
is lost with it. While Spark will recompute the lost partition automatically
using lineage, Nipype will fail for all tasks requiring the data.
Nevertheless, Spark will also fail if RDDs of node-local filenames are
shuffled, as the data associated to the filenames will not be shuffled with
the RDD and there will be no mechanism in place to fetch it.

When executed on Lustre, Nipype will checkpoint itself, ensuring resumption
from last checkpoint during re-execution. This is particularly important in
the case of compute-intensive applications, such as those found in
neuroimaging. Spark also provides a checkpointing mechanism, however, it
requires HDFS.

It is common in neuroimaging applications for users to want access to all
intermediate data. Such a feature is currently only possible when writing to
shared filesystem. It would also not be an option with Spark in-memory. To
enable this, burst buffers or heterogeneous storage managers (e.g. Triple-H
[25]) could be used to ensure fast processing and that all outputs (including
intermediate outputs) will be sent asynchronously to the shared filesystem.

It was expected that Spark would experience longer processing times,
particularly with the small datasets, due to Java serialization. This was not
found to be the case. Unlike Spark, Nipype performs a more thorough provenance
capture, potentially owing to longer processing times.

In this paper, we analyzed the effects of Big Data strategies on an map-only
artificial neuroimaging pipeline. This allowed us to examine the effects of
these strategies without being significantly obscured by other conditions.
Studying the effects of such strategies on map-reduce type workflows, in
addition to real neuroimaging pipelines, would allows us to gain further
insight on the added value of the Big Data performance stragies for
neuroimaging use cases remains to be done.

##  V Conclusion

Big Data performance optimization strategies help improve performance of
typical neuroimaging applications. Our experiments indicate that overall, in-
memory computing enables greater speedups than what can be obtained by using
page cache and tmpfs. While page cache and tmpfs do give memory-like
performance, they are likely to fill up faster than memory, leading to
increased performance penalties when compared to in-memory computing. We
conclude that extending Big Data processing engines to better support
neuroimaging applications, including developing their provenance, fault-
tolerance, and reproducibility features, is worthwhile.

Data locality plays an important role in application performance. Local disk
was found to perform better than the shared filesystem despite having lower
bandwidth, due to increased contention on the shared filesystem. Since local
disk typically has less contention than shared filesystems, it is recommended
to store data locally. However, using local storage without a distributed file
system may limit fault tolerance.

Although a more thorough analysis of lazy evaluation remains to be performed,
it is speculated that this may be the cause of the general performance
difference between Spark and Nipype. Furthermore, it was found that lazy
evaluation optimizations increase the likelihood of cache hits, thus improving
overall performance.

Even though Big Data strategies are beneficial to the processing of large
images, it is estimated that it would require running a functional MRI dataset
with 280 concurrent subjects for any noticeable impact using our Lustre
bandwidth estimate. Benchmarking Spark and Nipype using such a large fMRI
dataset would be a relevant follow-up experiment, to test this hypothesis. It
would also be useful to evaluate other types of applications, such as ones
containing data shuffling steps.

Finally, we plan to extend this study by including task scheduling strategies
in a multi-tenant environment. We expect to observe important differences
between Spark and Nipype, due to Spark’s use of overlay scheduling. The impact
of other Big Data technologies, such as distributed in-memory file
systems(e.g. Apache Ignite) and Lustre scalability issues [26], could also be
investigated.

##  VI Acknowledgments

We are thankful to the Dell EMC HPC and AI Innovation Lab in Roundrock, TX,
for providing the infrastructure and high-quality technical support.

## References

  * [1] J. Dean and S. Ghemawat, “MapReduce: simplified data processing on large clusters,” _Comm. of the ACM_ , vol. 51, no. 1, pp. 107–113, 2008.
  * [2] M. Zaharia, R. S. Xin, P. Wendell, T. Das, M. Armbrust, A. Dave, X. Meng, J. Rosen, S. Venkataraman, M. J. Franklin, _et al._ , “Apache Spark: a unified engine for big data processing,” _Comm. of the ACM_ , vol. 59, no. 11, pp. 56–65, 2016.
  * [3] M. Rocklin, “Dask: Parallel computation with blocked algorithms and task scheduling,” in _Proc. of the 14th Python in Science Conference_ , no. 130-136. Citeseer, 2015.
  * [4] E. Deelman, G. Singh, M.-H. Su, J. Blythe, Y. Gil, C. Kesselman, G. Mehta, K. Vahi, G. B. Berriman, J. Good, _et al._ , “Pegasus: A framework for mapping complex scientific workflows onto distributed systems,” _Scientific Programming_ , vol. 13, no. 3, pp. 219–237, 2005.
  * [5] T. Oinn, M. Addis, J. Ferris, D. Marvin, M. Senger, M. Greenwood, T. Carver, K. Glover, M. R. Pocock, A. Wipat, _et al._ , “Taverna: a tool for the composition and enactment of bioinformatics workflows,” _Bioinformatics_ , vol. 20, no. 17, pp. 3045–3054, 2004.
  * [6] J. Leskovec, A. Rajaraman, and J. D. Ullman, _Mining of massive datasets_. Cambridge university press, 2014\.
  * [7] K. Shvachko, H. Kuang, S. Radia, and R. Chansler, “The Hadoop distributed file system,” in _Mass storage systems and technologies (MSST), 2010 IEEE 26th symposium on_. Ieee, 2010, pp. 1–10.
  * [8] P. Schwan _et al._ , “Lustre: Building a file system for 1000-node clusters,” in _Proc. of the 2003 Linux symposium_ , 2003, pp. 380–386.
  * [9] R. Love, _Linux kernel development_. Pearson Education, 2010.
  * [10] J. D. Van Horn and A. W. Toga, “Human neuroimaging as a Big Data science,” _Brain imaging and behavior_ , vol. 8, no. 2, pp. 323–331, 2014\.
  * [11] K. Gorgolewski, C. D. Burns, C. Madison, D. Clark, Y. O. Halchenko, M. L. Waskom, and S. S. Ghosh, “Nipype: a flexible, lightweight and extensible neuroimaging data processing framework in Python,” _Frontiers in neuroinformatics_ , vol. 5, p. 13, 2011.
  * [12] P. Bellec, S. Lavoie-Courchesne, P. Dickinson, J. Lerch, A. Zijdenbos, and A. C. Evans, “The pipeline system for Octave and Matlab (PSOM): a lightweight scripting framework and execution engine for scientific workflows,” _Frontiers in neuroinformatics_ , vol. 6, p. 7, 2012.
  * [13] J. Freeman, N. Vladimirov, T. Kawashima, Y. Mu, N. J. Sofroniew, D. V. Bennett, J. Rosen, C.-T. Yang, L. L. Looger, and M. B. Ahrens, “Mapping brain activity at scale with cluster computing,” _Nature methods_ , vol. 11, no. 9, p. 941, 2014.
  * [14] M. Makkie, H. Huang, Y. Zhao, A. V. Vasilakos, and T. Liu, “Fast and scalable distributed deep convolutional autoencoder for fMRI Big Data analytics,” _Neurocomputing_ , vol. 325, pp. 20–30, 2019.
  * [15] P. Mehta, S. Dorkenwald, D. Zhao, T. Kaftan, A. Cheung, M. Balazinska, A. Rokem, A. Connolly, J. Vanderplas, and Y. AlSayyad, “Comparative evaluation of Big-Data systems on scientific image analytics workloads,” _Proc. of the VLDB Endowment_ , vol. 10, no. 11, pp. 1226–1237, 2017.
  * [16] T. Guedes, V. Silva, M. Mattoso, M. V. N. Bedo, and D. de Oliveira, “A practical roadmap for provenance capture and data analysis in spark-based scientific workflows,” in _Workshop on Workflows in Support of Large-Scale Science (WORKS)_ , Dallas, TX, 2018.
  * [17] S. Krishnan, M. Tatineni, and C. Baru, “myHadoop-Hadoop-on-Demand on traditional HPC resources,” _San Diego Supercomputer Center Technical Report TR-2011-2, University of California, San Diego_ , 2011.
  * [18] M. Zaharia, M. Chowdhury, M. J. Franklin, S. Shenker, and I. Stoica, “Spark: Cluster computing with working sets.” _HotCloud_ , vol. 10, no. 10-10, p. 95, 2010.
  * [19] V. K. Vavilapalli, A. C. Murthy, C. Douglas, S. Agarwal, M. Konar, R. Evans, T. Graves, J. Lowe, H. Shah, S. Seth, _et al._ , “Apache Hadoop YARN: Yet another resource negotiator,” in _Proc. of the 4th annual Symposium on Cloud Computing_. ACM, 2013, p. 5.
  * [20] B. Hindman, A. Konwinski, M. Zaharia, A. Ghodsi, A. D. Joseph, R. H. Katz, S. Shenker, and I. Stoica, “Mesos: A platform for fine-grained resource sharing in the data center.” in _NSDI_ , vol. 11, no. 2011, 2011, pp. 22–22.
  * [21] T. White, _Hadoop: The definitive guide_. O’Reilly Media, Inc., 2012.
  * [22] I. Paraskevakos, A. Luckow, G. Chantzialexiou, M. Khoshlessan, O. Beckstein, G. C. Fox, and S. Jha, “Task-parallel analysis of molecular dynamics trajectories,” _CoRR_ , vol. abs/1801.07630, 2018.
  * [23] K. Amunts, C. Lepage, L. Borgeat, H. Mohlberg, T. Dickscheid, M.-É. Rousseau, S. Bludau, P.-L. Bazin, L. B. Lewis, A.-M. Oros-Peusquens, _et al._ , “BigBrain: an ultrahigh-resolution 3D human brain model,” _Science_ , vol. 340, no. 6139, pp. 1472–1475, 2013.
  * [24] G. Ananthanarayanan, A. Ghodsi, S. Shenker, and I. Stoica, “Disk-locality in datacenter computing considered irrelevant.” in _HotOS_ , vol. 13, 2011, pp. 12–12.
  * [25] N. S. Islam, X. Lu, M. Wasi-ur Rahman, D. Shankar, and D. K. Panda, “Triple-h: A hybrid approach to accelerate hdfs on hpc clusters with heterogeneous storage architecture,” in _Cluster, Cloud and Grid Computing (CCGrid), 2015 15th IEEE/ACM International Symposium on_. IEEE, 2015, pp. 101–110.
  * [26] N. Chaimov, A. Malony, S. Canon, C. Iancu, K. Z. Ibrahim, and J. Srinivasan, “Scaling spark on hpc systems,” in _Proc. of the 25th ACM International Symposium on High-Performance Parallel and Distributed Computing_. ACM, 2016, pp. 97–110.


# Analytical models of supermassive black holes in galaxies surrounded by dark
matter halos

Zibo Shena,b<EMAIL_ADDRESS>Anzhong Wanga 111Corresponding author
<EMAIL_ADDRESS>Corresponding author Yungui Gongc
<EMAIL_ADDRESS>Shaoyu Yind<EMAIL_ADDRESS>a GCAP-CASPER, Physics
Department, Baylor University, Waco, Texas 76798-7316, USA
b Department of Astronomy, School of Astronomy and Space Science, University
of Science and Technology of China, Hefei, Anhui 230026, China
cDepartment of Physics, School of Physical Science and Technology, Ningbo
University, Ningbo, Zhejiang 315211, China
d Institute for Theoretical Physics and Cosmology, Zhejiang University of
Technology, Hangzhou 310023, China.

###### Abstract

In this Letter, three analytical models are constructed in closed forms, each
representing a supermassive black hole (SMBH) located at the center of a
galaxy surrounded by a dark matter halo. The latter is modeled by the Einstein
cloud. The density profile of the halo vanishes inside two Schwarzschild radii
of the SMBH and satisfies the weak, strong and dominant energy conditions. The
spacetime is asymptotically flat, and the difference among the three models
lies in the slopes of the density profiles in the far region from the galaxy
center.

##  I Introduction

The overwhelming evidence for dark matter (DM) is clear and compelling on a
wide range of scales [1, 2, 3, 4, 5]. The motion of galaxies within clusters
[6, 7], galaxy rotation curves [8, 9], bullet cluster systems [10], baryonic
acoustic oscillation [11] and cosmic microwave background [12] all point to
its existence. According to the latest Planck data [13], about $5/6$ of the
mass in the Universe is made of DM. However, despite of numerous efforts in
the search of DM in the past decades, our understanding is still based only on
their gravitational interactions. Therefore, uncovering the nature and origin
of DM is arguably one of the greatest challenges of modern physics and
cosmology [14, 15, 16].

According to the current galaxy formation scenario, every galaxy forms within
a DM halo, and at the center of almost every galaxy a supermassive black hole
(SMBH) is present [17]. As a matter of fact, two such SMBHs have been already
observed by the Event Horizon Telescope (EHT), one is the Sagittarius $A^{*}$
black hole (BH) at the centre of our own Milky Way galaxy [18], and the other
called $M87^{*}$ is located at the centre of the more distant Messier 87
galaxy [19]. With the arrival of the gravitational wave (GW) astronomy [20],
the studies of SMBHs in galaxies with DM halos have gained further momenta,
see, for example, [21, 22, 23, 24, 25] and references therein. This is because
such studies bear important consequences for compact objects orbiting around a
SMBH. The mass ratio $q\equiv m/M$ of such binaries is typically of
${\cal{O}}(10^{-3}\sim 10^{-4})$ for intemediate mass-ratio inspirals (IMRIs)
[26], and ${\cal{O}}(10^{-5}\sim 10^{-6})$ for extreme mass-ratio inspirals
(EMRIs) [27]. GWs emitted by IMRIs can be detected by the upcoming ground-
based detectors, such as the Einstein Telescope [28] and the Cosmic Explorer
[29], and the ones from EMRIs are ideal targets for space-based detectors,
such as LISA [30], TianQin [31], Taiji [32] and DESIGO [33].

When a compact object moves inside the DM halo, gravitational drags will be
induced on the inspiraling object, changing its orbit and gravitational
waveforms emitted by it. For IMRIs, such effects can be detected by LISA [34,
35, 36]. Similar effects for EMRIs are also expected and can be used to study
the nature of DM [37, 38, 39]. For binaries with a very small $q$, the compact
objects can be well described as test particles, moving under the influence of
SMBHs and DM halos. Therefore, in this Letter we shall focus ourselves on the
gravitational field produced by a SMBH in a galaxy with the presence of a DM
halo.

Using Newtonian gravity, Gondolo and Silk first studied the effects of a
massive BH sitting at the center of a galaxy on the distribution of the DM
halo, and found that a spike is always produced at a distance of $4r_{s}$,
where $r_{s}$ denotes the Schwarzschild radius of the SMBH [40]. The orbits of
the DM particles with radii less than $4r_{s}$ become unstable, so that the
density of the DM halo vanishes for $r\lesssim 4r_{s}$. Sadeghian, Ferrer and
Will later took the relativistic effects into full considerations and found a
spike is indeed developed [41]. However, they found that the orbits of the DM
particles now start to be unstable at $2r_{s}$, instead of $4r_{s}$, purely
due to the relativistic effects, which is further confirmed by Speeney et al
[42]. Due to the complexity of the problem, such studies have been carried out
mainly numerically, except one particular case [43], in which Cardoso et al
constructed an analytic model with a density profile that asymptotically
approaches the Hernquist density contribution [44]. This solution has
immediately attracted lots of attention, and its various properties and
applications have already been explored, including its Love numbers, quasi-
normal modes (QNMs), tails and grey-body factors [43, 45, 46, 47], shadows
[48], and the detectability of GWs emitted by EMRIs by LISA, TianQin and Taiji
[49], to name only a few of them. Despite of all its interesting properties
and aspects, the model suffers from a couple of defects, for example, the
density of the DM halo vanishes only inside the BH horizon, instead in the
region $r\leq 2r_{s}$, as shown explicitly in [41, 42]. In addition, near the
BH horizon, the DM halo does not satisfy the dominant energy condition [50].
Extensions of this model to different circumstances have been also studied,
see, for example, [51, 52, 53] and references therein.

In this Letter, we present three analytic solutions, each of which represents
a SMBH located at the center of a galaxy surrounded by a DM halo, in which the
energy density of the halo satisfies all the energy conditions, weak, strong
and dominant [50], and identically vanishes 1for $r\leq 2r_{s}$. The spacetime
is smoothly joint to the vacuum Schwarzschild BH solution across $r=2r_{s}$.
The difference among the three models lies in the slopes of the density
profiles in the far region from the galaxy center.

##  II General Setup

Let us consider a spherical Einstein cloud made of collisionless DM particles
all moving tangentially with zero total angular momentum. Then, the cloud can
be considered as a fluid with a negligible radial pressure [43]. The
corresponding energy-momentum tensor of the DM halo is given by
$T^{\mu}_{\nu}=(-\rho(r),0,P(r),P(r))$, here $\rho$ and $P$ denote
respectively the energy density and tangential pressures of the halo. The
spacetime is describe by
$ds^{2}=-f(r)dt^{2}+(1-2m(r)/r)^{-1}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})$, where $f(r)$ and $m(r)$ are solutions of the Einstein field
equations. Then, the Einstein field equations have three independent
equations, given respectively by

$\displaystyle\frac{f^{\prime}}{f}=\frac{{2m}}{r(r-2m)},\;\;m^{\prime}=4\pi
Gr^{2}\rho,\;\;P=\frac{m\rho}{2(r-2m)},$ (1)

where a prime denotes the derivative wrt $r$. Note that in the current case we
have four unknowns but only three independent equations. To determine them
uniquely, the equation of state can be imposed. In this Letter, instead we
consider the choice [54, 42]

$\rho(r)=\frac{\rho_{0}\left(1-{4M}/{r}\right)}{({r}/{a})^{\gamma}\left(1+({r}/{a})^{\alpha}\right)^{(\beta-\gamma)/\alpha}},$
(2)

where $a$ represents the characteristic scale of the DM halo, the parameters
$(\alpha,\beta,\gamma)$ control the slopes of the density profile in the near
($r\ll a$), far ($r\gg a$) and transition ($r\simeq a$) regions, and
$\rho_{0}$ is a characteristic density of the halo, while $M$ the BH mass.
Observationally, it is found that $0.005{\text{kpc}}\lesssim a\lesssim
50{\text{kpc}}$. Recall that $r_{s}\equiv 2M=9.608\times
10^{-11}(M/10^{6}\;M_{\odot})\;{\text{kpc}}$. So, in general we have $a\gg
r_{s}$. The above model includes several important cases [54], such as
$(\alpha,\beta,\gamma)=\\{(1,4,1),(1,3,1)\\}$, which correspond respectively
to the Hernquist [44] and Navarro-Frenk-White (NFW) [55] profiles without BHs
($M=0$). As mentioned above, when a central BH is present, the trajectories of
the DM particles become unstable for $r<2r_{s}$ and will eventually fall
inside the BH, whereby a spike is developed [40, 41, 42]. So, we must have
$\rho(r)=0$, $f(r)=f_{o}(1-r_{s}/r)$ and $m(r)=M$ for $r<2r_{s}$, where
$f_{o}$ is a constant, which will be determined by the matching conditions
across $r=2r_{s}$. Note that in the rest of the Letter we shall use both $M$
and $r_{s}$ to denote the BH mass and radius. From (1) we can see that the
Einstein field equations contain only the first-order derivatives of $f$ and
$m$. Therefore, as long as they are continuous, these equations will be
satisfied. For the choice (2), it is clear that $\rho$ and $P$ are already
continuous, while $f$ can be always made continuous by properly choosing
$f_{o}$. Therefore, in the following we shall impose the boundary condition
$m(2r_{s})=M$. With such choices of $f_{o}$ and $m(r)$, the Einstein field
equations (1) can be solved in the whole spacetime $r\in(0,\infty)$. In
addition, we shall also assume that the spacetime is asymptotically flat, so
that $f(\infty)=1$ and $m(\infty)=M+M_{h}$, where $M_{h}$ is the total mass of
the DM halo. Observations show that $M\ll M_{h}\ll a$ [55, 1, 42].

##  III Analytical Models

With the above setup, we are ready to integrate (1) for the choice of the
density profile (2), for which we find a large class of analytical solutions.
In the following, we shall present only three of them, as the expressions of
others are rather lengthy and shall present them in detail in the following-up
paper [56].

Model I: The first solution corresponds to the Hernquist choice of
$(\alpha,\beta,\gamma)$ and is given by

$\displaystyle m(r)$ $\displaystyle=$ $\displaystyle
M+M_{h}\frac{(r-4M)^{2}}{(r+a)^{2}},$ $\displaystyle f(r)$ $\displaystyle=$
$\displaystyle\frac{f_{I}}{r}\left(\frac{\Delta_{1}(r)}{\Delta_{1}}\right)^{\alpha_{1}}\left(\frac{\Delta_{2}(r)}{\Delta_{2}}\right)^{{(1-\alpha_{1})}/{2}}$
(3)
$\displaystyle\times\exp\Bigg{[}D_{o}\left(\tan^{-1}\left(\frac{a+r+k}{\Delta_{3}}\right)-\frac{\pi}{2}\right)\Bigg{]},~{}~{}~{}$

where $\Delta_{1}(r)\equiv a+r-x_{1}$,
$\Delta_{2}(r)\equiv(a+r+k)^{2}+|x_{2}|^{2}$,
$\Delta_{1,2}\equiv\Delta_{1,2}(2r_{s})$,
$\Delta_{3}\equiv\sqrt{|x_{2}|^{2}-k^{2}}$,
$f_{I}\equiv\Delta_{1}^{\alpha_{1}}\Delta_{2}^{(1-\alpha_{1})/2}$, and
$D_{o}\equiv({\alpha_{3}-k(1-\alpha_{1})})/{\Delta_{3}}$. Additionally,
$x_{1,2,3}$ are the three roots of ${\cal{D}}_{1}(x)\equiv
x^{3}-Ax^{2}+Bx-C=0$ with $x_{3}=x_{2}^{*}$ and $x_{1}=x_{1}^{*}$, $A\equiv
a+2(M+M_{h})$, $B\equiv 8\pi a^{4}\rho_{0}$, $C\equiv 2M_{h}(a+4M)^{2}$,
$M_{h}=2\pi\rho_{0}a^{4}/(a+4M)$, and

$\displaystyle\alpha_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{x_{1}^{2}}{x_{1}^{2}+|x_{2}|^{2}+2kx_{1}},\;\;\;\alpha_{2}\equiv\frac{|x_{2}|^{2}x_{1}}{x_{1}^{2}+|x_{2}|^{2}+2kx_{1}},$
$\displaystyle k$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left(\sqrt[3]{\frac{1}{2}\left(-q+\sqrt{\Delta}\right)}+\sqrt[3]{\frac{1}{2}\left(-q-\sqrt{\Delta}\right)}\right)+\frac{A}{3},$

where $p\equiv B-{A^{2}}/{3}$, $q\equiv-{2A^{3}}/{27}+{AB}/{3}-C$,
$\Delta\equiv q^{2}+{4p^{3}}/{27}$. It can be shown that $\Delta>0$ for
$10M<M_{h}<a/10$ [56], so that the two roots $x_{2,3}$ of ${\cal{D}}_{1}(x)=0$
are complex conjugate, while $x_{1}$ is real. It can be also shown that the
matching conditions at $r=2r_{s}$ and the asymptotically-flat conditions are
all satisfied.

On the other hand, defining ${\cal{G}}(r)\equiv r-2m(r)$, we find
${\cal{G}}^{\prime}(r)=1-4M_{h}(r-4M)(a+4M)/(r+a)^{3}>0$. Therefore,
${\cal{G}}(r)$ is a monotonically increasing function with its minimal value
$2M$ given at $r=2r_{s}$. Since $\rho+P=(1+m/2{\cal{G}})\rho$,
$\rho+2P=(1+m/{\cal{G}})\rho$, we can see that the conditions $\rho\geq
0,\;\rho+P\geq 0$ and $\rho+2P\geq 0$ always hold for $r\geq 2r_{s}$, so the
DM halo satisfies both of the weak and strong energy conditions [50]. We also
have $\rho-P=\rho[5-r/{\cal{G}}(r)]/4$. Since ${\cal{G}}(r)\geq 2M$, we find
that $m/r=M/r+M_{h}(r-4M)^{2}/[r(r+a)^{2}]<1/4+M_{h}r/(r+a)^{2}$, where
$r/(r+a)^{2}$ has its maximal value at $r=a$, so we have $m/r<(a+M_{h})/4a$,
and $r/{\cal{G}}(r)<\underline{}2a/(a-M_{h})$. Hence, we have
$\rho-P>(3a-5M_{h})\rho/[4(a-M_{h})]>0$ for $10M<M_{h}<a/10$. Thus, the halo
satisfies the dominant energy condition, too [50].

Model II: In this model, we choose $(\alpha,\beta,\gamma)=(1,5,1)$. Then, we
find

$\displaystyle
m(r)=M+M_{h}\Bigg{(}1+2\frac{(a+4M)^{3}}{(r+a)^{3}}-3\frac{(a+4M)^{2}}{(r+a)^{2}}\Bigg{)},~{}$
(5)

but now with $M_{h}\equiv 2\pi a^{5}\rho_{0}/[3(a+4M)^{2}]$. Inserting the
above expression into (1), we find
${f^{\prime}}/{f}={(r+a)^{3}}/[{a^{4}{\cal{D}}_{2}(y)}]-{1}/{r}$, where
${\cal{D}}_{2}(y)\equiv y^{4}-{\cal{A}}y^{3}+{\cal{B}}y-{\cal{C}}$ with
$y\equiv 1+r/a$, ${\cal{A}}\equiv 1+2(M+M_{h})/a$, ${\cal{B}}\equiv
6M_{h}(a+4M)^{2}/a^{3}$ and ${\cal{C}}\equiv 4M_{h}(a+4M)^{3}/a^{4}$. It can
be shown that ${\cal{D}}_{2}(y)=0$ has four roots, and two of them are real,
denoted by $y_{1}$ and $y_{2}$, and the other two are complex conjugate,
$y_{4}^{*}=y_{3}\equiv y_{3R}+iy_{3I}$, which are all given explicitly in [56]
in terms of $a,M$ and $M_{h}$. Then, we find that

$\displaystyle
f(r)=\frac{1}{y-1}\left(y-y_{1}\right)^{\beta_{1}}\left(y-y_{2}\right)^{\beta_{2}}\left[\left(y-y_{3R}\right)^{2}+{y_{3I}}^{2}\right]^{\beta_{3}}\exp\left\\{\hat{\cal{D}}_{o}\left[\tan^{-1}\left(\frac{y-y_{3R}}{y_{3I}}\right)-\frac{\pi}{2}\right]\right\\},$
(6)

where $\beta_{1}\equiv y_{1}^{3}/[(y_{1}-y_{2})D_{1}]$, $\beta_{2}\equiv
y_{2}^{3}/[(y_{2}-y_{1})D_{2}]$, $\beta_{3}\equiv(1-\beta_{1}-\beta_{2})/2$,
$\beta_{4}\equiv|y_{3}|^{2}((y_{1}+y_{2})y_{3}^{2}-2y_{3R}y_{1}y_{2})/D_{1}D_{2}$,
$\hat{\cal{D}}_{o}\equiv(\beta_{4}+2y_{3R}\beta_{3})/y_{3I}$, and
$D_{n}\equiv|y_{3}|^{2}+y_{n}^{2}-2y_{3R}y_{n}$. It can be shown that the
above solution satisfies the asymptotically-flat and matching conditions
mentioned above, by properly choosing the integration constant $f_{o}$ of
$f(r)$ in the region $r<2r_{s}$.

To show that all three energy conditions hold also in the current case, we
first note that
${\cal{G}}^{\prime\prime}(r)=12M_{h}(a+4M)^{2}(3r-a-16M)/(a+r)^{5}$, which
vanishes at $r_{\text{min}}=(a+16M)/3$. This corresponds to a minimum of
${\cal{G}}^{\prime}(r)$, where
${\cal{G}}^{\prime}(r_{\text{min}})=1-3^{4}M_{h}/(4^{3}(a+4M))>0$ for
$10M<M_{h}<a/10$. Therefore, ${\cal{G}}^{\prime}(r)>0$ for $r>2r_{s}$. Hence,
in this region ${\cal{G}}(r)$ is monotonically increasing and
${\cal{G}}(r)\geq r_{s}$. Then, the DM halo also satisfies the weak and strong
energy conditions. To check the dominant energy condition, we first note that
the maximum of $[m(r)/r]$ is 1/4 for $10M<M_{h}<a/10$. Then, we find that
$\rho-P=(5-r/{\cal{G}})\rho/4>(5-1/(1-0.5))\rho>0$. Hence, the dominant energy
condition also holds in the present case.

Model III: In this model, we choose $(\alpha,\beta,\gamma)=(1,\frac{5}{2},1)$.
Then we find

$\displaystyle m(r)$ $\displaystyle=$ $\displaystyle
M+{\cal{M}}\Bigg{[}\left(1+\frac{r}{a}\right)^{1/2}-2\left(1+\frac{4M}{a}\right)^{1/2}$
(7)
$\displaystyle+\left(1+\frac{4M}{a}\right)\left(1+\frac{r}{a}\right)^{-1/2}\Bigg{]},$

where ${\cal{M}}\equiv 8\pi Ga^{3}\rho_{0}$, and

$\displaystyle\frac{f^{\prime}(r)}{f(r)}$ $\displaystyle=$
$\displaystyle-\frac{1}{a(z^{2}-1)}+\frac{z}{a{\cal{D}}_{3}(z)},$ (8)

with $z\equiv(1+r/a)^{1/2}$, ${\cal{D}}_{3}(z)\equiv
z^{3}-uz^{2}+(2u\sqrt{1+4v}-2v-1)z-u(1+4v)$, where $v\equiv M/a,u\equiv
2{\cal{M}}/a$. Thus, ${\cal{D}}_{3}(z)=0$ has three roots, denoted by
$z_{1,2,3}$. Depending on the values of $u$ and $v$, the nature of these roots
are different. In particular, if $u<0.1709$, we find that all of them are
real. If $u>0.1716$, two of them are complex conjugate, and only one is real,
say, $z_{1}$. If $0.1709<u<0.1716$, we find that whether the three roots are
all real or not depends on the values of $v$. In the following, let us
consider the two cases separately.

Before doing so, it is interesting to note that
$m(r)\propto\sqrt{r}\rightarrow\infty$ as ${r}\rightarrow\infty$. Thus, now
the mass of the DM halo becomes infinitely large. In fact, this is always the
case for $\beta\leq 3$, including the NFW model [55], for which $\beta=3$ and
$m(r)\propto\ln r$. In reality, a halo is always restricted to a finite
region, so one needs to restrict the expression of (7) to be valid only up to
a maximal radius, say, $r_{\infty}$. Then, the total mass of the halo will be
finite. In our current case, the mass does not increase fast enough to close
the whole spacetime, so we still have an asymptotically-flat spacetime, as now
we still have $f(\infty)=1$ and $m(r)/r\rightarrow 0$ as
${r}\rightarrow\infty$, as to be shown below.

With the above in mind, let us first consider the case in which all three
roots are real. Then, we find that

$\displaystyle
f(r)=\frac{a}{r}\frac{\big{(}z-z_{1}\big{)}^{\gamma_{1}}\big{(}z-z_{2}\big{)}^{\gamma_{2}}}{\big{(}z-z_{3}\big{)}^{\gamma_{1}+\gamma_{2}-2}},$
(9)

where $\gamma_{1}\equiv 2z_{1}^{2}(z_{2}-z_{3})D_{3}/[z_{3}(z_{1}-z_{2})]$,
$\gamma_{2}\equiv 2z_{2}^{2}(z_{1}-z_{3})D_{3}/[z_{3}(z_{2}-z_{1})]$,
$D_{3}\equiv
z_{3}(z_{2}-z_{1})/[z_{1}(z_{2}^{2}-z_{3}^{2})+z_{2}(z_{3}^{2}-z_{1}^{2})+z_{3}(z_{1}^{2}-z_{2}^{2})]$,
with $z_{1,2,3}$ being the three real roots of ${\cal{D}}_{3}(z)=0$ [56]. The
above expressions shows that indeed we have $f(\infty)=1$. On the other hand,
from (7) we find that $m(2r_{s})=M$. Therefore, by properly choosing the
integration constant $f_{o}$ in the region $r<2r_{s}$, the solution can be
smoothly matched to that of the Schwarzschild vacuum solution valid in the
region $r\leq 2r_{s}$. We have also found that the corresponding DM halo
satisfies all the three energy conditions. However, the proof is rather
tedious and will be presented in detail in [56]. So, here we shall omit the
proof. The same will be done for the case in which only one of the three roots
is really.

When only one of the three roots is real, which will be chosen as $z_{1}$, we
find

$\displaystyle f(r)$ $\displaystyle=$
$\displaystyle\frac{a}{r}\big{(}z-z_{1}\big{)}^{\delta_{1}}\big{(}(z-z_{2R})^{2}+z_{2I}^{2}\big{)}^{\frac{2-\delta_{1}}{2}}$
(10)
$\displaystyle\times\exp\left\\{\frac{\tilde{\cal{D}}_{o}}{z_{2I}}\left[\tan^{-1}\left(\frac{z-z_{2R}}{z_{2I}}\right)-\frac{\pi}{2}\right]\right\\},~{}~{}~{}$

where $\delta_{1}=2az_{1}^{2}/\hat{\cal{D}}_{3}$, $z_{2}=z_{2R}+iz_{2I}$,
$\tilde{\cal{D}}_{o}=2a(K_{3}+2z_{2R}K_{2})z_{2R}-2a|z_{2}|^{2}K_{2}$,
$K_{2}=-z_{1}/\hat{\cal{D}}_{3}$, $K_{3}=|z_{2}|^{2}/\hat{\cal{D}}_{3}$, and
$\hat{\cal{D}}_{3}\equiv a[z_{1}^{2}+|z_{2}|^{2}-2z_{1}z_{2R}]$. Thus, we have
$f(\infty)=1$ and $m(r)/r\rightarrow 0$ as $r\rightarrow\infty$. As a result,
the spacetime is asymptotically flat, despite the fact
$m(r)\sim{\cal{M}}\sqrt{r/a}$ for $r\gg a$.

##  IV Conclusions

In this Letter, we have presented three analytic models, each of which
describe a SMBH sitting at the center of a galaxy and being surrounded by a DM
halo. The halo satisfies all the three energy conditions, weak, strong and
dominant [50], and vanishes within two Schwarzschild radii, $r\leq 2r_{s}$, of
the SMBH. The latter is consistent with previous investigations of
relativistic DM halo models in the presence of a SMBH [41, 42], in which it
was shown that the trajectories of DM particles become unstable when
$r<2r_{s}$ and these particles shall all fall inside the BH. As a result, we
have $\rho=0$ for $r<2r_{s}$. To our best knowledge, these are the only
analytic models found so far that have such desirable properties.

We expect that the construction of these analytic models will significantly
facilitate the studies of QNMs of SMBHs and GWs emitted by IMRIs and EMRIs,
which are the ideal sources of the next-generation GW detectors, both ground-
and space-based detectors [28, 29, 30, 31, 32, 33]. As the small compact
companion of the SMBH passes through the DM halo, it feels gravitational
dragging forces, whereby its orbits and gravitational waveforms emitted by it
will be changed. Such effects can be detected by LISA and other next
generation detectors [34, 35, 36, 37, 38, 39]. Hence. a completely different
window to study the nature of DM is provided.

Certainly, the applications of these models are not restricted to the studies
mentioned above but also equally applicable to other fields, such as shadows
of BHs, which can be the sources of EHT [18, 19], and the dynamics of galaxy
systems, including our own Milky Way galaxy [1, 2, 3, 4, 5, 44, 55].

## Acknowledgements

We thank Profs. Katherine Freese and Monica Valluri for valuable comments and
suggestions. Z.S. thanks the hospitality of Physics Department and Baylor
University during his visit. A.W. is partially supported by the U.S. Natural
Science Foundation (NSF) under the Grant No. PHY2308845. Y.G. is partially
supported by the National Key Research and Development Program of China under
Grant No. 2020YFC2201504.

## References

  * [1] G. Bertone, D. Hooper, and J. Silk, Particle Dark Matter: Evidence, Candidates and Constraints, Phys. Rep. 405, 279 (2005) [arXiv:hep-ph/0404175].
  * [2] K. Freese, Review of Observational Evidence for Dark Matter in the Universe and in upcoming searches for Dark Stars, EAS Publ. Ser. 36, 113 (2009) [arXiv:0812.4005].
  * [3] J. G. de Swart, G. Bertone, and J. van Dongen, How dark matter came to matter, Nature Astron. 1, 0059 (2017) [arXiv:1703.00013].
  * [4] R.H. Wechsler, J.L. Tinker, The Connection between Galaxies and their Dark Matter Halos, Annu. Rev. Astron. Astrophys. 56, 435 (2018) [arXiv:1804.03097 ].
  * [5] A. Arbey, F. Mahmoudi, Dark matter and the early Universe: a review, Prog. Part. Nucl. Phys. 119 (2021) 103865 [arXiv: 2104.11488].
  * [6] F. Zwicky, Die Rotverschiebung von extragalaktischen Nebeln, Helv. Phys. Acta, 6, 110 (1933)[arXiv: 1711.01693].
  * [7] F. Zwicky, On the masses of nebulae and of clusters of nebulae, Astrophys. J. 86, 217 (1937).
  * [8] V.C. Rubin, W.K. Jr. Ford, Rotation of the Andromeda nebula from a spectroscopic survey of emission regions, Astrophys. J. 159, 379 (1970).
  * [9] K.C. Freeman, On the disks of spiral and S0 galaxies, Astrophys. J. 160, 811 (1970).
  * [10] D. Clowe, et al., A direct empirical proof of the existence of dark matter, Astrophys. J. 648, L109 (2006) [arXiv:astro-ph/0608407].
  * [11] D. J. Eisenstein, et al., Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies, Astrophys. J. 633, 560 (2005) [arXiv:astro-ph/0501171].
  * [12] C. L. Bennett, et al., Nine-year Wilkinson Microwave anisotropy probe (WMAP) observa- tions: Final maps and results, Astrophys. J. Suppl. Ser. 208, 20 (2013) [arXiv:1212.5226].
  * [13] N. Aghanim, etal (Planck Collaboration), Planck 2018 results VI. Cosmological parameters, A$\&$A 641, A6 (2020) [arXiv:1807.06209].
  * [14] S. Cebrián, Review on dark matter searches, J. Phys.: Conf. Ser. 2502 (2023) 012004 [arXiv:2205.06833].
  * [15] C. Smarra, et al (European Pulsar Timing Array), Second Data Release from the European Pulsar Timing Array: Challenging the Ultralight Dark Matter Paradigm, Phys. Rev. Lett. 131, 171001 (2023)[arXiv:2306.16228].
  * [16] M. Misiaszeka and N. Rossib, Direct detection of dark matter: a critical review,[ arXiv:2310.20472].
  * [17] M. Valluri, D. Merritt, and E. Emsellem, Difficulties with Recovering the Masses of Supermassive Black Holes from Stellar Kinematical Data, Astrophys. J. 602, 66 (2004) [arXiv:astro-ph/0210379].
  * [18] K. Akiyama, et al (Event Horizon Telescope Collaboration), First Sagittarius $A^{*}$ Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way, Astrophys. J. Lett. 930 (2022) L12.
  * [19] K. Akiyama, et al (Event Horizon Telescope Collaboration), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875 (2019) L1[arXiv:1906.11238].
  * [20] R. Abbott, et al (LIGO Scientific, VIRGO and KAGRA Collaborations), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, [arXiv:2111.03606].
  * [21] A. El-Zant, Z. Roupas, J. Silk, Ejection of supermassive black holes and implications for merger rates in fuzzy dark matter haloes, MNRAS 499, 2575 (2020) [arXiv:2009.10167].
  * [22] L.E. Padilla, T. Rindler-Daller, P.R. Shapiro, T. Matos, J. A. Vázquez, Core-Halo Mass Relation in Scalar Field Dark Matter Models and its Consequences for the Formation of Supermassive Black Holes, Phys. Rev. D103, 063012 (2021) [arXiv:2010.12716].
  * [23] A. Marasco, G. Cresci, L. Posti, F. Fraternali, F. Mannucci, A. Marconi, F. Belfiore, S.M. Fall, A universal relation between the properties of supermassive black holes, galaxies, and dark matter haloes, MNRAS 507, 4274 (2021) [arXiv:2105.10508].
  * [24] A. Bansal, K. Ichiki, H. Tashiro, Y. Matsuoka, Evolution of the mass relation between supermassive black holes and dark matter halos across the cosmic time, MNRAS 523, 3840 (2023) [arXiv:2206.01443].
  * [25] H. Zhang, P. Behroozi, M. Volonteri, J. Silk, X. Fan, P.F. Hopkins , J. Yang, J. Aird, TRINITY I: self-consistently modelling the dark matter halo-galaxy-supermassive black hole connection from $z=0-10$, MNRAS 518, 2123 (2023) [arXiv:2105.10474].
  * [26] P. Amaro-Seoane, J.R. Gair, M. Freitag, M. Coleman Miller, I. Mandel, C.J. Cutler, and S. Babak, Intermediate and extreme mass-ratio inspirals—astrophysics, science applications and detection using LISA, Class. Quantum Grav. 24, R113 (2007)[arXiv:astro-ph/0703495].
  * [27] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Science with the space-based interferometer LISA. V. Extreme mass-ratio inspirals, Phys. Rev. D95, 103012 (2017)[arXiv:1703.09722].
  * [28] F. Amann, et al., Site-selection criteria for the Einstein Telescope, Rev. Sci. Instrum. 91, 094504 (2020)[arXiv:2003.03434]; B. Sathyaprakash et al., Corrigendum: Scientific objectives of Einstein telescope, Class. Quantum Grav. 30 (2013) 079501; M. Punturo, et al., The Einstein Telescope: a third-generation gravitational wave observatory, Class. Quantum Grav. 27, 194002 (2010); http://www.et-gw.eu.
  * [29] E.D. Hall, et al., Gravitational-wave physics with Cosmic Explorer: limits to low-frequency sensitivity, arXiv:2012.03608; B.P. Abbott, et al., Exploring the sensitivity of next generation gravitational wave detectors, Class. Quantum Grav. 34, 044001 (2017)[arXiv:1607.08697 ]; S. Dwyer, D. Sigg, S.W. Ballmer, L. Barsotti, N. Mavalvala, M. Evans, Gravitational wave detector with cosmological reach, Phys. Rev. D91, 082001 (2015)[arXiv:1410.0612]; https://cosmicexplorer.org.
  * [30] P. Amaro-Seoane, et al., Laser Interferometer Space Antenna, [arXiv:1702.00786v3].
  * [31] J. Luo, et al., TianQin: a space-borne gravitational wave detector, Class. Quantum Grav. 33, 035010 (2016)[arXiv:1512.02076].
  * [32] Z.-R. Luo, Z.-K. Guo, G. Jin, Y.-L. Wu, and W.-R. Hua, A brief analysis to Taiji: Science and technology, Results in Phys. 16 102918 (2020); W.-H. Ruan, Z.-K. Guo, R.-G. Cai, Y.-Z. Zhang, Taiji Program: Gravitational-Wave Sources, arXiv:1807.09495v2 [gr-qc].
  * [33] S. Kawamura, et al., Current status of space gravitational wave antenna DECIGO and B-DECIGO, arXiv:2006.13545; S. Sato, et al., The status of DECIGO, J. Phys.: Conf. Series 840, 012010 (2017).
  * [34] K. Eda, Y. Itoh, S. Kuroyanagi, J. Silk, Gravitational waves as a probe of dark matter minispikes, Phys. Rev. D91, 044045 (2015) [arXiv:1408.3534].
  * [35] B.J. Kavanagh, D.A. Nichols, G. Bertone, D. Gaggero, Detecting dark matter around black holes with gravitational waves: Effects of dark-matter dynamics on the gravitational waveform, Phys. Rev. D102, 083006 (2020) [arXiv:2002.12811].
  * [36] A. Coogan, G. Bertone, D. Gaggero, B.J. Kavanagh, D.A. Nichols, Measuring the dark matter environments of black hole binaries with gravitational waves, Phys. Rev. D105, 043009 (2022) [arXiv:2108.04154].
  * [37] O.A. Hannuksela, K.C.Y. Ng, T.G.F. Li, Extreme Dark Matter Tests with Extreme Mass Ratio Inspirals, Phys. Rev. D102, 103022 (2020) [arXiv:1906.11845].
  * [38] T.-G. Zi, J.-D. Zhang, H.-M. Fan, X.-T. Zhang, Y.-M. Hu, C. Shi, J. Mei, Science with the TianQin Observatory: Preliminary Results on Testing the No-hair Theorem with EMRI, Phys. Rev. D 104, 064008 (2021) [arXiv:2104.06047].
  * [39] H.-M. Fan, S. Zhong, Z.C. Liang, Z. Wu, J.-d. Zhang, Y.-M. Hu, Extreme-mass-ratio burst detection with TianQin, Phys. Rev. D106, 124028 (2023) [arXiv:2209.13387].
  * [40] P. Gondolo and J. Silk, Dark Matter Annihilation at the Galactic Center, Phys. Rev. Lett. 83, 1719 (1999).
  * [41] L. Sadeghian, F. Ferrer, and C. M. Will, Dark-matter distributions around massive black holes: A general relativistic analysis, Phys. Rev. D88, 063522 (2013)[arXiv:1305.2619].
  * [42] Ni. Speeney, A. Antonelli, V. Baibhav, and E. Berti, Impact of relativistic corrections on the detectability of dark-matter spikes with gravitational waves, Phys. Rev. D106, 044027 (2022) [arXiv:2204.12508v1].
  * [43] V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, A. Maselli, Black holes in galaxies: environmental impact on gravitational-wave generation and propagation, Phys. Rev. D105, L061501 (2022) [arXiv:2109.00005].
  * [44] L. Hernquist, An Analytical Model for Spherical Galaxies and Bulges, Astrophys. J. 356, 359 (1990).
  * [45] R.A. Konoplya, Black holes in galactic centers: quasinormal ringing, grey-body factors and Unruh temperature, Phys. Lett. B823 (2021) 136734 [arXiv:2109.01640].
  * [46] J. Liu, S. Chen, J. Jing, Tidal effects of dark matter halo around a galactic black hole, Chin. Phys. C46, 105104 (2022) [arXiv:2203.14039].
  * [47] K. Destounis, A. Kulathingal, K.D. Kokkotas, G.O. Papadopoulos, Gravitational-wave imprints of compact and galactic-scale environments in extreme-mass-ratio binaries, Phys. Rev. D107, 084027 (2023) [arXiv:2210.09357].
  * [48] S.V.M.C.B. Xavier, H.C.D. Lima Junior, L.C.B. Crispino, Shadows of black holes with dark matter halo, Phys. Rev. D107, 064040 (2023) [arXiv:2303.17666].
  * [49] N. Dai, Y. Gong, Y. Zhao, T. Jiang, Extreme mass ratio inspirals in galaxies with dark matter halos, [arXiv:2301.05088].
  * [50] S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time (Cambridge University Press, Cambridge, England, 1973).
  * [51] S. Nampalliwar, Saurabh K., K. Jusufi, Q. Wu, M. Jamil, P. Salucci, Modelling the Sgr $A^{*}$ Black Hole Immersed in a Dark Matter Spikez, Astrophys. J. 916, 116 (2021) [arXiv:2103.12439].
  * [52] R. A. Konoplya, A. Zhidenko, Solutions of the Einstein equations for a black hole surrounded by a galactic halo, Astrophys. J. 933 (2022) 166 [arXiv:2202.02205].
  * [53] D. Liu, Y. Yang, Z. Xu, Z.-W. Long, Modelling the black holes surrounded by a dark matter halo in the galactic center of M87, arXiv:2307.13553.
  * [54] H.S. Zhao, Analytical models for galactic nuclei, Mon. Not. R. Astron. Soc. 278, 488 (1996) [arXiv:astro-ph/9509122].
  * [55] J.F. Navarro, C.F. Frenk, S.D.M. White, a universal density profile from hierarchical clustering, Astrophys. J. 490, 493 (1997) [arXiv:astro-ph/9611107].
  * [56] Z. Shen, A. Wang, Y. Gong, and S.-Y. Yin, Supermassive black holes of galaxies surrounded by dark matter halos: Realistic and analytical models, in preparation (2023).