# Nonlinear coherent state of an exciton in a wide quantum dot

M. Bagheri Harouni, R. Roknizadeh, M. H. Naderi  Quantum Optics Group, Physics
Department, University of Isfahan<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>, ,

###### Abstract

In this paper, we derive the dynamical algebra of a particle confined in an
infinite spherical well by using the $f$-deformed oscillator approach. We
consider an exciton with definite angular momentum in a wide quantum dot
interacting with two laser beams. We show that under the weak confinement
condition, and quantization of the center-of-mass motion of exciton, the
stationary state of it can be considered as a special kind of nonlinear
coherent states which exhibits the quadrature squeezing.

##  1 Introduction

The conventional coherent states of the quantum harmonic oscillator, defined
by Glauber [1] as the right-hand eigenstates of non-hermitian annihilation
operator $\hat{a}([\hat{a},\hat{a}^{{\dagger}}]=1)$, have found many
interesting applications in different areas of physics such as quantum optics,
condensed matter physics, statistical physics and atomic physics [2]. These
states play an important role in the quantum theory of coherence, are
considered as the most classical ones among the pure quantum states, and laser
light can be supposed as a physical realization of them. Due to the vast
application of these states, there have been many attempts to generalize them
[3]. Among the all generalizations, nonlinear coherent states (NLCS) [4] have
been paid attention in recent years because they exhibit nonclassical features
such as quadrature squeezing and sub-poissonian statistics [5]. These states
are defined as the right-hand eigenstates of a deformed operator $\hat{A}$

$\hat{A}=\hat{a}f(\hat{n})\hskip
56.9055pt\hat{A}|\alpha,f\rangle=\alpha|\alpha,f\rangle,$ (1)

where the deformation function $f(\hat{n})$ is an operator-valued function of
the number operator $\hat{n}$. From (1) one can obtain an explicit form of
NLCS in the number state representation

$\displaystyle|\alpha,f\rangle=\mathcal{N}_{f}\sum_{n}\frac{\alpha^{n}}{\sqrt{n!}f(n)!}|n\rangle,$
$\displaystyle\mathcal{N}_{f}=\left[\sum_{n}\frac{|\alpha|^{2n}}{[f(n)!]^{2}n!}\right]^{-\frac{1}{2}}.$
(2)

A class of NLCS can be realized physically as the stationary state of the
center-of-mass motion of a laser driven trapped ion [6, 7]. Furthermore, it
has been proposed a theoretical scheme to show the possibility of generating
various families of NLCS [8] of the radiation field in a lossless coherently
pumped micromaser within the frame work of the intensity-dependent Jaynes-
Cummings model.
Recently, the influences of the spatial confinement [9] and the curvature of
physical space [10] on the algebraic structure of the coherent states of the
quantum harmonic oscillator have been investigated within the frame work of
nonlinear coherent states approach. It has been shown that if a quantum
harmonic oscillator be confined within a small region of order of its
characteristic length [9] or its physical space to be a sphere [10], then it
can be regarded as a deformed oscillator, i.e., an oscillator that its
creation and annihilation operators are deformed operator $\hat{A}$ and
$\hat{A}^{{\dagger}}$ given by Eq.(1).
On the other hand, we can consider nanostructures as systems whose physical
properties are related to the confinement effects. Thus, we expect that it is
possible to realize some natural deformations in these systems [9, 11]. In
addition, in nanostructures different kinds of quantum states can be prepared.
One of the most applicable of these states is exciton state. Exciton is an
elementary excitation in semiconductors interacting with light, electron in
conduction band which is bounded to hole in valance band that can easily move
through the sample. In one of the nano size systems, quantum dot (QD), due to
the confinement in three dimensions, energy bands reduce to quasi energy
levels. Therefore, in order to describe the interaction of QD with light we
can consider it as a few level atom [12]. These Exciton states can be used in
quantum information processes. It has been shown that excitons in coupled QDs
are ideal for preparation of entangled state in solid-state systems [13].
Entanglement of the exciton states in a single QD or in a QD molecule has been
demonstrated experimentally [14]. Entanglement of the coherent states of the
excitons in a system of two coupled QDs has been considered [15]. Recently,
coherent exciton states of excitonic nano-crystal-molecules has been
considered [16]. Theoretical approach for generating Dick states of excitons
in optically driven QD has been proposed in Ref.[17]. In a QD, the effects of
exciton-phonon interaction, exciton-impurity interaction and exciton-exciton
interaction play an important role. These effects are the main sources for the
decoherence of exciton states [18]. Furthermore, these effects cause the
exciton has the spontaneous recombination or scattered to other exciton modes
[19, 20].
In this paper we propose a theoretical scheme for generating excitonic NLCS.
We will show that under certain conditions the quantized motion of wave packet
of center-of-mass of exciton can be consider as a special kind of NLCSs. Our
scheme is based on the interaction of a quantum dot with two laser beams. By
using the approach considered in Ref.[6], we propose a theoretical scheme for
generation of NLCS of an exciton in a wide QD.
In section 2, we consider different confinement regimes in a QD, and the
explicit forms of the creation and annihilation operators for a particle
confined in an infinite well are derived by using the deformed quantum
oscillator approach. In section 3, we consider an exciton in a wide QD which
interacts with two laser beams. We shall show that under the weak confinement
condition, the stationary state of the exciton center-of-mass motion can be
considered as a NLCS.

##  2 Algebraic approach for a particle in an infinite spherical well

In nanostructures and confined systems, there are three different confinement
regimes. The criteria for this classification is based on the comparison
between excitation Bohr radius and the spatial dimensions of the system under
consideration. In the case of a QD, these regimes are defined as follows [21].
We first introduce three quantities $\Delta E_{c}$, $\Delta E_{v}$ and
$V_{exc}$ which, respectively denote: the electron energy due to the
confinement, the hole energy due to the confinement and Coulomb energy between
correlated electron-hole (exciton).
1) $V_{exc}>\Delta E_{c}-\Delta E_{v}$: In this case, the exciton energy is
much greater than the confinement energies of electron and hole. If we show
the system size by $L$ and the exciton Bohr radius by $a$, then in this regime
$L>a$. This regime corresponds to the weak confinement (in some literature the
weak confinement is characterized by the situation in which the electron and
the hole are not in the same matter, for example, hole be in QD and excited
electron in host matter. In this paper, by the weak confinement regime we mean
$L>a$ and the excitations in the same matter). In this regime due to the
confinement, the center-of-mass motion of the exciton is quantized and the
confinement do not affect electron and hole separately. Hence, the confinement
affect the exciton motion as a whole [22].
2) $V_{exc}<\Delta E_{c}\;,\Delta E_{v}$: This regime, in contrast to the
previous one, is associated with the cases where $L<a$. In this regime the
exciton is completely localized, and the confinement affects both the electron
and the hole independently and their states become quantized in conduction and
valance bands. This regime is called strong confinement.
3) $\Delta E_{c}>V_{exc}\;,\Delta E_{v}$: This condition is equivalent to the
situation $a_{c}<a<a_{v}$, where $a_{c}$ and $a_{v}$ are, respectively, the
Bohr radii of electron and hole. Here, due to the different effective masses
of electron and hole, the hole which has heavier effective mass is localized
and the electron motion will be quantized. This regime is called intermediate
confinement.
In the first case (weak confinement), in a wide QD, an exciton can move due to
its center-of-mass momentum, and because of the presence of the barriers, its
center-of-mass motion is quantized. Therefore, it moves as a whole between
energy levels of an infinite well. We consider a wide spherical QD whose
energy levels are equivalent to the energy levels of a spherical well

$E_{nl}=\frac{\hbar^{2}}{2M}\frac{\alpha_{nl}^{2}}{R^{2}},$ (3)

where $\alpha_{nl}$ is the n’th zero of the first kind Bessel function of
order $l$, $j_{l}(x)$. In this energy spectrum according to the azimuthal
symmetry around $z$ axis, we have a degenerate spectrum. As mentioned before,
in the weak confinement regime, the Coulomb potential plays an essential role
and its spectrum is given by

$E_{k}^{b}=\frac{\mu e^{4}}{2\hbar^{2}\varepsilon^{2}}\frac{1}{k^{2}}\hskip
14.22636pt,\hskip 14.22636pt\mu=\frac{m_{e}m_{h}}{m_{e}+m_{h}},$ (4)

where superscript $b$ shows binding energy related to the Coulomb interaction
and $\varepsilon$ shows dielectric constant of the system. As is usual, we
interpret the Coulomb part as an exciton and another degree of freedom (motion
between energy levels of the well) as the exciton center-of-mass motion.
Therefore, in a wide QD an exciton has two different kinds of degrees of
freedom: internal degrees of freedom due to the Coulomb potential and external
degrees of freedom related to the quantum confinement. Here we consider the
lowest exciton state, $1s$ exciton, because this exciton state has the largest
oscillator strength among other exciton state. Then the energy of the exciton
in a wide QD can be written as

$E_{nlm}=E_{g}-E_{1}^{b}+\frac{\hbar^{2}}{2MR^{2}}\alpha_{nl}^{2},$ (5)

where $E_{g}$ is the energy gap of QD, $E_{1}^{b}=E_{k}^{b}|_{k=1}$ is the
exciton binding energy, $M=m_{e}+m_{h}$ is the total mass of exciton, and $R$
is the radius of QD. Due to the relation of quantum numbers $l$ and $m$ with
the angular momentum and the selection rules for optical transitions, we can
fix $l$ and $m$ (by choosing a certain condition), and hence the energy of
exciton depends only on a single quantum number

$E_{n}=E_{g}-E_{1}^{b}+\frac{\hbar^{2}}{2MR^{2}}\alpha_{nl}^{2}.$ (6)

Therefore, we can prepare the conditions under which the exciton center-of-
mass motion has a one-dimensional degree of freedom. Due to the quantization
of the exciton center-of-mass motion, we can describe the exciton motion
between the energy levels by the action of a special kind of ladder operators.
In order to find these operators we use the $f$-deformed oscillator approach
[4].
As mentioned elsewhere [9], if the energy spectrum of the system is equally
spaced, such as harmonic oscillator, its creation and annihilation operators
satisfy the ordinary Weyl-Heisenberg algebra, otherwise we can interpret them
as the generators of a generalized Weyl-Heisenberg algebra.
The energy spectrum of a particle with mass $M$ confined in an infinite
spherical well can be written as (3). According to the conservation of angular
momentum, we assume that particle has been prepared with definite angular
momentum (for example by measuring its angular momentum). Then $l$ becomes
completely determined, i.e., in the energy spectrum the number $l$ is a
constant. By determining the number $l$ and considering the rotational
symmetry of the system around the $z$ axis, the angular part of the spectrum
becomes completely determined, and the radius part is described by (3). Now we
use a factorization method and write the Hamiltonian of the center-of-mass
motion of the system as follows

$\hat{H}=\frac{1}{2}(\hat{A}\hat{A}^{{\dagger}}+\hat{A}^{{\dagger}}\hat{A}),$
(7)

where $\hat{A}$ and $\hat{A}^{{\dagger}}$ are defined through the relation
(1). Therefore the spectrum of $\hat{H}$, after straightforward calculation,
is obtained as

$E_{n}=\frac{1}{2}[(n+1)f^{2}(n+1)+nf^{2}(n)].$ (8)

By comparing (8) with Eq.(3) we arrive at the following expression for the
corresponding deformation function $f_{1}(\hat{n})$

$f_{1}(n)=\sqrt{\frac{\hbar^{2}}{MR^{2}}\frac{(-1)^{n}}{n}\sum_{i=1}^{n}(-1)^{i}\alpha_{i-1l}^{2}}.$
(9)

