# Actor-Critic Methods using Physics-Informed Neural Networks: Control of a 1D
PDE Model for Fluid-Cooled Battery Packs

Amartya Mukherjee1 , Jun Liu1
1Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario
N2L 3G1, Canada Corresponding author. e-mail<EMAIL_ADDRESS>
###### Abstract

This paper proposes an actor-critic algorithm for controlling the temperature
of a battery pack using a cooling fluid. This is modeled by a coupled 1D
partial differential equation (PDE) with a controlled advection term that
determines the speed of the cooling fluid. The Hamilton-Jacobi-Bellman (HJB)
equation is a PDE that evaluates the optimality of the value function and
determines an optimal controller. We propose an algorithm that treats the
value network as a Physics-Informed Neural Network (PINN) to solve for the
continuous-time HJB equation rather than a discrete-time Bellman optimality
equation, and we derive an optimal controller for the environment that we
exploit to achieve optimal control. Our experiments show that a hybrid-policy
model that updates the value network using the HJB equation and updates the
policy network identically to PPO achieves the best results in the control of
this PDE system.

keywords: Actor-Critic Method, Physics-Informed Neural Network, Fluid Cooled
Battery Packs, Hamilton Jacobi Bellman Equation

##  1 Introduction

In recent years, there has been a growing interest in Reinforcement Learning
(RL) for continuous control problems. RL has shown promising results in
environments with unknown dynamics through a balance of exploration in the
environment and exploitation of the learned policies. Since the advent of
REINFORCE with Baseline, the value network in RL algorithms has shown to be
useful towards finding optimal policies as a critic network (Sutton & Barto
(2018)). This value network continues to be used in state-of-the-art RL
algorithms today.

Proximal Policy Optimization (PPO) is an actor-critic method introduced by
Schulman et al. (2017). It limits the update of the policy network to a trust
region at every iteration. This ensures that the objective function of the
policy network is a good approximation of the true objective function and
forces smooth and reliable updates to the value network as well.

In discrete-time RL, the value function estimates returns from a given state
as a sum of the returns over time steps. This value function is obtained by
solving the Bellman Optimality Equation. On the other hand, in continuous-time
RL, the value function estimates returns from a given state as an integral
over time. This value function is obtained by solving a partial differential
equation (PDE) known as the Hamilton-Jacobi-Bellman (HJB) equation Munos
(1999). Both equations are difficult to solve analytically and numerically,
and therefore the RL agent must explore the environment and make successive
estimations.

The introduction of physics-informed neural networks (PINNs) by Raissi et al.
(2019) has led to significant advancements in scientific machine learning.
PINNs leverage auto-differentiation to compute derivatives of neural networks
with respect to their inputs and model parameters exactly. This enables the
laws of physics (described by ODEs or PDEs) governing the dataset of interest
to act as a regularization term for the neural network. As a result, PINNs
outperform regular neural networks on such datasets by exploiting the
underlying physics of the data.

Control of PDEs is considered to be challenging compared to control of ODEs.
Works such as Vazquez & Krstic (2017) introduced the backstepping method for
the boundary control of reaction-advection-diffusion equations using kernels.
For PDE control problems where the control input is encoded in the PDE, the
HJB equation has been used (Sirignano & Spiliopoulos (2018),Kalise & Kunisch
(2017)). Works from control of ODEs have been used by writing the PDE as an
infinite-dimensional ODE.

To the best of our knowledge, this paper is the first to explore the
intersection between PINNs and RL in a PDE control problem. We discretize the
PDE as an ODE to derive an HJB equation. In order to force the convergence of
the value network in PPO towards the solution of the HJB equation, we utilize
PINNs to encode this PDE and train the value network. Upon deriving the HJB
equation, we also derive an optimal controller. We introduce two algorithms:
HJB value iteration and Hamilton-Jacobi-Bellman Proximal Policy Optimization
(HJBPPO), that train the value function using the HJB equation and use the
optimal controller. The HJBPPO algorithm shows superior performance compared
to PPO and HJB value iteration on the PackCooling environment.

##  2 Preliminaries

###  2.1 The 1D pack cooling problem

The 1D system for fluid-cooled battery packs was introduced by Kato & Moura
(2021) and is modeled by the following coupled PDE:

$\displaystyle u_{t}(x,t)=$
$\displaystyle-D(x,t)u_{xx}(x,t)+h(x,t,u(x,t))+\frac{1}{R(x,t)}(w-u)$ (1)
$w_{t}=-\sigma(t)w_{x}+\frac{1}{R(x,t)}(u-w),$ (2)

with the following boundary conditions:

$u_{x}(0,t)=u_{x}(1,t)=0$ (3) $w(0,t)=U(t)$ (4)

where $u(x,t)$ is the heat distribution across the battery pack, $w(x,t)$ is
the heat distribution across the cooling fluid, $D(x,t)$ is the thermal
diffusion constant across the battery pack, $R(x,t)$ is the heat resistance
between the battery pack and the cooling fluid, $h(x,t,u)$ is the internal
heat generation in the battery pack, $U(t)$ is the temperature of the cooling
fluid at the boundary, and $\sigma(t)$ is the transport speed of the cooling
fluid, which will be the controller in this paper.

The objective of the control problem in this paper is to determine $\sigma(t)$
such that $u(x,t)$ is as close to zero as possible. The transport speed
$\sigma(t)$ is strictly non-negative so the cooling fluid travels only in the
positive $x$-direction. We restrict $\sigma(t)$ to $[0,1]$.

###  2.2 Hamilton-Jacobi-Bellman equation

To achieve optimal control for the 1D PDE pack cooling problem, we will
utilize works from control theory for ODEs. Consider a controlled dynamical
system modeled by the following equation:

$\dot{x}=f(x,\sigma),\quad x(t_{0})=x_{0},$ (5)

where $x(t)$ is the state and $\sigma(t)$ is the control input. In control
theory, the optimal value function $V^{*}(x)$ is useful towards finding a
solution to control problems (Munos et al. (1999)):

$\displaystyle V^{*}(x)=\sup_{\sigma}\frac{1}{\Delta
t}\int_{t_{0}}^{\infty}\gamma^{\frac{t}{\Delta t}}L($ $\displaystyle
x(\tau;t_{0},x_{0},\sigma(\cdot)),\sigma(\tau))d\tau,$ (6)

where $L(x,\sigma)$ is the reward function, $\Delta t$ is the time step size
for numerical simulation, and $\gamma$ is the discount factor. The following
theorem introduces a criteria for assessing the optimality of the value
function (Liberzon (2012), Kamalapurkar et al. (2018)).

######  Theorem 2.1.

A function $V(x)$ is the optimal value function if and only if:

  1. 1.

$V\in C^{1}(\mathbb{R}^{n})$ and $V$ satisfies the Hamilton-Jacobi-Bellman
(HJB) Equation

$\displaystyle(\gamma-1)V(x)+\sup_{\sigma\in U}\\{L(x,\sigma)+\gamma\Delta
t\nabla_{x}V^{T}(x)f(x,\sigma)\\}=0$ (7)

for all $x\in\mathbb{R}^{n}$.

  2. 2.

For all $x\in\mathbb{R}^{n}$, there exists a controller $\sigma^{*}(\cdot)$
such that:

$\displaystyle(\gamma-1)V(x)+L(x,\sigma^{*}(x))+\gamma\Delta
t\nabla_{x}V^{T}(x)f(x,\sigma^{*}(x))$
$\displaystyle=(\gamma-1)V(x)+\sup_{\hat{\sigma}\in
U}\\{L(x,\hat{\sigma})+\gamma\Delta t\nabla_{x}V^{T}(x)f(x,\hat{\sigma})\\}.$
(8)

The proof of part 1 of this theorem is in Appendix B. The HJB equation will be
used in this paper to determine a new loss function for the value network
$V(x)$ in this pack cooling problem and an optimal controller $\sigma^{*}(t)$.

##  3 Related work

The HJB equation we intend to solve is a first-order quasi-linear PDE. The use
of HJB equations for continuous RL has sparked interest in recent years among
the RL community as well as the control theory community and has led to
promising works. Kim et al. (2021) introduced an HJB equation for Q Networks
and used it to derive a controller that is Lipschitz continuous in time. This
algorithm has shown improved performance over Deep Deterministic Policy
Gradient (DDPG) in three out of the four tested MuJoCo environments without
the need for an actor network. Wiltzer et al. (2022) introduced a
distributional HJB equation to train the FD-WGF Q-Learning algorithm. This
models return distributions more accurately compared to Quantile Regression TD
(QTD) for a particle-control task. Finite difference methods are used to solve
this HJB equation numerically. Furthermore, the authors mentioned the use of
auto-differentiation for increased accuracy of the distributional HJB equation
as a potential area for future research in their conclusion.

The use of neural networks to solve the HJB equation has been an area of
interest across multiple research projects. Jiang et al. (2016) uses a
structured Recurrent Neural Network to solve the HJB equation and achieve
optimal control for the Dubins car problem. Tassa & Erez (2007) uses the
Pineda architecture (Pineda (1987)) to estimate partial derivatives of the
value function with respect to its inputs. They used the iterative least
squares method to solve the HJB equation. This algorithm shows convergence in
several control problems without the need for an initial stable policy.

RL for PDE control is a challenging field that has been of interest to the
machine learning community lately. Farahmand et al. (2017) introduces the Deep
Fitted Q Iteration to solve a boundary control problem for a 2D convection-
diffusion equation. The model stabilizes the temperature in the environment
without encoding any knowledge of the governing PDE. Sirignano & Spiliopoulos
(2018) develops the DGM algorithm to solve PDEs. They use auto-differentiation
to compute first-order derivatives and Monte Carlo methods to estimate higher-
order derivatives. This algorithm was used to solve the HJB equation to
control a stochastic heat equation and achieved an error of 0.1%. Kalise &
Kunisch (2017) approximates the solution to the HJB equation using
polynomials. This was used to control a semilinear parabolic PDE.

PINNs have been used for the control of dynamical systems in recent works.
Antonelo et al. (2021) uses a PINN for model predictive control of a dynamical
system over a long time interval. The PINN takes the initial condition, the
control input, and the spatial and temporal coordinates as input and estimates
the trajectory of the dynamical system while repeatedly shifting the time
interval towards zero to allow for long-range interval predictions. Nicodemus
et al. (2022) uses a PINN-based model predictive control for the tracking
problem of a multi-link manipulator. Djeumou et al. (2022) uses a PINN to
incorporate partial knowledge about a dynamical system such as symmetry and
equilibrium points to estimate the trajectory of a controlled dynamical
system.

The use of a PINN to solve the HJB equation for the value network was done by
Nakamura-Zimmerer et al. (2020) in an optimal feedback control problem
setting. The paper achieves results similar to that of the true optimal
control function in high-dimensional problems.

##  4 HJB control of the pack cooling problem

In this section, we will connect the pack cooling PDE model with the HJB
equation to derive a new loss function for the value network $V(u,w)$ using
the HJB equation and an optimal controller. The HJB equation has been useful
in finding optimal controllers for systems modeled by ODEs. In Kalise &
Kunisch (2017), the controlled PDE system has been discretized in space to
form an ODE that can be used in the HJB equation. Similarly, to form the HJB
equation for this paper, we need to write equations 1 and 2 as an ODE.

###  4.1 ODE discretization of PDE

We can write equations 1 and 2 as an ODE by discretizing it in the $x$
variable. By letting $\Delta x=\frac{1}{N_{x}}$ where $N_{x}$ is the number of
points we choose to discretize the system along the x-axis, we arrive at a
$2N_{x}$ dimensional ODE:

$\dot{\hat{U}}=-DA\hat{U}+h(\hat{U})+\frac{1}{R}(\hat{W}-\hat{U})$ (9)
$\dot{\hat{W}}=-\sigma(t)B\hat{W}+\frac{1}{R}(\hat{U}-\hat{W}),$ (10)

where

$\hat{W}(t)=\begin{pmatrix}w(x_{1},t)\\\ \vdots\\\
w(x_{N_{x}},t)\end{pmatrix},\hat{U}(t)=\begin{pmatrix}u(x_{1},t)\\\ \vdots\\\
u(x_{N_{x}},t)\end{pmatrix},$

and $A\hat{U}$ is a second-order discretization of $u_{xx}$, e.g.,

$[A\hat{U}]_{k}=\frac{u(x_{k+1},t)-2u(x_{k},t)+u(x_{k-1},t)}{\Delta x^{2}},$

$B\hat{W}$ is a second-order discretization of $w_{x}$, e.g.,

$[B\hat{W}]_{k}=\frac{w(x_{k+1},t)-w(x_{k-1},t)}{2\Delta x}.$

###  4.2 Derivation of the optimal controller

The ODE system derived in section 4.1 can be used in the HJB equation to
determine a loss function and an optimal controller.

######  Theorem 4.1.

Let $u(\cdot,t),w(\cdot,t)\in L_{2}[0,1]$. With $\sigma(t)\in[0,1]$ and the
reward function $L(U_{t},W_{t},\sigma_{t})=-||U_{t+1}||_{2}^{2}\Delta x$, the
HJB equation for the 1D pack cooling problem is:

