KOBE-COSMO-21-06

# Indirect detection of gravitons through quantum entanglement

Sugumi Kanno  Department of Physics, Kyushu University, 744 Motooka, Nishi-ku,
Fukuoka 819-0395, Japan  Jiro Soda  Department of Physics, Kobe University,
Kobe 657-8501, Japan  Junsei Tokuda  Department of Physics, Kobe University,
Kobe 657-8501, Japan

###### Abstract

We propose an experiment that the entanglement between two macroscopic mirrors
suspended at the end of an equal-arm interferometer is destroyed by the noise
of gravitons through bremsstrahlung. By calculating the correlation function
of the noise, we obtain the decoherence time from the decoherence functional.
We estimate that the decoherence time induced by the noise of gravitons in
squeezed states stemming from inflation is approximately 20 seconds for 40 km
long arms and 40 kg mirrors. Our analysis shows that observation of the
decoherence time of quantum entanglement has the potential to detect gravitons
indirectly. This indirect detection of gravitons would give strong evidence of
quantum gravity.

Introduction.— It is widely believed that gravity should be quantized because
other forces are quantized. Nevertheless, there is no satisfactory quantum
theory of gravity. Some physicists think that gravity may not be quantized
after all Jacobson:1995ab . Thus, to find an experimental evidence of quantum
gravity is quite important. A clear consequence of quantum theory of gravity
would be the existence of gravitons. Hence, it is desired to come up with a
novel way to observe gravitons experimentally.

It has been known that gravitons are imperceptible. In fact, Dyson has
conjectured that no conceivable experiment in our universe can detect a single
graviton DYSON:2013jra ; Rothman:2006fp . If it is true, we need to seek
alternative ways for confirming the existence of gravitons. One possible way
is to focus on an inflationary scenario. It is believed that primordial
gravitational waves can be generated during inflation from quantum
fluctuations of geometry. If we succeed in observing the primordial
gravitational waves, it would imply a discovery of gravitons Kanno:2017dci ;
Kanno:2018cuk ; Kanno:2019gqw . However, even if the primordial gravitational
waves arrive at the interferometers, the direct detection of the gravitons is
difficult with current experimental techniques. Another possible way is to
focus on the statistical property of primordial gravitational waves. The state
of the primordial gravitational waves becomes squeezed during inflation
Grishchuk:1974ny ; Grishchuk:1990bj . Unfortunately, it is pointed out that
observing the squeezed states is also practically impossible Allen:1999xw .

Let us recall the history of discovery of atoms. Einstein used Brownian motion
to deduce the existence of atoms. In the same way, instead of direct detection
of gravitons, indirect search for gravitons might be possible. Recently, the
noise induced by gravitons is discussed in Parikh:2020nrd ; Parikh:2020kfh ;
Parikh:2020fhy ; Haba:2020jqs . An indirect detection of gravitons by making
use of the process of decoherence through the noise of gravitons is also
proposed in Riedel:2013yca ; Suzuki:2015nva ; Kanno:2020usf . However, no
feasible experiment has been proposed yet. The goal of this letter is to
propose a feasible experimental setup for indirect detection of gravitons in
the squeezed states stemming from inflation. In the following, we work in the
natural unit: $c=\hbar=1$.

Figure 1: The proposed setup: an equal-arm Michelson interferometer for a
single photon where there is a macroscopic suspended mirror at the end of each
arm

Experimental setups.— The proposed setup, shown in FIG. 1 consists of an
equal-arm Michelson interferometer which has a macroscopic suspended mirror at
the end of each arm. An incident photon (photon beams with very low intensity)
is injected from left and the beam splitter converts the photon into a
superposition state of being in both the upper and lower arms simultaneously
until it is detected by the oscillation of either mirror. Here we assume the
oscillation of the upper and lower mirrors ($2$ and $1$) as a semiclassical
state described by

$\displaystyle\vec{\xi}_{1}(t)=(\xi_{1},0,0)\,,\quad\xi_{1}=A\cos\omega t\,,$
$\displaystyle\vec{\xi}_{2}(t)=(0,\xi_{2},0)\,,\quad\xi_{2}=A\cos\omega t\,,$
(1)

where the amplitude of the oscillation $A$ and the angular frequency $\omega$
of both mirrors are set to the same value for simplicity. Let $H_{1}$ and
$H_{2}$ denote each Hilbert space of mirrors. Then the Hilbert space for the
combined system of the two mirrors is $H_{1}\otimes H_{2}$. Let $H_{1}$ be
spanned by a basis $\\{|0\rangle,\big{|}\vec{\xi}_{1}\rangle\\}$ and $H_{2}$
be spanned by $\\{|0\rangle,\big{|}\vec{\xi}_{2}\rangle\\}$ where $|0\rangle$
is the vacuum state in equilibrium and $|\vec{\xi}_{\rm i}\rangle$ (${\rm
i}=1,2$) represents a coherent state as an excited state induced by a single
photon. After an incident photon is injected, the state of both mirrors with
mass $m$ is described by the maximally entangled state where mirror 1
oscillates and mirror 2 is in equilibrium or mirror 1 in equilibrium and
mirror 2 oscillates

$\frac{1}{\sqrt{2}}\,\big{|}\vec{\xi}_{1}\rangle\otimes\big{|}0\rangle+\frac{1}{\sqrt{2}}\,\big{|}0\rangle\otimes\big{|}\vec{\xi}_{2}\rangle\,.$
(2)

Now, let us consider the influence of environmental quantum gravitational
fields (gravitons) on the system of the mirrors. We assume that the initial
total system is given by
$\big{|}0\rangle\otimes\big{|}0\rangle\otimes\big{|}{\rm vac}\rangle\,,$ where
$|{\rm vac}\rangle$ represents the vacuum state of gravitons. After the
incident photon is injected, the system of the mirrors gets entangled and the
total system becomes

$\big{|}\psi(t_{\rm
i})\rangle=\left\\{\frac{1}{\sqrt{2}}\big{|}\vec{\xi}_{1}\rangle\otimes\big{|}0\rangle+\frac{1}{\sqrt{2}}\big{|}0\rangle\otimes\big{|}\vec{\xi}_{2}\rangle\right\\}\otimes\big{|}{\rm
vac}\rangle\,.$ (3)

where $t_{\rm i}$ is the initial time. If we focus on the system of the
mirrors, the reduced density operator is obtained by tracing out the degree of
freedom of environmental gravitons such as

$\displaystyle\rho_{m}(t_{\rm i})$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{\rm grav}|\psi(t_{\rm i})\rangle\langle\psi(t_{\rm i})|$ (4)
$\displaystyle=$ $\displaystyle\rho_{11}(t_{\rm i})+\rho_{22}(t_{\rm
i})+\rho_{12}(t_{\rm i})+\rho_{21}(t_{\rm i})\,,\qquad$

where we defined
$\rho_{11}\equiv|\vec{\xi}_{1}\rangle|0\rangle\langle\vec{\xi}_{1}|\langle 0|$
etc. The interference term $\rho_{12}+\rho_{21}$ shows the initial entangled
state between the mirrors.
The total system evolves in time according to the Langevin equation of
geodesic deviation of the mirrors in the presence of gravitons Parikh:2020nrd
; Kanno:2020usf . Since each mirror interacts with the surrounding gravitons
as we will see in Eq. (13), the vacuum state of the gravitons changes due to
bremsstralung in accordance with the oscillation of the mirror on either side.
The total system is then found to be

$\displaystyle\big{|}\psi(t_{\rm f})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\big{|}\vec{\xi}_{1}\rangle\otimes\big{|}0\rangle\otimes\big{|}{\rm
vac};{\vec{\xi}_{1}}\rangle$ (5)
$\displaystyle+\frac{1}{\sqrt{2}}\big{|}0\rangle\otimes\big{|}\vec{\xi}_{2}\rangle\otimes\big{|}{\rm
vac};{\vec{\xi}_{2}}\rangle\,,$

where $t_{\rm f}$ is the final time and $|\rm{vac};\vec{\xi}_{i}\rangle$
represents the state of gravitons after the interaction with the oscillating
mirror ${\rm i}$. We ignored the force of radiation reaction of the mirror
because the gravitational back reaction is negligible and then the state of
each mirror remains unchanged. The point here is that the system of the
mirrors gets entangled with environmental gravitons due to the matter-gravity
interaction. In other words, the decoherence of the system occurs due to
bremsstralung of gravitons. The effect of the bremsstralung of gravitons is
expressed by the reduced density operator of the form

$\displaystyle\rho_{m}(t_{f})$ $\displaystyle=$ $\displaystyle\rho_{11}(t_{\rm
i})+\rho_{22}(t_{\rm i})$ (6) $\displaystyle+\exp(i\Phi)\rho_{12}(t_{\rm
i})+\exp(-i\Phi^{*})\rho_{21}(t_{\rm i})\,,\qquad$

where the influence of the gravitons on the reduced system is expressed by the
influence functional $\Phi$. The imaginary part of the influential functional
suppress the interference term and that is referred to as the decoherence
functional $\Gamma={\rm Im}\Phi$ Feynman:1963fq ; Breuer . If $\Gamma$ is
calculated, we can read off the decoherence time from it. In order to quantify
the change of entanglement over time, the entanglement negativity as a measure
of entanglement is computed as

$\displaystyle{\cal N}=\big{|}\rho_{12}(t_{\rm i})\big{|}\exp(-\Gamma)\ .$ (7)

By observing the decoherence time or the change of entanglement over time, we
could detect gravitons indirectly. In the next section, we see how the mirrors
interact with gravitons.

Action for the total system.— Let us consider the action for gravitational
waves in the Minkowski space first. The metric describing gravitational waves
in the transverse traceless gauge is expressed as

$\displaystyle ds^{2}=-dt^{2}+(\delta_{ij}+h_{ij})dx^{i}dx^{j}\,\,,$ (8)

where $t$ is the time, $x_{i}$ are spatial coordinates, $\delta_{ij}$ and
$h_{ij}$ are the Kronecker delta and the metric perturbations which satisfy
the transverse traceless conditions $h_{ij,\,j}=h_{ii}=0$. The indices $(i,j)$
run from $1$ to $3$. Substituting the metric Eq. (8) into the Einstein-Hilbert
action, we obtain the quadratic action for the metric perturbations

$\displaystyle S_{g}=\frac{M_{\rm p}^{2}}{8}\int{\rm
d}^{4}x\,\left[\,\dot{h}^{ij}\dot{h}_{ij}-h^{ij,k}\,h_{ij,k}\,\right]\,,$ (9)

where the reduced Planck mass is defined by $M_{\rm p}^{-2}=8\pi G$ and a dot
denotes the derivative with respect to the time. We can expand the metric
field $h_{ij}(x^{i},t)$ in terms of the Fourier modes

$\displaystyle h_{ij}(x^{i},t)=\frac{2}{M_{\rm p}\sqrt{V}}\sum_{{\bf k},B}\
h^{B}_{\bf{k}}(t)\,e^{i{\bf k}\cdot{\bf x}}\ e_{ij}^{B}(\bf k)\,,\qquad$ (10)

where we introduced the polarization tensor $e^{B}_{ij}({\bf k})$ normalized
as $e^{*B}_{ij}({\bf k})e^{C}_{ij}({\bf k})=\delta^{BC}$. Here, the index $B$
denotes the linear polarization modes $B=+,\times$. Note that we consider
finite volume $V=L_{x}L_{y}L_{z}$ and discretize the ${\bf k}$-mode with a
width ${\bf k}=\left(2\pi n_{x}/L_{x}\,,2\pi n_{y}/L_{y}\,,2\pi
n_{z}/L_{z}\right)$ where ${\bf n}=(n_{x},n_{y},n_{z})$ are integers. Note
that we used $k=|\bf k|$. We then see that a gravitational wave consists of an
infinite number of harmonic oscillators.
Next, we consider the action for two mirrors. Here, we regard the mirror as a
point particle effectively because the dynamical degree of the center of mass
is essential for discussing the noise of gravitons later. Then we consider how
the particle feels gravitational waves. In fact, a single particle does not
feel the gravitational waves because of the Einstein’s equivalence principle
at least classically. Thus, we need to consider geodesic deviation between the
two particles in order to see the effect of gravity on them. We introduce the
Fermi normal coordinates for calculating the geodesic deviation between the
time-like geodesics of two particles $\gamma_{\tau}$ and
$\gamma_{\tau^{\prime}}$. The Fermi normal coordinates are local inertial
coordinates that are adapted to a geodesic. We expand the coordinate along the
time-like and space-like geodesics that are orthogonal to each other at the
position of the beam splitter as shown in FIG. 1 where the space-like
geodesics are simplified as $(x,y)$-plane. The dynamics of the $\gamma_{\tau}$
and $\gamma_{\tau^{\prime}}$ is described by the position $\xi_{\rm i}(t),{\rm
i}=1,2$ respectively in the vicinity of the origin of the beam splitter. The
action for a particle along $\gamma_{\tau}$ is given by