Then, the ladder operators associated with the radial motion of a confined
particle in a spherical infinite well is given by

$\displaystyle\hat{A}=\hat{a}\sqrt{\frac{\hbar^{2}}{MR^{2}}\frac{(-1)^{n}}{n}\sum_{i=1}^{n}(-1)^{i}\alpha_{i-1l}^{2}},$
(10)
$\displaystyle\hat{A}^{{\dagger}}=\sqrt{\frac{\hbar^{2}}{MR^{2}}\frac{(-1)^{n}}{n}\sum_{i=1}^{n}(-1)^{i}\alpha_{i-1l}^{2}}\,\,\hat{a}^{{\dagger}}.$

These two deformed operators obey the following commutation relation

$[\,\hat{A}\,,\,\hat{A}^{{\dagger}}]=-nf_{1}^{2}(n)+\frac{\hbar^{2}}{MR^{2}}\alpha_{nl}^{2}.$
(11)

As is usual in the $f$-deformation approach, for a particular limit of the
corresponding deformation parameter, the deformed algebra should be reduced to
the conventional oscillator algebra. However, in this treatment we note that
there is no thing in common between the harmonic oscillator potential and an
infinite spherical well. Only in the limit $R\rightarrow\infty$, the system
reduces to a free particle which has continuous spectrum.
As a result, in this section we conclude that the radial motion of a particle
confined in a three-dimensional infinite spherical well can be interpreted by
an $f$-deformed Weyl-Heisenberg algebra.

##  3 Exciton dynamics in QD

Now we consider the formation of an exciton and its dynamics in a wide QD
during the exciton lifetime. As mentioned before, in this situation the
center-of-mass motion of the exciton is quantized. The exciton is created
during the interaction of a QD with light, and because of the angular momentum
conservation, the exciton has a well-defined angular momentum. The exciton is
a quasiparticle composed of an electron and a hole and thus the exciton spin
state can be in a singlet state or a triplet state. According to the optical
transition selection rules, the triplet state is optically inactive and is
called dark exciton [23]. By adding spin and angular momentum of absorbed
photons, the angular momentum of the exciton state can be determined. Hence,
the exciton behaves like a particle in a spherical well with the definite
angular momentum. According to the previous section, the center-of-mass motion
of the exciton in the QD and the barriers of QD can be described by an
oscillator-like Hamiltonian expressed in terms of the $f$-deformed
annihilation and creation operators given by Eq.(10)

$H_{well}=\frac{1}{2}(\hat{A}\hat{A}^{{\dagger}}+\hat{A}^{{\dagger}}\hat{A}),$
(12)

where we interpret the operator $\hat{A}$ $(\hat{A}^{{\dagger}})$ as the
operator whose action causes the transition of exciton center-of-mass motion
to a lower (an upper) energy state. In fact the Hamiltonian (12) is related to
the external degree of freedom of exciton. On the other hand, one can imagine
QD as a two-level system with the ground state $|g\rangle$ and the excited
state $|e\rangle$ (associated with the presence of exciton). Thus, for the
internal degree of freedom we can consider the following Hamiltonian

$H_{ex}=\hbar\omega_{ex}\hat{S}_{22},$ (13)

where $\hat{S}_{22}=|e\rangle\langle e|-|g\rangle\langle g|$ and
$\hbar\omega_{ex}=E_{g}-E_{1}^{b}$ is the exciton energy.
We consider a single exciton of frequency $\omega_{ex}$ confined in a wide QD
interacting with two laser fields, respectively, tuned to the internal degree
of freedom of the frequency $\omega_{ex}$ and to the non-equal spaced energy
levels of the infinite well. It is necessary that the second laser has special
conditions, because it should give rise to the transitions between energy
levels whose frequencies depend on intensity. The interacting system can be
described by the Hamiltonian

$\hat{H}=\hat{H}_{0}+\hat{H}_{int},$ (14)

where $\hat{H}_{0}=\hat{H}_{well}+\hat{H}_{ex}$ and

$H_{int}=g[E_{0}e^{-i(k_{0}\hat{x}-\omega_{ex}t)}+E_{1}e^{-i(k_{1}\hat{x}-(\omega_{ex}-\omega_{\overline{n}})t)}]\hat{S}_{12}+H.c.,$
(15)

in which $g$ is the coupling constant, $k_{0}$ and $k_{1}$ are the wave
vectors of the laser fields, $\hat{S}_{12}=|g\rangle\langle e|$ is the exciton
annihilation operator, and $\omega_{\overline{n}}$ is the frequency of exciton
transition between energy levels of QD due to the spatial confinement. Here,
we consider transition between specific side-band levels hence, we show the
frequency transition with definite dependence to $n$. We show this by a
c-number quantity $\overline{n}$.
The exciton has a finite lifetime that in systems with small dimension, is
increased [24]. The interaction with phonons is the main reason of damping of
the exciton [25]. Phonons in bulk matter have a continuous spectrum while in a
confined system such as QD their spectrum becomes discrete. Hence in a QD, the
resonant interaction between the exciton and phonons decreases and in this
system the exciton lifetime will increase. Therefore during the lifetime of an
exciton, its dynamics is under influence of a bath reservoir, and its damping
play an important role. We assume that during the presence of the exciton in
QD, it interacts with the reservoir and hence we can consider its steady
state. We consider an exciton in dark state. Experimental preparation methods
of such exciton has been described in [23]. In this situation lifetime of
exciton will increase and exciton has not spontaneously recombination
radiation. However, its interaction with phonons causes a finite lifetime for
it.
The operator of the center-of-mass motion position $\hat{x}$ of the exciton in
a spherical QD may be defined as

$\hat{x}=\frac{\kappa}{k_{ex}}(\hat{A}+\hat{A}^{{\dagger}}),$ (16)

where $\kappa$ being a parameter similar to the Lamb-Dick parameter in ion
trapped systems and is defined as the ratio of QD radius to the wavelength of
the driving laser (because of the spatial confinement of exciton, its wave
function width is determined by the barriers of QD), and we assume
$k_{0}\simeq k_{1}\simeq k_{ex}$ ($k_{ex}$ is the wavevector of the exciton).
The operators $\hat{A}$ and $\hat{A}^{{\dagger}}$ are defined in Eq.(10). The
interaction Hamiltonian (15) can be written as

$H_{int}=\hbar
e^{i\omega_{ex}t}\Omega_{1}\left[\frac{\Omega_{0}}{\Omega_{1}}+e^{-i\omega_{\overline{n}}t}\right]e^{i\kappa(\hat{A}+\hat{A}^{{\dagger}})}\hat{S}_{12}+H.c.,$
(17)

where $\Omega_{0}=\frac{gE_{0}}{\hbar}$ and $\Omega_{1}=\frac{gE_{1}}{\hbar}$
are the Rabi frequencies of the lasers, respectively, tuned to the electronic
transition of QD (internal degree of freedom) and the first center-of-mass
motion transition of exciton. Since the external degree of freedom is
definite, then $\omega_{\overline{n}}$ depends on a special value of $n$ such
that it can be consider as a c-number quantity. The frequency
$\omega_{\overline{n}}$ is depend on the number of quanta for each transition
and hence the laser tuned to the center-of-mass motion must be so strong that
causes transition. This allows us to treat the interaction of the confined
exciton in a wide QD with two lasers separately, by using a nonlinear Jaynes-
Cummings Hamiltonian [26] for each coupling. The interaction Hamiltonian in
the interaction picture can be written as

$H_{I}=\hbar\Omega_{1}\hat{S}_{12}\left[\frac{\Omega_{0}}{\Omega_{1}}+e^{i\omega_{\overline{n}}t}\right]\exp[i\kappa(e^{-i\omega_{\hat{n}}t}\hat{A}+\hat{A}^{{\dagger}}e^{i\omega_{\hat{n}}t})]+H.c.,$
(18)

where
$\omega_{\hat{n}}=\frac{1}{2\hbar}[(\hat{n}+2)f_{1}(\hat{n}+2)-\hat{n}f_{1}(\hat{n})]$.
By using the vibrational rotating wave approximation [6], applying the
disentangling formula introduced in [27], and using the fact that in the
present case the Lamb-Dick parameter is small, the interaction Hamiltonian
(18) is simplified to

$H_{I}^{(1)}=\hbar\Omega_{1}\hat{S}_{12}\left[F_{0}(\hat{n},\kappa)\frac{\Omega_{0}}{\Omega_{1}}+i\kappa
F_{1}(\hat{n},\kappa)\hat{a}\right]+H.c.,$ (19)

where the function $F_{i}(\hat{n},\kappa)\;(i=0,1)$ is defined by

$\displaystyle F_{i}(\hat{n},\kappa)$ $\displaystyle=$ $\displaystyle
e^{-\frac{\kappa^{2}}{2}((n+1+i)f_{1}^{2}(n+1+i)-(n+i)f_{1}^{2}(n+i))}\times$
$\displaystyle\sum_{l=0}^{n}\frac{\left(i\kappa\right)^{2l}}{l!(l+i)!}\frac{f_{1}(\hat{n})f_{1}(\hat{n}+i)}{[f_{1}(\hat{n}-l)!]^{2}}(\hat{a}^{{\dagger}})^{l}\hat{a}^{l}.$

It should be noted that this function in the limit $f_{1}(\hat{n})\rightarrow
1$ (which is equivalent to the harmonic confinement) is proportional to the
associated Laguerre polynomials

$F_{i}(\hat{n},\kappa)|_{f_{1}(\hat{n})\rightarrow
1}=\frac{e^{-\frac{\kappa^{2}}{2}}}{\hat{n}+i}L_{\hat{n}}^{i}\left(\kappa^{2}\right).$
(21)

Now we write the function $F_{i}(\hat{n},\kappa)$ (3)

$F_{i}(\hat{n},\kappa)=\frac{e^{-\frac{\kappa^{2}}{2}((n+1+i)f_{1}^{2}(n+1+i)-(n+i)f_{1}^{2}(n+i))}}{\hat{n}+i}f_{1}(\hat{n})!f_{1}(\hat{n}+1)!L_{f,\hat{n}}^{i}\left(\kappa^{2}\right),$
(22)

where the function $L_{f,\hat{n}}^{i}(x)$ is defined as

$L_{f,\hat{n}}^{i}(x)=\sum_{l=0}^{n}\frac{1}{[f_{1}(\hat{n}-l)!]^{2}}\frac{(\hat{n}+i)!}{(\hat{n}-l)!l!(l+i)!}(-x)^{l}.$
(23)

This function is similar to the associated Laguerre function.
The time evolution of the system under consideration is characterized by the
master equation

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}[\hat{H}_{I}^{(1)},\hat{\rho}]+\mathfrak{L}\hat{\rho},$
(24)

where $\mathfrak{L}\hat{\rho}$ defines damping of the system due to the
different kinds of interactions which lead to annihilation of exciton. We
assume a bosonic reservoir that causes damping of exciton system. Due to the
properties of dark exciton, the rate of spontaneous recombination and hence
spontaneous emission is decrease. On the other hand, interactions of exciton-
phonon and exciton-impurities cause the exciton to be damped. In fact in low
temperatures it is possible to ignore the phonon effects and by assuming a
pure system we neglect the impurity effects. Hence we can write

$\mathfrak{L}\hat{\rho}=\frac{\Gamma}{2}(2\hat{b}\hat{\rho}\hat{b}^{{\dagger}}-\hat{b}^{{\dagger}}\hat{b}\hat{\rho}-\hat{\rho}\hat{b}^{{\dagger}}\hat{b}),$
(25)