$\displaystyle(\gamma-1)V-||u(\cdot,t+\Delta t)||^{2}$ $\displaystyle+\langle
V_{u}(u(\cdot,t),w(\cdot,t)),u_{t}(\cdot,t)\rangle$
$\displaystyle+\frac{1}{R}\langle
V_{w}(u(\cdot,t),w(\cdot,t)),u(\cdot,t)-w(\cdot,t)\rangle$
$\displaystyle+\max(0,-\langle
V_{w}(u(\cdot,t),w(\cdot,t)),w_{x}(\cdot,t)\rangle)=0$ (11)

where $||\cdot||$ is the $L_{2}[0,1]$ norm and $\langle\cdot,\cdot\rangle$ is
the $L_{2}[0,1]$ inner product.

The proof of this theorem is in Appendix C. Theorem 2.1 shows that there
exists a controller that satisfies equation 8. This allows us to determine an
optimal controller, as shown in the following corollary:

######  Corollary 4.2.

Let $w(\cdot,t)\in L_{2}[0,1]$. With $\sigma(t)\in[0,1]$ and the reward
function $L(U_{t},W_{t},\sigma_{t})=-||U_{t+1}||_{2}^{2}\Delta x$, provided
the optimal value function $V^{*}(u,w)$ with $V^{*}_{w}(\cdot,t)\in
L_{2}[0,1]$, the optimal controller for the 1D pack cooling problem is:

$\sigma^{*}(t)=\begin{cases}1,&\langle
V^{*}_{w}(u(\cdot,t),w(\cdot,t)),w_{x}(\cdot,t)\rangle<0,\\\
0,&\text{otherwise},\end{cases}$ (12)

where $\langle\cdot,\cdot\rangle$ is the $L_{2}[0,1]$ inner product.

The proof of this corollary is in Appendix D. These results will be used in
our algorithms to achieve optimal control of the pack cooling problem.

##  5 Algorithm

For the control of the PDE, we introduce two algorithms. The first algorithm,
called HJB Value Iteration, uses only a value network and exploits the HJB
equation and optimal controller derived in Theorem 4.1 and Corollary 4.2. The
second algorithm, called HJBPPO, is a hybrid-policy model that uses policy
network updates from PPO and value network updates from HJB Value Iteration.

To define these algorithms, we first define two loss functions. The first loss
function is derived from the proof of theorem 4.1.

$\displaystyle MSE_{f}=\frac{1}{T}\sum_{t=0}^{T-1}($
$\displaystyle(\gamma-1)V(\hat{U}_{t},\hat{W}_{t})$
$\displaystyle-||\hat{U}_{t+1}||_{2}^{2}\Delta x$
$\displaystyle+\nabla_{U}V^{T}(\hat{U}_{t},\hat{W}_{t})\dot{\hat{U}}_{t}\Delta
t$ $\displaystyle+\frac{1}{R}\nabla_{W}V^{T}(\hat{U}_{t}-\hat{W}_{t})\Delta t$
$\displaystyle+\max(0,-\nabla_{W}V^{T}B\hat{W}))^{2}\Delta t$ (13)

The second loss function provides an initial condition. At
$u(x,T)=0,w(x,T)=-R(x,t)=-2$, we have: $u(x,T)=0$ and $u_{t}(x,T)=0$. As a
result, we have $L(0,-R(x,t))=0$ and $L_{t}(0,-R(x,t))=0$. This shows us that
$u(x,T)=0,w(x,T)=-R(x,t)$ is considered a stable point that maximizes the
reward. Thus, we choose to let $V(0,-R(x,t))=0$ be the Dirichlet boundary
condition for the HJB equation. This leads to the second loss function:

$MSE_{u}=(V(0,-R(x,t)))^{2}.$ (14)

Since the value function achieves its global maximum at $u(x,T)=0,w(x,T)=-2$,
this means the derivatives of $V$ must be zero along all directions. Thus, we
choose to let $\frac{\partial V}{\partial n}=0$ at $u(x,T)=0,w(x,T)=-R(x,t)$
along every normal be the Neumann boundary condition for the HJB equation.
This leads to the third loss function:

$MSE_{n}=||\nabla_{U}V(0,-R(x,t))||_{2}^{2}+||\nabla_{W}V(0,-R(x,t))||_{2}^{2}$
(15)

We derived an optimal controller in corollary 4.2. Gym environments recommend
that actions be in the range $[-1,1]$. We can use the proof of the optimal
controller in Appendix D to derive a way of selecting actions:

$a_{t}=-\text{sign}(\nabla_{W}V^{T}B\hat{W}(t))$ (16)

The algorithms introduced in this paper will focus on minimizing both of the
loss functions defined above and using the optimal controller.

###  5.1 HJB value iteration

The HJB Value Iteration trains the loss function without the need for an
actor-network. We treat the value network as a PINN, using auto-
differentiation to estimate gradient vectors to compute the loss in equation 5
and the control in equation 16. At every time step, it uses the controller
given in equation 16. It updates the value network using the loss functions as
shown above. The method is provided in Algorithm 1.

Algorithm 1 HJB Value Iteration

1: Initiate value network parameter $\phi$

2: Run the control as given in equation (16) in the environment for $T$
timesteps and observe samples $\\{(s_{t},a_{t},R_{t},s_{t+1})\\}_{t=1}^{T}$.

3: Compute the value network loss as: $J(\phi)=MSE_{f}+MSE_{u}+MSE_{n}$
described in equations (5), (14), and (15)

4: Update $\phi\leftarrow\phi-\alpha_{2}\nabla_{\phi}J(\phi)$

5: Run steps 2–4 for multiple iterations

###  5.2 HJBPPO

HJBPPO is an algorithm that combines policy optimization from PPO with HJB
value iteration. This is implemented by modifying the PPO implementation by
Barhate (2021).

To facilitate exploration of the environment and exploitation of the models,
we introduce an action selection method that uses the policy network and
equation 16 with equal probability, as shown in Algorithm 2. Upon running the
policy $\pi_{\theta}$, we sample from a distribution $N(\mu,s)$ where $\mu$ is
the output from the policy network. We initiate $s$ to $0.3$ and decrease it
by $0.01$ every $1000$ episodes until it reaches $0.1$. After sampling an
action from the normal distribution, we clip it between $-1$ and $1$.

Algorithm 2 HJBPPO action selection

1: Retrieve state $s_{t}$, policy network parameter $\theta$ and value network
parameter $\phi$

2: Sample $i\in\\{0,1\\}$

3: if $i=0$ then  Select controller based on equation (16)

4: else Run policy $\pi_{\theta}$

5: end

This action selection method ensures that we select actions that are not only
in $\\{-1,1\\}$ but also in $[-1,1]$. It introduces a new method of
exploration of the environment by choosing from two different methods of
action selection. Actions selected using equation 16 are also stored in the
memory buffer and are used to train the policy network $\pi_{\theta}$. The
method is provided in Algorithm 3.

Algorithm 3 HJBPPO

1: Initiate policy network parameter $\theta$ and value network parameter
$\phi$

2: Run action selection as given in algorithm 2 in the environment for $T$
timesteps and observe samples $\\{(s_{t},a_{t},R_{t},s_{t+1})\\}_{t=1}^{T}$.

3: Compute the advantage $A_{t}$

4: Compute
$r_{t}(\theta)=\frac{\pi_{\theta}(a_{t}|s_{t})}{\pi_{\theta_{\text{old}}}(a_{t}|s_{t})}$

5: Compute the objective function of the policy network:  $\displaystyle
L(\theta)=\frac{1}{T}\sum_{t=0}^{T-1}$
$\displaystyle\min[r_{t}(\theta)A_{t},\text{clip}(r_{t}(\theta),1-\epsilon,1+\epsilon)A_{t}],$

6: Update $\theta\leftarrow\theta+\alpha_{1}\nabla_{\theta}L(\theta)$

7: Compute the value network loss as: $J(\phi)=MSE_{f}+MSE_{u}+MSE_{n}$
described in equations (5), (14), and (15)

8: Update $\phi\leftarrow\phi-\alpha_{2}\nabla_{\phi}J(\phi)$

9: Run steps 2–8 for multiple iterations

We will train PPO, HJB value iteration, and HJBPPO on the PackCooling
environment and compare these algorithms.

##  6 Results

###  6.1 Training

To ensure the reproducibility of our results, we have posted our code in the
following link:

https://github.com/amartyamukherjee/PPO-PackCooling. We posted our
hyperparameters in Appendix E. The details of the implementation of the
PackCooling gym environment are posted in Appendix A. The code was run using
Kaggle CPUs. Each algorithm was trained for a million timesteps. Training each
algorithm took approximately 5 hours.

###  6.2 Reward Curves

Figure 1: Reward curves of PPO (red), HJB value iteration (blue), and HJBPPO
(green) averaged over 5 seeds. Shaded area indicates 0.2 standard deviations.

The reward curves have been plotted in Figure 1, comparing PPO, HJB value
iteration, and HJBPPO. Each algorithm was run for 5 different seeds. We
plotted the mean reward over each seed and over 20 consecutive episodes, and
shaded the area 0.2 standard deviations from the mean. HJB value iteration
shows the worst performance, as its rewards decrease past PPO after training
for multiple episodes. PPO shows a rapid increase in average rewards after the
first episode and a slow increase in average rewards afterward. HJBPPO shows
the best performance in the graph, achieving the highest average reward in
each episode and an increase in average rewards after training for multiple
episodes.

The significantly higher average reward in HJBPPO in the first episode shows
that the action selection method described in Algorithm 2 provides a robust
strategy to explore the environment and train the models. The higher average
rewards are due to the exploitation of the dynamics of the environment as done
by the HJB equation.

###  6.3 Trajectories

The plots of the trajectories have been posted in Appendix F. After training
for a million timesteps, we tested our models on the PackCooling environment
and produced the plots. These plots were generated using the rendering feature
explained in section A.4.

The trajectory of HJB value iteration shows the worst results. $\sigma(t)$
returns $1.0$ only once. It achieves a cumulative reward of $-7294.51$. Thus,
the input of the cooling fluid from the boundary is minimal. As a result of
the internal heat generation in the battery pack, $u(x,t)$ reaches high values
of roughly $5$ at $t=10$, and as a result, $w(x,t)$ also reaches high values
of roughly $4$. This shows that the training of the value function in HJB
value iteration is inadequate and we have not arrived at an optimal controller
for the pack cooling problem. This is because exploration of the environment
was at a minimum, as we only exploited equation 16 at each time step.

The trajectory of PPO shows that the values $\sigma(t)$ takes at every
timestep have a large variance with its neighboring timesteps. It achieves a
cumulative reward of $-3970.02$. Control of the temperature of the battery
pack has been achieved as $u(x,t)$ takes values between $-2$ and $2$ at
$t=10$.

The trajectory of u(x,t) with HJBPPO shows it takes values between $-2$ and
$2$ at $t=10$. The values $\sigma(t)$ takes at every timestep have a lower
variance with its neighboring timesteps compared to PPO. It achieves a
cumulative reward of $-881.55$. For $t\in[4,6]$, $u(x,t)$ shows an increasing
trend towards $u=2$. In response, the controller $\sigma(t)$ took values
closer to $1.0$ to allow for greater input of cooling fluid from the boundary
so that $u(x,t)$ decreases towards zero. Due to higher average rewards as
shown in Figure 1, this shows that a model that exploits the dynamics of the
environment to return a controller shows improved performance compared to a
model that returns noisy control centered at $\sigma=0.5$.

##  7 Conclusion

In this paper, we have introduced two algorithms that use PINNs to solve the
pack cooling problem. This paper combines PINNs with RL in a PDE control
setting. In the HJB value iteration algorithm, the HJB equation is used to
introduce a loss function and a controller using a value network. The HJBPPO
algorithm is a hybrid-policy model that combines the training of the value
network from HJB value iteration and the training of the policy network from
PPO. HJBPPO shows an overall improvement in performance compared to PPO due to
its ability to exploit the physics of the environment to improve the learning
curve of the agent.

##  8 Future research

Despite showing an overall improvement in the reward curves, the HJBPPO
algorithm leaves room for improved RL algorithms using PINNs.

In this paper, we computed the HJB equation by expressing the PDE as an ODE by
discretizing in $x$. This was possible because the pack cooling problem was
modeled by 1D PDEs. Currently existing works such as (Sirignano & Spiliopoulos
(2018)) and (Kalise & Kunisch (2017)) solve the HJB equation for 1D PDEs by
discretizing it in $x$. It will be interesting to see how HJB control can be
extended to higher dimensional PDEs.

The goal of PINNs is to solve PDEs without the need for numerical methods. In
this paper, we solved the pack cooling problem numerically using the Crank-
Nicolson method and the method of characteristics. An area for further
research may be the use of PINNs to solve for the HJB equation and the PDE
that governs the dynamics of the system.

In the PackCooling environment, the HJBPPO algorithm showed an improvement
compared to PPO. But this is due to the fact that we knew the dynamics of the
system, thus allowing for the physics of the environment to be exploited. The
environments give all the details of the state needed to choose an action. One
limitation of HJBPPO is that it may not perform well in partially observable
environments because the estimate of the dynamics of the system may be
inaccurate. Deep Transformer Q Network (DTQN) was introduced by Esslinger et
al. (2022) and achieves state-of-the-art results in many partially observable
environments. A potential area for further research may be the introduction of
an HJB equation that facilitates partial observability. The DTQN algorithm may
be improvised by incorporating this HJB equation using PINNs.