$\displaystyle
S_{p}=-m\int_{\gamma_{\tau}}d\tau=-m\int_{\gamma_{\tau}}dt\sqrt{-g_{\mu\nu}(t,\xi^{i}(t))\dot{\xi}^{\mu}\dot{\xi}^{\nu}}\,,$
(11)

where $\tau$ is proper time and the position of the particle is represented by
$\xi^{\mu}=(t,\xi^{i}(t))$. The metric $g_{\mu\nu}$ up to the second order of
arbitrary position $x^{i}$ in the Fermi normal coordinates is computed as

$\displaystyle ds^{2}$ $\displaystyle\simeq$
$\displaystyle\left(-1-R_{0i0j}\,x^{i}x^{j}\right)dt^{2}-\frac{4}{3}R_{0jik}\,x^{j}x^{k}\,dtdx^{i}$
(12)
$\displaystyle+\left(\delta_{ij}-\frac{1}{3}R_{ikj\ell}\,x^{k}x^{\ell}\right)dx^{i}dx^{j}\
.$

Here the Riemann tensor is evaluated at the origin $x^{i}=0$ in the Fermi
normal coordinate system. Because the Riemann tensor $R_{0i0j}$ is gauge
invariant at the leading order in the metric fluctuation $h_{ij}$, we can
evaluate it in the transverse traceless gauge and obtain
$R_{0i0j}(0,t)=-\ddot{h}_{ij}(0,t)/2$. Substituting the metric (12) into the
action (11), we can read off the interaction between the particles and
gravitational waves as

$\displaystyle S_{\rm int}$ $\displaystyle\simeq$ $\displaystyle\int
dt\left[\frac{m}{4}\ddot{h}_{11}(0,t)\,{\xi^{1}_{1}}^{2}+\frac{m}{4}\ddot{h}_{22}(0,t)\,{\xi^{2}_{2}}^{2}\right]\,.$
(13)

where $x,y$ components are expressed by $1,2$ respectively. We find that a
qubic derivative interaction between gravitational waves and geodesic
deviation appeared in the action. We can then derive the Langevin equation of
geodesic deviation of the mirrors in the presence of gravitons from this
action Parikh:2020nrd ; Kanno:2020usf . Note that the geodesic deviation of
particles in the graviton background is studied in DeLorenci:2014vwa ;
Quinones:2017wka . Note also that a particle motion in the graviton background
is discussed in Oniga:2015lro ; Oniga:2017pyq  where a different form of
interaction is used.

Quantization and the noise induced by gravitons.— Now we canonically quantize
the total system. For simplicity, we work in the interaction picture below. We
promote the free metric field $h^{A}({\bf k},t)$ to the operator
$\hat{h}^{A}({\bf k},t)$ in terms of the creation and annihilation operators
such as

$\displaystyle\hat{h}^{A}({\bf k},t)=\hat{a}_{A}({\bf
k})u_{k}(t)+\hat{a}^{\dagger}_{A}(-{\bf k})u^{*}_{k}(t)\ \,,$ (14)

where the creation and annihilation operators satisfy the standard commutation
relations $[\,{\hat{a}}_{A}({\bf k})\,,\,\hat{a}^{\dagger}_{A^{\prime}}({\bf
k^{\prime}})\,]=\delta_{{\bf k},{\bf{k}^{\prime}}}\delta_{AA^{\prime}}\,,$ and
$u_{k}(t)$ denotes a mode function properly normalized as
$\dot{u}_{k}(t)u_{k}^{*}(t)-u_{k}(t)\dot{u}_{k}^{*}(t)=-i\,.$ The Minkowski
vacuum $\ket{0}$ is defined by $\hat{a}_{A}({\bf k})\ket{0}=0$, with choosing
the mode function as $u_{k}(t)=\frac{1}{\sqrt{2k}}e^{-ikt}\equiv u_{k}^{\rm
M}(t)\ .$ From the Langevin equation of geodesic deviation of the mirrors in
the presence of gravitons, the noise of gravitons is identified as
Kanno:2020usf

$\displaystyle\hat{N}_{ij}(t)\equiv\frac{1}{M_{\rm
p}\sqrt{V}}\sum_{A}\sum_{{\bf k}\leq\Omega_{\rm m}}k^{2}e^{A}_{ij}({\bf
k})\hat{h}^{A}({\bf k},t)\,,\qquad$ (15)

where $\sum_{{\bf k}\leq\Omega_{\rm m}}$ represents the mode sum with the UV
cutoff. This noise of gravitons always exists if the gravitational waves are
quantized.

The decoherence functional and the decoherence time.— Since we found that
gravitons are quantum fluctuations of gravitational waves in the form of the
noise Eq. (15), we calculate the decoherence functional and then find the
necessary experimental setups for obtaining a measurable time of decoherence
due to the noise of gravitons. For this purpose, we consider squeezed states
$|\zeta\rangle$ below. For instance, the state of the primordial gravitational
waves becomes squeezed during inflation Grishchuk:1974ny ; Grishchuk:1990bj .
The decoherence functional reads Breuer ; Kanno:2020usf

$\displaystyle\Gamma(t_{\rm f})$ $\displaystyle\approx$
$\displaystyle\frac{m^{2}}{8}\int^{t_{\rm
f}}_{0}\mathrm{d}t\,\Delta(\xi^{i}\xi^{j})(t)\int^{t_{\rm
f}}_{0}\mathrm{d}t^{\prime}\,\Delta(\xi^{k}\xi^{\ell})(t^{\prime})$ (16)
$\displaystyle\times\Big{\langle}\left\\{\hat{N}_{ij}(t),\,\hat{N}_{k\ell}(t^{\prime})\right\\}\Big{\rangle}\,.$

Here
$\Delta(\xi^{i}\xi^{j})(t)=\xi^{i}_{1}(t)\xi^{j}_{1}(t)-\xi^{i}_{2}(t)\xi^{j}_{2}(t)$
denotes a difference of $\xi^{i}(t)\xi^{j}(t)$ in the superposition. The
anticommutator correlation function of $N_{ij}(t)$ in the squeezed state can
be computed in the infinite volume limit $L_{x},L_{y},L_{z}\to\infty$ as

$\displaystyle\Big{\langle}\zeta\Big{|}\left\\{\hat{N}_{ij}(t),\,\hat{N}_{k\ell}(t^{\prime})\right\\}\Big{|}\zeta\Big{\rangle}$
$\displaystyle=\left(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk}-\frac{2}{3}\delta_{ij}\delta_{k\ell}\right)\frac{F(t-t^{\prime})}{10\pi^{2}M_{\rm
p}^{2}}\,,\quad$ (17)

where we defined the anticommutator symbol $\\{\cdot\,,\cdot\\}$ as
$\\{\hat{X},\hat{Y}\\}\equiv(\hat{X}\hat{Y}+\hat{Y}\hat{X})/2$ and

$\displaystyle
F(t-t^{\prime})\equiv\int^{\Omega_{m}}_{0}\mathrm{d}k\,k^{6}\mathrm{Re}\left[u^{\rm
sq}_{k}(t){u^{\rm sq}_{k}}^{*}(t^{\prime})\right]\ .$ (18)

Here, the mode function in the squeezed state is given by

$\displaystyle u^{\rm sq}_{k}(t)\equiv u^{\rm M}_{k}(t)\cosh
r_{k}-e^{-i\varphi_{k}}{u_{k}^{\rm M}}^{*}(t)\sinh r_{k}\,.$ (19)

In general, the squeezing parameter $r_{k}$ and the phase $\varphi_{k}$ depend
on $k$ but we assume $\varphi_{k}=\varphi$ for simplicity and $r_{k}\gg 1$
which corresponds to the end of inflation. In the conventional inflationary
scenario, we obtain $\sinh 2r_{k}\simeq\cosh 2r_{k}\simeq\left(k_{\rm
c}/k\right)^{4}$ where $k_{\rm c}=2\pi f_{\rm c}$ and $f_{\rm c}$ is the
cutoff frequency of primordial gravitational waves. The bounds on the cutoff
frequency from the CMB is $f_{\rm c}\lesssim 10^{9}$ Hz Maggiore:1999vm . For
alternatives to inflation, the bounds become $f_{\rm c}\lesssim 4.3\times
10^{10}$ Hz Maggiore:1999vm . We have a maximum value of the noise correlation
(17) for $\varphi=\pi$. By a change of variables $y=k/\Omega_{\rm
m},x=\Omega_{\rm m}(t-t^{\prime})$, we can evaluate the integral as

$\displaystyle F(x)=k_{\rm c}^{4}\Omega_{m}^{2}\frac{x\sin x+\cos x-1}{x^{2}}\
.$ (20)

In our setup, the decoherence functional (16) is

$\displaystyle\Gamma=\frac{m^{2}}{120\pi^{2}M_{\rm p}^{2}}\int_{0}^{t_{\rm
f}}dt\Delta\xi_{1}^{2}(t)\int_{0}^{t_{\rm
f}}dt^{\prime}\Delta\xi_{1}^{2}(t^{\prime})F(t-t^{\prime})$
$\displaystyle+\frac{m^{2}}{120\pi^{2}M_{\rm p}^{2}}\int_{0}^{t_{\rm
f}}dt\Delta\xi_{2}^{2}(t)\int_{0}^{t_{\rm
f}}dt^{\prime}\Delta\xi_{2}^{2}(t^{\prime})F(t-t^{\prime})$
$\displaystyle-\frac{m^{2}}{120\pi^{2}M_{\rm p}^{2}}\int_{0}^{t_{\rm
f}}dt\Delta\xi_{1}^{2}(t)\int_{0}^{t_{\rm
f}}dt^{\prime}\Delta\xi_{2}^{2}(t^{\prime})F(t-t^{\prime})\,,$ (21)

where $\Delta\xi_{1}^{2},\Delta\xi_{1}^{2}$ are

$\displaystyle\Delta\xi_{1}^{2}(t)$ $\displaystyle=$
$\displaystyle\left(L+A\cos\omega t\right)^{2}-L^{2}\,,$ (22)
$\displaystyle\Delta\xi_{2}^{2}(t)$ $\displaystyle=$ $\displaystyle
L^{2}-\left(L+A\cos\omega t\right)^{2}\,.$ (23)

Here, we stress that the decoherece process we considered is non-Markovian,
that is, we took into account the correlated noise of gravitons. In the
literature Blencowe:2012mp ; Anastopoulos:2013zya , however, the decoherence
rate is computed in the Markovian approximation.
Let us evaluate the double integrals in Eq. (21). Since the amplitude of
oscillation induced by a photon is negligible compared with the arm length
$A\ll L$, $\Delta\xi_{i}^{2}\sim 2LA\cos\omega t$ and we obtain

$\displaystyle\Gamma=\frac{m^{2}}{10\pi^{2}M_{\rm p}^{2}}\left(AL\right)^{2}I\
,$ (24)

where

$\displaystyle I$ $\displaystyle=$ $\displaystyle\ \frac{1}{\Omega_{\rm
m}^{2}}\int_{0}^{t_{\rm f}}dx\cos\beta x$ (25)
$\displaystyle\times\int_{0}^{t_{\rm f}}dx^{\prime}\cos\beta
x^{\prime}F(x-x^{\prime})\,.$

Here we defined the dimensionless time and frequency $x=\Omega_{\rm m}t$,
$\beta=\omega/\Omega_{\rm m}$, respectively. After the integration, we find
the following terms become dominant for large $x$

$\displaystyle\frac{I(x)}{k_{c}^{4}}$ $\displaystyle\simeq$
$\displaystyle-\frac{1}{2}\beta x\,{\rm Si}\left((1+\beta)x\right)$ (26)
$\displaystyle+\frac{1}{2}\beta x\,{\rm Si}\left((1-\beta)x\right)+\beta
x\,\rm{Si}(\beta x)\,,$

where the sine integral is defined by ${\rm Si}(x)=\int^{x}_{0}\frac{\sin
t}{t}dt$ whose value at $x=\infty$ is $\pi/2$. Thus, the decoherence
functional is found to be

$\displaystyle\Gamma=\frac{4\pi^{3}}{5}\left(\frac{m}{M_{\rm
p}}\right)^{2}\left(Lf_{\rm c}\right)^{4}\left(\frac{A}{L}\right)^{2}N\ ,$
(27)