where $\Gamma$ is the energy relaxation rate, $\hat{b}$ and
$\hat{b}^{{\dagger}}$ are the annihilation and creation operators of the
reservoir. Due to the confinement and dark state properties, spontaneous
recombination of exciton decreases and hence the lifetime of exciton becomes
so long that we can consider the stationary solution of Eq.(24). We assume a
finite lifetime for exciton, and during this time we neglect damping effects.
The stationary solution of the master equation (24) in the time scales of our
interest is

$\hat{\rho}=|e\rangle|\psi\rangle\langle\psi|\langle e|,$ (26)

where $|e\rangle$ is the electronic excited state correspond to the presence
of exciton and $|\psi\rangle$ is the center-of-mass motion state of the
exciton, which can be considered as a right-hand eigenstate of the deformed
operator $\hat{A}=\frac{F_{1}(\hat{n},\kappa)}{F_{0}(\hat{n},\kappa)}\hat{a}$

$\frac{F_{1}(\hat{n},\kappa)}{F_{0}(\hat{n},\kappa)}\hat{a}|\psi\rangle=\frac{i\Omega_{0}}{\Omega_{1}\kappa}|\psi\rangle.$
(27)

According to Eq.(22) the corresponding deformation function reads as

$\displaystyle f(\hat{n})$ $\displaystyle=$
$\displaystyle\frac{F_{1}(\hat{n}-1,\kappa)}{F_{0}(\hat{n}-1,\kappa)}$
$\displaystyle=$
$\displaystyle\frac{f_{1}(\hat{n})L_{f,\hat{n}-1}^{1}(\kappa^{2})}{nL_{f,\hat{n}-1}^{0}(\kappa^{2})}e^{-\frac{\kappa^{2}}{2}\left((n+1)f_{1}^{2}(n+1)-(n-1)f_{1}^{2}(n-1)\right)}.$

Hence, we can express the state $|\psi\rangle$ in the Fock space
representation as

$|\psi\rangle=\mathcal{N}_{f}\sum_{n}\frac{\chi^{n}}{\sqrt{n!}f(n)!}|n\rangle,$
(29)

where $\chi=\frac{i\Omega_{0}}{\kappa\Omega_{1}}$. According to the definition
(1), it is evident that the state $|\psi\rangle$ can be regarded as a special
kind of NLCS. As is seen from equation (27), the eigenvalues of the deformed
operator $\hat{A}$ depends on some physical parameters such as the Rabi
frequencies, the parameter $\kappa$ and radius of QD.
As is clear from equation (3), the deformation function $f(\hat{n})$ depends
on the quantum number $\hat{n}$ and physical parameters such as QD radius and
$\kappa$ which characterizes the confinement regime. In the limit
$f_{1}(\hat{n})\rightarrow 1$, (harmonic confinement), which corresponds, for
example, to a QD in lens shape [28], the function $L_{f,\hat{n}}^{i}$ reduces
to the ordinary associated Laguerre polynomials, its argument tends to
$\kappa^{2}$ and therefore, the deformation function (3) takes the following
form

$f(\hat{n})=e^{-\kappa^{2}}L_{\hat{n}}^{1}(\kappa^{2})[(\hat{n}+1)L_{\hat{n}}^{0}(\kappa^{2})]^{-1}.$
(30)

This is the deformation function that appears in the center-of-mass motion of
a trapped ion confined in a harmonic trap [6].
In order to investigate the nonclassical behavior of the NLCS $|\psi\rangle$
we consider the quadrature squeezing of the center-of-mass motion. For this
purpose, we define the deformed quadratures operators as follows

$\hat{X}_{1}=\frac{1}{2}(\hat{A}e^{i\phi}+\hat{A}^{{\dagger}}e^{-i\phi}),\hskip
28.45274pt\hat{X}_{2}=\frac{1}{2i}(\hat{A}e^{i\phi}-\hat{A}^{{\dagger}}e^{-i\phi}).$
(31)

In the limiting case $f(\hat{n})\rightarrow 1$, these two operators reduce to
the conventional (non-deformed) quadrature operators [29]. The commutation
relation of $\hat{X}_{1}$ and $\hat{X}_{2}$ is

$[\hat{X}_{1},\hat{X}_{2}]=\frac{i}{2}[(\hat{n}+1)f^{2}(\hat{n}+1)-\hat{n}f^{2}(\hat{n})].$
(32)

The variances
$\langle(\Delta\hat{X}_{i})^{2}\rangle\equiv\langle\hat{X}_{i}^{2}\rangle-\langle\hat{X}_{i}\rangle^{2}(i=1,2)$
satisfy the uncertainty relation

$\langle(\Delta\hat{X}_{1})^{2}\rangle\langle(\Delta\hat{X}_{2})^{2}\rangle\geq\frac{1}{16}(\langle(\hat{n}+1)f^{2}(\hat{n}+1)-\hat{n}f^{2}(\hat{n})\rangle)$
(33)

A quantum state is said to be squeezed when one of the quadratures components
$\hat{X}_{1}$ and $\hat{X}_{2}$ satisfies the relation

$\langle(\Delta\hat{X}_{i})^{2}\rangle<\frac{1}{4}\langle(\hat{n}+1)f^{2}(\hat{n}+1)-\hat{n}f^{2}(\hat{n})\rangle\hskip
28.45274pti=1\;or\;2$ (34)

The degree of squeezing can be measured by the squeezing parameter
$s_{i}(i=1,2)$ defined by

$s_{i}=4\langle(\Delta\hat{X}_{i})^{2}\rangle-\frac{1}{4}\langle(\hat{n}+1)f^{2}(\hat{n}+1)-\hat{n}f^{2}(\hat{n})\rangle.$
(35)

Then the condition for squeezing in the quadrature component can be simply
written as $s_{i}<0$. In Fig.(1) we plot the squeezing parameter $s_{1}$
versus the parameter $\frac{R}{a_{B}}$ defined as the ratio of the QD radius
to the Bohr radius of exciton for two different values of ratio
$\frac{\Omega_{1}}{\Omega_{0}}$. As is clear from Fig.(1) for small values of
the parameter $\frac{R}{a_{B}}$ the state shows quadrature squeezing and by
increasing this parameter the quadrature squeezing disappears.

##  4 Conclusion

In this paper, we first considered a particle confined in a spherical infinite
well and we found the explicit forms of its creation and annihilation
operators by using the $f$-deformed oscillator approach. Then we considered an
exciton in a wide QD interacts with two laser beams. We showed that under the
weak confinement condition, the exciton is influenced as a whole and its
center-of-mass motion will be quantized. Within the framework of the
$f$-deformed oscillator approach, we found that under certain circumstances of
exciton-laser interaction the stationary state of the exciton center-of-mass
is a nonlinear coherent state which exhibits the quadrature squeezing.

Acknowledgment The authors wish to thank the Office of Graduate Studies of the
University of Isfahan and Iranian Nanotechnology initiative for their support.

## References

## References

  * [1] R. J. Glauber, Phys. Rev. 130 (1963) 2529; R. J. Glauber, Phys. Rev. 131 (1963) 2766; R. J. Glauber, Phys. Rev. Lett. 10, 84 (1963).
  * [2] J. R. Klauder and B. S. Skagerstam, Coherent states, Applications in Physics and Mathematical Physics (Singapore: World Scientific, 1985).
  * [3] P. A. Perelomov, Generalized Coherent States and Their Applications, (Berlin: Springer, 1986).
S. T. Ali, J-P Antoine and J-P Gazeau, Coherent States, Wavelets and Their
Generalization, (New York: Springer, 2000).

  * [4] V. I. Man’ko, G. Marmo, E. C. G. Sudarshan and F. Zaccaria, Phys. Scr. 55, 528 (1997).
  * [5] V. I. Man’ko, G. Marmo, E. C. G. Sudarshan and F. Zaccaria, in: N. M. Atakishiev(Ed.), Proc. IV Wigner Symp. (Guadalajara, Mexico, July 1995), World Scientific, Singapore, 1996, P. 427; S. Mancini, Phys. Lett. A 233, 291 (1997); S. Sivakumar, J. Phys. A: Math. Gen. 33, 2289 (2000); B. Roy, Phys. Lett. A 249, 25 (1998); H. C. Fu and R. Sasaki, J. Phys. A: Math. Gen. 29, 5637 (1996); R. Roknizadeh and M. K. Tavassoli, J. Phys. A: Math. Gen. 37, 5649 (2004); M. H. Naderi, M. Soltaolkotabi and R. Roknizadeh, J. Phys. A: Math. Gen. 37, 3225 (2004).
  * [6] R. L. deMatos Filho and W. Vogel, Phys. Rev. A 54, 4560 (1996).
  * [7] V. I. Man’ko, G. Marmo, A. Porzio, S. Solimeno and F. Zaccaria, Phys. Rev. A 62, 053407 (2000).
  * [8] M. H. Naderi, M. Soltanolkotabi and R. Roknizadeh, Eur. Phys. J. D. 32, 397 (2005).
  * [9] M. Bagheri Harouni, R. Roknizadeh and M. H. Naderi,
  * [10] A. Mahdifar, R. Roknizadeh and M. H. Naderi, J. Phys. A: Math. Gen. 39, 7003 (2006).
  * [11] Y. X. Liu, C. P. Sun, S. X. Yu and D. L. Zhou, Phys. Rev. A 63,023802 (2001).
  * [12] S. Schmitt-Rink, D. A. B. Miller and D. S. Chemla, Phys. Rev. B 35, 8113 (1987).
  * [13] L. Quiroga and N. F. Johnson, Phys. Rev. Lett. 83, 2270 (1999).
  * [14] G. Chen, N. H. Bonadeo, D. G. Steel, D. Gammon, D. S. Katzer, D. Park and L. J. Sham, Science 289, 1906 (2000);
M. Bayer, P. Hawrylak, K. Hinzer, S. Fafard, D. Korkushiski, Z. R. Wasilewski,
O. Stern and A. Forchel, Science 291, 451 (2001).

  * [15] Y. X. Liu, S. K. Özdemir, M. Koashi and N. Imoto, Phys. Rev. A 65, 042326 (2002).
  * [16] S-K Hong, S. W. Nam and K-H Yeon, Phys. Rev. B 76, 115330 (2007).
  * [17] X. Zou, K. Pahlke and W. Mathis, Phys. Rev. A 68, 034306 (2003).
  * [18] Ka-Di Zhu, Z-J Wu, X-Z Yuam and Zhen, Phys. Rev. B 71, 135312 (2005).
  * [19] F. Tassone and Y. Yamamoto, Phys. Rev. B 59, 1003 (1999).
  * [20] I. V. Bondarev, S. A. Maksimenko, G. Ya. Slepyan, I. L. Krestnikov and A. Hoffmann, Phys. Rev. B 68, 073310 (2003).
  * [21] E. Hanamura, Phys. Rev B 37, 1273 (1988).
  * [22] A. D’Andrea and R. Del Sole, Phys. Rev B 41, 1413 (1990);
S. Jazirl, G. Bastard and R. Bennaceur, Semicond. Sci. Tech. 8, 670 (1993).

  * [23] M. Nirmal, D. J. Norris, M. Kuno, M. G. Bawendi, Al. L. Efros and M. Rosen, Phys. Rev. Lett. 75, 3728 (1995).
  * [24] M. Sugawara, Phys. Rev. B 51, 10743(1995);
M. Califano, A. Franceschetti, and A. Zunger, Phys. Rev. B 75, 115401 (2007).

  * [25] P. Machikowski and L. Jacak, Phys. Rev. B 71, 115309 (2005).
  * [26] W. Vogel and R. L. de Matos Filho, Phys. Rev. A 52,4214(1995); D. M. Meekhof, C. Monroe, B. E. King, W. M. Itano and D. J. Wineland, Phys. Rev. Lett. 76, 1796(1996).
  * [27] R. P. Feynman, Phys. Rev. 84, 108 (1951).
  * [28] A. Wojs, P. Hawrylak, S. Fafard and L. Jacak, Phys. Rev. B 54, 5604 (1996).
  * [29] M. O. Scully and M. S. Zubairy, Quantum optics, (Cambridge University Press, Cambridge, (1997)).