## References

  * Antonelo et al. (2021) Antonelo, E. A., Camponogara, E., Seman, L. O., de Souza, E. R., Jordanou, J. P., and Hubner, J. F.  Physics-informed neural nets for control of dynamical systems.  _arXiv preprint arXiv:2104.02556_ , 2021.
  * Barhate (2021) Barhate, N.  Minimal pytorch implementation of proximal policy optimization.  https://github.com/nikhilbarhate99/PPO-PyTorch, 2021.
  * Djeumou et al. (2022) Djeumou, F., Neary, C., Goubault, E., Putot, S., and Topcu, U.  Neural networks with physics-informed architectures and constraints for dynamical systems modeling.  In _Learning for Dynamics and Control Conference_ , pp. 263–277. PMLR, 2022.
  * Esslinger et al. (2022) Esslinger, K., Platt, R., and Amato, C.  Deep transformer q-networks for partially observable reinforcement learning.  _Preprint arXiv:2206.01078_ , 2022.
  * Farahmand et al. (2017) Farahmand, A.-m., Nabi, S., and Nikovski, D. N.  Deep reinforcement learning for partial differential equation control.  In _2017 American Control Conference (ACC)_ , pp. 3120–3127, 2017\.  doi: 10.23919/ACC.2017.7963427.
  * Jiang et al. (2016) Jiang, F., Chou, G., Chen, M., and Tomlin, C. J.  Using neural networks for fast reachable set computations.  _Preprint arXiv:611.03158_ , 2016.
  * Kalise & Kunisch (2017) Kalise, D. and Kunisch, K.  Polynomial approximation of high-dimensional hamilton–jacobi–bellman equations and applications to feedback control of semilinear parabolic pdes.  _SIAM Journal on Scientific Computing_ , 40, 02 2017.  doi: 10.1137/17M1116635.
  * Kamalapurkar et al. (2018) Kamalapurkar, R., Rosenfeld, J., Dixon, W. E., and Walters, P.  _Reinforcement Learning for Optimal Feedback Control: A Lyapunov-Based Approach_.  Springer Cham, 2018.
  * Kato & Moura (2021) Kato, D. and Moura, S. J.  1d pde model for thermal dynamics in fluid-cooled battery packs: Numerical methods and sensor placement.  In _2021 American Control Conference (ACC)_ , pp. 3102–3107, 2021\.  doi: 10.23919/ACC50511.2021.9483248.
  * Kim et al. (2021) Kim, J., Shin, J., and Yang, I.  Hamilton–jacobi deep q-learning for deterministic continuous-time systems with lipschitz continuous controls.  _Journal of Machine Learning Research_ , 22:1–34, 2021\.
  * Liberzon (2012) Liberzon, D.  _Calculus of variations and optimal control theory: a concise introduction_.  Princeton University Press, 2012.
  * Munos (1999) Munos, R.  A Study of Reinforcement Learning in the Continuous Case by the Means of Viscosity Solutions.  _Machine Learning_ , 40:265–299, 1999.
  * Munos et al. (1999) Munos, R., Baird, L., and Moore, A.  Gradient descent approaches to neural-net-based solutions of the hamilton-jacobi-bellman equation.  In _IJCNN’99. International Joint Conference on Neural Networks. Proceedings (Cat. No.99CH36339)_ , volume 3, pp. 2152–2157 vol.3, 1999.
  * Nakamura-Zimmerer et al. (2020) Nakamura-Zimmerer, T., Gong, Q., and Kang, W.  A causality-free neural network method for high-dimensional hamilton-jacobi-bellman equations.  In _2020 American Control Conference (ACC)_ , pp. 787–793, 2020\.
  * Nicodemus et al. (2022) Nicodemus, J., Kneifl, J., Fehr, J., and Unger, B.  Physics-informed neural networks-based model predictive control for multi-link manipulators.  _IFAC-PapersOnLine_ , 55(20):331–336, 2022.  ISSN 2405-8963.  10th Vienna International Conference on Mathematical Modelling MATHMOD 2022.
  * Pineda (1987) Pineda, F.  Generalization of back propagation to recurrent and higher order neural networks.  In _Neural information processing systems_ , 1987.
  * Raffin et al. (2022) Raffin, A., Hill, A., Gleave, A., Kanervisto, A., Ernestus, M., and Dormann, N.  Reinforcement learning tips and tricks.  https://stable-baselines3.readthedocs.io/en/master/guide/rl_tips.html, 2022\.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019\.
  * Schulman et al. (2017) Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O.  Proximal policy optimization algorithms.  _Preprint arXiv:1707.06347_ , 2017.
  * Sirignano & Spiliopoulos (2018) Sirignano, J. and Spiliopoulos, K.  Dgm: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, 2018\.  ISSN 0021-9991.
  * Sutton & Barto (2018) Sutton, R. S. and Barto, A. G.  _Reinforcement Learning: An Introduction (2nd edition)_.  Cambridge, MA : The MIT Press, 2018.
  * Tassa & Erez (2007) Tassa, Y. and Erez, T.  Least squares solutions of the hjb equation with neural network value-function approximators.  _IEEE Transactions on Neural Networks_ , 18(4):1031–1041, 2007.
  * Vazquez & Krstic (2017) Vazquez, R. and Krstic, M.  Boundary control of coupled reaction-advection-diffusion systems with spatially-varying coefficients.  _IEEE Transactions on Automatic Control_ , 62(4):2026–2033, 2017.  doi: 10.1109/TAC.2016.2590506.
  * Wiltzer et al. (2022) Wiltzer, H. E., Meger, D., and Bellemare, M. G.  Distributional Hamilton-jacobi-Bellman equations for continuous-time reinforcement learning.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162, pp. 23832–23856. PMLR, 2022.

##  Appendix A The PackCooling gym environment

To run numerical experiments in this paper, we had to create a gym environment
that solves the pack cooling model numerically and takes a controller as input
every time step.

In gym environments, it is recommended to set the action range to $[-1,1]$
(Raffin et al. (2022)). In our control problem, $\sigma(t)$ returns values in
$[0,1]$. Thus, we let $\sigma(t)=\frac{a_{t}+1}{2}$.

###  A.1 Internal heat generation

In this paper, we set the internal heat generation in the battery pack,
$h(x,t,u(x,t))$ as:

$h(x,t,u(x,t))=e^{0.1u(x,t)}.$ (17)

The objective of the internal heat generation was to increase the temperature
distribution across the battery pack at every time step, which serves as a
motivation for choosing a strictly positive exponential function. This makes
the desired state of the controller, $u=0$, unstable.

###  A.2 Initial conditions

For $u(x,t)$ we chose initial conditions that satisfy the boundary conditions
3. A set of solutions that satisfy the boundary conditions are:

$u(x,0)=\sum_{n=0}^{\infty}C_{n}\cos(\pi nx),$ (18)

since the derivative of $\cos(\pi nx)$ is zero at $x=0$ and $x=1$ for all
integers $n$. Thus, we set our initial conditions to the following:

$u(x,0)=\sum_{n=0}^{N}C_{n}\cos(\pi nx),$ (19)

where $N$ is the number of Fourier nodes. Note that $N\leq N_{x}$ where
$N_{x}$ is the number of points we discretize the numerical solution into. In
our experiments, we set $N=9$. The Fourier coefficients $C_{n}$ are sampled
from the uniform distribution $U(-2,2)$. Adding randomness in our initial
conditions leads to more robustness in our controller.

For $w(x,t)$, we set the initial condition to $w(x,0)=U(0)$ where $U$ is
defined in the boundary condition in equation 4. In our experiments, we set
$U(t)$ to a constant: $U(t)=-5.0$.

This means if we set the controller to a constant $\sigma(t)=0$, then $u(x,t)$
will increase to high values due to the forcing term $h(x,t,u(x,t))$. And if
we set the controller to a constant $\sigma(t)=1$, then $u(x,t)$ will decrease
to low values due to the cooling fluid. So the objective of the controller is
to choose $\sigma(t)$ so that the temperature distribution in the battery pack
is as close to zero as possible.

###  A.3 Numerical solution

In our numerical solutions, we discretize the heat distributions $u(x,t)$ and
$w(x,t)$ in $x$ and $t$ as: $\Delta x=0.01,\Delta t=0.01$, thus satisfying
$\Delta x\geq\Delta t\geq\sigma(t)\Delta t$. We set $D(x,t)$ and $R(x,t)$ from
equations 1 and 2 to constants as $D(x,t)=0.01,R(x,t)=2.0$.

The numerical solution to the PDE was implemented as explained in Kato & Moura
(2021). We solve for a system

$A^{+}z(t+\Delta t)-\frac{1}{2}h(z(t+\Delta t))=Az(t)+\frac{1}{2}h(z(t)),$
(20)

where $z(t)$ is the concatenation of $u(x,t)$ and $w(x,t)$ discretized in
space. The matrices $A^{+}$ and $A$ are derived using a combination of the
Crank-Nicolson method, characteristics, and interpolation. Finally, to solve
this non-linear system, we run $10$ iterations of the Newton-Raphson method.

###  A.4 Rendering

In the rendering of the environment, we plot the trajectories of
$\sigma(t),u(x,t),w(x,t)$ in the environment. We store the numerical solutions
of the three variables in a buffer every episode and reset the buffer upon
resetting the environment.

Examples of rendered plots are shown in Appendix F. We return a line plot of
$\sigma(t)$ along $t$, and a mesh plot of $u(x,t)$ and $w(x,t)$ along $x$ and
$t$.

##  Appendix B Proof of theorem 2.1 part 1

Theorem A function $V(x)$ is the optimal value function if and only if $V\in
C^{1}(\mathbb{R}^{n})$ and $V$ satisfies the HJB Equation

$\displaystyle(\gamma-1)V(x)+\sup_{\sigma\in U}\\{L(x,\sigma)+\gamma\Delta
t\nabla_{x}V^{T}(x)f(x,\sigma)\\}=0$ (21)

for all $x\in\mathbb{R}^{n}$, assuming for all $x\in\mathbb{R}^{n}$, there
exists a controller $\sigma^{*}(\cdot)$ such that:

$\displaystyle(\gamma-1)V(x)+L(x,\sigma^{*}(x))+\gamma\Delta
t\nabla_{x}V^{T}(x)f(x,\sigma^{*}(x))$
$\displaystyle=(\gamma-1)V(x)+\sup_{\hat{\sigma}\in
U}\\{L(x,\hat{\sigma})+\gamma\Delta t\nabla_{x}V^{T}(x)f(x,\hat{\sigma})\\}.$
(22)

Proof We start with equation 6:

$V^{*}(x(t_{0}))=\sup_{\sigma}\frac{1}{\Delta
t}\int_{t_{0}}^{\infty}\gamma^{\frac{t}{\Delta
t}}L(x(\tau),\sigma(\tau))d\tau$ (23)

Separating the integral term gives:

$\displaystyle V^{*}(x(t_{0}))$
$\displaystyle=\sup_{\sigma}\left[\frac{1}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta
t}\gamma^{\frac{t}{\Delta t}}L(x(\tau),\sigma(\tau))d\tau+\frac{\gamma}{\Delta
t}\int_{t_{0}+\Delta t}^{\infty}\gamma^{\frac{t}{\Delta
t}}L(x(\tau),\sigma(\tau))d\tau\right]$ (24)
$\displaystyle=\sup_{\sigma}\left[\frac{1}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta
t}\gamma^{\frac{t}{\Delta t}}L(x(\tau),\sigma(\tau))d\tau+\gamma
V^{*}(x(t_{0}+\Delta t))\right]$ (25)

As $\Delta t\to 0$, we can rely on the following two assumptions:

$\frac{1}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta t}\gamma^{\frac{t}{\Delta
t}}L(x(\tau),\sigma(\tau))d\tau\approx\frac{1}{\Delta t}(\Delta
tL(x(t_{0}),\sigma(t_{0})))=L(x(t_{0}),\sigma(t_{0}))$ (26)
$V^{*}(x(t_{0}+\Delta t))\approx V^{*}(x(t_{0}))+\Delta
t\nabla_{x}V^{*}(x(t_{0}))^{T}f(x(t_{0}))$ (27)

This gives us the following equation:

$V^{*}(x(t_{0}))=\sup_{\sigma}\left[L(x(t_{0}),\sigma(t_{0}))+\gamma
V^{*}(x(t_{0}))+\gamma\Delta t\nabla_{x}V^{*}(x(t_{0}))^{T}f(x(t_{0}))\right]$
(28)
$0=(\gamma-1)V^{*}(x(t_{0}))+\sup_{\sigma}\left[L(x(t_{0}),\sigma(t_{0}))+\gamma\Delta
t\nabla_{x}V^{*}(x(t_{0}))^{T}f(x(t_{0}))\right]$ (29)

##  Appendix C Proof of theorem 4.1

Theorem Let $u(\cdot,t),w(\cdot,t)\in L_{2}[0,1]$. With $\sigma(t)\in[0,1]$
and the reward function $L(U_{t},W_{t},\sigma_{t})=-||U_{t+1}||_{2}^{2}\Delta
x$, the HJB equation for the 1D pack cooling problem is:

$\displaystyle(1-\gamma)V(u,w)-||u(\cdot,t+\Delta t)||^{2}+\langle
V_{u},u_{t}\rangle+\frac{1}{R}\langle V_{w},u-w\rangle+\max(0,-\langle
V_{w},w_{x}\rangle)=0$ (30)

where $||\cdot||$ is the $L_{2}[0,1]$ norm and $\langle\cdot,\cdot\rangle$ is
the $L_{2}[0,1]$ inner product.