where we defined $N=\omega t_{f}$. Remarkably, the resultant decoherence
functional is independent of the UV cutoff $\Omega_{\rm m}$, which indicates
the result is reliable. If we consider conventional inflation and take the
parameters, $\omega=1$ kHz, $L=40$ km, $m=40$ kg, $f_{\rm c}=10^{9}$ Hz,
$N=2\times 10^{4}$, we obtain $\Gamma\simeq 1$ where the amplitude of
oscillation is supposed to be about 10 times zero-point fluctuations
$A=10/\sqrt{2m\omega}$ Gely:2021fhv . Then the decoherence time ${t_{\rm f}}$
is approximately $20$ s. In this case, the time evolution of negativity is
plotted in FIG. 2. If we increase the arm length $L$ by $10^{2}$, the
decoherence time becomes $2$ ms.
For alternatives to inflation ($f_{\rm c}=4.3\times 10^{10}$ Hz
Maggiore:1999vm ), the decoherence time becomes $6$ $\mu$s with the same
parameters.

Figure 2: Time evolution of negativity normalized by the initial value ${\cal
N}_{\rm i}$ for $\omega=1$ kHz, $L=40$ km, $m=40$ kg, $f_{\rm c}=10^{9}$ Hz.
The negativity decays with the decoherence time $20$ s.

Conclusion.— We performed a detailed study of the experimental setup for
detecting gravitons indirectly by observing the decoherence time of the
entanglement between two macroscopic mirrors suspended at the end of an equal-
arm interferometer. In principle the proposed setup can be realized in near
future. We found that longer arms or heavier mirrors make the decoherence time
shorter. We estimated that the decoherence time induced by the noise of
gravitons in the squeezed state stemming from inflation is approximately 20
seconds for 40 km long arms and 40 kg mirrors. The arm length is 1 order of
magnitude longer than that of LIGO. To generate entanglement between 40
kilogram mirrors appears to be within reach of current technology Yu:2020ece .
The other values used in the estimation is achieved in laboratory to date.
There remains an issue to suppress other sources of the decoherence. While
very demanding, the indirect detection of gravitons appears to be within reach
in near future.

Acknowledgments.— We would like to thank Nobuyuki Matsumoto for useful
discussion and suggestions. S. K. was supported by the Japan Society for the
Promotion of Science (JSPS) KAKENHI Grant Number JP18H05862. J. S. was in part
supported by JSPS KAKENHI Grant Numbers JP17H02894, JP17K18778, JP20H01902. J.
T. was supported by JSPS Postdoctoral Fellowship No. 202000912.

## References

  * (1) T. Jacobson, Phys. Rev. Lett. 75, 1260-1263 (1995) [arXiv:gr-qc/9504004 [gr-qc]].
  * (2) F. Dyson, Int. J. Mod. Phys. A 28, 1330041 (2013)
  * (3) T. Rothman and S. Boughn, Found. Phys. 36, 1801-1825 (2006) [arXiv:gr-qc/0601043 [gr-qc]].
  * (4) S. Kanno and J. Soda, Phys. Rev. D 96, no.8, 083501 (2017) [arXiv:1705.06199 [hep-th]].
  * (5) S. Kanno and J. Soda, Phys. Rev. D 99, no.8, 084010 (2019) [arXiv:1810.07604 [hep-th]].
  * (6) S. Kanno, Phys. Rev. D 100, no.12, 123536 (2019) [arXiv:1905.06800 [hep-th]].
  * (7) L. P. Grishchuk, Zh. Eksp. Teor. Fiz. 67, 825-838 (1974)
  * (8) L. P. Grishchuk and Y. V. Sidorov, Phys. Rev. D 42, 3413-3421 (1990) doi:10.1103/PhysRevD.42.3413
  * (9) B. Allen, E. E. Flanagan and M. A. Papa, Phys. Rev. D 61, 024024 (2000) [arXiv:gr-qc/9906054 [gr-qc]].
  * (10) M. Parikh, F. Wilczek and G. Zahariade, [arXiv:2005.07211 [hep-th]].
  * (11) M. Parikh, F. Wilczek and G. Zahariade, [arXiv:2010.08205 [hep-th]].
  * (12) M. Parikh, F. Wilczek and G. Zahariade, [arXiv:2010.08208 [hep-th]].
  * (13) Z. Haba, Eur. Phys. J. C 81, no.1, 40 (2021) [arXiv:2009.12306 [gr-qc]].
  * (14) C. J. Riedel, [arXiv:1310.6347 [quant-ph]].
  * (15) F. Suzuki and F. Queisser, J. Phys. Conf. Ser. 626, no.1, 012039 (2015) [arXiv:1502.01386 [gr-qc]].
  * (16) S. Kanno, J. Soda and J. Tokuda, Phys. Rev. D 103, no.4, 044017 (2021) [arXiv:2007.09838 [hep-th]].
  * (17) R. P. Feynman and F. L. Vernon, Jr., Annals Phys. 24, 118-173 (1963)
  * (18) Heinz-Peter Breuer and Francesco Petruccione Phys. Rev. A 63, 032102
  * (19) V. A. De Lorenci and L. H. Ford, Phys. Rev. D 91, no.4, 044038 (2015) [arXiv:1412.4685 [gr-qc]].
  * (20) D. A. Quinones, T. Oniga, B. T. H. Varcoe and C. H. T. Wang, Phys. Rev. D 96, no.4, 044018 (2017) [arXiv:1702.03905 [gr-qc]].
  * (21) T. Oniga and C. H. T. Wang, Phys. Rev. D 93, no.4, 044027 (2016) [arXiv:1511.06678 [quant-ph]].
  * (22) T. Oniga and C. H. T. Wang, Phys. Rev. D 96, no.8, 084014 (2017) [arXiv:1701.04122 [gr-qc]].
  * (23) M. Maggiore, Phys. Rept. 331, 283-367 (2000) [arXiv:gr-qc/9909001 [gr-qc]].
  * (24) M. P. Blencowe, Phys. Rev. Lett. 111, no.2, 021302 (2013) [arXiv:1211.4751 [quant-ph]].
  * (25) C. Anastopoulos and B. L. Hu, Class. Quant. Grav. 30, 165007 (2013) [arXiv:1305.5231 [gr-qc]].
  * (26) M. F. Gely and G. A. Steele, [arXiv:2103.12729 [quant-ph]].
  * (27) H. Yu et al. [LIGO Scientific], Nature 583, no.7814, 43-47 (2020) [arXiv:2002.01519 [quant-ph]].


XXXX-XXXX

# Axion cloud evaporation during inspiral of black hole binaries - the effects
of backreaction and radiation -

Takuya Takahashi<EMAIL_ADDRESS>Department of Physics,
Kyoto University, Kyoto 606-8502, Japan  Hidetoshi Omiya  Takahiro Tanaka
Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University, Kyoto 606-8502, Japan

###### Abstract

Ultralight scalar fields such as axions can form clouds around rotating black
holes (BHs) by the superradiant instability. It is important to consider the
evolution of clouds associated with BH binaries for the detectability of the
presence of clouds through gravitational wave signals and observations of the
mass and spin parameters of BHs. The impact on the axion cloud due to the
tidal perturbation from the companion in a binary system is first studied in
Ref Baumann:2019ztm . Here, we re-examine this issue taking into account the
following points. First, we study the influence of higher-multipole moments.
Second, we consider the backreaction due to the angular momentum transfer
between the cloud and the orbital motion. This angular momentum transfer
further causes the backreaction to the hyperfine split through the change in
geometry. Finally, we calculate the particle number flux to the infinity
induced by the tidal interaction. As a result, we found that the scalar field
is not reabsorbed by the BH. Instead, the scalar particles are radiated away
to evaporate during the inspiral, irrespectively of the direction of the
orbital motion, for almost equal mass binaries.

##  1 Introduction

Ultralight bosons such as axions are one of the best-motivated particles in
the context of an extension of the Standard Model. The existence of plenty of
axion-like particles is predicted by string theory Arvanitaki:2009fg ;
Svrcek:2006yi , and such particles are also a good candidate for the dark
matter Dine:1982ah ; Preskill:1982cy ; Abbott:1982af ; Hui:2016ltb .
Interestingly, they can make an observable signature on cosmology and
astrophysical phenomena. As one of them, we consider the observability of the
existence of an ultralight scalar field associated with black holes (BHs).

A massive scalar field around a BH can be bound by gravitational interaction,
and it can extract the rotation energy from a spinning BH through the
superradiance Brito:2020oca . When the Compton wavelength of the scalar field
is comparable to the radius of the BH, this energy extraction mechanism works
most efficiently, and thus the field forms a macroscopic condensate around the
BH. For astrophysical BHs, a scalar field with the mass in the range
$10^{-20}\sim 10^{-10}\mathrm{eV}$ has an appropriate Compton wavelength to
form this condensate Detweiler:1980uk ; Dolan:2007mj . In the following, we
refer to such a scalar field simply as an axion, and its condensate formed
around a BH as an axion cloud.

One observational signature of the existence of an axion cloud is the
appearance of a forbidden region in the parameter distribution of the mass and
spin of BHs, owing to the spin-down due to the superradiance Arvanitaki:2010sy
; Brito:2014wla ; Brito:2017zvb . In addition, gravitational waves (GWs)
emitted from an axion cloud and the modification of the gravitational waveform
from a binary coalescence accompanied by an axion cloud are considered to be
possible targets for observations Arvanitaki:2014wva ; Arvanitaki:2016qwi ;
Yoshino:2013ofa ; Choudhary:2020pxy ; Chung:2021roh ; DeLuca:2021ite ;
Baumann:2018vus ; Liu:2021llm . From observations, we are now getting new
constraints on the axion mass Tsukada:2018mbp ; Ng:2020ruv ;
LIGOScientific:2021jlr .

The GW events observed so far are those originating from binary coalescences
Abbott:2016blz ; TheLIGOScientific:2016qqj ; LIGOScientific:2021djp , and most
of BHs whose spin is measured are belonging to a binary system. On the other
hand, the tidal interaction from the companion in a binary greatly affects the
evolution of the axion cloud during the inspiral phase. Axion clouds are
initially dominated by the superradiant growing mode, but when the binary
orbital frequency coincides with the difference in the phase velocities with
another mode, a resonant level transition is induced Baumann:2018vus . The
disruption of the cloud due to tidal interaction from the static perturber is
also studied by numerical simulations Cardoso:2020hca .

If axions are transferred to a mode with a large decay rate, axions would be
reabsorbed by the BH and the cloud is depleted before the merger Berti:2019wnn
; Takahashi:2021eso . Furthermore, axions in the decaying mode typically have
an angular momentum with the opposite sign to that of the central BH, and the
reabsorption may accelerate the spin-down of the BH. As a result, the
forbidden region in the mass and spin distribution of BHs might be modified.
Therefore, to clarify the detectability of axions and to make the obtained
constraints more robust, we need to understand the dynamics of axion clouds in
binary systems more precisely.

The level transition of an axion cloud due to the tidal interaction is
formulated in Ref. Baumann:2019ztm . We develop this work further in the
following respects. In the previous works, the main focus was only on the
effect of the leading quadrupole moment of the tidal potential. However, since
the resonance caused by the effect of higher-multipole moments occurs earlier
during the binary inspiral, it may dominate the evolution of clouds. We also
consider the backreaction due to the angular momentum transfer between the
cloud and the orbital motion. We describe the transition between different
axion energy eigenstates taking into account the backreaction to the orbital
evolution using a non-linear model Hamiltonian, and derive the persistent rate
of axions in a strongly non-linear regime. Furthermore, we propose a model to
take into account the modification to the hyperfine split due to the
backreaction of the angular momentum transfer, and show this modification can
drastically change the final fate of the transition of axion clouds. In
addition, we examine a new channel that turns out to deplete the axion clouds;
axion emission induced by the tidal perturbation. Once axions are transferred
to an upper level by the higher-multipoles of the tidal field, they can be
easily further excited to an unbound state. Then, axions evaporate away to
infinity.

In this paper, we analyze axion clouds in the non-relativistic regime,
assuming that the binary’s orbit is quasi-circular and the orbital angular
mementum is aligned to the BH spin, for simplicity. Also, we neglect the axion
self-interaction, which may play an important role in the evolution of axion
clouds Baryakhtar:2020gao ; Omiya:2020vji . Hence, we consider the simplest
model of a massive scalar field and assume that there is nothing that prevents
the growth of clouds until the superradiance condition is saturated by the
spin-down of the central BH. We believe that this work will give a foundation
of the analysis of more complicated cases.

This paper is organized as follows. In Sec. 2, we review how the level
transition of axion clouds is described and show the importance of higher
multipole moments of the tidal potential. In Sec. 3, we analyze the effect of
the backreaction of the angular momentum transfer in detail. In Sec. 4, we
calculate the particle number flux to the infinity. Finally, In Sec. 5, we
summarize our results and give a brief discussion about their implications. In
the rest of this paper, we use the unit with $G=c=\hbar=1$.

##  2 Level transition of axion clouds

In this section, we first review the level transition of axion clouds induced
by the tidal perturbation due to the binary companion, formulated in Ref.
Baumann:2019ztm . Then, we show that higher-multipole moments of the tidal
perturbation actually play an important role in the evolution of clouds during
the inspiral phase.