Figure 1: Plots of squeezing versus ration $\frac{r}{a_{B}}$. Solid line
correspond to $\frac{\Omega_{0}}{\Omega_{1}}=0.5$, dash line correspond to
$\frac{\Omega_{0}}{\Omega_{1}}=0.2$. In both plots Lamb-Dick parameter is
chosen as $\kappa=0.3$.


# Black holes surrounded by generic matter distributions:
polar perturbations and energy flux

Nicholas Speeney<EMAIL_ADDRESS>William H. Miller III Department of Physics
and Astronomy, Johns Hopkins University,
3400 N. Charles Street, Baltimore, Maryland, 21218, USA  Emanuele Berti
<EMAIL_ADDRESS>William H. Miller III Department of Physics and Astronomy, Johns
Hopkins University,
3400 N. Charles Street, Baltimore, Maryland, 21218, USA  Vitor Cardoso
<EMAIL_ADDRESS>Niels Bohr International Academy, Niels Bohr
Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  CENTRA, Departamento de
Física, Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Avenida
Rovisco Pais 1, 1049-001 Lisboa, Portugal  Yukawa Institute for Theoretical
Physics, Kyoto University, Kyoto  Andrea Maselli<EMAIL_ADDRESS>Gran
Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy  INFN, Laboratori
Nazionali del Gran Sasso, I-67100 Assergi, Italy

###### Abstract

We develop a numerical approach to compute polar parity perturbations within
fully relativistic models of black hole systems embedded in generic,
spherically symmetric, anisotropic fluids. We apply this framework to study
gravitational wave generation and propagation from extreme mass-ratio
inspirals in the presence of several astrophysically relevant dark matter
models, namely the Hernquist, Navarro-Frenk-White, and Einasto profiles. We
also study dark matter spike profiles obtained from a fully relativistic
calculation of the adiabatic growth of a BH within the Hernquist profile, and
provide a closed-form analytic fit of these profiles. Our analysis completes
prior numerical work in the axial sector, yielding a fully numerical pipeline
to study black hole environmental effects. We study the dependence of the
fluxes on the DM halo mass and compactness. We find that, unlike the axial
case, polar fluxes are not adequately described by simple gravitational-
redshift effects, thus offering an exciting avenue for the study of black hole
environments with gravitational waves.

##  I Introduction

Astrophysical black hole (BH) binaries are, in general, nonvacuum systems. The
environment surrounding a BH binary is expected to play an important role in
the generation and propagation of gravitational waves (GWs) by modifying the
orbital trajectory of the binary and the GW phase Barausse  _et al._ (2014);
Yunes  _et al._ (2011); Cardoso and Maselli (2020). These modifications are
particularly interesting because the properties of the environment are encoded
in the gravitational radiation. The relation between the GW signal and the
properties of the environment has been studied in many cases of interest,
including accretion disks Yunes  _et al._ (2011); Speri  _et al._ (2022), dark
matter (DM) clouds and overdensities Eda  _et al._ (2013, 2015); Bošković  _et
al._ (2018); Annulli  _et al._ (2020); Coogan  _et al._ (2022); Traykova  _et
al._ (2021); Cole  _et al._ (2022); Speeney  _et al._ (2022), bosonic clouds
Macedo  _et al._ (2013); Baumann  _et al._ (2022); Barsanti  _et al._ (2022a,
b); Vicente and Cardoso (2022); Duque  _et al._ (2023), etcetera.

All of these prior studies are based on simplifying assumptions. Very often,
the effect of the environment is modeled by including Newtonian corrections to
the GW quadrupole formula (see e.g. Eda  _et al._ (2013, 2015)). In the few
cases in which the environment is treated relativistically Speeney  _et al._
(2022); Vicente and Cardoso (2022), the orbital trajectory and GW emission are
modeled within the post-Newtonian (PN) approximation, a perturbative expansion
in the ratio $v/c$ of the binary’s orbital velocity to the speed of light.

Extreme mass-ratio inspirals (EMRIs) of stellar-mass objects with mass $m_{p}$
around massive BHs, with mass ratios $q=m_{p}/M_{\rm
BH}\sim\mathcal{O}(10^{-5}-10^{-8})$, can spend $\sim 10^{5}$ orbital cycles
in the low-frequency regime accessible to space-based detectors such as the
Laser Interferometer Space Antenna (LISA) Amaro-Seoane  _et al._ (2017), and
they are expected to lead to the best constraints on environmental effects
Cardoso and Maselli (2020).

The approximations commonly used in the treatment of environmental effects are
not accurate enough for EMRIs. These binaries are not well modeled within the
PN framework, which is more amenable to comparable mass systems Yunes and
Berti (2008); Zhang  _et al._ (2011); Sago  _et al._ (2016). In addition,
these studies do not consistently include the backreaction effect of the
environment on the background solution (which is also expected to modify the
orbital dynamics and the GW signal) in a fully relativistic manner. It is
therefore imperative to correctly model the environment and the orbital
dynamics as precisely as possible to fully capture secular effects which may
accumulate throughout the orbit.

Recent work has established a fully relativistic approach to address the need
for accurate modeling of EMRIs embedded in astrophysical environments. This
program was initiated in Ref. Cardoso  _et al._ (2022a), which employed an
extension of the Einstein cluster model to find relativistic solutions for the
background spacetime and to compute observable quantities. Later work
established a complete framework to study the coupled gravitational and fluid
perturbations of either parity (axial or polar) for generic, spherically
symmetric density distributions affecting EMRI inspirals Cardoso  _et al._
(2022b). While the framework developed in these works is generic, the analysis
focused on the Hernquist distribution Hernquist (1990), which yields closed-
form analytic results for the background metric. Further progress was made in
Ref. Figueiredo  _et al._ (2023), which developed a fully numerical pipeline
to treat generic density distributions in the axial sector.

In this paper we complete this research program by extending the analysis of
Ref. Figueiredo  _et al._ (2023) to include polar perturbations. Polar
perturbations in the presence of generic density profiles can be reduced to 5
frequency-domain ordinary differential equations that couple the variables
describing perturbations of the fluid with the variables describing
gravitational perturbations. This is different from the axial case, in which
the gravitational and fluid variables are decoupled, and therefore they can be
analyzed separately. The framework we develop here is valid for any
spherically symmetric distribution surrounding the central BH, but we focus on
a set of specific DM models (including the Hernquist, Navarro-Frenk-White, and
Einasto profiles) which are commonly used to describe the DM environment
affecting the evolution of the EMRI. We compute the polar GW energy flux
multipoles for various DM models and configurations of the halo, and then we
compare these calculations with previously computed axial multipoles.
Throughout the paper we use geometrical units ($G=c=1$).

##  II Background and Density Profiles

In this section we review the formalism used to construct relativistic BH
solutions embedded within a DM halo. We consider a spherically symmetric,
static spacetime, described by the line element

$\displaystyle ds^{2}=$ $\displaystyle g^{(0)}_{\mu\nu}dx^{\mu}dx^{\nu}$
$\displaystyle=$
$\displaystyle-a(r)dt^{2}+\frac{dr^{2}}{1-2m(r)/r}+r^{2}d\Omega^{2}\ ,$ (1)

which solves the nonvacuum Einstein equations

$G_{\mu\nu}=8\pi T^{\text{env}}_{\mu\nu}\ .$ (2)

Following the Einstein cluster prescription Einstein (1939); Geralico  _et
al._ (2012), the properties of the environment are described by an anisotropic
stress-energy tensor
$(T^{\text{env}})^{\mu}{{}_{\nu}}=\text{diag}(-\rho,0,P_{t},P_{t})$, where
$\rho(r)$ and $P_{t}(r)$ are the density and the tangential pressure of the
matter distribution. The continuity equation determines the mass profile
through the relation $m^{\prime}(r)=4\pi\rho(r)$. The metric functions $a(r)$
and $P_{t}(r)$ are found from the $rr$ component of Eqs. (2) and from the
Bianchi identities, respectively:

$\frac{a^{\prime}(r)}{a(r)}=\frac{2m(r)/r}{r-2m(r)}\quad\ ,\ \quad
P_{t}(r)=\frac{m(r)/2}{r-2m(r)}\rho(r)\ .$ (3)

The solutions of Eqs. (3) and the function $m(r)$, for a given profile
$\rho(r)$, fully specify the background spacetime. These equations can be
integrated numerically Figueiredo  _et al._ (2023).

The form of $m(r)$ and $a(r)$ also determines the properties of particle
geodesics in the spacetime. The energy and the angular momentum of a particle
on a circular orbit with radius $r_{p}$ are given by

$E_{p}=\left[\frac{r-2m(r)}{r-3m(r)}a(r)\right]^{1/2}_{r=r_{p}}\
,~{}~{}~{}L_{p}=\left[\frac{m(r)}{r-3m(r)}\right]^{1/2}_{r=r_{p}}\ ,$ (4)

while the orbital frequency is $\Omega_{p}=a(r_{p})L_{p}/r_{p}^{2}E_{p}$.
Hereafter, without loss of generality, we assume that the orbits are
equatorial and set $\theta=\pi/2$.

With the background spacetime solution in hand, we can now examine the orbital
dynamics of EMRIs and assess the relevance of the environment on the resulting
GW emission. At the leading dissipative order, the secondary can be treated as
a particle of mass $m_{p}$ inducing small perturbations of the background
metric and of the stress-energy tensor, such that

$\displaystyle g_{\mu\nu}=g_{\mu\nu}^{(0)}+g_{\mu\nu}^{(1)}\quad\ ,\quad
T_{\mu\nu}^{\text{env}}=T_{\mu\nu}^{(0)\text{env}}+T_{\mu\nu}^{(1)\text{env}}\
,$ (5)

and

$G_{\mu\nu}^{(1)}=8\pi T_{\mu\nu}^{(1)\text{env}}+8\pi T_{\mu\nu}^{p}\,,$ (6)

where $T_{\mu\nu}^{p}$ is the stress-energy tensor of the particle. The metric
and matter perturbations can be decomposed in terms of axial and polar tensor
spherical harmonics, yielding decoupled radial equations due to the symmetries
of the background spacetime Regge and Wheeler (1957); Zerilli (1970); Lindblom
and Detweiler (1983).

The axial and polar equations have been derived in Ref. Cardoso  _et al._
(2022b), and for brevity we do not reproduce them here. The axial
perturbations have been extensively investigated in Ref. Figueiredo  _et al._
(2023) for both analytical and numerical backgrounds, assuming different DM
profiles. The polar sector has been studied in Ref. Cardoso  _et al._ (2022b)
for the Hernquist model, which has the advantage that the background spacetime
can be written down in a closed analytic form.

###  II.1 Dark matter density profiles

In this work we focus on polar perturbations, with the goal of extending
previous calculations through a fully numerical pipeline capable of describing
generic DM profiles.

Figure 1: Mass profiles for some representative DM configurations considered
in this work. The NFW model is truncated at the cutoff radius $r_{c}=5a_{0}$.
The Hernquist, NFW and Einasto profiles are obtained by the “cutoff factor”
recipe $\rho(r)\rightarrow(1-2M_{\text{BH}}/r)\rho(r)$. The “DM spike” profile
is obtained with the fully relativistic calculation of the adiabatic growth of
a BH within a Hernquist profile from Refs. Sadeghian  _et al._ (2013); Speeney
_et al._ (2022), as described in the main text. Going from the top to the
bottom panels, the DM halo mass grows relative to the central BH mass
($M_{\text{Halo}}/M_{\text{BH}}=1,\,10,\,100$), while the compactness of the
halo decreases going from the left to the right panels
($M_{\text{Halo}}/a_{0}=0.1,\,0.01,\,0.001$). The vertical lines correspond to
a representative orbital radius $r_{p}=7.9456M_{\text{BH}}$ that will be used
in Table 1 below to compare GW energy fluxes.