Proof In the pack cooling problem, the value function is a function of
$\hat{U}(t)$ and $\hat{W}(t)$, which is $u(x,t)$ and $w(x,t)$ discretized in
$x$. The HJB equation from equation 7 can thus be written as follows:

$(\gamma-1)V(\hat{U}(t),\hat{W}(t))+\sup_{\sigma\in[0,1]}\\{L(\hat{U}(t),\hat{W}(t),\sigma)+\gamma\Delta
t\nabla_{\hat{U}}V^{T}(\hat{U}(t),\hat{W}(t))\dot{\hat{U}}+\gamma\Delta
t\nabla_{\hat{W}}V^{T}(\hat{U}(t),\hat{W}(t))\dot{\hat{W}}\\}=0$ (31)

We know that $L(\hat{U}(t),\hat{W}(t),\sigma)=-||\hat{U}(t+\Delta
t)||_{2}^{2}\Delta x$, where $\Delta x$ and $\Delta t$ are the discretizations
in $x$ and $t$ used in the environment respectively. $\hat{U}(t+\Delta t)$ is
determined by $\hat{U}(t)$ and $\dot{\hat{U}}(t)$. We know from equation 9
that neither $\hat{U}(t)$ or $\dot{\hat{U}}(t)$ depend on $\sigma(t)$. Thus,
we can bring $L(\hat{U}(t),\hat{W}(t),\sigma)$ and
$\nabla_{\hat{U}}V^{T}(\hat{U}(t),\hat{W}(t))\dot{\hat{U}}$ outside the
supremum.

$(\gamma-1)V(\hat{U}(t),\hat{W}(t))+L(\hat{U}(t),\hat{W}(t),\sigma)+\gamma\Delta
t\nabla_{\hat{U}}V^{T}(\hat{U}(t),\hat{W}(t))\dot{\hat{U}}+\sup_{\sigma\in[0,1]}\\{\gamma\Delta
t\nabla_{\hat{W}}V^{T}(\hat{U}(t),\hat{W}(t))\dot{\hat{W}}\\}=0$ (32)

We expand $\dot{\hat{W}}$ from equation 10.

$\displaystyle\sup_{\sigma\in[0,1]}\\{\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})\dot{\hat{W}}\\}$
$\displaystyle=\sup_{\sigma\in[0,1]}\\{\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})(-\sigma(t)B\hat{W}+\frac{1}{R}(\hat{U}-\hat{W}))\\}$
(33)
$\displaystyle=\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})(\frac{1}{R}(\hat{U}-\hat{W}))+\sup_{\sigma\in[0,1]}\\{\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})(-\sigma(t)B\hat{W})\\}$
(34)

For the expression to achieve its supremum, we let $\sigma(t)$ be defined as
the following:

$\sigma(t)=\begin{cases}1&\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})B\hat{W}<0\\\
0&\text{otherwise}\end{cases}$ (35)

We get:

$\sup_{\sigma\in[0,1]}\\{\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})(-\sigma(t)B\hat{W})\\}=\max(0,-\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})B\hat{W})$
(36)

Thus, we get the following HJB equation in discretized $x$:

$(\gamma-1)V(\hat{U},\hat{W})+L(\hat{U},\hat{W},\sigma)+\gamma\Delta
t\nabla_{\hat{U}}V^{T}(\hat{U},\hat{W})\dot{\hat{U}}+\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})(\frac{1}{R}(\hat{U}-\hat{W}))+\max(0,-\gamma\Delta
t\nabla_{\hat{W}}V^{T}(\hat{U},\hat{W})B\hat{W})=0$ (37)

Let $\Delta x=\Delta t$ to satisfy the condition for numerical stability. As
we let $\Delta x\to 0$, we see a Riemann sum forming.

$(\gamma-1)V-\int_{0}^{1}(u(\cdot,t+\Delta
t))^{2}dx+\int_{0}^{1}V_{u}u_{t}dx+\int_{0}^{1}V_{w}(\frac{1}{R}(u-w))dx+\max(0,-\int_{0}^{1}V_{w}w_{x}dx)=0$
(38)

We can write this as inner products and norms in the $L_{2}[0,1]$ space.

$(\gamma-1)V-||u(\cdot,t+\Delta t)||^{2}+\langle V_{u},u_{t}\rangle+\langle
V_{w},\frac{1}{R}(u-w)\rangle+\max(0,-\langle V_{w},w_{x}\rangle)=0$ (39)

##  Appendix D Proof of corollary 4.2

Corollary Let $u(\cdot,t),w(\cdot,t)\in L_{2}[0,1]$. With $\sigma(t)\in[0,1]$
and the reward function $L(U_{t},W_{t},\sigma_{t})=-||U_{t+1}||_{2}^{2}\Delta
x$, provided the optimal value function $V^{*}(u,w)$ the optimal controller
for the 1D pack cooling problem is:

$\sigma^{*}(t)=\begin{cases}1,&\langle
V^{*}_{w}(\cdot,t),w_{x}(\cdot,t)\rangle<0,\\\
0,&\text{otherwise},\end{cases}$ (40)

where $||\cdot||$ is the $L_{2}[0,1]$ norm and $\langle\cdot,\cdot\rangle$ is
the $L_{2}[0,1]$ inner product.

Proof The optimal controller was derived in the proof of the HJB equation in
equation 35. As we let $\Delta x\to 0$, we get:

$\sigma(t)=\begin{cases}1&\int_{0}^{1}V_{w}w_{x}dx<0\\\
0&\text{otherwise}\end{cases}$ (41)

We can write this as inner products and norms in the $L_{2}[0,1]$ space.

$\sigma(t)=\begin{cases}1&\langle V_{w},w_{x}\rangle<0\\\
0&\text{otherwise}\end{cases}$ (42)

##  Appendix E Hyperparameters

Hyperparameter | Value
---|---
Horizon (T) | 1024
Actor learning rate | 3e-04
Critic learning rate | 1e-03
Num. epochs | 10
Minibatch size | 64
Discount ($\gamma$) | 0.99
GAE parameter ($\lambda$) | 0.95
Table 1: Hyperparameters

##  Appendix F Plots of trajectories

Figure 2: Trajectory of PPO. Cumulative reward: $-3970.02$ Figure 3:
Trajectory of HJB value iteration. Cumulative reward: $-7294.51$ Figure 4:
Trajectory of HJBPPO. Cummulative reward: $-881.55$


# Compact Binary Merger Rate in Dark-Matter Spikes

Saeed Fakhry<EMAIL_ADDRESS>Department of Physics, Shahid Beheshti
University, Evin, Tehran 19839, Iran  PDAT Laboratory, Department of Physics,
K.N. Toosi University of Technology, P.O. Box 15875-4416, Tehran, Iran  Zahra
Salehnia<EMAIL_ADDRESS>Department of Physics, K.N. Toosi
University of Technology, P.O. Box 15875-4416, Tehran, Iran  PDAT Laboratory,
Department of Physics, K.N. Toosi University of Technology, P.O. Box
15875-4416, Tehran, Iran  Azin Shirmohammadi<EMAIL_ADDRESS>Department of Physics, K.N. Toosi University of Technology, P.O. Box
15875-4416, Tehran, Iran  PDAT Laboratory, Department of Physics, K.N. Toosi
University of Technology, P.O. Box 15875-4416, Tehran, Iran  Mina Ghodsi
Yengejeh<EMAIL_ADDRESS>Department of Physics, Shahid Beheshti
University, Evin, Tehran 19839, Iran  PDAT Laboratory, Department of Physics,
K.N. Toosi University of Technology, P.O. Box 15875-4416, Tehran, Iran  Javad
T. Firouzjaee<EMAIL_ADDRESS>Department of Physics, K.N. Toosi
University of Technology, P.O. Box 15875-4416, Tehran, Iran  School of
Physics, Institute for Research in Fundamental Sciences (IPM), P.O. Box
19395-5531, Tehran, Iran  PDAT Laboratory, Department of Physics, K.N. Toosi
University of Technology, P.O. Box 15875-4416, Tehran, Iran

###### Abstract

Nowadays, the existence of supermassive black holes (SMBHs) in the center of
galactic halos is almost confirmed. An extremely dense region referred to as
dark-matter spike is expected to form around central SMBHs as they grow and
evolve adiabatically. In this work, we calculate the merger rate of compact
binaries in dark-matter spikes while considering halo models with spherical
and ellipsoidal collapses. Our findings exhibit that ellipsoidal-collapse dark
matter halo models can potentially yield the enhancement of the merger rate of
compact binaries. Finally, our results confirm that the merger rate of
primordial black hole binaries is consistent with the results estimated by the
LIGO-Virgo detectors, while such results can not be realized for primordial
black hole-neutron star binaries.

Dark-Matter Spike — Primordial Black Hole — Neutron Star — Ellipsoidal
Collapse

##  1 Introduction

Over the last few decades, gravitational waves (GWs) have been studied as an
interesting cosmological observable. Regarding this, a large part of
astrophysical and cosmological phenomena have been evaluated through the study
of GWs and their direct detections. Meanwhile, compact binary mergers are
always considered potential sources for the propagation of GWs (Mandel &
Broekgaarden, 2022). During recent years, dozens of compact binary mergers
have been recorded by the LIGO-Virgo detectors (e.g., Abbott et al. (2016a, b,
c, 2020a, 2020b)). In this regard, three main classes of compact binaries as
binary black holes (BBHs), black hole-neutron star (BH-NS) binaries, and
binary neutron stars (BNSs) can potentially be captured by the LIGO-Virgo
detectors. Meanwhile, most of the recorded GW events belong to BBH merger
events in the mass range $10\,M_{\odot}\mbox{-}\,100\,M_{\odot}$ (Abbott et
al., 2019, 2021a, 2021b). There have been several investigations into the
origin of such BHs (e.g., Raidal et al. (2019); Bouffanais et al. (2021); Nitz
& Wang (2021a, b); Mandel & Farmer (2022)). However, our information on this
matter is still insufficient. There is a possibility that they have formed
during stellar collapse (likely via different channels) or that they have
formed at the beginning of the Universe as a result of gravitational collapse
of density fluctuations.

Surprisingly, most of the BBH merger events reported by the LIGO-Virgo
collaboration are consistent with the primordial black hole (PBH) scenario. In
the cosmological perturbation theory, PBHs are predicted to form due to
nonlinear peaks in initial density fluctuations during horizon reentry (see,
e.g., Zel’dovich & Novikov (1967); Hawking (1971); Carr & Hawking (1974)).
Several studies show that the existence of a critical state in primordial
density fluctuations seems necessary for the formation of PBHs, which is
achieved through exceeding a threshold value (e.g., Shibata & Sasaki (1999);
Polnarev & Musco (2007); Musco & Miller (2013); Young et al. (2014);
Bloomfield et al. (2015); Allahyari et al. (2017)). According to this
argument, exceeding a certain threshold value is equivalent to the direct
collapse of primordial density fluctuations and consequently the formation of
PBHs. In addition, PBHs are characterized by their broad range of masses,
which makes them distinct from astrophysical black holes (Sasaki et al.,
2018). Meanwhile, the standard model of cosmology tries to explain the nature
of two basic components of the dark sector of the Universe: dark matter and
dark energy. In this regard, dark energy is considered a type of energy that
governs the accelerating expansion of the late-time Universe, which
corresponds to the cosmological constant in general relativity (see, e.g.,
Sherwin et al. (2011); Yang & Xu (2014); Huterer & Shafer (2018); Di Valentino
et al. (2020); Ghodsi Y. et al. (2022); Ghodsi Yengejeh et al. (2023)). Dark
matter is also attributed to an invisible matter that is thought to make up
approximately $85\%$ of the matter in the Universe. Nowadays, it is believed
that PBHs can be regarded as one of the potential candidates for dark matter
(e.g., Bird et al. (2016); Sasaki et al. (2016); Kashlinsky (2016); Blinnikov
et al. (2016); Clesse & García-Bellido (2017); Bellomo et al. (2018);
Villanueva-Domingo et al. (2021); Carr & Kuhnel (2021)). However, other
candidates for dark matter are being discussed seriously (e.g., Spergel &
Steinhardt (2000); Alcock et al. (2000); Liebling & Palenzuela (2012);
Roszkowski et al. (2018); Braine et al. (2020); Di Giovanni et al. (2020,
2021)). Efficient observational methods have been used to constrain the
abundance of PBHs in different mass ranges, which in turn are compelling
contexts for studying the early Universe at small scales (e.g., Carr et al.
(2017); Lehmann et al. (2018); Carr et al. (2021)). However, it is possible to
obtain strong constraints on the contribution of PBHs to dark matter by
calculating their merger rate and validating the results via GW data (Bird et
al., 2016; Clesse & García-Bellido, 2017).