###  2.1 Hamiltonian

The equation of motion of a real scalar field (axion) of mass $\mu$ around a
rotating BH is

$\left(g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}-\mu^{2}\right)\phi=0~{},$ (1)

where $g_{\mu\nu}$ is the Kerr metric. Let the BH mass be $M$ and the angular
momentum be $J=aM$. In the non-relativistic limit, it is appropriate to start
with the ansatz

$\phi=\frac{1}{\sqrt{2\mu}}\left(e^{-i\mu t}\psi+e^{i\mu
t}\psi^{\ast}\right)~{}.$ (2)

and $\psi$ is a complex scalar function whose timescale of variation is much
larger than $\mu^{-1}$. With this ansatz, the equation of motion (1) is
rewritten as

$i\frac{\partial}{\partial
t}\psi=\left(-\frac{1}{2\mu}\nabla^{2}-\frac{\alpha}{r}\right)\psi~{},$ (3)

at the leading order in $r^{-1}$. Here, we have introduced the gravitational
fine structure constant

$\alpha=M\mu~{}\simeq
0.2\left(\frac{M}{30M_{\odot}}\right)\left(\frac{\mu}{10^{-12}\mathrm{eV}}\right)~{}.$
(4)

Eq.(3) is formally equivalent to the Schr$\mathrm{\ddot{o}}$dinger equation
and the bound state solution is given by

$\psi_{nlm}=e^{-i(\omega_{nlm}-\mu)t}R_{nl}(r)Y_{lm}(\theta,\varphi)~{},$ (5)

where $R_{nl}(r)$ represents the radial wave function of the hydrogen atom
(Appendix. A) and $Y_{lm}(\theta,\varphi)$ is the spherical harmonics. The
radius corresponding to the Bohr radius for a hydrogen atom is
$r_{0}=M/\alpha^{2}$, and this approximation is well justified for $\alpha\ll
1$. In the non-relativistic approximation, the eigenfrequency is given by
Baumann:2019eav ; Baumann:2018vus ; Detweiler:1980uk

$\omega_{nlm}=(\omega_{R})_{nlm}+i(\omega_{I})_{nlm}~{},$ (6)

with

$\displaystyle(\omega_{R})_{nlm}$
$\displaystyle=\mu\left(1-\frac{\alpha^{2}}{2n^{2}}-\frac{\alpha^{4}}{8n^{4}}+\frac{(2l-3n+1)\alpha^{4}}{n^{4}(l+1/2)}+\frac{2m\alpha^{5}}{n^{3}l(l+1/2)(l+1)}\frac{a}{M}\right)~{},$
(7) $\displaystyle(\omega_{I})_{nlm}$
$\displaystyle=2(r_{+}/M)C_{nlm}(a,\alpha)(m\Omega_{H}-\omega)\alpha^{4l+5}~{},$
(8)

where $r_{+}$ is the horizon radius, $\Omega_{H}=a/2Mr_{+}$ is the angular
velocity of the BH horizon and explicit form of $C_{nlm}(a,\alpha)$ is given
in Ref Detweiler:1980uk . Note that the labels of the eigenstates that we can
take are restricted to $n\geq 1,l\leq n-1,|m|\leq l$. As you can see from
Eq.(8), the mode satisfying $\omega_{R}<m\Omega_{H}$ is a growing mode. We can
consider the situation in which axions occupy the fastest growing mode
initially. In the following, we assume that there is no mechanism to terminate
the growth of the cloud. Then, the growth of the cloud will saturate, when the
the BH spin is reduced to satisfy the condition for the critical spin

$\frac{a}{M}=\frac{a_{\mathrm{crit}}}{M}\equiv\frac{4m\alpha}{m^{2}+4\alpha^{2}}~{}.$
(9)

When a BH with an axion cloud forms a binary system, the tidal perturbation
from the binary companion of mass $M_{\ast}$ at
$(R_{\ast},\Theta_{\ast},\Phi_{\ast})$ can be described by the additional
potential Baumann:2018vus

$V_{\ast}(t)=-M_{\ast}\mu\sum_{l_{\ast}=2}^{\infty}\sum_{|m_{\ast}|\leq
l_{\ast}}\frac{4\pi}{2l_{\ast}+1}Y_{l_{\ast}m_{\ast}}^{\ast}(\Theta_{\ast},\Phi_{\ast})Y_{l_{\ast}m_{\ast}}(\theta,\varphi)\left(\frac{r^{l_{\ast}}}{R_{\ast}^{l_{\ast}+1}}\theta(R_{\ast}-r)+\frac{R_{\ast}^{l_{\ast}}}{r^{l_{\ast}+1}}\theta(r-R_{\ast})\right)$
(10)

in Eq.(3). The orbital angular velocity of the companion can be defined by
$\dot{\Phi}_{\ast}(t)=\pm\Omega(t)$, and the upper (lower) sign represents the
case of co-rotating (counter-rotating) orbits relatively to the central BH
spin. For simplicity, we focus on the case in which the binary orbit is on the
equatorial plane, $\Theta=\pi/2$.

We denote the $i$-th eigenstate of an axion as $\ket{i}$, and a general bound
state can be written by a superposition of eigenstates, i.e.,
$\ket{\psi(t)}=\sum_{i}\tilde{c}_{i}(t)\ket{i}$. Since the tidal interaction
acts as an periodically oscillating tiny external field and becomes important
only when the orbital angular velocity coincides with the difference between
the phase velocities of the two levels before and after transition. Therefore,
it is sufficient to consider a two level system, neglecting all the other
levels. $|\tilde{c}_{i}(t)|^{2}$ represents the normalized particle number
occupying the level $i$, and its time evolution is given by

$\displaystyle
i\frac{\mathrm{d}\tilde{c}_{i}}{\mathrm{d}t}=\sum_{j=1,2}\tilde{\mathcal{H}}_{ij}\tilde{c}_{j}~{},$
(11) $\displaystyle\tilde{\mathcal{H}}=\left(\begin{array}[]{cc}-\frac{\Delta
E}{2}&\eta e^{i\Delta m\Phi_{\ast}(t)}\\\ \eta e^{-i\Delta
m\Phi_{\ast}(t)}&\frac{\Delta E}{2}\end{array}\right)~{},$ (14)

where $\Delta E=(\omega_{R})_{2}-(\omega_{R})_{1}$ and $\Delta m=m_{2}-m_{1}$.
Writing

$V_{\ast}(t)=\sum_{l_{\ast}m_{\ast}}V_{\ast,\l_{\ast}m_{\ast}}e^{-im\Phi_{\ast}(t)}~{},$
(15)

we parameterize the amplitude of perturbation by
$\eta\equiv|\braket{1}{V_{\ast,l_{\ast}m_{\ast}}}{2}|$, whose explicit
expression is

$\displaystyle\frac{\eta}{\Omega}=$
$\displaystyle\left|\frac{q}{(1+q)^{(1+l_{\ast})/3}}\frac{(M\Omega)^{(2l_{\ast}-1)/3}}{\alpha^{2l_{\ast}-1}}\frac{4\pi}{2l_{\ast}+1}Y_{l_{\ast}m_{\ast}}\left(\Theta_{\ast},\Phi_{\ast}\right)\
I_{\tilde{r}}^{\mathrm{in}}\times I_{A}\right.$
$\displaystyle\left.+q(1+q)^{l_{\ast}/3}\frac{\alpha^{2l_{\ast}+3}}{(M\Omega)^{(2l_{\ast}+3)/3}}\frac{4\pi}{2l_{\ast}+1}Y_{l_{\ast}m_{\ast}}\left(\Theta_{\ast},\Phi_{\ast}\right)\
I_{\tilde{r}}^{\mathrm{out}}\times I_{A}\right|~{},$ (16)

with

$\displaystyle
I_{\tilde{r}}^{\mathrm{in}}=\int_{0}^{R_{\ast}/r_{0}}d\tilde{r}\
\tilde{r}^{2+l_{\ast}}\tilde{R}_{n_{1}l_{1}}\tilde{R}_{n_{2}l_{2}}~{},\quad
I_{\tilde{r}}^{\mathrm{out}}=\int_{R_{\ast}/r_{0}}^{\infty}d\tilde{r}\
\tilde{r}^{1-l_{\ast}}\tilde{R}_{n_{1}l_{1}}\tilde{R}_{n_{2}l_{2}}~{},$ (17)
$\displaystyle I_{A}=\int\sin\theta d\theta d\varphi\
Y_{l_{1}m_{1}}^{\ast}Y_{l_{\ast}m_{\ast}}Y_{l_{2}m_{2}}~{},$ (18)

where $\tilde{r}=r/r_{0}$, $\tilde{R}_{nl}=r_{0}^{3/2}R_{nl}$. Here we defined
the mass ratio by

$q=\frac{M_{\ast}}{M}~{}.$ (19)

Performing the unitary transformation with the transformation matrix
$\mathcal{U}=\mathrm{diag}(e^{i\Delta m\Phi_{\ast}(t)/2}$, $e^{-i\Delta
m\Phi_{\ast}(t)/2})$, we can remove the rapidly oscillating terms in Eq.(14).
In the new frame, the coefficients and Hamiltonian are, respectively,
transformed to $c(t)=\mathcal{U}^{-1}\tilde{c}(t)$ and
$\mathcal{H}=\mathcal{U}^{\dagger}\mathcal{H}\mathcal{U}-i\mathcal{U}^{\dagger}\dot{\mathcal{U}}$.
As a result, we can describe the level transition due to the tidal interaction
by the following equation with the Hamiltonian;

$\displaystyle
i\frac{\mathrm{d}c_{i}}{\mathrm{d}t}=\sum_{j=1,2}\mathcal{H}_{ij}c_{j}~{},$
(20) $\displaystyle\mathcal{H}=\left(\begin{array}[]{cc}\pm\frac{\Delta
m}{2}(\Omega(t)-\Omega_{0})&\eta\\\ \eta&\mp\frac{\Delta
m}{2}(\Omega(t)-\Omega_{0})\end{array}\right)~{},$ (23)

where we defined the resonance frequency by

$\Omega_{0}=\pm\frac{\Delta E}{\Delta m}~{}.$ (24)

###  2.2 Linear orbital evolution

For the binary BH system, the orbital angular velocity gradually increases due
to the GW radiation reaction. During the inspiral phase, the time evolution of
the orbital angular velocity is given by PhysRev.131.435 ; Blanchet:2013haa

$\displaystyle\frac{\mathrm{d}\Omega}{\mathrm{d}t}=\gamma\left(\frac{\Omega}{\Omega_{0}}\right)^{11/3}~{},$
(25)

with

$\displaystyle\frac{\gamma}{\Omega_{0}^{2}}=\frac{96}{5}\frac{q}{(1+q)^{1/3}}(M\Omega_{0})^{5/3}~{}.$
(26)

If we approximate $\Omega(t)$ linearly around the resonance frequency;

$\Omega(t)=\Omega_{0}+\gamma t~{},$ (27)

the Hamiltonian (23) becomes

$\mathcal{H}=\left(\begin{array}[]{cc}\pm\frac{\Delta m}{2}\gamma t&\eta\\\
\eta&\mp\frac{\Delta m}{2}\gamma t\end{array}\right)~{}.$ (28)

The problem defined by the the Hamiltonian (28) is known as the Landau Zener
problem 10011873546 ; Zener:1932ws  and the dynamics of this system is
characterized by the parameter

$z\equiv\frac{\eta^{2}}{|\Delta m|\gamma}~{},$ (29)

which measures the adiabaticity of the transition. This set of equations have
an exact solution, and we can calculate the persistent rate
$\Gamma=|c_{1}(+\infty)|^{2}/c_{1}(-\infty)|^{2}$ analytically. With the
initial condition $c_{1}(-\infty)=1$ and $c_{2}(-\infty)=0$, the persistent
rate is given by Baumann:2019ztm

$\Gamma=\exp(-2\pi z)~{}.$ (30)

This means that almost all axions are transferred to another mode in the
adiabatic limit, $z\gg 1$.

###  2.3 Higher-multipole moments

In Ref. Baumann:2019ztm , the main focus is on the effect of leading
quadrupole tidal perturbation; $l_{\ast}=2$ in Eq.(10). In this case, only the
transition with $\Delta m\leq 2$ is allowed because of the selection rules
reflecting the conservation of the angular momentum; $-m_{1}+m_{\ast}+m_{2}=0$
and $|m_{\ast}|\leq l_{\ast}$. If we take into account the effect of higher-
multipole moments $l_{\ast}>2$, a wider class of transitions can occur. From
Eq.(24), we can see that, if we have the same order of the difference between
the energy levels, $\Delta E$, the larger $\Delta m$ and $l_{\ast}$, the
smaller resonance frequency $\Omega_{0}$ becomes. Therefore, if the magnitude
of the perturbation of higher multipole moments is sufficiently large, since
their resonances occur earlier, they might dominate the transition. In the
following, we examine to which mode the axions that initially occupy
$\ket{n_{1}l_{1}m_{1}}=\ket{211}$ make the first transition effectively.