We consider a selection of DM distributions proposed to model cold DM profiles
within galaxies. Among them, the Hernquist Hernquist (1990) and Navarro-Frenk-
White (NFW) Navarro  _et al._ (1997) distributions, which have been widely
studied and compared to both observations and N-body simulations, can be
parametrically described as

$\rho(r)=\rho_{0}(r/a_{0})^{-\gamma}[1+(r/a_{0})^{\alpha}]^{(\gamma-\beta)/\alpha}\
,$ (7)

where $a_{0}$ is the scale radius of the halo and $\rho_{0}$ is a scale factor
proportional to the density through the relation
$\rho_{0}=2^{(\beta-\gamma)/\alpha}\,\rho(a_{0})$. The exponents $\gamma$ and
$\beta$ control the slope of the distribution at small and large radii,
respectively. The values of $a_{0}$ and $\alpha$ determine the location and
“sharpness” of the change in slope. By fixing $(\alpha,\beta,\gamma)=(1,4,1)$
and $(\alpha,\beta,\gamma)=(1,3,1)$ we recover the Hernquist profile and the
NFW profile, respectively. Note that the total mass of the NFW distribution is
logarithmically divergent. Therefore, we define a cutoff radius $r_{c}$ such
that $M_{\text{Halo}}(r>r_{c})=0$ (as in, for example, Ref. Figueiredo  _et
al._ (2023)).

We also consider a third, observationally driven model: the Einasto profile,
commonly used to study the galaxy luminosity distribution Einasto (1965, 1969)
and such that

$\rho(r)=\rho_{e}\exp\left\\{-d_{n}[(r/r_{e})^{1/n}-1]\right\\}\ ,$ (8)

with $n=6$ and $d_{n}=53/3$ Graham  _et al._ (2006); Prada  _et al._ (2006).
Here $\rho_{e}$ is the density at radius $r_{e}$, which defines a volume
containing half of the halo mass. Throughout this work, we set $r_{e}=a_{0}$.

The DM density profile arising from the (adiabatic) accretion growth of a
nonrotating BH sitting at the core is expected to develop an overdense cusp
followed by a sharp cutoff close to the horizon Gondolo and Silk (1999). A
general relativistic analysis predicts that the density falls off at
$r=4M_{\text{BH}}$ Sadeghian  _et al._ (2013). In our numerical profiles, we
mimic this behavior by considering two different approaches. The first
approach consists of following the strategy adopted in Ref. Figueiredo  _et
al._ (2023), i.e., we multiply the density of the DM models discussed above by
a “cutoff factor” such that $\rho(r)\rightarrow(1-2M_{\text{BH}}/r)\rho(r)$.
In the second approach we consider the fully relativistic “DM spike” models
investigated by some of us in previous work Speeney  _et al._ (2022). In that
work, following Ref. Sadeghian  _et al._ (2013), we derived semianalytic fits
valid for $r\lesssim a_{0}$ for spiked profiles based on the Hernquist and NFW
distributions, and used these fits to assess the impact of the overdensities
on the GW emission of binary sources.

From now on we will focus on the Hernquist subclass of the DM spikes, i.e., we
consider the adiabatic growth of a spike starting from Eq. (7) with
$(\alpha,\beta,\gamma)=(1,4,1)$. This is because the Hernquist and NFW DM
spikes are very similar at small values of $r$ (see Fig. 1) and they mostly
differ at large values of $r$, which are not of interest for our present
goals.

We further improve the fits of the Hernquist DM spike profile from Ref.
Speeney  _et al._ (2022) and extend their domain of validity to large radii
(i.e., to the domain $r\gg a_{0}$). The result is of the form

$\rho(r)=\frac{M_{\text{Halo}}}{M_{\text{BH}}^{2}}\left(\int_{4M_{\text{BH}}}^{r_{c}}4\pi
r^{2}\bar{\rho}(r)dr\right)^{-1}\bar{\rho}(r)$ (9)

with

$\bar{\rho}(r)=\left(1-\frac{4M_{\text{BH}}}{r}\right)^{\alpha}\left(\frac{rM_{\text{BH}}}{M_{\text{Halo}}a_{0}}\right)^{\beta}\left(1+\frac{rM_{\text{BH}}}{M_{\text{Halo}}a_{0}}\right)^{\gamma}\
,$ (10)

where $\alpha=2.366$, $\beta=-2.320$, and $\gamma=-1.370$. The functional
dependence of $\bar{\rho}(r)$ on $M_{\text{Halo}}$ and $a_{0}$ was determined
empirically by studying the behavior of an ensemble of DM profiles computed
numerically for selected values of $M_{\text{Halo}}$ and $a_{0}$ Sadeghian
_et al._ (2013). The values of the exponents $(\alpha,~{}\beta,~{}\gamma)$
were then found by fitting with Mathematica a numerical spike with
$M_{\text{Halo}}=10^{4}M_{\text{BH}}$ and $M_{\text{Halo}}/a_{0}=0.001$.

In the calculation of the GW fluxes (see Sec. III below) we enforce a cutoff
radius of $r_{c}=100M_{\text{Halo}}a_{0}/M_{\text{BH}}$ such that
$M_{\text{Halo}}(r>r_{c})=0$ for the DM spikes. Equation (9) reproduces the
numerical results with $2\%$ accuracy in the range $r_{\text{ISCO}}<r<r_{c}$
for the values of $(M_{\text{Halo}}/M_{\text{BH}},M_{\text{Halo}}/a_{0})$
which we investigated in this work. The accuracy of the fit generally improves
as both $M_{\text{Halo}}$ and $a_{0}$ increase.

In Fig. 1 we plot some of the DM density distributions described so far for
representative choices of the halo parameters. The Hernquist and NFW profiles
are always rather similar, except at the highest values of the halo mass and
smallest values of the compactness (bottom right panel), where the density of
the Hernquist model is larger than the NFW density by about one order of
magnitude. When the halo is very compact (left panels) the Einasto profiles
have the highest DM density close to the BH, orders of magnitude larger than
the DM spike profiles. The trend reverses as we move to the less compact
profiles on the right: for the largest halo masses and for the smallest values
of the compactness (bottom right panel) the DM spike contains by far the
largest amount of DM at the orbital radii of interest for EMRIs.

##  III Numerical Procedure

We outline here the numerical procedure followed to compute the polar
gravitational perturbations.

We first make a specific choice of $\rho(r)$, and then we solve for the
background metric component $a(r)$ and for the mass function $m(r)$. This is
accomplished by numerically integrating $m^{\prime}(r)=4\pi r^{2}\rho(r)$ from
$r=2M_{\text{BH}}$ to the cutoff radius $r_{\infty}=10^{11}M_{\text{BH}}$,
which corresponds to our numerical approximation of spatial infinity. We then
plug the numerical solution for $m(r)$ into Eq. (3) to solve for $a(r)$ by
integrating “backwards” from $r_{\infty}$, with the following boundary
condition:

$\displaystyle m(r_{\infty})=$ $\displaystyle M_{\text{BH}}+M_{\text{Halo}}\
,$ $\displaystyle a(r_{\infty})=$ $\displaystyle
1-\frac{2m(r_{\infty})}{r_{\infty}}\ .$ (11)

Here, $M_{\text{Halo}}$ corresponds to the total mass of the environment
surrounding the BH.

The numerical solutions for $a(r)$ and $m(r)$ are used to calculate the
tangential pressure $P_{t}(r)$ from Eq. (3), as well as the energy, angular
momentum, and orbital frequency of particles in circular orbits using Eqs.
(4).

Our next task is to solve the five coupled first-order ODEs governing the
polar sector. These ODEs can be cast in the compact form

$\frac{d\vec{\psi}}{dr}-\boldsymbol{\alpha}\vec{\psi}=\vec{S}\quad\
,\quad\vec{\psi}=\\{K,H_{0},H_{1},W,\delta\rho\\}\ ,$ (12)

where $\\{K,\,H_{0},\,H_{1}\\}$ and $\\{W,\,\delta\rho\\}$ are the
perturbation variables related to the metric and to the fluid, respectively.
The components of the matrix $\boldsymbol{\alpha}$ and of the source vector
$\vec{S}$ are given explicitly in Ref. Cardoso  _et al._ (2022b). We require
the gravitational perturbation variables to satisfy ingoing- and outgoing-wave
boundary conditions at the horizon and at spatial infinity, respectively.
Following Ref. Cardoso  _et al._ (2022b) we assume, for simplicity, that the
fluid variables vanish at the two boundaries, i.e.
$W(2M_{\text{BH}})=\delta\rho(2M_{\text{BH}})=W(r_{\infty})=\delta\rho(r_{\infty})=0$.
The gravitational perturbations at the horizon can be expanded as follows:

$\displaystyle K^{(\text{in})}(r)$ $\displaystyle=e^{-i\omega
r*}\sum_{i=0}^{N_{\text{in}}}K_{H}^{(i)}(r-2M_{\text{BH}})^{i}\ ,$ (13)
$\displaystyle H_{0}^{(\text{in})}(r)$ $\displaystyle=e^{-i\omega
r*}\sum_{i=0}^{N_{\text{in}}}H_{0,H}^{(i)}(r-2M_{\text{BH}})^{i}\ ,$
$\displaystyle H_{1}^{(\text{in})}(r)$ $\displaystyle=e^{-i\omega
r*}\sum_{i=0}^{N_{\text{in}}}H_{1,H}^{(i)}(r-2M_{\text{BH}})^{i}\ ,$

where in our numerical calculations we set $N_{\text{in}}=5$. We similarly
expand the background functions at the horizon,

$\displaystyle m(r)$
$\displaystyle=M_{\text{BH}}+\sum_{i=2}^{N_{\text{in}}}m_{H}^{(i)}(r-2M_{\text{BH}})^{i}\
,$ (14) $\displaystyle a(r)$
$\displaystyle=\sum_{i=1}^{N_{\text{in}}}a_{H}^{(i)}(r-2M_{\text{BH}})^{i}\ ,$

and we use the numerically computed profiles for $m(r)$ and $a(r)$ to
determine the coefficients $a_{H}^{(i)},\,m_{H}^{(i)}$. Given these
coefficients, we can expand Eqs. (12) in powers of $r-2M_{\text{BH}}$ and
solve order-by-order to find the coefficients $K_{H}^{(i)}$, $H_{0,H}^{(i)}$,
$H_{1,H}^{(i)}$. This completely specifies the boundary conditions at the
horizon, modulo the value of $K_{H}^{(0)}$.