Furthermore, the random distribution of BHs in the Universe allows them to
form binaries with NSs. Interestingly, BH-NS mergers can provide important
information about multi-messenger astronomy (Ruiz et al., 2021). In addition
to emitting of GWs, they can also transmit electromagnetic signals during the
merger process (Barbieri et al., 2020). In such events, residual matter from
the NS is usually accreted by the BH and results in a luminous event
(Fragione, 2021). Moreover, recording and processing data from BH-NS binary
mergers by GW detectors can include individual information about the NS
nuclear equation of state and accretion processes of BHs, as well as
constraining their spin and abundance (see, e.g., Hinderer et al. (2019);
Zappa et al. (2019); Fragione et al. (2021); Tiwari et al. (2021)). In recent
years, two merger events have been captured by the LIGO-Virgo detectors, which
are attributed to BH-NS binaries and their component masses are
($8.9^{+1.2}_{-1.5}M_{\odot},1.9^{+0.3}_{-0.2}M_{\odot}$) and
($5.7^{+1.8}_{-2.1}M_{\odot},1.5^{+0.7}_{-0.3}M_{\odot}$), respectively Abbott
et al. (2021c). There are many uncertainties surrounding the formation of BH-
NS binaries and their merging. Nevertheless, describing this class of merger
events using the evolution of the field binaries can be a viable approach. It
is expected that a significant number of BH-NS merger events can be detected
by GW detectors in the upcoming years. Due to these circumstances, it seems
imperative to fully understand how compact objects involved in these events
are formed. In Sasaki et al. (2022), taking into account spherical-collapse
dark matter halo models and the framework of the PBH scenario, the merger rate
of BH-NS binaries has been calculated. Also, by comparing the results with the
GW observations, it has been claimed that the BH components participating in
such events can not have primordial origins. However, by considering more
realistic halo models (e.g., those with ellipsoidal-collapse), a more accurate
picture of galactic halos can be obtained (e.g., Fakhry et al. (2021, 2022a,
2022b); Fakhry & Del Popolo (2022)). With this argument, in Fakhry et al.
(2022c), we have indicated that by considering ellipsoidal-collapse dark
matter halo models and validating the results with relevant events estimated
by the LIGO-Virgo detectors, the BHs contributing to the BH-NS events are most
likely PBHs.

On the other hand, there is convincing evidence of the presence of
supermassive black holes (SMBHs) in the center of galactic halos (see, e.g.,
Volonteri et al. (2021); Prokhorenko & Sazonov (2021); Shapiro & Heggie
(2022); Diana et al. (2022)). This can be deduced from the Keplerian behavior
of the velocity dispersion of the stars in the inner regions of galactic halos
(Ghez et al., 1998). It is believed that central SMBHs are capable of
amplifying the density of surrounding dark matter particles to a certain
extent. In this regard, in Gondolo & Silk (1999), it has been proposed that a
dense region known as the dark-matter spike will form around the central SMBH
if it evolves adiabatically and initially has a power-law cusp. Given that the
dark matter density is very high in the spike regions, it is expected that the
number density of PBHs is also significant in such regions. It is noteworthy
that the structure of dark-matter spikes is determined by the growth of the
central SMBH and dark matter halo model. In Nishikawa et al. (2019), the
merger rate of PBH-PBH binaries has been calculated in dark-matter spikes
while accounting for spherical-collapse dark matter halo models. However, as
mentioned earlier, more realistic halo models (e.g., those with ellipsoidal
collapse) can potentially affect the final results.

In this work, we propose to employ the ellipsoidal-collapse dark matter halo
models to calculate the merger rate of compact binaries within dark-matter
spikes. In this regard, the outline of this work is as follows. In Sec. 2, we
introduce a convenient model for dark-matter spikes and discuss some crucial
quantities such as spike density profile, SMBH mass function, and
concentration parameter. Next, in Sec. 3, we calculate the merger rate of
compact binaries in the framework of ellipsoidal-collapse dark matter halo
models and compare it with that obtained from spherical-collapse dark matter
halo models. We also compare the relevant results of the present analysis with
the corresponding data estimated by the LIGO-Virgo detectors and constrain the
value of the power-law index. Finally, we scrutinize the results and summarize
the findings in Sec. 4.

##  2 Model of dark-matter spike

###  2.1 The density profile

According to the standard model of cosmology, dark matter halos are known as
nonlinear structures that have been formed hierarchically and distributed in
the Universe as a result of the collapse of linear cosmological fluctuations.
Based on indirect observations from the rotation curve of galaxies, it can be
argued that dark matter particles should not be uniformly distributed in
galactic halos. In this regard, special attention should be paid to the SMBHs
structured at the galactic center. Although ordinary black holes might be born
from a stellar collapse, SMBHs are difficult to accommodate in standard
astrophysical scenarios at high redshifts. It is suggested that SMBH mass can
be related to the mass of the dark matter halo, implying that SMBHs could
coevolve with their host halos (Volonteri et al., 2021; Bansal et al., 2022).
Accordingly, self-interacting dark matter halo models predict early seeds for
supermassive black holes through the gravothermal catastrophe (Choquette et
al., 2019; Feng et al., 2021, 2022). Consequently, the SMBH is expected to be
surrounded by a highly dense spike of dark matter at the center of the
galactic halo.

Assume $M_{\rm SMBH}$ is the mass of SMBH at the galactic center, which takes
a density profile of the form $\rho(r)\simeq\rho_{0}(r_{0}/r)^{\gamma}$, where
$\rho_{0}$ and $r_{0}$ are characteristic parameters of the halo, and $\gamma$
represents the power-law index. In light of this reasoning, one can expect to
form a dark-matter spike whose radius is specified by the following relation
(Gondolo & Silk, 1999)

$r_{\rm sp}=a_{\gamma}r_{0}\left(\frac{M_{\rm
SMBH}}{\rho_{0}r_{0}^{3}}\right)^{1/(3-\gamma)},$ (1)

where $a_{\gamma}$ is determined through numerical suggestions for each power-
law index $\gamma$.

Accordingly, the density profile of the dark-matter spike for $r$ in the range
of $4r_{\rm s}<r<r_{\rm sp}$ can be obtained as follows (Gondolo & Silk, 1999;
Nishikawa et al., 2019)

$\rho_{\rm sp}(r)=\rho_{0}\left(\frac{r_{0}}{r_{\rm
sp}}\right)^{\gamma}\left(1-\frac{4r_{\rm s}}{r}\right)^{3}\left(\frac{r_{\rm
sp}}{r}\right)^{\gamma_{\rm sp}},$ (2)

where $\gamma_{\rm sp}=(9-2\gamma)/(4-\gamma)$. Furthermore, $r_{\rm s}$
indicates the Schwarzchild radius of SMBH which is in the form of

$r_{\rm s}=\frac{2GM_{\rm SMBH}}{c^{2}}\simeq 2.95\,{\rm
km}\,\left(\frac{M_{\rm SMBH}}{M_{\odot}}\right),$ (3)

in which $G$ is the gravitational constant and $c$ is the velocity of light in
a vacuum. The results of numerical simulations and analytical approaches
exhibit that the density profile in small radii has a power-law form (Merritt
et al., 2006; Stadel et al., 2009; Navarro et al., 2010). However, there are
different predictions about the value of the power-law index $\gamma$. In this
work, we consider a range for the power-law index as $0<\gamma\leq 2$. On the
other hand, to describe dark matter distribution in galactic halos, a
convenient density profile was presented by Navarro, Frenk, and White (NFW)
(Navarro et al., 1996), which is specified by the following formula

$\rho_{\rm
NFW}(r)=\frac{\rho_{0}}{\left(r/r_{0}\right)\left(1+r/r_{0}\right)^{2}}.$ (4)

In Fig. 1, we have depicted the difference in the behavior of the density
profile corresponding to the dark-matter spike with the NFW density profile
while considering the mass of SMBH as $M_{\rm SMBH}=10^{6}M_{\odot}$. As can
be seen from the figure, the density profile in the regions specific to the
dark-matter spike is extremely higher than the NFW density profile. Therefore,
it seems interesting to calculate the merger rate of compact binaries in dark-
matter spikes. In addition, it should be noted that $r=r_{\rm sp}$ defines the
radius within which the merger rate of compact objects must be calculated as
dark-matter spike profiles cross the NFW density profile.

Figure 1: Comparison of the NFW density profile with the profile of dark-
matter spike with $\gamma=1$ and $2$ structured around the SMBH with a mass
$M_{\rm SMBH}=10^{6}M_{\odot}$.

###  2.2 The $M_{\rm SMBH}\mbox{-}\sigma$ relation

Nowadays, there are several lines of evidence indicating that the growth of
SMBHs and the evolution of their host halos are intricately linked (e.g.,
McLure & Dunlop (2002); Tremaine et al. (2002); Häring & Rix (2004); Graham et
al. (2007); Hopkins et al. (2007)). Regarding this, it is proposed that the
mass of SMBH can be strongly correlated with the velocity dispersion of dark
matter particles in a galactic halo, $\sigma$, which is called the $M_{\rm
SMBH}\mbox{-}\sigma$ relation. In other words, the halo characteristic
parameters, i.e., $\rho_{0}$ and $r_{0}$, can be related to $M_{\rm SMBH}$ by
employing such a relation. A convenient form of the $M_{\rm
SMBH}\mbox{-}\sigma$ relation has been obtained as (Ferrarese & Merritt, 2000;
Gebhardt et al., 2000)

$\log\left(\frac{M_{\rm
SMBH}}{M_{\odot}}\right)=a+b\log\left(\frac{\sigma}{200\,{\rm
kms^{-1}}}\right),$ (5)

where $a$ and $b$ are determined via empirical situation. In Gültekin et al.
(2009), values $a=8.12\pm 0.08$ and $b=4.24\pm 0.41$ have been suggested,
which fits reasonably well with the various types of galactic halos (Kormendy
& Ho, 2013). The NFW profile is assumed to manage the density profile outside
of spike regions, i.e. $r\gg r_{\rm sp}$, up to the virial radius $r_{\rm
vir}>r_{0}$. In fact, $r_{\rm vir}$ is attributed to a radius that includes a
volume in which the average halo density reaches $200$ to $500$ times the
critical density of the Universe. Under such assumptions, the total mass
enclosed by a volume of radius $r$ is determined as follows

$M(r)=4\pi\rho_{0}r_{0}\int_{0}^{r}\frac{rdr}{\left(1+r/r_{0}\right)^{2}}=4\pi\rho_{0}r_{0}^{3}g(r/r_{0}),$
(6)

where $g(x)=\log(1+x)-x/(1+x)$. It needs to be mentioned that the contribution
of dark-matter spike and central SMBH are negligible in comparison to the
total mass of dark matter halo. Moreover, the concentration parameter
determines the central density of dark matter halos, which is defined as
$C\equiv r_{\rm vir}/r_{0}$. Hence, the virial mass takes the following form

$M_{\rm vir}=4\pi\rho_{0}r_{0}^{3}g(C).$ (7)

Also, the circular velocity of dark matter particles reaches the maximum value
at a distance $r_{\rm m}=C_{\rm m}r_{0}=2.16\,r_{0}$, and corresponds to the
one-dimensional velocity dispersion of dark matter particles, namely

$\sigma^{2}=\frac{GM(C_{\rm m}r_{0})}{C_{\rm m}r_{0}}=4\pi
G\rho_{0}r_{0}^{2}\frac{g(C_{\rm m})}{C_{\rm m}}.$ (8)

As a result, a relation between $\rho_{0}$, $r_{0}$, and $M_{\rm SMBH}$ can be
established via Eqs. (5) and (8). Based on the results of $N$-body
simulations, the concentration parameter is a decreasing function of halo mass
and is a function of redshift at constant mass (e.g. Prada et al. (2012);
Dutton & Macciò (2014); Ludlow et al. (2016); Okoli & Afshordi (2016)), which
is consistent with the dynamics expected from the evolution of dark matter
haloes. In this work, to calculate the merger rate of compact binaries in the
present-time Universe, we utilize the concentration parameter presented in
Ludlow et al. (2016) for spherical-collapse dark matter halo models, and we
employ the corresponding one obtained in Okoli & Afshordi (2016) for
ellipsoidal-collapse dark matter halo models.

###  2.3 The mass function of SMBHs

One of the most fundamental challenges of extragalactic astronomy is to
understand how SMBHs grow and evolve. Accordingly, the SMBH mass function
provides comprehensive information on the mass of SMBHs and their evolution at
the center of galactic halos. Therefore, the SMBH mass function can be
considered a powerful and available tool to investigate the growth of SMBHs
and constrain related theoretical models. On the other hand, the SMBH mass
function might play a significant role in the structuring of upcoming surveys
because it provides an estimate of the mass classification of SMBHs Kelly &
Merloni (2012). It should be noted that obtaining an accurate mass function
for SMBHs is a relatively difficult task. For this reason, the current
estimates of the SMBH mass function include many theoretical uncertainties,
which in turn may affect the accuracy of calculating the merger rate of
compact binaries in dark-matter spikes. A reasonable approach for managing
these uncertainties is to compare the results from several different empirical
SMBH mass functions.

In Benson et al. (2007), by employing the Galactica code, a sample of $8839$
SDSS galaxies was employed to extrapolate the luminosity functions of spheroid
and disc galaxies, and a mass function of SMBHs was obtained as follows

$\phi(M_{\rm SMBH})=10^{9}\left(\frac{\phi_{0}M_{\rm
SMBH}^{\alpha}}{M_{*}^{\alpha+1}}\right)\exp\left[-\left(\frac{M_{\rm
SMBH}}{M_{*}}\right)^{\beta}\right],$ (9)

in which $\alpha=-0.65$, $\beta=0.6$, $\phi_{0}=2.9\times 10^{-3}\,h^{3}{\rm
Mpc^{-3}}$, and $M_{*}=4.07\times 10^{7}\,h^{-2}M_{\odot}$.