For the co-rotating case, the transition with $\Delta E/\Delta m>0$ occurs. In
particular, the transition between two levels separated by the hyperfine
split, i.e., $\ket{211}\to\ket{21\mbox{-1}}$, is possible. As one can see from
Eq.(7), this energy gap $|\Delta E|$ is very tiny, since it is suppressed by a
factor of $\mathcal{O}(\alpha^{4})$ compared to the gap between different $n$
modes, which is $\mathcal{O}(\alpha^{2})$. For this reason, even if we take
higher-multipole moments into account, binaries evolution hits the frequency
for the transition due to the quadrupole moment first.

For the counter-rotating case, the transition with $\Delta E/\Delta m<0$
occurs. In this case, only the transition to the level with a higher $n$
occurs. Since $\Delta E$ is the same order for all transitions that are
associated with the change of $n$, the transition with a large $|\Delta m|$
works earlier. The question that we need to answer is which is the largest
$l_{\ast}$ that starts to contribute to the transition effectively, because
the amplitude of perturbation $\eta$ is suppressed like $\propto
R_{\ast}^{-(l_{\ast}+1)}$. As an example, we show the evolution of the
particle number $|c_{1}|^{2}$ for $\alpha=0.1,q=1$ in Fig. 1. This evolution
takes into account the higher-multipole moments, and $|c_{1}|^{2}$ is reduced
by multiplying the persistent rate (30) at each resonance frequency for all
transitions satisfying the selection rules. For this set of parameters, the
$(l_{\ast}=)7$-th moment is the first one that gives a non-negligible
contribution to the transition, and a part of axions in $\ket{211}$ are
transferred to $\ket{76-6}$. Successively, the axions remaining in $\ket{211}$
are gradually transferred to the levels with a high $n$ (e.g.,
$\ket{86-6},\ket{96-6},\cdots$).

|
---|---
Figure 1: Evolution of the normalized particle number occupying $\ket{211}$,
$|c_{1}|^{2}$, for the co-rotating case (left) and the counter-rotating case
(right) with $\alpha=0.1,q=1$. The blue solid curves show the transition
taking into account the higher-multipole moments. The yellow dashed curves
show the transition due to the quadrupole moment. For the co-rotating case,
they are the same. While, for counter-rotating case, $(l_{\ast}=)7$-th moment
first gives a non-negligible contribution to the transition. The right side of
the arrows shows the transition destinations for each step.

Next, we discuss the dependence on the parameter set $(\alpha,q)$. For a given
$l_{\ast}$ and an initial state $(n_{1},l_{1},m_{1})$, we consider the
combinations of $(n_{2},l_{2},m_{2})$ for the transition destination that
minimizes the resonance frequency, i.e. the energy difference $\Delta E$
divided by $\Delta m$. First, to maximize $\Delta m$, $m_{\ast}$ is determined
to be $\pm l_{\ast}$, then we have $m_{2}=m_{1}\pm l_{\ast}$. From Eq.(7),
$\omega_{R}$ is minimized for the minimum $n$ and $l$ that we can take. Hence,
$\Delta E$ is minimized for $n_{2}=|m_{2}|+1$ and $l_{2}=|m_{2}|$, assuming
that the hyperfine split is negligible. Considering the transition from
$\ket{211}$ and calculating the persistent rate from the adiabaticity $z$, we
determine the region in which the $l_{\ast}$-th moment can contribute to the
transition in the parameter space of $(\alpha,q)$ as shown in Fig. 2 111For
the counter-rotating case, we consider the transition
$\ket{211}\to\ket{31\mbox{-1}}$ with $l_{\ast}=2$ moment.. We set the
threshold as $z=0.01$. In the parameter region below the curves, the
transition rate, which is determined by $z$, is so small that one would be
able to neglect its contribution to the cloud depletion. As one can see from
this figure, it is important to take into account the higher-multipole
moments, at least for almost equal mass binaries.

|
---|---
Figure 2: Boundary curves above which the $l_{\ast}$-th moment can give a non-
negligible contribution to the transitions for the co-rotating case (left) and
for the counter-rotating case (right). The solid curves show the boundaries
where the adiabaticity $z$ is equal to 0.01 for the transition with the
smallest resonance frequency for the respective specified values of
$l_{\ast}$. In the parameter region above the curve, the transition due to the
$l_{\ast}$-th moment effectively works.

Here, we confirmed that higher-multipole moments play an important role in the
level transition during the binary inspiral. However, as we will see in the
next section, we will find that the backreaction due to the angular momentum
transfer cannot be ignored for both co-rotating and counter-rotating cases.
Therefore, we should understand that the discussion presented here is not the
final conclusion about the transition destination.

##  3 Backreaction

In this section, we investigate the backreaction caused by the angular
momentum transfer associated with the level transition. When higher-multipole
moments are considered, the transitions occur at smaller resonance frequencies
compared with the case in which only the quadrupole moment is considered.
Therefore, since the change rate of the orbital angular momentum is smaller,
the effect of the backreaction becomes more important. We also investigate the
effect of the backreaction on the hyperfine split due to the angular momentum
transfer between the axion cloud and the orbital motion. In addition, we
analytically derive a formula for the persistent rate taking into account
these backreaction effects.

###  3.1 Non-linear Hamiltonian

For the system that we are considering, the conservation of angular momentum
reads

$\frac{\mathrm{d}}{\mathrm{d}t}\left(J_{\mathrm{binary}}\pm
J_{\mathrm{c}}\right)=-{\cal T}_{\mathrm{GW}}~{},$ (31)

where the right hand side is the torque caused by the radiation reaction due
to the gravitational wave emission. We choose the signs of $J_{\rm binary}$
and ${\cal T}_{GW}$ to be positive, irrespectively of the direction of the
orbital rotation. The upper (lower) sign corresponds to the co-rotating
(counter-rotating) case. The angular momentum of the cloud can be written as

$J_{\mathrm{c}}=J_{\mathrm{c},0}\left[m_{1}|c_{1}(t)|^{2}+m_{2}|c_{2}(t)|^{2}\right]~{},$
(32)

where we defined $J_{c,0}$ so that the initial angular momentum of the cloud
is given by $m_{1}J_{c,0}$. When we consider the transition between the modes
with a small decay rate, we can ignore the time variation of the BH spin
during the transition. Then, we have the time evolution of the orbital angular
velocity around the resonance frequency as Baumann:2019ztm

$\frac{\mathrm{d}\Omega}{\mathrm{d}t}=\gamma\pm
3R_{J}\Omega_{0}\frac{\mathrm{d}}{\mathrm{d}t}\left[m_{1}|c_{1}|^{2}+m_{2}|c_{2}|^{2}\right]~{},$
(33)

where

$R_{J}=\frac{J_{c,0}}{M^{2}}\frac{(1+q)^{1/3}}{q}(M\Omega_{0})^{1/3}~{}.$ (34)

Integrating Eq.(33) with the initial condition $\Omega(t)\to\Omega_{0}+\gamma
t$ at $t\to-\infty$, we have

$\Omega(t)=\Omega_{0}+\gamma t\pm
3R_{J}\Omega_{0}\left[m_{1}(|c_{1}|^{2}-1)+m_{2}|c_{2}|^{2}\right]~{}.$ (35)

Therefore, from Eq.(23), the level transition with backreaction is described
by the following non-linear Hamiltonian

$\mathcal{H}=\left(\begin{array}[]{cc}\theta(t)&\eta\\\
\eta&-\theta(t)\end{array}\right)~{},$ (36) $\theta(t)=\pm\frac{\Delta
m}{2}\gamma t+\frac{3|\Delta m|^{2}}{2}R_{J}\Omega_{0}|c_{2}|^{2}~{},$ (37)

where we used $|c_{1}|^{2}+|c_{2}|^{2}=1$. We now introduce a parameter which
represents the strength of the non-linearity,

$g\equiv\frac{3|\Delta m|^{2}R_{J}\Omega_{0}}{2\eta}~{}.$ (38)

As an example, we show the value of the non-linearity $g$ for some transitions
with $q=1$ in Fig. 3. When higher-multipole moments are concerned, the effect
of the non-linearity can be especially large. For this reason, it is important
to investigate the transition including this backreaction.

|
---|---
Figure 3: Non-linearity $g$ of the transition from $\ket{211}$ for the co-
rotating case (left) and the counter-rotating case (right) with $q=1$. Each
line corresponds to a different transition destination.

As shown in Ref. Baumann:2019ztm , there are two cases in the way how the
backreaction affects the evolution of the orbital motion. They are classified
by the relative sign of the terms in Eq.(37). Case(1) is the case where the
relative sign is positive. This occurs for the ”co-rotating and $\Delta m>0$”
case or for the ”counter-rotating and $\Delta m<0$” case. In case(1), the
cloud receives the angular momentum from the orbital motion of the binary, and
the orbit sinks faster. Hence, the binary inspiral passes through the
resonance band more quickly, which makes the persistent rate larger. Case(2)
is the case where the relative sign is negative. This occurs for the ”co-
rotating and $\Delta m<0$” case or for the ”counter-rotating and $\Delta m>0$”
case. In case(2), the cloud gives the angular momentum to the orbital motion,
and the orbital evolution tends to stagnate. Hence, the binary inspiral passes
through the resonance band more slowly, which makes the transition more
adiabatic and the persistent rate smaller.

The first transition during the binary inspiral for the counter-rotating case
belongs to the case(1). On the other hand, the first transition between
hyperfine split for the co-rotating case seems to belong to the case(2). We
will discuss this point in detail in Sec. 3.2.

###  3.2 Backreaction to the hyperfine split through the change in the
geometry

We consider a cloud that occupies only the level with the maximum superradiant
growth rate. Usually, we can ignore the backreaction to the geometry when we
discuss the evolution of the cloud because of the low energy density of the
cloud Brito:2014wla . However, focusing on the very small hyperfine split,
which is proportional to the angular momentum, we may not be able to ignore
the backreaction due to the change of the angular momentum of the saturated
cloud, because it is expected to be comparable to that of the central BH. To
take this effect into account, we model the frequency shift related to the
coupling with the background angular momentum, which corresponds to the fifth
term in Eq.(7), as

$\delta\omega^{h}(t)=\mu\frac{2m\alpha^{5}}{n^{3}l(l+1/2)(l+1)}\frac{J}{M^{2}}\left(1+{x}\frac{J_{\mathrm{c}}(t)}{J}\right)~{}.$
(39)

We parameterized the efficiency of the frequency shift due to the angular
momentum of the cloud itself by ${x}$, which should be $\mathcal{O}(1)$. In
fact, we evaluate the value of $x$ in Appendix. B for the flat background
222If we consider the maximally grown cloud, it will have a larger angular
momentum than that of the BH for $\alpha\ll 1$. Hence, the frequency shift due
to the angular momentum of the cloud becomes dominant, and the resonance
frequency deviates from the one calculated from Eq.(7). This means that we
should consider the deviation from the Kerr metric when we need the precise
value of the resonance frequency. However, it would not change our results
qualitatively, and we leave this problem as a future work. .

In this section our interest is in the first transition for the co-rotating
case, $\ket{211}\to\ket{21\mbox{-1}}$, i.e., $\Delta m<0$. So far, we have
considered the resonance frequency $\Omega_{0}$ to be a constant. However, in
this case, since the level gap varies with time, we consider the resonance
frequency to be time-dependent 333It is only in this section that the
resonance frequency to be time-dependent.. We define the initial resonance
frequency by

$\Omega_{0}^{(-\infty)}=\mu\frac{2\alpha^{5}}{n^{3}l(l+1/2)(l+1)}\frac{J}{M^{2}}\left(1+x\frac{J_{\mathrm{c}}(-\infty)}{J}\right)~{},$
(40)

where $J_{\mathrm{c}}(-\infty)=m_{1}J_{\mathrm{c},0}$. Then, the time-
dependent resonant frequency becomes

$\Omega_{0}(t)=\Omega_{0}^{(-\infty)}+\frac{xJ_{\mathrm{c},0}}{J+xm_{1}J_{\mathrm{c},0}}\Omega_{0}^{(-\infty)}\left[m_{1}(|c_{1}(t)|^{2}-1)+m_{2}|c_{2}(t)|^{2}\right]~{}.$
(41)

This implies that the resonance band shrinks as the angular momentum is
transferred from the cloud to the orbital motion. Thus, we have the diagonal
term of the Hamiltonian (23) that describes the transition between the levels
separated by the hyperfine split as