$\ell~{}m$ |  Vacuum |  |  Hernquist |  NFW $r_{c}=5a_{0}$ |  Einasto $r_{e}=a_{0}$ |  DM Spike
---|---|---|---|---|---|---
|  |  $\tfrac{M_{\text{Halo}}}{M_{\rm BH}}$$\tfrac{M_{\text{Halo}}}{a_{0}}$ |  0.1 0.001 |  0.1 0.001 |  0.1 0.001 |  0.1 0.001
$2~{}2$ | 1.7068e-4 |  1 100 |  2.4341e-4 1.7038e-4 1.4133e-4 1.7037e-4 |  2.1997e-4 1.7036e-4 1.4396e-4 1.7036e-4 |  8.0265e-4 1.6972e-4 5.8136e-4 1.6833e-4 |  1.8033e-4 1.6997e-4 1.6533e-4 1.7037e-4
$3~{}1$ | 2.2443e-9 |  1 100 |  5.8115e-9 2.1823e-9 2.3212e-9 2.2010e-9 |  4.2975e-10 2.1788e-9 2.1489e-9 2.2013e-9 |  2.3454e-6 2.4672e-9 7.1498e-8 2.1772e-9 |  5.4913e-9 2.4933e-9 3.1469e-9 2.2568e-9
$3~{}3$ | 2.5488e-5 |  1 100 |  3.5564e-5 2.5429e-5 2.1075e-5 2.5448e-5 |  3.1495e-5 2.5429e-5 2.1446e-5 2.5434e-5 |  1.4468e-4 2.5301e-5 7.6403e-6 2.5131e-5 |  2.5779e-5 2.5323e-5 2.4470e-5 2.5415e-5
$4~{}2$ | 2.5855e-9 |  1 100 |  3.1319e-8 2.4966e-9 1.1507e-9 2.5389e-9 |  2.5378e-8 2.5028e-9 1.5796e-9 2.5384e-9 |  1.9289e-7 2.0099e-9 3.3474e-9 2.5045e-9 |  4.6899e-9 2.0384e-9 8.3514e-10 2.4449e-9
$4~{}4$ | 4.7327e-6 |  1 100 |  6.0451e-6 4.7169e-6 3.8716e-6 4.7200e-6 |  5.2628e-6 4.7177e-6 3.9592e-6 4.7204e-6 |  3.6665e-5 4.6776e-6 1.2052e-6 4.6639e-6 |  4.3577e-6 4.6794e-6 4.4617e-6 4.7165e-6
Table 1: Selected polar multipolar components $\dot{E}^{\infty}_{(\ell,m)}$ of
the energy flux at spatial infinity in units of
$q^{2}=m_{p}^{2}/M_{\text{BH}}^{2}$ for a particle orbiting at
$r_{p}=7.9456M_{\text{BH}}$. In each case, we set the width of the Gaussian to
$\sigma=1/50$. For nonvacuum fluxes, the four entries in each cell span two
values of the DM halo mass and compactness:
$M_{\text{Halo}}/M_{\text{BH}}=1,\,100$, $M_{\text{Halo}}/a_{0}=0.1,\,0.001$.

We perform a similar procedure at spatial infinity, where we expand the
gravitational perturbations as

$\displaystyle K^{(\text{out})}(r)$ $\displaystyle=e^{i\omega
r*}\sum_{i=0}^{N_{\text{out}}}\frac{K_{\infty}^{(i)}}{r^{i}}\ ,$ (15)
$\displaystyle H_{0}^{(\text{out})}(r)$ $\displaystyle=e^{i\omega
r*}\sum_{i=0}^{N_{\text{out}}}\frac{H_{0,\infty}^{(i)}}{r^{i}}\ ,$
$\displaystyle H_{1}^{(\text{out})}(r)$ $\displaystyle=e^{i\omega
r*}\sum_{i=0}^{N_{\text{out}}}\frac{H_{1,\infty}^{(i)}}{r^{i}}\ ,$

and the metric functions as

$m(r)=\sum_{i=0}^{N_{\text{out}}}\frac{m_{\infty}^{(i)}}{r^{i}}\quad\ ,\quad
a(r)=\sum_{i=0}^{N_{\text{out}}}\frac{a_{\infty}^{(i)}}{r^{i}}\,.$ (16)

Once again we set $N_{\text{out}}=5$. By matching Eqs. (16) with the numerical
profiles of $a(r)$ and $m(r)$ we can determine the values of
$a_{\infty}^{(i)}$ and $m_{\infty}^{(i)}$. Then we expand Eqs. (12) in powers
of $1/r$ and solve order-by-order for the coefficients $K_{\infty}^{(i)}$,
$H_{0,\infty}^{(i)}$, $H_{1,\infty}^{(i)}$. Fixing $K_{\infty}^{(0)}=1$
completely specifies the boundary values for the gravitational variables at
spatial infinity.

We solve the five ODEs in the frequency domain for $K(r)$, $H_{0}(r)$,
$H_{1}(r)$, $W(r)$, and $\delta\rho(r)$ using standard Mathematica numerical
routines for differential equations. We assume the boundary conditions listed
in Eq. (13) at $r=2M_{\text{BH}}$, and those listed in Eq. (15) at
$r=r_{\infty}$. The source term in Eq. (12) depends on the particle location
through a delta function $\delta(r-r_{p})$, which we approximate by a Gaussian
$\delta(r-r_{p})\sim\text{exp}[-(r-r_{p})^{2}/2\sigma^{2}]/\sqrt{2\pi}\sigma$,
as in Refs. Figueiredo  _et al._ (2023); Cardoso  _et al._ (2022b). We also
fix the radial and transverse sound speeds which appear in the matrix
$\boldsymbol{\alpha}$ of Eq. (12) to be $(c_{s,t},c_{s,r})=(0,9/10)$ (see
Cardoso  _et al._ (2022b) for a discussion). The gravitational perturbations
are then extracted at $r_{\text{obs}}=\max[3\times 10^{3}/\omega,\,2a_{0}]$.
We make this choice for all the cases we consider, except for profiles with
$M_{\text{Halo}}=100$, $a_{0}=10^{5}M_{\text{BH}}$, which require a larger
extraction radius $r_{\text{obs}}=2\times 10^{4}/\omega$ to achieve sufficient
numerical accuracy.

For a given orbital radius $r_{p}$, the solution of Eqs. (12) is found by a
shooting method, with the shooting parameter $K_{H}^{(0)}$ determined by
requiring that the perturbations and their derivatives must be continuous at
some radius $r_{\textnormal{obs}}$:

$\displaystyle\lim_{r\to r_{\text{obs}}}\bigg{[}K_{\text{num}}(r)$
$\displaystyle K^{\prime}_{\text{exp}}(r)$
$\displaystyle-K^{\prime}_{\text{num}}(r)K_{\text{exp}}(r)\bigg{]}=0\ ,$ (17)

In this equation, the quantity $K_{\text{num}}$ refers to the numerically
computed $K(r)$, while $K_{\text{exp}}$ refers to the expansion of Eq. (15).

Figure 2: Solid lines: $\ell=m=2$ mode GW energy fluxes, in units of
$q^{2}=m_{p}^{2}/M_{\text{BH}}^{2}$, in the presence of a DM spike matter
background. The frequency range on the x axis corresponds to particle orbital
radii ranging from $r_{p}=50M_{\text{BH}}$ to $r_{p}=6M_{\text{BH}}$. In the
left panel we compare DM spike halos with the same value of
$M_{\text{Halo}}/M_{\text{BH}}=100$ but different compactness
$M_{\text{Halo}}/a_{0}$, while in the right panel we compare DM spike halos
with fixed compactness $M_{\text{Halo}}/a_{0}=0.01$ but different values of
$M_{\text{Halo}}/M_{\text{BH}}$. For reference, the dashed line shows the
$\ell=m=2$ multipolar component of the vacuum energy flux.

Figure 3: Solid lines: relative difference between the $\ell=m=2$ GW energy
fluxes, in units of $m_{p}^{2}/M_{\text{Halo}}^{2}$, in vacuum ($\dot{E}_{v}$)
and in the presence of a matter background ($\dot{E}_{m}$) for different total
halo masses, fixing the compactness to $M_{\text{Halo}}/a_{0}=1/100$. This
frequency range corresponds to particle orbital radii ranging from
$r_{p}=50M_{\text{BH}}$ to $r_{p}=6M_{\text{BH}}$. Dashed lines: relative
difference between $\dot{E}_{m}$ and the vacuum fluxes redshifted according to
Eq. (19). Dot-dashed lines: relative difference between $\dot{E}_{m}$ and the
post-Newtonian fluxes of Eq. (20).

Figure 4: Solid lines: relative difference between the $\ell=m=2$ GW energy
fluxes in vacuum ($\dot{E}_{v}$) and in the presence of matter ($\dot{E}_{m}$)
for different compactness, fixing the halo mass to
$M_{\text{Halo}}=100M_{\text{BH}}$. This frequency range corresponds to
particle orbital radii ranging from $r_{p}=50M_{\text{BH}}$ to
$r_{p}=6M_{\text{BH}}$. Dashed lines: relative difference between
$\dot{E}_{m}$ and the vacuum fluxes redshifted according to Eq. (19). Dot-
dashed lines: relative difference between $\dot{E}_{m}$ and the post-Newtonian
fluxes of Eq. (20).

We can finally compute the GW fluxes at infinity for different multipolar
components $(\ell,\,m)$ from the relation

$\dot{E}^{\infty}_{(\ell,m)}/q^{2}=\lim_{r\to
r_{\text{obs}}}\frac{1}{32\pi}\frac{(\ell+2)!}{(\ell-2)!}|K(r)|^{2}\,,~{}~{}(\ell+m~{}\text{even})\,,$
(18)

where $q=m_{p}/M_{\text{BH}}$. The codes used to calculate the background
quantities and the polar fluxes are publicly available online SGR .

##  IV Results

Examples of the GW fluxes $\dot{E}^{\infty}_{(\ell,m)}$ computed in this way
are shown in Table 1. All fluxes are evaluated at a reference particle orbital
radius $r_{p}=7.9456M_{\text{BH}}$ to facilitate comparison with the vacuum
fluxes computed in Ref. Martel (2004). For nonvacuum fluxes and for a given DM
profile, the four entries in each cell of the table span two values of the DM
halo mass and compactness: $M_{\text{Halo}}/M_{\text{BH}}=1,\,100$ and
$M_{\text{Halo}}/a_{0}=0.1,\,0.001$. In other words, we consider the four most
extreme cases shown in Fig. 1: the top (bottom) entries in each cell
correspond to a small (large) halo mass, while the left (right) entries
correspond to a large (small) compactness.

Let us first focus on halos with small compactness,
$M_{\text{Halo}}/a_{0}=0.001$ (right entries in each cell of Table 1). In this
case, changes in the halo mass $M_{\text{Halo}}/M_{\text{BH}}$ have nearly no
effect: in fact the fluxes are very similar for any choice of the DM halo
profile, and they differ by at most a few percent from the corresponding
vacuum fluxes. This makes sense when we look at the DM distributions plotted
in the right panels of Fig. 1: if the compactness is small, most of the DM
density is at large radii. The particle evolution is most sensitive to the DM
density inside the orbital radius (marked by vertical lines in Fig. 1), which
is always very low for “dilute” DM profiles with
$M_{\text{Halo}}/a_{0}=0.001$, with peak values of
$M_{\text{BH}}^{2}\rho\lesssim 10^{-6}$. Therefore, the DM distribution close
to the BH has an almost negligible effect on the GW energy flux and on the
orbital evolution of the particle.

Consider next halos with large compactness, $M_{\text{Halo}}/a_{0}=0.1$ (left
entries in each cell of Table 1). We now observe that the multipolar
contributions to the GW energy fluxes are typically larger for small values of
the halo mass ($M_{\text{Halo}}/M_{\text{BH}}=1$) than for large values of the
halo mass ($M_{\text{Halo}}/M_{\text{BH}}=100$). Once again, this can be
understood by looking at the DM distributions plotted in the left panels of
Fig. 1. Consider for reference the peak of the Einasto DM density distribution
(dotted green line): it decreases from $M_{\text{BH}}^{2}\rho\sim 10^{-3}$ for
$M_{\text{Halo}}/M_{\text{BH}}=1$, to $M_{\text{BH}}^{2}\rho\sim 10^{-4}$ for
$M_{\text{Halo}}/M_{\text{BH}}=100$. This is because at fixed compactness, the
more massive halos are more diluted, and extend out to larger orbital radii,
so the amount of DM within the particle orbit decreases as the total halo mass
$M_{\text{Halo}}/M_{\text{BH}}$ increases. This explains why, at fixed
compactness, the multipolar contributions to the GW energy fluxes are
typically smaller for larger values of the halo mass. There are some
exceptions to this general rule, such as the $(\ell,\,m)=(3,\,1)$ multipole
for the NFW profile and the $(\ell,\,m)=(4,\,4)$ multipole for the DM spike,
but as we will discuss below these are largely subdominant and, ultimately,
they do not contribute much to the energy fluxes and to the orbital evolution
of the particle.