In addition, in Vika et al. (2009), a convenient mass function for SMBHs has
been obtained via the Millenium Galaxy Catalogue (Liske et al., 2003) for
$1743$ galaxies. This mass function is based on the experimental relation
between the mass SMBH and the luminosity of the host spheroid, which has the
following form

$\phi(M_{\rm SMBH})=\phi_{*}\left(\frac{M_{\rm
SMBH}}{M_{*}}\right)^{\alpha+1}\exp\left[1-\left(\frac{M_{\rm
SMBH}}{M_{*}}\right)\right],$ (10)

where $\log\phi_{*}=-3.15\,h^{3}{\rm Mpc^{-3}dex^{-1}}$, $\log(M_{*}/M)=8.71$,
and $\alpha=-1.20$. This mass function is valid for the masse range
$10^{6}M_{\odot}<M_{\rm SMBH}<10^{10}M_{\odot}$.

Also, in Shankar et al. (2004), another suitable mass function was derived for
SMBHs according to the observational relation between the SMBH mass and the
halo velocity dispersion and using kinematic and photometric data, which takes
the following formula

$\phi(M_{\rm SMBH})=\phi_{*}\left(\frac{M_{\rm
SMBH}}{M_{*}}\right)^{\alpha+1}\exp\left[1-\left(\frac{M_{\rm
SMBH}}{M_{*}}\right)^{\beta}\right],$ (11)

where $\phi_{*}=7.7\times 10^{-3}\,{\rm Mpc^{-3}}$, $M_{*}=6.4\times
10^{7}\,M_{\odot}$, $\beta=0.49$, and $\alpha=-1.11$. This mass function is
valid for the mass range $10^{6}M_{\odot}\leq M_{\rm SMBH}\leq 5\times
10^{9}M_{\odot}$.

##  3 Compact binary merger rate

Assume in dark-matter spike, a compact object with mass $m_{1}$ suddenly
encounters another one with mass $m_{2}$ on a hyperbolic orbit, and their
relative velocity at large separation is $v_{\rm rel}=|v_{1}-v_{2}|$. Hence,
based on two-body scattering, highly significant gravitational radiation emits
at the periastron $r_{\rm a}$. Keplerian mechanics states that such a system
is gravitationally bound when the emitted gravitational energy dominates the
kinetic energy of the system. Under these conditions, a maximum value for
periastron can be obtained as follows

$r_{\rm
mp}=\left[\frac{85\pi}{6\sqrt{2}}\frac{G^{7/2}m_{1}m_{2}(m_{1}+m_{2})^{3/2}}{c^{5}v_{\rm
rel}^{2}}\right]^{2/7}.$ (12)

Furthermore, in the Newtonian limit, the impact parameter is determined by the
periastron as follows:

$b^{2}(r_{\rm p})=\frac{2G(m_{1}+m_{2})r_{\rm p}}{v_{\rm rel}^{2}}+r_{\rm
p}^{2}.$ (13)

For the regions of dark-matter spikes that are gravitationally active, a
strong limit of gravitational focusing, i.e., $r_{\rm p}\ll b$, can be
considered in such a way that the relevant distortions of surrounding compact
objects on the formed binaries can be ignored. Thus, the cross-section for the
binary formation can be obtained via the following equation

$\xi(m_{1},m_{2},v_{\rm rel})=\pi b^{2}(r_{\rm mp})\simeq\frac{2\pi
G(m_{1}+m_{2})r_{\rm mp}}{v_{\rm rel}^{2}}.$ (14)

Hence, by Substituting Eq. (12) into Eq. (14), the cross section for the
binary formation can be derived as

$\xi\simeq
2\pi\left(\frac{85\pi}{6\sqrt{2}}\right)^{2/7}\frac{G^{2}(m_{1}+m_{2})^{10/7}(m_{1}m_{2})^{2/7}}{c^{10/7}v_{\rm
rel}^{18/7}}.$ (15)

Figure 2: The merger rate of compact binaries in a single dark-matter spike as
a function of SMBH mass for different values of $\gamma$. The top panels
demonstrate this relation for PBH-PBH events in spherical- and ellipsoidal-
collapse dark matter halo models, while the bottom panels display the
corresponding results for PBH-NS events.

Therefore, the merger rate of compact binaries within each region of the dark-
matter spike is determined as follows

$N_{\rm sp}=4\pi\int_{4r_{\rm s}}^{r_{\rm
sp}}\mathcal{T}(\rho,m_{1},m_{2})\langle\xi v_{\rm rel}\rangle\,r^{2}dr,$ (16)

where for the PBH-PBH events:

$\mathcal{T}=\left\\{\frac{1}{2}\frac{\left[f_{\rm PBH}\,\rho_{\rm
sp}(r)\right]^{2}}{m_{1}m_{2}}\right\\},$ (17)

and for the PBH-NS events:

$\mathcal{T}=\left(\frac{f_{\rm PBH}\,\rho_{\rm
sp}(r)}{m_{1}}\right)\left(\frac{\rho_{\rm NS}(r)}{m_{2}}\right).$ (18)

In the above relation, $0<f_{\rm PBH}\leq 1$ represents the fraction of PBHs
that specifies their contribution to dark matter, and the angle bracket
denotes an average over the relative velocity distribution at the vicinity of
the central SMBH. Furthermore, $\rho_{\rm NS}(r)$ is the NS density profile
that we define through the spherically symmetric form:

$\rho_{\rm NS}=\rho^{*}_{\rm NS}\exp\left(-\frac{r}{r^{*}_{\rm NS}}\right),$
(19)

where $r^{*}_{\rm NS}$ and $\rho^{*}_{\rm NS}$ are characteristic radius and
density of NSs, respectively, which need to be determined. For the
characteristic radius of NSs, an approximative value has been proposed as
$r^{*}_{\rm NS}\simeq 0.1\,r_{\rm s}$ (Sasaki et al., 2022), which we use in
our calculations.

Figure 3: The merger rate of PBH-PBH binaries per unit time and volume as a
function of SMBH mass for different values of $\gamma$. The top panels show
this relation for spherical-collapse dark matter halo models while considering
three different SMBH mass functions, whereas the bottom panels exhibit the
corresponding results for ellipsoidal-collapse dark matter halo models.

Furthermore, the characteristic density of NSs must be obtained by normalizing
the distribution of NSs to their estimated population in an arbitrary galaxy.
To accomplish this, we utilize the time-independent form of the initial
Salpeter stellar mass function, which is in the form $\chi(m_{*})\approx
m_{*}^{-2.35}$. Our main assumption is based on the fact that the entire
population of stars in the mass range of $8M_{\odot}\mbox{-}20M_{\odot}$ will
eventually yield a supernova explosion, and their outcome will be an NS.
Consequently, the number of NSs in a single galaxy with stellar mass $M_{*}$
is given by

$n_{\rm NS}=M_{*}\int_{m^{\rm min}_{*}}^{m^{\rm max}_{*}}\chi(m_{*})dm_{*},$
(20)

where $\chi(m_{*})m_{*}$ is normalized to unity. It should be noted that to
characterize the galactic stellar mass $M_{*}$, the stellar mass–halo mass
relation $M_{*}(M_{\rm halo})$ must be determined. For this purpose, the
stellar mass–halo mass relation obtained in Behroozi et al. (2013) can be
used, the basic assumption of which is the presence of the maximum number of
NSs at the center of the galactic halo. One must note that in the present
analysis, for the relative velocity near the central SMBH, we use the circular
velocity $v(r)=\sqrt{GM_{\rm SMBH}/r}$ at each radius bounded by the dark-
matter spike. This is a reasonable choice since the total mass enclosed by the
region of dark-matter spike is negligible versus the mass of the central SMBH.
We consider the mass of involving PBHs in PBH-PBH events as $M_{\rm
PBH}=30\,M_{\odot}$ and fix the masses of PBHs and NSs participating in PBH-NS
events as $M_{\rm PBH}=5\,M_{\odot}$ and $M_{\rm NS}=1.4\,M_{\odot}$. Also, in
the present analysis, we consider the contribution of PBHs in dark matter to
be $f_{\rm PBH}=1$. It is obvious from Eqs. (16), (17) and (18) that the
merger rate of PBH-PBH binaries is straightly proportional to the $f_{\rm
PBH}^{2}$, while it changes directly with $f_{\rm PBH}$ for the PBH-NS events.

In Fig. 2, we have plotted the merger rate of compact binaries within a single
spike as a function of SMBH mass for several values of power-law index
$\gamma$ while accounting for dark matter halo models with spherical and
ellipsoidal collapses. As can be seen from the figure, the merger rate of PBH-
PBH binaries changes inversely with the mass of the SMBH for both spherical-
and ellipsoidal-collapse dark matter halo models. It is also evident that the
merger rate of PBH-PBH binaries is directly proportional to the value of the
power-law index $\gamma$. However, comparing PBH-PBH events for spherical- and
ellipsoidal-collapse dark matter halo models, it can be inferred that the
merger rate per spike for ellipsoidal-collapse dark matter halo models is
higher than the corresponding one derived from spherical-collapse dark matter
halo models.

Figure 4: The merger rate of PBH-NS binaries per unit time and volume as a
function of SMBH mass for different values of $\gamma$. The top panels
indicate this relation for spherical-collapse dark matter halo models while
accounting for three different SMBH mass functions, whereas the bottom panels
show the corresponding results for ellipsoidal-collapse dark matter halo
models.

However, for the merger rate of PBH-NS binaries, the situation is slightly
different. Interestingly, unlike the previous case, the merger rate of PBH-NS
binaries per spike reaches the maximum value for the SMBH mass of $M_{\rm
SMBH}\simeq 10^{7}M_{\odot}$ for both spherical- and ellipsoidal-collapse dark
matter halo models. Also, it can be deduced from the calculation that the
merger rate of PBH-NS binaries per spike for a power-law index of $\gamma=1.7$
has a maximum value that is completely different from the relevant results of
PBH-PBH events. In other words, the merger rate of PBH-NS binaries increases
monotonically up to $\gamma=1.7$ and decreases for the larger values. In the
following, we will discuss this in more detail. In addition, it can be
inferred from the results that the merger rate of PBH-NS binaries for
ellipsoidal-collapse dark matter halo models is much higher than the
corresponding one derived from spherical-collapse dark matter halo models.

On the other hand, the cumulative merger rate of compact binaries is
considered the main quantity to be recorded and processed through the LIGO-
Virgo detectors. Therefore, the overall merger rate of compact binaries per
unit volume and per unit time needs to be specified. To perform this task, one
has to convolve the mas function of SMBH, $\phi(M_{\rm SMBH})$, with the
merger rate of compact binaries per spike, $N_{\rm sp}(M_{\rm SMBH})$:

$\mathcal{R}=\int_{M_{\rm min}}^{M_{\rm max}}N_{\rm sp}(M_{\rm
SMBH})\phi(M_{\rm SMBH})dM_{\rm SMBH}.$ (21)

According to Eqs. (9), (10) and (11), the mentioned mass functions have a
decreasing exponential term with respect to the mass of SMBHs. Therefore, it
can be concluded that $M_{\rm max}$ does not have a significant effect on the
final result. In contrast, the introduced mass functions indicate that the
maximum abundance belongs to the smallest central black holes in the Universe.
Consequently, $M_{\rm min}$ can have a significant contribution to the merger
rate of compact binaries in dark-matter spikes.

In Fig. 3, we have depicted the merger rate of PBH-PBH binaries per unit time
and volume as a function of SMBH mass for several values of power-law index
$\gamma$ while considering dark matter halo models with spherical and
ellipsoidal collapses. We have provided the results for three mass functions
Shankar et al. (2004), Benson et al. (2007), and Vika et al. (2009) to justify
the possible uncertainties in the present analysis as much as possible. The
merger rate of PBH-PBH binaries in both dark matter halo models with spherical
and ellipsoidal collapses decreases monotonically with increasing the mass of
SMBHs. Given the classification of the mass functions of SMBHs for their
abundance, this result seems reasonable. As it is clear from the figures, the
merger rate of PBH-PBH binaries for all three mass functions and both dark
matter halo models reaches the maximum value as $M_{\rm
SMBH}=10^{6}M_{\odot}$. Also, the direct proportionality of the merger rate of
PBH-PBH binaries to the values of the power-law index is evident. Moreover,
the merger rate of PBH-PBH binaries in ellipsoidal-collapse dark matter halo
models is slightly higher than that obtained from spherical-collapse dark
matter halo models. In the best case, which can be realized at the minimum
value of the power-law index, e.g. $\gamma=0.05$, the amplification of the
overall merger rate is $62\%$. Additionally, it can be concluded that the
merger rate of PBH-PBH binaries yields the highest, middle, and lowest values
while considering Shankar et al. (2004), Vika et al. (2009), and Benson et al.
(2007) mass functions respectively.

Figure 5: The overall merger rate of compact binaries as a function of power-
law index $\gamma$ for both spherical- and ellipsoidal-collapse dark matter
halo models. The top panels demonstrate this relation for PBH-PBH events for
three different SMBH mass functions, while the bottom panels show the
corresponding results for PBH-NS events. The shaded cyan bands represent the
BH-BH mergers estimated by the LIGO-Virgo detectors during the latest
observing run, i.e., $(17.9\mbox{-}44)\,{\rm Gpc^{-3}Yr^{-1}}$.