$\displaystyle\theta^{h}(t)$ $\displaystyle=-\frac{|\Delta
m|}{2}\left(\Omega(t)-\Omega_{0}(t)\right)$ $\displaystyle=-\frac{|\Delta
m|}{2}\gamma t+\frac{|\Delta
m|^{2}}{2}\Omega_{0}^{(-\infty)}\left(3R_{J}-\frac{xJ_{\mathrm{c},0}}{J+xm_{1}J_{\mathrm{c},0}}\right)|c_{2}|^{2}~{}.$
(42)

When we ignore the backreaction on the geometry, i.e., for $x=0$, this
transition belongs to the case(2). However, if
$3R_{J}<xJ_{\mathrm{c},0}/(J+xm_{1}J_{\mathrm{c},0})$, the relative sign
reverses and this transition belongs to the case(1). In Fig. 4, we show the
parameter region that satisfies
$3R_{J}<xJ_{\mathrm{c},0}/(J+xm_{1}J_{\mathrm{c},0})$. For example, in the
case of equal mass binaries, this condition is satisfied and the transition
effectively belongs to the case(1). In other words, since the effect of the
shrinking the resonance frequency works more efficiently than the effect of
slowing down the evolution of the orbit, the resonance band is passed through
quickly and the persistent rate becomes larger.

In this hyperfine transition, the non-linearity $g$ given in Eq.(38) should be
replaced with the effective one,

$g_{\mathrm{eff}}\equiv\frac{|\Delta
m|^{2}}{2\eta}\left|3R_{J}-\frac{xJ_{\mathrm{c},0}}{J+xm_{1}J_{\mathrm{c},0}}\right|\Omega_{0}^{(-\infty)}~{}.$
(43)

In particular, for the transition $\ket{211}\to\ket{21\mbox{-1}}$ with
$3R_{J}<xJ_{\mathrm{c},0}/(J+xm_{1}J_{\mathrm{c},0})$, the effective non-
linearity can be approximately evaluated as

$g_{\mathrm{eff}}\simeq 3\times
10^{3}\frac{1+q}{q}\left(\frac{0.1}{\alpha}\right)^{3}\frac{xJ_{{\rm
c},0}/M^{2}}{(J/M^{2}+xJ_{{\rm c},0}/M^{2})^{2}}~{}.$ (44)

Therefore, even for equal mass binaries, this transition is in the strongly
non-linear regime. This means that the persistent rate is large and most of
the axions remain in the state $\ket{211}$. The resonances that the binary
experiences after the first hyperfine resonance are the ones associated with
the change of $n$ and $\Delta m>0$. Hence, they belong to the case(1). As a
result, when we consider the binary system with $q=\mathcal{O}(1)$, almost all
transitions from $\ket{211}$ belong to the case(1) and are in the strongly
non-linear regime.

Figure 4: Parameter region that satisfies
$3R_{J}<xJ_{\mathrm{c},0}/(J+xm_{1}J_{\mathrm{c},0})$ for
$\ket{211}\to\ket{21\mbox{-1}}$ with
$J/M^{2}=a_{\mathrm{crit}}/M,J_{\mathrm{c},0}/M^{2}=1-J/M^{2},x=7/32$. In the
shaded region, this hyperfine transition belongs to the case(1) due to the
shrink of the level gap according to the angular momentum transfer.

###  3.3 Persistent rate for the strongly non-linear regime

Focusing on the binary system with $q=\mathcal{O}(1)$, we investigate the
level transition in the strongly non-linear regime $g\gg 1$ for the case(1).
In particular, we derive the self-consistent equation for the persistent rate
by referring to Ref. PhysRevA.66.023404 .

For the case(1), $\theta(t)$ in the diagonal elements of the Hamiltonian
Eq.(36) is

$\theta(t)=\frac{|\Delta m|}{2}\gamma t+\frac{3|\Delta
m|^{2}}{2}R_{J}\Omega_{0}|c_{2}|^{2}~{}.$ (45)

First, we transform the coefficients as

$\displaystyle c_{1}=\bar{c}_{1}e^{-i\int_{0}^{t}dt^{\prime}\theta}~{},$ (46)
$\displaystyle c_{2}=\bar{c}_{2}e^{+i\int_{0}^{t}dt^{\prime}\theta}~{}.$ (47)

By this transformation, the diagonal elements in Eq.(36) vanish, and the time
evolution of $\bar{c}_{2}$ is determined by

$\frac{\mathrm{d}\bar{c}_{2}}{\mathrm{d}t}=\frac{\eta}{i}\bar{c}_{1}e^{-2i\int_{0}^{t}dt^{\prime}\theta}~{}.$
(48)

Integrating this equation with the initial condition $\bar{c}_{2}(-\infty)=0$,
we formally obtain

$\bar{c}_{2}(t)=\frac{\eta}{i}\int_{-\infty}^{t}dt^{\prime}\bar{c}_{1}(t^{\prime})\exp\left[-2i\int_{0}^{t^{\prime}}dt^{\prime\prime}\theta(t^{\prime\prime})\right]~{}.$
(49)

Because $g\gg 1$ and the exponent in the above expression shows rapid
oscillations, we evaluate the above integral using the stationary phase
approximation. We denote the second derivative coefficient at the stationary
point $t_{0}$ by

$\displaystyle|\Delta m|\bar{\gamma}$
$\displaystyle\equiv\left.\frac{\mathrm{d^{2}}}{\mathrm{d}t^{2}}\left(2\int_{0}^{t}dt^{\prime}\theta(t^{\prime})\right)\right|_{t_{0}}$
$\displaystyle=|\Delta m|\gamma+3|\Delta
m|^{2}R_{J}\Omega_{0}\left.\frac{\mathrm{d}|c_{2}|^{2}}{\mathrm{d}t}\right|_{t_{0}}~{}.$
(50)

Then, we have an expression for the particle number as

$\displaystyle|c_{2}(t)|^{2}=$
$\displaystyle\eta^{2}|\bar{c}_{1}(t_{0})|^{2}\left(\int_{-\infty}^{t}dt^{\prime}\exp\left[\frac{i}{2}|\Delta
m|\bar{\gamma}(t^{\prime}-t_{0})^{2}\right]\right)$ (51)
$\displaystyle\qquad\qquad\qquad\times\left(\int_{-\infty}^{t}dt^{\prime}\exp\left[-\frac{i}{2}|\Delta
m|\bar{\gamma}(t^{\prime}-t_{0})^{2}\right]\right)~{}.$ (52)

From this expression, we get the form of Fresnel integral and can calculate
the following quantities;

$\left.\frac{\mathrm{d}|c_{2}|^{2}}{\mathrm{d}t}\right|_{t_{0}}=\eta^{2}|\bar{c}_{1}(t_{0})|^{2}\sqrt{\frac{\pi}{|\Delta
m|\bar{\gamma}}}~{},$ (53)
$|c_{2}(+\infty)|^{2}=\eta^{2}|\bar{c}_{1}(t_{0})|^{2}\frac{2\pi}{|\Delta
m|\bar{\gamma}}~{}.$ (54)

Using the above results, we evaluate the persistent rate
$\Gamma=|c_{1}(+\infty)|^{2}=1-|c_{2}(+\infty)|^{2}$ self-consistently. Here,
$\Gamma$ is nearly one because of the strong non-linearity, so we approximate
$|\bar{c}_{1}(t_{0})|^{2}\simeq 1$. From Eq.(54), we get

$\frac{1}{1-\Gamma}=\frac{|\Delta m|\bar{\gamma}}{2\pi\eta^{2}}~{}.$ (55)

Substituting Eq.(3.3) into this equation and using Eq.(53), finally we obtain
a closed form of the equation for $\Gamma$ as

$\frac{1}{1-\Gamma}=\frac{1}{2\pi z}+\frac{g}{\sqrt{2}\pi}\sqrt{1-\Gamma}~{},$
(56)

where the defenitions of $z$ and $g$ are given by Eq.(29) and Eq.(38),
respectively. Solving this equation numerically, we can estimate the
persistent rate including the backreaction without solving non-linear
differential equations. We compare the solution of Eq.(56) with the numerical
result for the persistent rate obtained by solving Eq.(36) directly in Fig. 5.
We can see that Eq.(56) gives a good approximation of the persistent rate with
backreaction for $g\gg 1$.

|
---|---
Figure 5: Persistent rate $\Gamma$ including backreaction whose strength is
characterized by the non-linearity $g$ with $z=5$ (left) and $z=0.1$ (right).
The red solid line shows the solution of the self-consistent equation (56) for
strong non-linear regime. The blue dots show the numerical results obtained by
solving differential equations described by the non-linear Hamiltonian (36).

Using this formula, we investigate how the evolution of the particle number
that we have shown in Fig. 1 is modified. We show the results obtained by
replacing the persistent rate Eq.(30) with the one evaluated using Eq.(56) in
Fig. 6 444In the case(1) with strong non-linearity, transition becomes non-
adiabatic and particle number changes with oscillations. Since the orbital
frequency also oscillates at the transition Baumann:2019ztm , strictly
speaking, it may not be appropriate to treat the orbital frequency as a time
coordinate. 555In this figure, we use the persisitent rate obtained by solving
Eq.(36) directly for transitions with $g<10$. . For the counter-rotating case,
the qualitative results do not change. By contrast, for the co-rotating case,
the axion cloud passes through the resonance corresponding to the transition
$\ket{211}\to\ket{21\mbox{-1}}$ with little transition. Later, the axions are
transferred to higher levels as in the counter-rotating case.

When the mass ratio is small, i.e., $q\ll 1$, the first transition for the co-
rotating case $\ket{211}\to\ket{21\mbox{-1}}$ belongs to the case(2), even if
we take into account the reduction of the hyperfine split due to the angular
momentum loss of the cloud. In this case, we cannot apply the method we
introduced in this section straightforwardly because of the large change in
the number of particles occupying the initial state. In addition, for such a
transition, since the decay rate of the transition destination is large
compared with the timescale of the evolution of the binary orbital frequency,
we may not be able to treat the level transition as usual. We leave a study on
the extreme mass ratio inspiral case as future work.

|
---|---
Figure 6: Evolution of the normalized particle number $|c_{1}|^{2}$ occupying
$\ket{211}$, including the backreaction for the co-rotating case (left) and
the co-rotating case (right) with $\alpha=0.1,q=1$. The red solid curve shows
the transition taking into account the backreaction and the higher-multipole
moments. The blue and yellow curves are the same as in Fig. 1. Note that, for
the co-rotating case, since the transition between hyperfine split is passed
through with little transition, the blue and yellow lines are almost zero in
this range.

##  4 Axion emission

So far, we discussed the first transition destination and found that the
axions in $\ket{211}$ are transferred to higher levels for both directions of
the orbital rotation, at least for $q=\mathcal{O}(1)$. At this point, since
the decay (growth) rates of the transition destinations are generally very
small compared with the timescale of the binary evolution, axions are not
reabsorbed by the BH before the merger. However, axions once excited to higher
levels can be easily further excited by the tidal interaction to an unbounded
state. In this section, we evaluate the particle number flux to the infinity.

###  4.1 Particle number flux

We evaluate the particle number flux perturbatively. We solve the following
approximate equations,

$\displaystyle\left(i\partial_{t}+\frac{1}{2\mu}\nabla^{2}+\frac{\alpha}{r}\right)\psi^{b}$
$\displaystyle=0~{},$ (57)
$\displaystyle\left(i\partial_{t}+\frac{1}{2\mu}\nabla^{2}+\frac{\alpha}{r}\right)\psi^{r}$
$\displaystyle=V_{\ast}\psi^{b}~{},$ (58)

where

$\psi^{b}(x)=e^{-i(\omega_{0}-\mu)t}R_{n_{0}l_{0}}(r)Y_{l_{0}m_{0}}(\theta,\varphi)$
(59)

is the background bound state solution and $\psi^{r}$ is the radiative part.
Here, we normalize the number of particles that occupy the background state as

$\int d^{3}x|\psi^{b}(x)|^{2}=1~{}.$ (60)

We construct the solution for the radiative part perturbatively using Green’s
function as

$\psi^{r}(x)=\int
d^{4}x^{\prime}G(x,x^{\prime})V_{\ast}\psi^{b}(x^{\prime})~{}.$ (61)

Here, the Green’s function is defined to satisfy

$\left(i\partial_{t}+\frac{1}{2\mu}\nabla^{2}+\frac{\alpha}{r}\right)G(x,x^{\prime})=\delta^{4}(x-x^{\prime})~{}.$
(62)

We can decompose the Green’s function as

$G(x,x^{\prime})=\sum_{lm}\int\frac{d\omega}{2\pi}e^{-i(\omega-\mu)(t-t^{\prime})}Y_{lm}(\theta,\varphi)Y_{lm}^{\ast}(\theta^{\prime},\varphi^{\prime})G_{l}^{\omega}(r,r^{\prime})~{}.$
(63)