In Fig. 2 we focus on the fiducial “DM spike” model, which was computed from a
fully relativistic calculation of the adiabatic growth of a Hernquist profile
Sadeghian  _et al._ (2013); Speeney  _et al._ (2022). The trends visible in
this figure are consistent with Table 1. The left panel shows that the
quadrupolar energy flux decreases as the DM halo become less compact at fixed
values of $M_{\text{Halo}}/M_{\text{BH}}$: this is easy to understand, because
the strong-field dynamics of the particle is dominated by the BH when the DM
profile is more dilute (i.e., when the DM density close to the BH decreases).
The right panel shows that the quadrupolar energy flux is smaller for larger
values of $M_{\text{Halo}}/M_{\text{BH}}$ at fixed compactness, again because
more dilute DM profiles have less DM within the particle orbit. These are the
same trends we observed at $r_{p}=7.9456M_{\text{BH}}$ in Table 1, and they
apply to all values of $r_{p}/M_{\text{BH}}$ (or $M_{\text{BH}}\Omega_{p}$).

However, Fig. 2 allows us to better appreciate the dependence of these effects
on the particle radius $r_{p}/M_{\text{BH}}$. For example, the right panel
shows that the effect of gravitational redshift is still important. When we
fix $M_{\text{Halo}}/a_{0}=0.01$, $a_{0}$ changes by two orders of magnitude
as we vary $M_{\text{Halo}}/M_{\text{BH}}$, so the emitted GWs propagate
through very different density profiles. For $r_{p}=50M_{\text{BH}}$
($M_{\text{BH}}\Omega_{p}\sim 0.003$) the difference between the three cases
is more pronounced, because the particle is (roughly) located halfway through
the DM distribution when $M_{\text{Halo}}/M_{\text{BH}}=1$ (so that
$a_{0}=100M_{\text{Halo}}=100M_{\text{BH}}$). On the contrary, when
$r_{p}=6M_{\text{BH}}$ ($M_{\text{BH}}\Omega_{p}\sim 0.05$), the particle is
located so deep within the DM distribution that any differences induced by
gravitational redshift saturate, and the three fluxes become nearly
indistinguishable.

The solid lines in the top-left panel of Fig. 3 show the relative difference
between the $\ell=m=2$ multipolar contribution to the energy flux in the
frequency domain for the Hernquist distribution and the corresponding vacuum
result for different DM halos with fixed compactness
$M_{\textnormal{Halo}}/a_{0}=0.01$. The plot spans values of the orbital
frequencies such that the circular orbit radius of the particle varies between
$r_{p}=50M_{\text{BH}}$ and $r_{p}=6M_{\text{BH}}$. Results are in excellent
agreement with those in Fig. 3 of Ref. Cardoso  _et al._ (2022b), which made
use of closed-form analytic expressions for the metric functions $m(r)$ and
$a(r)$. In fact, the energy fluxes shown here are smoother compared to those
of Ref. Cardoso  _et al._ (2022b), mainly because of the larger value of
$r_{\text{out}}$ considered in this paper to increase numerical accuracy.

In the remaining panels of Fig. 3 we compare vacuum and matter fluxes for the
NFW, Einasto, and DM spike profiles. The solid lines in different panels shows
trends consistent with the DM profiles plotted in Fig. 1. For example, by
comparing the solid lines in the two top panels we observe minimal differences
between the Hernquist and NFW fluxes, which is hardly surprising given the
similarity of the corresponding DM profiles. The GW flux in the presence of DM
differs more from the vacuum case when we consider the DM spike and (even more
so) the Einasto profiles. The reason can be easily understood by a glance at
the middle column of Fig. 1: the DM density in the regime of interest is
highest for the Einasto profile, followed by the DM spike profile and (at some
distance) by the Hernquist and NFW profiles.

The dashed lines in each panel of Fig. 3 show the relative difference between
the numerical values of $\dot{E}^{\infty}_{(2,2)}$ and the redshifted vacuum
fluxes obtained through the scaling Cardoso  _et al._ (2022b)

$\Omega_{p}\rightarrow\Omega_{p}/\gamma\quad\
,\quad\omega\rightarrow\omega/\gamma\quad\ ,\quad m_{p}\rightarrow\gamma
m_{p}\ ,$ (19)

where the redshift factor is given by $\gamma=1-M_{\text{Halo}}/a_{0}$ for
Hernquist, NFW, and DM spike distributions, and by
$\gamma=1-M_{\text{Halo}}/r_{e}$ for the Einasto model. Interestingly, the
nonvacuum fluxes are not well captured by the scaling (19), even for small
values of the halo mass $M_{\text{Halo}}$. This is a key difference with
respect to the energy fluxes in the axial sector: in that case, this simple
redshift scaling turned out to significantly reduce the difference between the
vacuum energy fluxes and those in the presence of matter Cardoso  _et al._
(2022b); Figueiredo  _et al._ (2023).

Finally, the dot-dashed lines in each panel of Fig. 3 show the relative
difference between the numerical fluxes and the corrections to the vacuum PN
GW fluxes computed in Ref. Speeney  _et al._ (2022), i.e.,

$\dot{E}^{\infty}_{\text{GW}}=\frac{32}{5}(M_{\text{BH}}\Omega_{p})^{10/3}\left(1+\frac{4}{3}\epsilon
q(r)\right)\ .$ (20)

Equation (20) is obtained by expanding the energy balance equation in powers
of the small quantity $\epsilon q(r)=M_{\text{DM}}(r)/M_{\text{BH}}$, where
$M_{\text{DM}}(r)$ is the mass of the halo enclosed in a sphere of radius $r$.
Therefore, not surprisingly, the approximation works best when
$M_{\text{Halo}}/M_{\text{BH}}$ is smallest. Note that at variance with Ref.
Speeney  _et al._ (2022), in Eq. (20) we do not include dynamical friction
terms, because our framework does not accommodate these effects at leading
order in the perturbations. In fact, the agreement with our numerical fluxes
gets worse when we compare them against PN expressions which include dynamical
friction.

Figure 4 is similar to Fig. 3, except now we fix the ratio
$M_{\text{Halo}}/M_{\text{BH}}=100$ and vary the halo compactness:
$M_{\text{Halo}}/a_{0}=0.1,0.01,0.001$. This plot confirms the general trends
observed earlier. The Hernquist and NFW energy fluxes behave in very similar
ways. The Einasto and (to a lesser extent) the DM spike profile produce much
larger deviations from the vacuum fluxes. Finally, and most importantly, the
dephasing induced by DM generally increases with the compactness of the halo.

Figures 3 and 4 show that Eq. (20) is not sufficient to reproduce the
numerical fluxes at the orbital frequencies of interest for EMRIs for any of
the DM profiles we consider. The agreement between Eq. (20) and
$\dot{E}^{\infty}_{2,2}$ improves at large radii and small orbital
frequencies, as expected in any PN treatment, but it gets worse for more
astrophysically realistic (i.e., less compact) halos.

Overall, these results lead to two main conclusions: (i) the polar fluxes
cannot be adequately described by PN expansions, and they cannot be simply
attributed to gravitational redshift effects; (ii) the degeneracies due to the
redshift scaling (19) found in the axial sector are generally broken when we
consider polar perturbations. This observation paves the way for the
interesting possibility of constraining the halo properties through GW
observations with space-based detectors.

Figure 5: Multipolar components of the GW energy flux in the presence of a DM
spike background as a function of the particle’s orbital radius
$r_{p}/M_{\text{BH}}$. We fix the halo mass to
$M_{\text{Halo}}/M_{\text{BH}}=100$ and the compactness to
$M_{\text{Halo}}/a_{0}=0.001$. Solid (dashed) lines refer to polar (axial)
modes.

Finally, in Fig. 5 we show the multipolar components of the GW energy flux as
a function of the orbital radius $r_{p}/M_{\text{BH}}$ for a DM spike model
with $M_{\rm Halo}/a_{0}=0.001$ and $M_{\rm Halo}=100M_{\rm BH}$. For
completeness we also include the $(2,\,1)$ and $(3,\,2)$ axial contributions,
that were first computed in Ref. Figueiredo  _et al._ (2023). Similarly to the
vacuum case, we observe a clear hierarchy in the amplitude of the GW
multipoles. The total energy flux is dominated by the $(2,\,2)$, $(3,\,3)$ and
$(4,\,4)$ modes. Among the subdominant components, the axial $(2,\,1)$
multipole dominates over the polar $(5,\,5)$ multipole at large values of
$r_{p}/M_{\text{BH}}$, while the opposite is true for
$r_{p}/M_{\text{BH}}\lesssim 8.5$. We observe an analogous behavior for the
axial $(3,\,2)$ and polar $(6,\,6)$ modes, with a “crossing radius” around
$r_{p}/M_{\rm BH}\sim 7$. Similar considerations apply to the Hernquist, NFW,
and Einasto models, which are not shown here for brevity.

##  V Conclusions

This work extends the fully relativistic modeling of EMRIs embedded in a DM
environment to include a numerical treatment of the polar sector, thus
completing the program initiated in Refs. Cardoso  _et al._ (2022b);
Figueiredo  _et al._ (2023).

We have used this formalism to study the GW fluxes produced by different DM
density distribution models, and we have found that the differences can be of
order unity or larger with respect to the Schwarzschild fluxes, depending on
our assumptions on the DM density profile.

Some general trends can be identified from Table 1. For halos with large
compactness (say, $M_{\text{Halo}}/a_{0}=0.1$), the multipolar contributions
to the GW energy fluxes are typically large for small values of the halo mass
($M_{\text{Halo}}/M_{\text{BH}}=1$), because at fixed compactness the more
massive halos are more diluted, extending out to larger orbital radii, so the
amount of DM within the particle orbit decreases as
$M_{\text{Halo}}/M_{\text{BH}}$ increases. In the more realistic case of halos
with small compactness, $M_{\text{Halo}}/a_{0}\ll 1$, the fluxes are very
similar for any choice of the DM halo profile and for any value of the halo
mass $M_{\text{Halo}}/M_{\text{BH}}$, differing by at most a few percent from
the corresponding vacuum fluxes. This is because most of the DM density is at
large radii, and therefore the DM distribution close to the BH has an almost
negligible effect on the GW energy flux and on the orbital evolution of the
particle.

In Fig. 2 we focus on the fiducial “DM spike” model, computed from a fully
relativistic calculation of the adiabatic growth of a BH within the Hernquist
profile Sadeghian  _et al._ (2013); Speeney  _et al._ (2022). The quadrupolar
energy flux decreases as the DM halo become less compact at fixed values of
$M_{\text{Halo}}/M_{\text{BH}}$, because the strong-field dynamics of the
particle is dominated by the BH when the DM profile is more dilute (i.e., when
the DM density close to the BH decreases). The energy flux is smaller for
large values of $M_{\text{Halo}}/M_{\text{BH}}$ at fixed compactness for a
similar reason: more dilute DM profiles have less DM within the particle
orbit. Realistic DM profiles are probably in the regime
$M_{\text{Halo}}/a_{0}\ll 1$, and we plan to explore more carefully this
“dilute DM” limit in future work.