Similarly to the merger rate per spike discussed earlier, Fig. 4 exhibits that
the overall merger rate of PBH-NS binaries differs from that of PBH-PBH
binaries. As a common point of whole applied dark matter models and mass
functions, it is deduced that the merger rate of PBH-NS binaries has a maximum
at $\gamma=1.7$ and experiences a decreasing behavior for values greater than
that (see Fig. 5). It is also evident that the merger rate of PBH-NS binaries
for ellipsoidal-collapse dark matter halo models is in turn higher than that
obtained for spherical-collapse dark matter halo models. Such a relative
advantage could stem from the fact that, one can hope to improve the
theoretical predictions of the merger rate of compact binaries in galactic
halos by using more realistic dark matter halo models. Interestingly, for PBH-
NS binaries, the inverse proportionality of the overall merger rate with the
mass of SMBHs is not monotonic and has a plateau in some cases. If spherical-
collapse dark matter halo models are reliable, the results obtained for
Shankar et al. (2004) and Vika et al. (2009) mass functions for $\gamma<1.7$
decrease monotonically with the mass of PBHs, while for $\gamma>1.7$, it has a
plateau around $M_{\rm SMBH}=10^{7}M_{\odot}$. However, the result obtained
from the Benson et al. (2007) mass function has a plateau around the same SMBH
mass for all values of $\gamma$. On the other hand, if the ellipsoidal-
collapse dark matter halo models are plausible, the existence of the plateau
is evident for all the mass functions used and for all values of the power-law
index. These results suggest that dark-matter spikes structured around central
SMBHs with mass $M_{\rm SMBH}=10^{7}M_{\odot}$ may contain a non-standard
abundance of NSs that disturbs the monotonic dependence of the merger rate on
the mass SMBHs.

In Fig. 5, to quantitatively compare the results obtained from ellipsoidal-
collapse dark matter halo models with those obtained from spherical-collapse
dark matter halo models, we have displayed the merger rate of compact binaries
in terms of power-law index $\gamma$ while taking into account Shankar et al.
(2004), Vika et al. (2009), and Benson et al. (2007) mass functions. To assess
the theoretical predictions via the experimental data, we have also included
the relevant mergers estimated by the GW detectors in the case of PBH-PBH
events, while such an action has not been taken in the case of PBH-NS events.
This is because the prediction of the present analysis from the merger rate of
PBH-PBH binaries is capable of justifying the data recorded by GW detectors,
i.e., $(17.9\mbox{-}44)\,{\rm Gpc^{-3}Yr^{-1}}$ (Abbott et al., 2021b),
whereas the merger rate of PBH-NS binaries cannot in any way perform this
task. Specifically, the estimate of the LIGO-Virgo detectors from the total
merger rate of BH-NS binaries is presented as $(7.8\mbox{-}140)\,{\rm
Gpc^{-3}Yr^{-1}}$ (Abbott et al., 2021d), while the relevant prediction of
both dark matter halo models in the current analysis is extremely far from
this range.

As it is evident from the figures related to PBH-PBH events, and of course as
mentioned earlier, the merger rate of PBH-PBH binaries while accounting for
ellipsoidal-collapse dark matter halo models is higher than that extracted
from spherical-collapse dark matter halo models. In addition, such enhancement
in the merger rate is greater at lower power-law index values than that at
higher ones. In this regard, the merger rate of PBH-PBH binaries in spherical-
collapse dark matter halo models while considering Shankar et al. (2004), Vika
et al. (2009), and Benson et al. (2007) mass functions will be consistent with
the BH-BH mergers estimated by GW detectors if the value of power-law index
lies in the interval $\gamma=(1.05\mbox{-}1.15)$, $\gamma=(1.10\mbox{-}1.20)$,
and $\gamma=(1.40\mbox{-}1.50)$, respectively. However, the corresponding
results obtained from ellipsoidal-collapse dark matter halo models can
potentially modify these values as $\gamma=(0.95\mbox{-}1.05)$,
$\gamma=(1.0\mbox{-}1.20)$, and $\gamma=(1.30\mbox{-}1.40)$, respectively.

On the other hand, the results of PBH-NS events indicate that despite the
mismatch of the outcome of the present analysis with GW data, the effect of
ellipsoidal-collapse dark matter halo models in the amplification of the
merger rate of such binaries is significant. It is also obvious that the
merger rate of PBH-NS binaries in both dark matter halo models increases
monotonical with the power-law index reaches a maximum at $\gamma=1.7$, and
decreases for higher values of $\gamma$.

##  4 Conclusions

In this work, we have calculated the merger rate of compact binaries in dark-
matter spikes, which are expected to be structured around SMBHs at the center
of galactic halos. For this purpose, we have initially described theoretical
models, which suit dark matter spikes. We have also discussed crucial
quantities for dark-matter spikes, such as the density profile, concentration
parameter, and $M_{\rm SMBH}\mbox{-}\sigma$ relation, which can specify the
distribution of dark matter particles in the region of spikes. On the other
hand, the strong correlation between the growth of central SMBHs and halo
parameters suggests that another quantity called the mass function of SMBHs
also plays a prominent role in the present analysis. However, the
insufficiency of our knowledge of the exact distribution of dark matter
particles in the central regions of galactic halos and the abundance of SMBHs
in the Universe may lead to uncertainties in the results. Hence, to manage
this uncertainty, we consider three empirical SMBH mass functions to compare
their results.

In the following, relying on the PBH scenario and with the assumption that
PBHs are capable of contributing to the structure of dark matter, we have
discussed the encountering conditions of compact objects such as PBHs and NSs
in dark-matter spikes. Consequently, we have calculated the merger rate of
compact binaries in a single dark-matter spike for spherical- and ellipsoidal-
collapse dark matter halo models. Our results confirm that the merger rate of
PBH-PBH binaries within each spike for ellipsoidal-collapse dark matter halo
models is slightly higher than that derived from spherical-collapse dark
matter halo models. It is concluded that the merger rate of PBH-PBH binaries
changes directly with the value of the power-law index $\gamma$. The results
obtained from the analysis of PBH-NS events show that the maximum value of the
merger rate in each spike takes place around the central SMBH with a mass
$M_{\rm SMBH}=10^{7}\,M_{\odot}$. Moreover, it turned out that the merger rate
of PBH-NS binaries per spike has a maximum value at $\gamma=1.7$. Also, the
results indicate that the merger rate of PBH-NS binaries per spike for
ellipsoidal-collapse dark matter halo models is much higher than that derived
from spherical-collapse dark matter halo models.

As the main measurable factor in GW detectors, we have calculated the total
merger rate of compact binaries in dark-matter spikes around central SMBHs
with masses of $M_{\rm SMBH}=(10^{6}\mbox{-}10^{9})M_{\odot}$. As mentioned
earlier, to account for possible uncertainties in our analysis, we have used
three different mass functions for calculating the total merger rate of
compact binaries. Our findings indicate that the merger rate of PBH-PBH
binaries in spherical- and ellipsoidal-collapse dark matter halo models has
the maximum value at $M_{\rm SMBH}=10^{6}M_{\odot}$ and decreases
monotonically with increasing the mass of SMBHs. Also, the results exhibit
that the overall merger rate of PBH-PBH binaries for ellipsoidal-collapse dark
matter halo models is slightly higher than that obtained from spherical-
collapse dark matter halo models, in such a way that in the best case, e.g. at
about $\gamma=0.05$, the amplification of the overall merger rate is about
$62\%$. This means that in the limit of $r_{\rm sp}\rightarrow 4r_{\rm s}$,
the effect of ellipsoidal-collapse dark matter halo models on the enhancement
of the merger rate of PBH-PBH binaries increases. Moreover, among the
considered mass functions, the highest, middle, and lowest values of the
merger rate of PBH-PBH binaries have been extracted from Shankar et al.
(2004), Vika et al. (2009), and Benson et al. (2007) mass functions,
respectively. However, the results of the merger rate of PBH-NS binaries were
obtained slightly different. Interestingly, the inverse proportionality of the
overall merger rate PBH-NS binaries with the mass of SMBHs is not monotonic
and has a plateau at $M_{\rm SMBH}=10^{7}M_{\odot}$ in almost all considered
models. This may be due to a non-standard abundance of NSs in dark-matter
spikes clustered around the central SMBH with of mass $M_{\rm
SMBH}=10^{7}M_{\odot}$, which should be validated with informative
observational data.

Finally, we have calculated the merger rate of compact binaries according to
the power-law index $\gamma$ for ellipsoidal-collapse dark matter halo models
and compared them with the corresponding results of spherical-collapse dark
matter halo models. To compare our findings with the experimental data of GWs,
we have also included the LIGO-Virgo sensitivity band for the merger rate of
BH-BH binaries. This task was not possible to do for PBH-NS events, as our
results were far from those estimated via the LIGO-Virgo detectors. Our
results show that the inclusion of ellipsoidal-collapse dark matter halo
models in the calculations of the merger rate of PBH-PBH binaries can reduce
the range of the power-law index obtained from spherical-collapse dark matter
halo models by $0.1$. This result comes from the fact that the dark-matter
spikes in ellipsoidal-collapse halo models are denser (and naturally smaller
in radius) than those in spherical-collapse halo models. Additionally, the
results show that, unlike PBH-PBH events, the effect of ellipsoidal-collapse
dark matter halo models in the amplification of the merger rate of PBH-NS
binaries is significant. A maximum at $\gamma=1.7$ is also confirmed for such
events.

Although the current analysis gives interesting and reasonable results, the
conditions considered in it are not necessarily satisfied and the problem
always includes uncertainties. For example, the inaccuracy of the $M_{\rm
SMBH}\mbox{-}\sigma$ relation, uncertainties arising from the SMBH mass
functions, the incomplete knowledge of the exact distribution of dark matter
in the spike regions and the possible disruption of the formed binaries via
the surrounding compact objects before their mergers may affect our final
results. Hopefully, one can achieve a more detailed insight into the merger
event of compact binaries in the vicinity of SMBHs by improving relevant
instruments and increasing our understanding of these issues.