The radial part of the Green’s function $G_{l}^{\omega}(r,r^{\prime})$ can be
constructed using the Coulomb wavefunction $R_{kl}(r)$ that satisfies

$\frac{1}{2\mu}\left(\frac{1}{r^{2}}\frac{{\rm d}}{{\rm d}r}r^{2}\frac{{\rm
d}}{{\rm
d}r}-\frac{l(l+1)}{r^{2}}+k^{2}+\frac{2\mu\alpha}{r}\right)R_{kl}(r)=0~{},$
(64)

where $k=\sqrt{2\mu(\omega-\mu)}$. We have

$\displaystyle
G_{l}^{\omega}(r,r^{\prime})=\frac{2\mu}{W_{kl}}\left(R^{0}_{kl}(r)R^{+}_{kl}(r^{\prime})\theta(r^{\prime}-r)+R^{0}_{kl}(r^{\prime})R^{+}_{kl}(r)\theta(r-r^{\prime})\right)~{},$
(65)

with

$\displaystyle W_{kl}=r^{2}\left(R^{0}_{kl}\frac{{\rm d}R^{+}_{kl}}{{\rm
d}r}-R^{+}_{kl}\frac{{\rm d}R^{0}_{kl}}{{\rm d}r}\right)~{}.$ (66)

$R_{kl}^{0}$ is regular at the origin and $R_{kl}^{+}(r)$ is the outgoing wave
at infinity (see Appendix A for the explicit form). Now, assuming the constant
orbital angular velocity of the binary, i.e., $\Phi_{\ast}(t)=\pm\Omega t$, we
consider the tidal field given by

$V_{\ast}=\sum_{l_{\ast}m_{\ast}}V_{\ast,l_{\ast}m_{\ast}}e^{\mp
im_{\ast}\Omega t}~{}.$ (67)

Then, we obtain the perturbative solution,

$\displaystyle\psi^{r}$ $\displaystyle=\sum_{lm}\sum_{l_{\ast}m_{\ast}}\int
d^{4}x^{\prime}\int\frac{d\omega}{2\pi}e^{-i(\omega-\mu)(t-t^{\prime})}Y_{lm}(\theta,\varphi)Y_{lm}^{\ast}(\theta^{\prime},\varphi^{\prime})$
$\displaystyle\qquad\qquad\times
V_{\ast,l_{\ast}m_{\ast}}(r^{\prime},\theta^{\prime},\varphi^{\prime})e^{\mp
im_{\ast}\Omega
t^{\prime}}G_{l}^{\omega}(r,r^{\prime})e^{-i(\omega_{0}-\mu)t^{\prime}}R_{n_{0}l_{0}}(r^{\prime})Y_{l_{0}m_{0}}(\theta^{\prime},\varphi^{\prime})$
$\displaystyle\to-\sum_{lm}\sum_{l_{\ast}m_{\ast}}\int d\omega\
\delta(\omega-\omega_{0}\mp
m_{\ast}\Omega)e^{-i(\omega-\mu)t}\frac{e^{ikr}}{r}Y_{lm}(\theta,\varphi)\frac{2\mu}{k}$
$\displaystyle\qquad\qquad\times\int
dr^{\prime}r^{\prime}{}^{2}d\theta^{\prime}\sin\theta^{\prime}d\varphi^{\prime}R^{0}_{kl}(r^{\prime})Y^{\ast}_{lm}(\theta^{\prime},\varphi^{\prime})V_{\ast,l_{\ast}m_{\ast}}(r^{\prime},\theta^{\prime},\varphi^{\prime})R_{n_{0}l_{0}}(r^{\prime})Y_{l_{0}m_{0}}(\theta^{\prime},\varphi^{\prime})\,,$
(68)

where we take the limit $r\to\infty$ and choose the normalization of
$R^{0}_{kl}$ as given in Appendix A. Therefore, we can calculate the
normalized particle number flux to infinity as

$\displaystyle F$ $\displaystyle=\int d\theta\sin\theta
d\varphi\frac{r^{2}}{2\mu
i}\left(\psi^{r\ast}\partial_{r}\psi^{r}-\psi^{r}\partial_{r}\psi^{r\ast}\right)$
$\displaystyle=\sum_{lm}\sum_{l_{\ast}m_{\ast}}\frac{4\mu}{\sqrt{2\mu((\omega_{0}\pm
m_{\ast}\Omega)-\mu)}}\left|\int
d^{3}x^{\prime}R^{0}_{kl}Y^{\ast}_{lm}V_{\ast,l_{\ast}m_{\ast}}R_{n_{0}l_{0}}Y_{l_{0}m_{0}}\right|^{2}~{}.$
(69)

###  4.2 Results

Here, we show how efficiently this emission works in comparison with the
timescale of the binary evolution;

$T_{\mathrm{binary}}\equiv\frac{\Omega}{\gamma}=\frac{5}{96}M\frac{(1+q)^{1/3}}{q}(M\Omega)^{-8/3}~{}.$
(70)

First, for simplicity, we give an analytic estimate of the flux from the cloud
composed of $\ket{211}$ state in the limit $k\to 0$. This limit corresponds to
the case where the binary orbital frequency is just at the threshold for the
excitation to an unbound state. To estimate the flux, we consider the
component that gives the maximum amplitude of the tidal field, i.e., the flux
due to the $(l_{\ast},m_{\ast})=(2,\pm 2)$ perturbation. Hence,

$M\Omega\simeq\frac{1}{|\Delta
m|}\frac{\alpha^{3}}{2n^{2}}=\frac{\alpha^{3}}{16}~{}.$ (71)

At this frequency, the timescale of the binary’s evolution is given by

$T_{\mathrm{binary}}\simeq 8\times
10^{9}\frac{(1+q)^{1/3}}{q}M\left(\frac{0.1}{\alpha}\right)^{8}~{}.$ (72)

In the limit $k\to 0$, the radiative wave function is given by

$R_{kl}^{0}\to\sqrt{k}\sqrt{\pi/r}J_{2l+1}(\sqrt{8\mu\alpha r})~{},$ (73)

where $J_{\nu}(x)$ is the Bessel function. Then, we can roughly estimate the
flux as 666Here, we consider only the inner part of perturbation in Eq.(10)
and underestimate the actual value. However, the flux is still so large that
it should not affect the conclusion.

$\displaystyle F$ $\displaystyle\simeq 2\times
10^{-5}\frac{1}{M}\frac{q^{2}}{(1+q)^{2}}\left(\frac{\alpha}{0.1}\right)^{3}\quad\textrm{(co-
rotating)}~{},$ (74) $\displaystyle F$ $\displaystyle\simeq 4\times
10^{-6}\frac{1}{M}\frac{q^{2}}{(1+q)^{2}}\left(\frac{\alpha}{0.1}\right)^{3}\quad\textrm{(counter-
rotating)}~{}.$ (75)

The particle number flux divided by the number of particles contained in the
cloud is sufficiently large compared with the timescale of the binary
evolution, except for the extreme mass ratio inspiral case with $q\ll 1$. It
means that if the orbital frequency reaches the threshold for the transition
to an unbound state, bounded axions are all radiated away.

In practice, as we showed in the previous section, axions are first
transferred to the levels with some higher $l$. Hence, We show the numerical
results of the comparison between the flux from higher levels and the
timescale of the orbital evolution due to radiation reaction in Fig. 7 for
$q=1$. The orbital frequency is set to the resonance frequency for the
transition from $\ket{211}$ to each level. As one can see from this figure,
axions transferred to higher levels will be radiated away to the infinity
rather immediately. Therefore, we conclude that the axion cloud disappears
just after the first transition, through the outward radiation without being
reabsorbed by the BH for equal mass binaries.

|
---|---
Figure 7: Comparison between the normalized particle number flux and the
timescale of the binary’s evolution for both the co-rotating case (left) and
the counter-rotating case (right). They are evaluated at the resonance
frequency for the transition from $\ket{211}$ to each level for $q=1$.

##  5 Summary and Discussion

In this paper, we studied the dynamics of the axion cloud associated with a BH
in a binary system during the inspiral phase, for the purpose to give a
reliable constraint on the properties of axions from observations of GWs and
BHs. The presence of the cloud should affect the gravitational waveform,
especially if it remains until a late inspiral phase close to the merger time.
In addition, the level transition of the cloud due to the binary tidal
interaction may change the effective forbidden region in the mass and spin
parameters of BHs observed in binary systems. Therefore, in this paper we
clarified in which case the cloud disappears, and when it disappears, when and
how the process proceeds. We mainly focused on equal mass binaries, which are
the main target of the ground based gravitational wave detectors, and extended
the previous studies in the following manner.

First, we studied level transitions taking into account the higher-multipole
moments of the tidal perturbation without the backreaction to the orbital
evolution. Though the magnitude of higher-multipole perturbation is smaller
than that of the leading quadrupole moment, they can be sufficiently strong to
make a transition by a non-negligible amount. The resonance frequency is
roughly inversely proportional to the difference of the azimuthal number
between two levels, unless the difference of the energies between two levels
is governed by the hyperfine split. Therefore, higher-multipole moments can
become at work earlier during the inspiral of the binary orbit. We confirmed
that this picture is correct and showed the parameter dependence.

Second, we investigated the effects of backreaction due to the angular
momentum transfer associated with the level transition, precisely. We
described the level transition with this backreaction by introducing a non-
linear model Hamiltonian and analyzed it. In particular, for the transition
between two levels gapped by the hyperfine split, which is proportional to the
angular momentum of the BH and the cloud, we considered the backreaction from
the transfer of the angular momentum between the cloud and the orbital motion.
This effect makes the duration staying at the resonance band shorter and the
persistent rate of the initial level larger. In addition, we derived the
analytic closed formula for the persistent rate in a strongly non-linear
regime. Using this formula, we analyzed the transitions with backreaction and
found that axions are transferred to the levels with a larger $l$ gradually
for both co-rotating and counter-rotating cases.

Finally, we considered the axion emission to the infinity from the cloud
induced by the tidal interaction. At the first transition, axions are
transferred to higher levels, which do neither decay nor grow rapidly by
themselves. Hence, axions are not reabsorbed by the BH. However, they can be
easily excited to a continuous unbound state. We calculated the normalized
particle number flux and compared it with the timescale of the evolution of
the binary orbital frequency. As a result, the flux turned out to be large
enough to deplete the axion clouds completely through the radiation to the
infinity.

From the above results, we can conclude that axion clouds in equal mass binary
systems will disappear during the early inspiral phase by the transition to
unbound states caused by the tidal perturbation. This fact may imply that it
is difficult to detect the modulation of the gravitational waveform caused by
the presence of axion clouds using ground-based gravitational wave detectors.
On the other hand, because axions are not reabsorbed by the BH, the tidal
interaction due to the companion in binaries does not affect the evolution of
the BH spin. It will imply that the current constraints on the possible axion
mass from the distribution of the mass and the spin of BHs will not be
altered, even if we use the estimates of these parameters derived from
gravitational wave observations of binary black holes. In this paper, we have
concentrated on the theoretical aspects, and we leave more careful examination
of possible implications to observations, along with the case of extreme mass
ratio inspirals, for our future study.

## Acknowledgment

This work is supported by JSPS Grant-in-Aid for Scientific Research JP17H06358
(and also JP17H06357), as a part of the innovative research area,
“Gravitational wave physics and astronomy: Genesis”, and also by JP20K03928.
T. Takahashi is supported by the establishment of university fellowships
towards the creation of science technology innovation.

##  Appendix A Wavefunction for a hydrogen atom

In this appendix, we summarize the wavefunction for a hydrogen atom that we
use in this paper. We consider the Schr$\mathrm{\ddot{o}}$dinger type equation
with the Coulomb potential

$\left(i\partial_{t}+\frac{1}{2\mu}\nabla^{2}+\frac{\alpha}{r}\right)\psi=0~{}.$
(76)

The general solution of this equation can be written by the superposition of
discrete unbound states and continuous unbound states as

$\psi=\sum_{nlm}c_{nlm}\psi_{nlm}+\int dk\sum_{lm}c_{lm}(k)\psi_{klm}~{}.$
(77)

The bound state solution can be written as

$\psi_{nlm}=e^{-i(\omega_{nlm}-\mu)t}R_{nl}(r)Y_{lm}(\theta,\varphi)~{}.$ (78)

In particular, the radial wavefunction that is regular at the origin is

$R_{nl}(r)=\sqrt{\left(\frac{2}{nr_{0}}\right)^{3}\frac{(n-l-1)!}{2n(n+l)!}}e^{-\tilde{r}/n}\left(\frac{2\tilde{r}}{n}\right)^{l}L_{n-l-1}^{2l+1}\left(\frac{2\tilde{r}}{n}\right)~{},$
(79)

where