As shown in Figs. 3 and 4, the numerical energy fluxes in the polar case
cannot be captured by simply taking into account the gravitational redshift of
the vacuum fluxes (at variance with the axial case), nor by leading-order
corrections to the mass profile combined with PN estimates of the flux.

Our analysis has several limitations that we plan to address in future
studies.

First of all, it is important to accurately model dynamical friction, which
may accumulate large dephasings over the $\mathcal{O}(10^{5}-10^{6})$ orbital
cycles typical of LISA EMRIs Eda  _et al._ (2015); Cardoso and Maselli (2020);
Speeney  _et al._ (2022); Coogan  _et al._ (2022); Cardoso  _et al._ (2019);
Traykova  _et al._ (2021, 2023). Moreover, our model does not properly take
into account matter fluxes. While EMRIs can be correctly described by
stationary configurations with matter profiles that vanish at the horizon and
at spatial infinity, intermediate-mass binaries will be sensitive to feedback
effects. This backreaction will affect the DM profile during the evolution of
the orbit Kavanagh  _et al._ (2020); Coogan  _et al._ (2022); Nichols  _et
al._ (2023). A proper treatment of backreaction requires the inclusion of
higher-order terms in the perturbative expansion. If simple estimates similar
to Eq. (20) hold at very small orbital frequencies, DM halos could produce
effects of phenomenological relevance for pulsar timing arrays Ghoshal and
Strumia (2023).

On the data analysis side, prior works have proposed to use EMRI observations
to distinguish various environmental signatures Cole  _et al._ (2022, 2023);
Barausse  _et al._ (2014, 2015); Speeney  _et al._ (2022); Vicente and Cardoso
(2022). The majority of these works uses a combination of fully relativistic
and PN estimates of environmental effects, but so far there is no consistent
description of the evolution of the BH spacetime embedded in DM during the
inspiral. Our framework can accommodate relativistic effects in the generation
and propagation of GWs. In the future, we plan to develop waveform models
based on these flux calculations and to use them in a Bayesian framework to
estimate the detectability of DM overdensities.

###### Acknowledgements.

N.S. and E.B. are supported by NSF Grants No. AST-2006538, PHY-2207502,
PHY-090003 and PHY-20043, by NASA Grants No. 20-LPS20-0011 and 21-ATP21-0010,
by the John Templeton Foundation Grant 62840, and by the Simons Foundation.
A.M. and E.B. acknowledge support from the Indo-US Science and Technology
Forum through the Indo-US Centre for Gravitational-Physics and Astronomy,
grant IUSSTF/JC-142/2019. This work was supported in part by the Italian
Ministry of Foreign Affairs and International Cooperation”, grant number
PGR01167. We acknowledge support by VILLUM Foundation (grant no. VIL37766) and
the DNRF Chair program (grant no. DNRF162) by the Danish National Research
Foundation. V.C. is a Villum Investigator and a DNRF Chair. V.C. acknowledges
financial support provided under the European Union’s H2020 ERC Advanced Grant
“Black holes: gravitational engines of discovery” grant agreement no.
Gravitas–101052587. Views and opinions expressed are however those of the
author only and do not necessarily reflect those of the European Union or the
European Research Council. Neither the European Union nor the granting
authority can be held responsible for them. This project has received funding
from the European Union’s Horizon 2020 research and innovation programme under
the Marie Sklodowska-Curie grant agreement No 101007855 and No 101131233. This
work was carried out at the Advanced Research Computing at Hopkins (ARCH) core
facility (rockfish.jhu.edu), which is supported by the NSF Grant No.
OAC-1920103.

## References

  * Barausse  _et al._ (2014) E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Yunes  _et al._ (2011) N. Yunes, B. Kocsis, A. Loeb, and Z. Haiman, Phys. Rev. Lett. 107, 171103 (2011), arXiv:1103.4609 [astro-ph.CO] .
  * Cardoso and Maselli (2020) V. Cardoso and A. Maselli, Astron. Astrophys. 644, A147 (2020), arXiv:1909.05870 [astro-ph.HE] .
  * Speri  _et al._ (2022) L. Speri, A. Antonelli, L. Sberna, S. Babak, E. Barausse, J. R. Gair, and M. L. Katz, (2022), arXiv:2207.10086 [gr-qc] .
  * Eda  _et al._ (2013) K. Eda, Y. Itoh, S. Kuroyanagi, and J. Silk, Phys. Rev. Lett. 110, 221101 (2013), arXiv:1301.5971 [gr-qc] .
  * Eda  _et al._ (2015) K. Eda, Y. Itoh, S. Kuroyanagi, and J. Silk, Phys. Rev. D 91, 044045 (2015), arXiv:1408.3534 [gr-qc] .
  * Bošković  _et al._ (2018) M. Bošković, F. Duque, M. C. Ferreira, F. S. Miguel, and V. Cardoso, Phys. Rev. D 98, 024037 (2018), arXiv:1806.07331 [gr-qc] .
  * Annulli  _et al._ (2020) L. Annulli, V. Cardoso, and R. Vicente, Phys. Rev. D 102, 063022 (2020), arXiv:2009.00012 [gr-qc] .
  * Coogan  _et al._ (2022) A. Coogan, G. Bertone, D. Gaggero, B. J. Kavanagh, and D. A. Nichols, Phys. Rev. D 105, 043009 (2022), arXiv:2108.04154 [gr-qc] .
  * Traykova  _et al._ (2021) D. Traykova, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, Phys. Rev. D 104, 103014 (2021), arXiv:2106.08280 [gr-qc] .
  * Cole  _et al._ (2022) P. S. Cole, G. Bertone, A. Coogan, D. Gaggero, T. Karydas, B. J. Kavanagh, T. F. M. Spieksma, and G. M. Tomaselli, (2022), arXiv:2211.01362 [gr-qc] .
  * Speeney  _et al._ (2022) N. Speeney, A. Antonelli, V. Baibhav, and E. Berti, Phys. Rev. D 106, 044027 (2022), arXiv:2204.12508 [gr-qc] .
  * Macedo  _et al._ (2013) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Astrophys. J. 774, 48 (2013), arXiv:1302.2646 [gr-qc] .
  * Baumann  _et al._ (2022) D. Baumann, G. Bertone, J. Stout, and G. M. Tomaselli, Phys. Rev. D 105, 115036 (2022), arXiv:2112.14777 [gr-qc] .
  * Barsanti  _et al._ (2022a) S. Barsanti, A. Maselli, T. P. Sotiriou, and L. Gualtieri, (2022a), arXiv:2212.03888 [gr-qc] .
  * Barsanti  _et al._ (2022b) S. Barsanti, N. Franchini, L. Gualtieri, A. Maselli, and T. P. Sotiriou, Phys. Rev. D 106, 044029 (2022b), arXiv:2203.05003 [gr-qc] .
  * Vicente and Cardoso (2022) R. Vicente and V. Cardoso, (2022), arXiv:2201.08854 [gr-qc] .
  * Duque  _et al._ (2023) F. Duque, C. F. B. Macedo, R. Vicente, and V. Cardoso, (2023), arXiv:2312.06767 [gr-qc] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Yunes and Berti (2008) N. Yunes and E. Berti, Phys. Rev. D 77, 124006 (2008), [Erratum: Phys.Rev.D 83, 109901 (2011)], arXiv:0803.1853 [gr-qc] .
  * Zhang  _et al._ (2011) Z. Zhang, N. Yunes, and E. Berti, Phys. Rev. D 84, 024029 (2011), arXiv:1103.6041 [gr-qc] .
  * Sago  _et al._ (2016) N. Sago, R. Fujita, and H. Nakano, Phys. Rev. D 93, 104023 (2016), arXiv:1601.02174 [gr-qc] .
  * Cardoso  _et al._ (2022a) V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Phys. Rev. D 105, L061501 (2022a), arXiv:2109.00005 [gr-qc] .
  * Cardoso  _et al._ (2022b) V. Cardoso, K. Destounis, F. Duque, R. Panosso Macedo, and A. Maselli, Phys. Rev. Lett. 129, 241103 (2022b), arXiv:2210.01133 [gr-qc] .
  * Hernquist (1990) L. Hernquist, Astrophys. J. 356, 359 (1990).
  * Figueiredo  _et al._ (2023) E. Figueiredo, A. Maselli, and V. Cardoso, Phys. Rev. D 107, 104033 (2023), arXiv:2303.08183 [gr-qc] .
  * Einstein (1939) A. Einstein, Annals Math. 40, 922 (1939).
  * Geralico  _et al._ (2012) A. Geralico, F. Pompi, and R. Ruffini, Int. J. Mod. Phys. Conf. Ser. 12, 146 (2012).
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
  * Zerilli (1970) F. J. Zerilli, Phys. Rev. D 2, 2141 (1970).
  * Lindblom and Detweiler (1983) L. Lindblom and S. L. Detweiler, Astrophys. J. Suppl. 53, 73 (1983).
  * Sadeghian  _et al._ (2013) L. Sadeghian, F. Ferrer, and C. M. Will, Phys. Rev. D 88, 063522 (2013), arXiv:1305.2619 [astro-ph.GA] .
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 .
  * Einasto (1965) J. Einasto, Trudy Astrofizicheskogo Instituta Alma-Ata 5, 87 (1965).
  * Einasto (1969) J. Einasto, Astrofizika 5, 137 (1969).
  * Graham  _et al._ (2006) A. W. Graham, D. Merritt, B. Moore, J. Diemand, and B. Terzic, Astron. J. 132, 2685 (2006), arXiv:astro-ph/0509417 .
  * Prada  _et al._ (2006) F. Prada, A. A. Klypin, E. Simonneau, J. Betancort-Rijo, S. Patiri, S. Gottlober, and M. A. Sanchez-Conde, Astrophys. J. 645, 1001 (2006), arXiv:astro-ph/0506432 .
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999), arXiv:astro-ph/9906391 .
  * (39) “sgrep repo,” (github.com/masellia/SGREP/).
  * Martel (2004) K. Martel, Phys. Rev. D 69, 044025 (2004), arXiv:gr-qc/0311017 .
  * Cardoso  _et al._ (2019) V. Cardoso, M. Kimura, A. Maselli, E. Berti, C. F. B. Macedo, and R. McManus, Phys. Rev. D 99, 104077 (2019), arXiv:1901.01265 [gr-qc] .
  * Traykova  _et al._ (2023) D. Traykova, R. Vicente, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, Phys. Rev. D 108, L121502 (2023), arXiv:2305.10492 [gr-qc] .
  * Kavanagh  _et al._ (2020) B. J. Kavanagh, D. A. Nichols, G. Bertone, and D. Gaggero, Phys. Rev. D 102, 083006 (2020), arXiv:2002.12811 [gr-qc] .
  * Nichols  _et al._ (2023) D. A. Nichols, B. A. Wade, and A. M. Grant, Phys. Rev. D 108, 124062 (2023), arXiv:2309.06498 [gr-qc] .
  * Ghoshal and Strumia (2023) A. Ghoshal and A. Strumia, (2023), arXiv:2306.17158 [astro-ph.CO] .
  * Cole  _et al._ (2023) P. S. Cole, A. Coogan, B. J. Kavanagh, and G. Bertone, Phys. Rev. D 107, 083006 (2023), arXiv:2207.07576 [astro-ph.CO] .
  * Barausse  _et al._ (2015) E. Barausse, V. Cardoso, and P. Pani, J. Phys. Conf. Ser. 610, 012044 (2015), arXiv:1404.7140 [astro-ph.CO] .