## References

  * Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2016b) —. 2016b, Phys. Rev. Lett., 116, 241103, doi: 10.1103/PhysRevLett.116.241103
  * Abbott et al. (2016c) —. 2016c, Phys. Rev. Lett., 116, 221101, doi: 10.1103/PhysRevLett.116.221101
  * Abbott et al. (2019) —. 2019, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2020b) —. 2020b, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abbott et al. (2021a) —. 2021a, Phys. Rev. X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021b, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021c) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021c, Astrophys. J. Lett., 915, L5, doi: 10.3847/2041-8213/ac082e
  * Abbott et al. (2021d) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021d, arXiv e-prints, arXiv:2111.03634.  https://arxiv.org/abs/2111.03634
  * Alcock et al. (2000) Alcock, C., Allsman, R. A., Alves, D. R., et al. 2000, Astrophys. J., 542, 281, doi: 10.1086/309512
  * Allahyari et al. (2017) Allahyari, A., Firouzjaee, J. T., & Abolhasani, A. A. 2017, J. Cosmol. Astropart. Phys., 2017, 041, doi: 10.1088/1475-7516/2017/06/041
  * Bansal et al. (2022) Bansal, A., Ichiki, K., Tashiro, H., & Matsuoka, Y. 2022, arXiv e-prints, arXiv:2206.01443.  https://arxiv.org/abs/2206.01443
  * Barbieri et al. (2020) Barbieri, C., Salafia, O. S., Perego, A., Colpi, M., & Ghirlanda, G. 2020, Eur. Phys. J. A, 56, 8, doi: 10.1140/epja/s10050-019-00013-x
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, Astrophys. J., 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Bellomo et al. (2018) Bellomo, N., Bernal, J. L., Raccanelli, A., & Verde, L. 2018, J. Cosmol. Astropart. Phys., 2018, 004, doi: 10.1088/1475-7516/2018/01/004
  * Benson et al. (2007) Benson, A. J., Džanović, D., Frenk, C. S., & Sharples, R. 2007, Mon. Notices Royal Astron. Soc., 379, 841, doi: 10.1111/j.1365-2966.2007.11923.x
  * Bird et al. (2016) Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, Phys. Rev. Lett., 116, 201301, doi: 10.1103/PhysRevLett.116.201301
  * Blinnikov et al. (2016) Blinnikov, S., Dolgov, A., Porayko, N. K., & Postnov, K. 2016, J. Cosmol. Astropart. Phys., 2016, 036, doi: 10.1088/1475-7516/2016/11/036
  * Bloomfield et al. (2015) Bloomfield, J., Bulhosa, D., & Face, S. 2015, arXiv e-prints, arXiv:1504.02071.  https://arxiv.org/abs/1504.02071
  * Bouffanais et al. (2021) Bouffanais, Y., Mapelli, M., Santoliquido, F., et al. 2021, Mon. Not. Roy. Astron. Soc, 507, 5224, doi: 10.1093/mnras/stab2438
  * Braine et al. (2020) Braine, T., Cervantes, R., Crisosto, N., et al. 2020, Phys. Rev. Lett., 124, 101303, doi: 10.1103/PhysRevLett.124.101303
  * Carr et al. (2021) Carr, B., Kohri, K., Sendouda, Y., & Yokoyama, J. 2021, Rept. Prog. Phys., 84, 116902, doi: 10.1088/1361-6633/ac1e31
  * Carr & Kuhnel (2021) Carr, B., & Kuhnel, F. 2021, arXiv e-prints, arXiv:2110.02821.  https://arxiv.org/abs/2110.02821
  * Carr et al. (2017) Carr, B., Raidal, M., Tenkanen, T., Vaskonen, V., & Veermäe, H. 2017, Phys. Rev. D, 96, 023514, doi: 10.1103/PhysRevD.96.023514
  * Carr & Hawking (1974) Carr, B. J., & Hawking, S. W. 1974, Mon. Not. Roy. Astron. Soc, 168, 399, doi: 10.1093/mnras/168.2.399
  * Choquette et al. (2019) Choquette, J., Cline, J. M., & Cornell, J. M. 2019, J. Cosmol. Astropart. Phys., 2019, 036, doi: 10.1088/1475-7516/2019/07/036
  * Clesse & García-Bellido (2017) Clesse, S., & García-Bellido, J. 2017, Phys. Dark Univ., 15, 142, doi: 10.1016/j.dark.2016.10.002
  * Di Giovanni et al. (2020) Di Giovanni, F., Fakhry, S., Sanchis-Gual, N., Degollado, J. C., & Font, J. A. 2020, Phys. Rev. D, 102, 084063, doi: 10.1103/PhysRevD.102.084063
  * Di Giovanni et al. (2021) —. 2021, Class. Quantum. Grav., 38, 194001, doi: 10.1088/1361-6382/ac1b45
  * Di Valentino et al. (2020) Di Valentino, E., Melchiorri, A., Mena, O., & Vagnozzi, S. 2020, Phys. Dark Universe, 30, 100666, doi: 10.1016/j.dark.2020.100666
  * Diana et al. (2022) Diana, A., Caccianiga, A., Ighina, L., et al. 2022, Mon. Notices Royal Astron. Soc., 511, 5436, doi: 10.1093/mnras/stac364
  * Dutton & Macciò (2014) Dutton, A. A., & Macciò, A. V. 2014, Mon. Not. Roy. Astron. Soc, 441, 3359, doi: 10.1093/mnras/stu742
  * Fakhry & Del Popolo (2022) Fakhry, S., & Del Popolo, A. 2022, arXiv e-prints, arXiv:2212.08646.  https://arxiv.org/abs/2212.08646
  * Fakhry et al. (2021) Fakhry, S., Firouzjaee, J. T., & Farhoudi, M. 2021, Phys. Rev. D, 103, 123014, doi: 10.1103/PhysRevD.103.123014
  * Fakhry et al. (2022a) Fakhry, S., Naseri, M., Firouzjaee, J. T., & Farhoudi, M. 2022a, Phys. Rev. D, 105, 043525, doi: 10.1103/PhysRevD.105.043525
  * Fakhry et al. (2022b) Fakhry, S., Sajad Tabasi, S., & Firouzjaee, J. T. 2022b, arXiv e-prints, arXiv:2210.13558.  https://arxiv.org/abs/2210.13558
  * Fakhry et al. (2022c) Fakhry, S., Salehnia, Z., Shirmohammadi, A., & Firouzjaee, J. T. 2022c, Astrophys. J., 941, 36, doi: 10.3847/1538-4357/aca523
  * Feng et al. (2022) Feng, W.-X., Parisi, A., Chen, C.-S., & Lin, F.-L. 2022, J. Cosmol. Astropart. Phys., 2022, 032, doi: 10.1088/1475-7516/2022/08/032
  * Feng et al. (2021) Feng, W.-X., Yu, H.-B., & Zhong, Y.-M. 2021, Astrophys. J. Lett., 914, L26, doi: 10.3847/2041-8213/ac04b0
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, Astrophys. J. Lett., 539, L9, doi: 10.1086/312838
  * Fragione (2021) Fragione, G. 2021, Astrophys. J. Lett., 923, L2, doi: 10.3847/2041-8213/ac3bcd
  * Fragione et al. (2021) Fragione, G., Loeb, A., & Rasio, F. A. 2021, Astrophys. J. Lett., 918, L38, doi: 10.3847/2041-8213/ac225a
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, Astrophys. J. Lett., 539, L13, doi: 10.1086/312840
  * Ghez et al. (1998) Ghez, A. M., Klein, B. L., Morris, M., & Becklin, E. E. 1998, Astrophys. J., 509, 678, doi: 10.1086/306528
  * Ghodsi Y. et al. (2022) Ghodsi Y., M., Behnamfard, A., Fakhry, S., & Firouzjaee, J. T. 2022, Phys. Dark Universe, 35, 100918, doi: 10.1016/j.dark.2021.100918
  * Ghodsi Yengejeh et al. (2023) Ghodsi Yengejeh, M., Fakhry, S., Firouzjaee, J. T., & Fathi, H. 2023, Phys. Dark Universe, 39, 101144, doi: 10.1016/j.dark.2022.101144
  * Gondolo & Silk (1999) Gondolo, P., & Silk, J. 1999, Phys. Rev. Lett., 83, 1719, doi: 10.1103/PhysRevLett.83.1719
  * Graham et al. (2007) Graham, A. W., Driver, S. P., Allen, P. D., & Liske, J. 2007, Mon. Notices Royal Astron. Soc., 378, 198, doi: 10.1111/j.1365-2966.2007.11770.x
  * Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, Astrophys. J., 698, 198, doi: 10.1088/0004-637X/698/1/198
  * Häring & Rix (2004) Häring, N., & Rix, H.-W. 2004, Astrophys. J. Lett., 604, L89, doi: 10.1086/383567
  * Hawking (1971) Hawking, S. 1971, Mon. Not. Roy. Astron. Soc, 152, 75, doi: 10.1093/mnras/152.1.75
  * Hinderer et al. (2019) Hinderer, T., Nissanke, S., Foucart, F., et al. 2019, Phys. Rev. D, 100, 063021, doi: 10.1103/PhysRevD.100.063021
  * Hopkins et al. (2007) Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Krause, E. 2007, Astrophys. J., 669, 67, doi: 10.1086/521601
  * Huterer & Shafer (2018) Huterer, D., & Shafer, D. L. 2018, Rep. Prog. Phys., 81, 016901, doi: 10.1088/1361-6633/aa997e
  * Kashlinsky (2016) Kashlinsky, A. 2016, Astrophys. J. Lett., 823, L25, doi: 10.3847/2041-8205/823/2/L25
  * Kelly & Merloni (2012) Kelly, B. C., & Merloni, A. 2012, Adv. Astron., 2012, 970858, doi: 10.1155/2012/970858
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, Annu. Rev. Astron. Astrophys., 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Lehmann et al. (2018) Lehmann, B. V., Profumo, S., & Yant, J. 2018, J. Cosmol. Astropart. Phys., 2018, 007, doi: 10.1088/1475-7516/2018/04/007
  * Liebling & Palenzuela (2012) Liebling, S. L., & Palenzuela, C. 2012, Living Rev. Relativ., 15, 6, doi: 10.12942/lrr-2012-6
  * Liske et al. (2003) Liske, J., Lemon, D. J., Driver, S. P., Cross, N. J. G., & Couch, W. J. 2003, Mon. Notices Royal Astron. Soc., 344, 307, doi: 10.1046/j.1365-8711.2003.06826.x
  * Ludlow et al. (2016) Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, Mon. Not. Roy. Astron. Soc, 460, 1214, doi: 10.1093/mnras/stw1046
  * Mandel & Broekgaarden (2022) Mandel, I., & Broekgaarden, F. S. 2022, Living Rev. ReL., 25, 1, doi: 10.1007/s41114-021-00034-3
  * Mandel & Farmer (2022) Mandel, I., & Farmer, A. 2022, Phys. Rept., 955, 1, doi: 10.1016/j.physrep.2022.01.003
  * McLure & Dunlop (2002) McLure, R. J., & Dunlop, J. S. 2002, Mon. Notices Royal Astron. Soc., 331, 795, doi: 10.1046/j.1365-8711.2002.05236.x
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, Astron. J., 132, 2685, doi: 10.1086/508988
  * Musco & Miller (2013) Musco, I., & Miller, J. C. 2013, Class. Quant. Grav., 30, 145009, doi: 10.1088/0264-9381/30/14/145009
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, Astrophys. J., 462, 563, doi: 10.1086/177173
  * Navarro et al. (2010) Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, Mon. Notices Royal Astron. Soc., 402, 21, doi: 10.1111/j.1365-2966.2009.15878.x
  * Nishikawa et al. (2019) Nishikawa, H., Kovetz, E. D., Kamionkowski, M., & Silk, J. 2019, Phys. Rev. D, 99, 043533, doi: 10.1103/PhysRevD.99.043533
  * Nitz & Wang (2021a) Nitz, A. H., & Wang, Y.-F. 2021a, Phys. Rev. Lett., 126, 021103, doi: 10.1103/PhysRevLett.126.021103
  * Nitz & Wang (2021b) —. 2021b, Phys. Rev. Lett., 127, 151101, doi: 10.1103/PhysRevLett.127.151101
  * Okoli & Afshordi (2016) Okoli, C., & Afshordi, N. 2016, Mon. Not. Roy. Astron. Soc, 456, 3068, doi: 10.1093/mnras/stv2905
  * Polnarev & Musco (2007) Polnarev, A. G., & Musco, I. 2007, Class. Quant. Grav., 24, 1405, doi: 10.1088/0264-9381/24/6/003
  * Prada et al. (2012) Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, Mon. Not. Roy. Astron. Soc, 423, 3018, doi: 10.1111/j.1365-2966.2012.21007.x
  * Prokhorenko & Sazonov (2021) Prokhorenko, S. A., & Sazonov, S. Y. 2021, Astron. Lett., 47, 515, doi: 10.1134/S106377372108003X
  * Raidal et al. (2019) Raidal, M., Spethmann, C., Vaskonen, V., & Veermäe, H. 2019, J. Cosmol. Astropart. Phys., 2019, 018, doi: 10.1088/1475-7516/2019/02/018
  * Roszkowski et al. (2018) Roszkowski, L., Sessolo, E. M., & Trojanowski, S. 2018, Rep. Prog. Phys., 81, 066201, doi: 10.1088/1361-6633/aab913
  * Ruiz et al. (2021) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2021, Front. Astron. Space Sci., 8, 39, doi: 10.3389/fspas.2021.656907
  * Sasaki et al. (2016) Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. 2016, Phys. Rev. Lett., 117, 061101, doi: 10.1103/PhysRevLett.117.061101
  * Sasaki et al. (2018) —. 2018, Class. Quant. Grav., 35, 063001, doi: 10.1088/1361-6382/aaa7b4
  * Sasaki et al. (2022) Sasaki, M., Takhistov, V., Vardanyan, V., & Zhang, Y.-l. 2022, Astrophys. J., 931, 2, doi: 10.3847/1538-4357/ac66da
  * Shankar et al. (2004) Shankar, F., Salucci, P., Granato, G. L., De Zotti, G., & Danese, L. 2004, Mon. Notices Royal Astron. Soc., 354, 1020, doi: 10.1111/j.1365-2966.2004.08261.x
  * Shapiro & Heggie (2022) Shapiro, S. L., & Heggie, D. C. 2022, Phys. Rev. D, 106, 043018, doi: 10.1103/PhysRevD.106.043018
  * Sherwin et al. (2011) Sherwin, B. D., Dunkley, J., Das, S., et al. 2011, Phys. Rev. Lett., 107, 021302, doi: 10.1103/PhysRevLett.107.021302
  * Shibata & Sasaki (1999) Shibata, M., & Sasaki, M. 1999, Phys. Rev. D, 60, 084002, doi: 10.1103/PhysRevD.60.084002
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760, doi: 10.1103/PhysRevLett.84.3760
  * Stadel et al. (2009) Stadel, J., Potter, D., Moore, B., et al. 2009, Mon. Notices Royal Astron. Soc., 398, L21, doi: 10.1111/j.1745-3933.2009.00699.x
  * Tiwari et al. (2021) Tiwari, S., Ebersold, M., & Hamilton, E. Z. 2021, Phys. Rev. D, 104, 123024, doi: 10.1103/PhysRevD.104.123024
  * Tremaine et al. (2002) Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, Astrophys. J., 574, 740, doi: 10.1086/341002
  * Vika et al. (2009) Vika, M., Driver, S. P., Graham, A. W., & Liske, J. 2009, Mon. Notices Royal Astron. Soc., 400, 1451, doi: 10.1111/j.1365-2966.2009.15544.x
  * Villanueva-Domingo et al. (2021) Villanueva-Domingo, P., Mena, O., & Palomares-Ruiz, S. 2021, Front. Astron. Space Sci., 8, 87, doi: 10.3389/fspas.2021.681084
  * Volonteri et al. (2021) Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nat. Rev. Phys., 3, 732, doi: 10.1038/s42254-021-00364-9
  * Yang & Xu (2014) Yang, W., & Xu, L. 2014, J. Cosmol. Astropart. Phys., 2014, 034, doi: 10.1088/1475-7516/2014/08/034
  * Young et al. (2014) Young, S., Byrnes, C. T., & Sasaki, M. 2014, J. Cosmol. Astropart. Phys., 2014, 045, doi: 10.1088/1475-7516/2014/07/045
  * Zappa et al. (2019) Zappa, F., Bernuzzi, S., Pannarale, F., Mapelli, M., & Giacobbo, N. 2019, Phys. Rev. Lett., 123, 041102, doi: 10.1103/PhysRevLett.123.041102
  * Zel’dovich & Novikov (1967) Zel’dovich, Y. B., & Novikov, I. D. 1967, Soviet Astron. AJ., 10, 602