$\tilde{r}=\frac{r}{r_{0}}~{},\quad r_{0}=\frac{1}{\mu\alpha}~{}.$ (80)

$L_{p}^{q}(x)$ is the generalized Laguerre polynomial, defined by

$L_{p}^{q}(x)=\sum_{n=0}^{p}(-1)^{n}\left(\begin{array}[]{c}p+q\\\
p-n\end{array}\right)\frac{x^{n}}{n!}~{}.$ (81)

The unbound state solution can be written by

$\psi_{klm}=e^{-i(\omega(k)-\mu)t}R_{kl}(r)Y_{lm}(\theta,\varphi)~{}.$ (82)

The radial wavefunctions which are regular at the origin and are outgoing at
the infinity are, respectively,

$\displaystyle
R^{0}_{kl}(r)=\frac{F_{l}(-1/kr_{0},kr)}{r}\to\left\\{\begin{array}[]{ll}C_{l}(-1/kr_{0})k(kr)^{l}&(r\to
0)\\\ \frac{\sin{kr}}{r}&(r\to\infty)\end{array}\right.~{},$ (85)
$\displaystyle
R^{+}_{kl}(r)=\frac{H_{l}^{+}(-1/ka_{0},kr)}{r}\to\frac{e^{ikr}}{r}\ \
(r\to\infty)~{},$ (86)

where

$\displaystyle
F_{l}(\eta,\rho)=C_{l}(\eta)\rho^{l+1}e^{-i\rho}{}_{1}F_{1}(l+1-i\eta,2l+2,2i\rho)~{},$
(87) $\displaystyle
C_{l}(\eta)=\frac{2^{l}e^{-\pi\eta/2}|\Gamma(l+1+i\eta)|}{(2l+1)!}~{}.$ (88)

${}_{1}F_{1}$ is the confluent hypergeometric function, and see
mathematicalfunction  for $H_{l}^{+}$, whose explicit form is not necessary in
this paper.

##  Appendix B Frequency shift due to the angular momentum of the axion cloud

In this appendix, we derive the frequency shift due to the angular momentum of
the axion cloud itself. This shift is sourced by the $t$-$\varphi$ component
of the metric perturbation from the background as Baumann:2018vus

$\left(i\partial_{t}+\frac{1}{2\mu}\nabla^{2}+\frac{\alpha}{r}\right)\psi=i\delta
g^{t\varphi}\partial_{\varphi}\psi~{}.$ (89)

Writing the frequency $\omega=\bar{\omega}+\delta\omega$, we have

$\delta\omega=\frac{-m\int d^{3}x~{}\delta
g^{t\varphi}|\psi_{nlm}(x)|^{2}}{\int d^{3}x|\psi_{nlm}(x)|^{2}}~{}.$ (90)

If we consider the leading term of the Kerr metric $\delta
g^{t\varphi}=-2Ma/r^{3}$, we can reproduce the fifth tern of Eq.(7). Here, we
evaluate the metric perturbation due to the self-gravity of the axion cloud.
Approximating the background by Minkowski spacetime, we calculate the
stationary $t$-$\varphi$ component of the metric perturbation by

$\Delta\delta g_{ti}\simeq-16\pi T_{ti}~{},$ (91)

where $\Delta$ is the Laplacian in flat space and $T_{ti}$ is the $t$-$i$
component of the energy momentum tensor. Here, we write $x^{i}=(x,y,z)$.
Solving Eq.(91), we have

$\delta g_{t\varphi}=-16\pi\frac{\partial x^{i}}{\partial\varphi}\int
d^{3}x^{\prime}\frac{-1}{4\pi|\bm{x}-\bm{x}^{\prime}|}T_{ti}~{}.$ (92)

On the other hand, the $t$-$\varphi$ component of the energy momentum tensor
is

$T_{t\varphi}=-\frac{2\omega m}{\mu}|\psi_{nlm}(x)|^{2}~{}$ (93)

and we can evaluate the angular momentum of the axion cloud as

$J_{\mathrm{c}}\simeq-\int d^{3}xT_{t\varphi}=\frac{2\omega m}{\mu}\int
d^{3}x|\psi_{nlm}(x)|^{2}~{}.$ (94)

As a result, using $\delta g^{t\varphi}=\delta
g_{t\varphi}/r^{2}\sin^{2}\theta$, we have the frequency shift as

$\displaystyle\frac{\delta\omega}{\mu}$
$\displaystyle=\frac{4m}{\mu}J_{\mathrm{c}}\frac{\int
d^{3}xd^{3}x^{\prime}\frac{1}{r\sin\theta}\frac{1}{r^{\prime}\sin\theta^{\prime}}\frac{\cos(\varphi-\varphi^{\prime})}{|\bm{x}-\bm{x}^{\prime}|}|\psi_{nlm}(x)|^{2}|\psi_{nlm}(x^{\prime})|^{2}}{\left(\int
d^{3}x|\psi_{nlm}(x)|^{2}\right)^{2}}$
$\displaystyle=4m\alpha^{5}\frac{J_{\mathrm{c}}}{M^{2}}\sum_{LM}\frac{4\pi}{2L+1}$
$\displaystyle\quad\times\int
d\tilde{r}d\tilde{r}^{\prime}r\tilde{r}\tilde{R}_{nl}^{2}(\tilde{r})\tilde{R}_{nl}^{2}(\tilde{r}^{\prime})\left(\frac{\tilde{r}^{\prime
L}}{\tilde{r}^{L+1}}\theta(\tilde{r}-\tilde{r}^{\prime})+\frac{\tilde{r}^{L}}{\tilde{r}^{\prime
L+1}}\theta(\tilde{r}^{\prime}-\tilde{r})\right)$
$\displaystyle\quad\times\int d\theta d\varphi
d\theta^{\prime}d\varphi^{\prime}\cos(\varphi-\varphi^{\prime})Y_{LM}(\Omega)|Y_{lm}(\Omega)|^{2}Y^{\ast}_{LM}(\Omega^{\prime})|Y_{lm}(\Omega^{\prime})|^{2}~{}.$
(95)

In particular, for the cloud with $\ket{nlm}=\ket{21\pm 1}$, we have

$\frac{\delta\omega}{\mu}=\frac{7}{384}m\alpha^{5}\frac{J_{\mathrm{c}}}{M^{2}}~{}.$
(96)

In this case, the parameter we defined in Sec. 3.2 become
$x=(7/384)/(1/12)=7/32$. This parameter $x$ depends on the choice of the pair
of levels, the one that produces the metric perturbation and the one whose
energy eigenvalue is affected by it. However, since $x$ remains to be
$\mathcal{O}(1)$ anyway, we adopt this value as an approximation for our
analysis.

Note that the frequency shift due to $t$-$t$ component of the self-gravity is
also calculated in the appendix of Ref. Baryakhtar:2020gao . This effect comes
from the order of $\mu\alpha^{3}$. Hence, when we consider the transition
between different $n$ modes, which is the order of $\mu\alpha^{2}$, we can
neglect it.

## References

  * (1) Daniel Baumann, Horng Sheng Chia, Rafael A. Porto, and John Stout, Phys. Rev. D, 101(8), 083019 (2020), arXiv:1912.04932.
  * (2) Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, Phys. Rev., D81, 123530 (2010), arXiv:0905.4720.
  * (3) Peter Svrcek and Edward Witten, JHEP, 06, 051 (2006), arXiv:hep-th/0605206.
  * (4) Michael Dine and Willy Fischler, Phys. Lett., 120B, 137–141 (1983).
  * (5) John Preskill, Mark B. Wise, and Frank Wilczek, Phys. Lett., 120B, 127–132 (1983).
  * (6) L. F. Abbott and P. Sikivie, Phys. Lett., 120B, 133–136 (1983).
  * (7) Lam Hui, Jeremiah P. Ostriker, Scott Tremaine, and Edward Witten, Phys. Rev. D, 95(4), 043541 (2017), arXiv:1610.08297.
  * (8) Richard Brito, Vitor Cardoso, and Paolo Pani, Lect. Notes Phys., 971, pp.1–237 (2020), arXiv:1501.06570.
  * (9) Steven L. Detweiler, Phys. Rev., D22, 2323–2326 (1980).
  * (10) Sam R. Dolan, Phys. Rev., D76, 084001 (2007), arXiv:0705.2880.
  * (11) Asimina Arvanitaki and Sergei Dubovsky, Phys. Rev., D83, 044026 (2011), arXiv:1004.3558.
  * (12) Richard Brito, Vitor Cardoso, and Paolo Pani, Class. Quant. Grav., 32(13), 134001 (2015), arXiv:1411.0686.
  * (13) Richard Brito, Shrobana Ghosh, Enrico Barausse, Emanuele Berti, Vitor Cardoso, Irina Dvorkin, Antoine Klein, and Paolo Pani, Phys. Rev., D96(6), 064050 (2017), arXiv:1706.06311.
  * (14) Asimina Arvanitaki, Masha Baryakhtar, and Xinlu Huang, Phys. Rev., D91(8), 084011 (2015), arXiv:1411.2263.
  * (15) Asimina Arvanitaki, Masha Baryakhtar, Savas Dimopoulos, Sergei Dubovsky, and Robert Lasenby, Phys. Rev. D, 95(4), 043001 (2017), arXiv:1604.03958.
  * (16) Hirotaka Yoshino and Hideo Kodama, PTEP, 2014, 043E02 (2014), arXiv:1312.2326.
  * (17) Sunil Choudhary, Nicolas Sanchis-Gual, Anshu Gupta, Juan Carlos Degollado, Sukanta Bose, and José A. Font, Phys. Rev. D, 103(4), 044032 (2021), arXiv:2010.00935.
  * (18) Adrian Ka-Wai Chung, Joseph Gais, Mark Ho-Yeuk Cheung, and Tjonnie G. F. Li, Phys. Rev. D, 104(8), 084028 (2021), arXiv:2107.05492.
  * (19) Valerio De Luca and Paolo Pani, JCAP, 08, 032 (2021), arXiv:2106.14428.
  * (20) Daniel Baumann, Horng Sheng Chia, and Rafael A. Porto, Phys. Rev. D, 99(4), 044001 (2019), arXiv:1804.03208.
  * (21) Tao Liu and Kun-Feng Lyu (7 2021), arXiv:2107.09971.
  * (22) Leo Tsukada, Thomas Callister, Andrew Matas, and Patrick Meyers, Phys. Rev. D, 99(10), 103015 (2019), arXiv:1812.09622.
  * (23) Ken K. Y. Ng, Salvatore Vitale, Otto A. Hannuksela, and Tjonnie G. F. Li (11 2020), arXiv:2011.06010.
  * (24) R. Abbott et al. (11 2021), arXiv:2111.15507.
  * (25) B. P. Abbott et al., Phys. Rev. Lett., 116(6), 061102 (2016), arXiv:1602.03837.
  * (26) B. P. Abbott et al., Phys. Rev. D, 93(12), 122003 (2016), arXiv:1602.03839.
  * (27) R. Abbott et al. (11 2021), arXiv:2111.03606.
  * (28) Vitor Cardoso, Francisco Duque, and Taishi Ikeda, Phys. Rev. D, 101(6), 064054 (2020), arXiv:2001.01729.
  * (29) Emanuele Berti, Richard Brito, Caio F. B. Macedo, Guilherme Raposo, and Joao Luis Rosa, Phys. Rev. D, 99(10), 104039 (2019), arXiv:1904.03131.
  * (30) Takuya Takahashi and Takahiro Tanaka, JCAP, 10, 031 (2021), arXiv:2106.08836.
  * (31) Masha Baryakhtar, Marios Galanis, Robert Lasenby, and Olivier Simon, Phys. Rev. D, 103(9), 095019 (2021), arXiv:2011.11646.
  * (32) Hidetoshi Omiya, Takuya Takahashi, and Takahiro Tanaka, PTEP, 2021(4), 043E02 (2021), arXiv:2012.03473.
  * (33) Daniel Baumann, Horng Sheng Chia, John Stout, and Lotte ter Haar, JCAP, 12, 006 (2019), arXiv:1908.10370.
  * (34) P. C. Peters and J. Mathews, Phys. Rev., 131, 435–440 (Jul 1963).
  * (35) Luc Blanchet, Living Rev. Rel., 17, 2 (2014), arXiv:1310.1528.
  * (36) L. D. LANDAU, Z. Sowjetunion, 2, 46–51 (1932).
  * (37) Clarence Zener, Proc. Roy. Soc. Lond. A, 137, 696–702 (1932).
  * (38) Jie Liu, Libin Fu, Bi-Yiao Ou, Shi-Gang Chen, Dae-Il Choi, Biao Wu, and Qian Niu, Phys. Rev. A, 66, 023404 (Aug 2002).
  * (39) NIST Digital Library of Mathematical Functions,  https://dlmf.nist.gov/ ().