# Mass transfer and boson cloud depletion in a binary black hole system

Yao Guo  Wenjie Zhong  Yiqiu Ma<EMAIL_ADDRESS>Daiqin Su
<EMAIL_ADDRESS>MOE Key Laboratory of Fundamental Physical Quantities
Measurement, Hubei Key Laboratory of Gravitation and Quantum Physics, PGMF,
Institute for Quantum Science and Engineering, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China

###### Abstract

Ultralight boson is one of the potential candidates for dark matter. If
exists, it can be generated by a rapidly rotating black hole via
superradiance, extracting the energy and angular momentum of the black hole
and forming a boson cloud. The boson cloud can be affected by the presence of
a companion star, generating fruitful dynamical effects and producing
characteristic gravitational wave signals. We study the dynamics of the boson
cloud in a binary black hole system, in particular, we develop a framework to
study the mass transfer between two black holes. It is found that bosons
occupying the growing modes of the central black hole can jump to the decaying
modes of the companion black hole, resulting in cloud depletion. This
mechanism of cloud depletion is different from that induced by the resonant
perturbation from the companion.

##  I Introduction

The detection of gravitational waves by current ground-based gravitational
wave detectors Abbott et al. (2016) (LIGO, Virgo, etc.) opens up a new avenue
to explore the astrophysical processes that involve strong gravitational
field. Future space-borne gravitational wave detectors (LISA, TQ Luo et al.
(2016), etc.) extend the detection frequency range and allow for exploration
of more fruitful physics, e.g., supermassive black holes. One of the
interesting target sources for gravitational wave detectors is the ultralight
boson, including axion and pseudoscalar axion-like particles. Axion provides a
solution to the strong CP problem Peccei and Quinn (1977); Weinberg (1978);
Wilczek (1978) and has the potential to demystifying the baryon asymmetry of
the universe, while more general axion-like particles are predicted by
symmetry breaking in string theory Svrcek and Witten (2006). These ultralight
bosons are also potential candidates for dark matter Preskill et al. (1983);
Abbott and Sikivie (1983); Dine and Fischler (1983); Li et al. (2014); Hui et
al. (2017). The ultralight bosons, if exists, can be produced by a rapidly
rotating black hole (BH) through superradiance instabilities Zel’Dovich (1971,
1972); Misner (1972); Starobinskiǐ (1973); Starobinskiǐ and Churilov (1974);
East and Pretorius (2017); East (2017); Brito et al. (2020); Cromb et al.
(2020), carrying away the mass and angular momentum of black hole and forming
a boson cloud/condensate Arvanitaki and Dubovsky (2011).

Since the presence of ultralight boson reduces the mass and spin of the black
hole, the precision measurement of black hole mass and spin via gravitational
wave detectors provides a powerful constraint on the properties of bosons
Brito et al. (2017); Arvanitaki et al. (2015, 2010). In addition, the boson
cloud can radiate continuous gravitational waves due to its asymmetric
distribution and self-annihilation Arvanitaki and Dubovsky (2011), and thus
can be detected by future gravitational wave detectors.

If a companion star is present, the boson cloud around a black hole is
distorted by a time-dependent tidal field. Under certain conditions, the time-
dependent tidal field can induce transitions between the boson’s various
energy levels, in particular, between the growing modes and decaying modes
Baumann et al. (2019). This would generate fruitful dynamical effects. The
transition of boson can transfer energy and angular momentum from the boson
cloud to the companion star, modifying its orbital evolution. In the extreme
case, the energy loss of the companion star due to the emission of
gravitational waves is balanced by the energy gain from the boson cloud, thus
forming the so called “floating orbit” Zhang and Yang (2019). When the boson
occupies the decaying mode, it may decay into the black hole, resulting in the
depletion of the boson cloud Baumann et al. (2019); Berti et al. (2019). In
some cases, the boson cloud could be completely cleaned up by the companion
star. The bosons that absorbed by the black hole increase the mass of the
black hole while reduce its spin. The reduction of the black hole spin can
turn some of the growing modes into decaying modes and then accelerate the
cloud depletion Takahashi et al. (2023). The companion can also induce
transitions between bound and unbound orbits of the boson, thus ionize it
Baumann et al. (2022). These fruitful dynamical processes of the boson cloud
produce characteristic gravitational wave signals and consequently can be
detected by current and future gravitational wave detectors Siemonsen et al.
(2023). This provides a way to infer the presence of the ultralight boson and
constrain its various properties.

When the black hole and its companion is sufficiently close, the boson can
escape to the companion, resulting in a redistribution of the mass of the
boson cloud. A gravitational “molecule”, an analog to the hydrogen molecule,
can be formed under certain conditions Wong (2020); Ikeda et al. (2021). If
the companion is also a rotating black hole, then the escaped boson may occupy
the decaying modes of the companion black hole and may decay into it. If
happens, then this provides a new channel for boson cloud depletion. In this
work, we develop a framework to study the mass transfer between two black
holes in a binary system, assuming that the two black holes have the same mass
and spin, and their spin orientation is parallel. We make analogy with the
hydrogen molecule ion and calculate the wave functions of the molecular orbits
and the corresponding energy eigenvalues using the variational method. We then
obtain the probability for the boson jumping to the decaying mode of the
companion black hole, and show that this leads to a strong cloud depletion
which almost completely cleans up the boson cloud.

This paper is organized as follows. In Sec. II, we briefly review the boson
cloud around a black hole and its depletion into the black hole due to the
time-dependent tidal field produced by the companion. In Sec. III, by making
analogy with hydrogen molecular ion, we use the variational method to derive
the wave functions and energy eigenvalues of the boson in the binary black
hole system. In Sec. IV, we use the adiabatic approximation to study the time
evolution of the wave function of the boson, evaluate the the probability of
jumping to the decaying mode of the companion black hole, and calculate the
time evolution of the cloud mass due to the decay into the companion black
hole. We finally conclude in Sec. V.

##  II Boson cloud around a black hole

###  II.1 Gravitational atom

A rapidly rotating black hole can radiate ultralight bosons via superradiance
instabilities. These bosons condensate in some of their orbits, forming a
boson cloud around the black hole. When the mass of the boson is small, the
size of the cloud can be much larger than the gravitational radius of the
rotating black hole. In this limit, the Newtonian approximation is sufficient
to describe the dynamics of the boson cloud. The orbits of the boson is
determined by a Schrödinger-like equation, similar to that of an electron in a
hydrogen atom. The eigenstates of the boson are denoted as
$\ket{\varphi_{n\ell m}}$ or $\varphi_{n\ell m}$, and the eigenfrequencies are
given by Baumann et al. (2019)

$\displaystyle\omega_{n\ell
m}\approx\mu\bigg{(}1-\frac{\alpha^{2}}{2n^{2}}\bigg{)},$ (1)

where $\mu$ is the mass of the boson, $\alpha\equiv GM\mu/\hbar c$ is the
dimensionless “fine-structure constant”. The radial profile of the wave
function peaks at

$\displaystyle
r_{c,n}\approx\bigg{(}\frac{n^{2}}{\alpha^{2}}\bigg{)}r_{g}=n^{2}r_{b},$ (2)

where $r_{g}$ is defined as the gravitational radius of the black hole,
$r_{g}\equiv GM/c^{2}$, and $r_{b}=r_{g}/\alpha^{2}$ is defined as the Bohr
radius.

However, there is a crucial difference between the electron in the hydrogen
atom and the boson around a black hole: the orbits of the electron are stable
while the orbits of the boson are not stable due to the presence of the black
hole horizon. This is characterized by the imaginary part of the
eigenfrequency, $\omega_{n\ell m}\rightarrow\omega_{n\ell m}+i\Gamma_{n\ell
m}$. In the limit $\alpha\ll 1$, $\Gamma_{n\ell m}$ can be approximated as
Detweiler (1980)

$\displaystyle\Gamma_{n\ell m}=\frac{2r_{+}}{M}C_{n\ell
m}(\alpha)(m\Omega_{\rm H}-\omega_{n\ell m})\alpha^{4\ell+5},$ (3)

where $C_{n\ell m}(\alpha)$ is positive and given by

$\displaystyle C_{n\ell m}(\alpha)$ $\displaystyle=$
$\displaystyle\frac{2^{4\ell+1}(n+\ell)!}{n^{2\ell+4}(n-\ell-1)!}\bigg{[}\frac{\ell!}{(2\ell)!(2\ell+1)!}\bigg{]}^{2}$
(4)
$\displaystyle\times\prod_{j=1}^{\ell}\bigg{[}j^{2}(1-\tilde{a}^{2})+(\tilde{a}m-2\tilde{r}_{+}\alpha)^{2}\bigg{]},$

with $\tilde{a}=a/M$ and $\tilde{r}_{+}=r_{+}/M$. Here $r_{+}$ is the size of
the event horizon, $a$ is the spin and $\Omega_{\rm H}$ is the angular
velocity of the rotating black hole. The orbits with positive $\Gamma_{n\ell
m}$ are growing modes, for which the number of boson grows exponentially;
while the orbits with negative $\Gamma_{n\ell m}$ are decaying modes, for
which the number of boson decays exponentially. Starting from a rapidly
rotating black hole, bosons are radiated due to the superradiance
instabilities and then occupy the growing modes. The radiated bosons carry
away angular momentum, slowing down the rotation of the black hole. At an
equilibrium point the black hole rotates slow enough such that no bosons can
be further radiated, and a quasi-stationary boson cloud forms around the black
hole.

###  II.2 Hyperfine and Bohr resonance

When the black hole with a boson cloud is part of a binary system, the
gravitational field of the companion star distorts the cloud, resulting in
transition of bosons between growing modes and decaying modes. The bosons that
jump to decaying modes can return to the black hole, reducing the total mass
of the cloud and transferring angular momentum to the companion star.
Therefore, the existence of boson cloud would affect the orbital evolution of
the companion star and the gravitational waveforms from the binary system.

The companion star induces a time-dependent perturbation to the Kerr metric,
which then introduces a time-dependent shift of gravitational potential to the
Schrödinger equation that dominates the dynamics of bosons. Under certain
conditions, the time-dependent perturbation can induce resonant transitions
between growing modes and decaying modes Baumann et al. (2019). There are two
types of resonances, the hyperfine (or Rabi) resonance and the Bohr resonance.
The resonance occurs at a specific orbital separation. For the hyperfine
resonance, the orbital separation is given by

$\displaystyle R_{\rm
res}^{(h)}=144^{1/3}\alpha^{-4}(1+q)^{1/3}\tilde{a}^{-2/3}r_{g},$ (5)

while for the Bohr resonance, the orbital separation is given by

$\displaystyle R_{\rm
res}^{(b)}=\bigg{(}\frac{144}{5}\bigg{)}^{2/3}\alpha^{-2}(1+q)^{1/3}r_{g}.$
(6)

In general, the orbital separation of the hyperfine resonance is much larger
than that of the Bohr resonance. This is because the hyperfine energy gap is
much narrower than the Bohr energy gap. This implies that the companion star
has to be closer to the black hole in order to excite the Bohr resonance.

##  III Boson orbits in a binary black hole system

We are concerned with the dynamics of the boson cloud when a black hole and
its companion are sufficiently close so that mass transfer between them cannot
be ignored. The process of mass transfer is generally very complicated and
numerical simulation is needed to fully uncover the cloud dynamics. It has
been shown by using effective field theory techniques Wong (2020) and
numerical calculation Ikeda et al. (2021) that a gravitational molecule can
form in a binary black hole system. Here, we use a simple model to describe
the mass transfer and the time evolution of a BH-cloud-companion system. The
model is based on the analogy between the hydrogen molecule ion $H_{2}^{+}$
and the BH-cloud-companion system, and it can capture the main characteristics
of mass transfer and cloud depletion.

For simplicity, we assume (1) the companion is also a black hole, therefore
forming a binary black hole system; (2) two black holes have the same mass and
spin; (3) their spin orientation is parallel and perpendicular to the orbital
plane. The schematic of the configuration is shown in Fig. 1. These
assumptions allow us to focus on the effect of mass transfer, and they can be
relaxed to incorporate more complicated effects. Under these assumptions, the
orbits of the boson in the BH-cloud-BH system is analogous to that of the
electron in the hydrogen molecule ion $H_{2}^{+}$, except that the two black
holes rotate with each other and that the existence of a horizon results in
boson absorption.

Figure 1:  Configuration of the BH-cloud-BH system and the coordinate system
used to describe the boson. The origin of the coordinate system is located at
the central black hole, and the $z$-axis is parallel to the spin of the
central black hole and the $x$-axis is pointing towards the companion black
hole. Spherical coordinates $(r_{1},\theta_{1},\varphi_{1})$ are used to
represent the position of the boson relative to the central black hole, and
$(r_{2},\theta_{2},\varphi_{2})$ are used to denote the position of the boson
relative to the companion black hole. Here $a_{1}$ and $a_{2}$ represent the
spin of the central and companion black holes, respectively, and $R$ denote
the orbital separation.

###  III.1 Orbits of hydrogen molecule ion

The hydrogen molecule ion $H_{2}^{+}$ consists of two protons and a single
electron. The electron moves in the potential produced by the two protons with
a fixed distance. The potential is time independent, so the electron wave
functions can be derived by solving the stationary Schrödinger equation.
Approximate energy eigenvalues and eigenfunctions can be solved using the
variational method. The energy eigenvalues are the stationary points of the
expectation value of the Hamiltonian. To approximate the energy eigenvalue,
one starts from a trial wave function and then vary the parameters in the
trial wave function to find the stationary point of the Hamiltonian.

For the hydrogen molecule ion, the initial trial wave functions can be
selected by using the symmetric properties of the system. For example, one can
choose a linear superposition of the two ground-state wave functions as a
trial function to find the ground state of the hydrogen molecule ion. The
excited states can also be found in a similar way. In this paper, we are
mainly interested in the excited states with $n=2$ because these states are
relevant to the states that dominate the evolution of the boson cloud. It can
be shown that there exists two molecular $\sigma$ orbits obtained from the
atomic orbit $2p_{x}$ Cohen-Tannoudji et al. (1986) (see Appendix B for
details),

$\displaystyle\ket{\sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}N_{1}}\big{(}\ket{\varphi_{2p_{x}}^{1}}-\ket{\varphi_{2p_{x}}^{2}}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2N_{1}}\big{(}\ket{\varphi_{2,1,1}^{1}}+\ket{\varphi_{2,1,-1}^{1}}-\ket{\varphi_{2,1,1}^{2}}-\ket{\varphi_{2,1,-1}^{2}}\big{)},$
$\displaystyle\ket{\sigma^{*}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}N_{2}}\big{(}\ket{\varphi_{2p_{x}}^{1}}+\ket{\varphi_{2p_{x}}^{2}}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2N_{2}}\big{(}\ket{\varphi_{2,1,1}^{1}}+\ket{\varphi_{2,1,-1}^{1}}+\ket{\varphi_{2,1,1}^{2}}+\ket{\varphi_{2,1,-1}^{2}}\big{)};$

and another two molecular $\pi$ orbits obtained from the atomic orbit
$2p_{y}$,

$\displaystyle\ket{\pi}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}N_{3}}\big{(}\ket{\varphi_{2p_{y}}^{1}}+\ket{\varphi_{2p_{y}}^{2}}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2iN_{3}}\big{(}\ket{\varphi_{2,1,1}^{1}}-\ket{\varphi_{2,1,-1}^{1}}+\ket{\varphi_{2,1,1}^{2}}-\ket{\varphi_{2,1,-1}^{2}}\big{)},$
$\displaystyle\ket{\pi^{*}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}N_{4}}\big{(}\ket{\varphi_{2p_{y}}^{1}}-\ket{\varphi_{2p_{y}}^{2}}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2iN_{4}}\big{(}\ket{\varphi_{2,1,1}^{1}}-\ket{\varphi_{2,1,-1}^{1}}-\ket{\varphi_{2,1,1}^{2}}+\ket{\varphi_{2,1,-1}^{2}}\big{)},$

where $N_{1},N_{2},N_{3}$ and $N_{4}$ are introduced to normalize these
states,

$\displaystyle N_{1}$ $\displaystyle=$
$\displaystyle\sqrt{1-\langle\varphi_{2p_{x}}^{1}\ket{\varphi_{2p_{x}}^{2}}},~{}~{}~{}~{}N_{2}=\sqrt{1+\langle\varphi_{2p_{x}}^{1}\ket{\varphi_{2p_{x}}^{2}}},$
$\displaystyle N_{3}$ $\displaystyle=$
$\displaystyle\sqrt{1+\langle\varphi_{2p_{y}}^{1}\ket{\varphi_{2p_{y}}^{2}}},~{}~{}~{}~{}N_{4}=\sqrt{1-\langle\varphi_{2p_{y}}^{1}\ket{\varphi_{2p_{y}}^{2}}}.$

It is evident that $N_{i}$ depends on the distance between two protons. In the
following, we will use the notation $i=1,2,3,4$ to denote the four molecular
orbits $\sigma,\sigma^{*},\pi$ and $\pi^{*}$, respectively.

###  III.2 Born-Oppenheimer approximation

We now turn to the problem of solving the orbits of the boson in a binary
black hole system with the two black holes have the same mass and spin. This
system is analogous to the hydrogen molecule ion $H_{2}^{+}$ except that the
two black holes rotate with each other and their separation shrinks due to the
emission of gravitational waves. The bosons experience a time-dependent field
rather than a static one. An important question is how to take into account
the effects of rotation and orbital shrinking.

In Ref. Baumann et al. (2019), the rotation of the companion star relative to
the central black hole produces a time-varying tidal field that perturbs the
central black hole, which then induces hyperfine and Bohr mixing of growing
and decaying modes, resulting in boson cloud depletion. Whilst we are
concerned with the transfer of bosons between the two black holes, the
framework developed in Ref. Baumann et al. (2019) does not apply here. As a
first approximation, we neglect the effects of rotation and orbital shrinking,
and treat the two black holes as stationary with a fixed orbital separation.
This is known as the Born-Oppenheimer approximation, in which the binary black
hole system with a boson cloud is analogous to the hydrogen molecule ion at
any given time. Before performing the detailed calculation, we compare various
timescales to validate the use of Born-Oppenheimer approximation.

The velocity of a boson can be estimated using the Virial theorem and the
uncertainty principle. The boson’s velocity $v_{a}$ satisfies

$\displaystyle\frac{1}{2}\mu
v_{a}^{2}\sim\frac{GM\mu}{r_{c}},~{}~{}~{}~{}~{}~{}\mu
v_{a}\sim\frac{\hbar}{r_{c}},$ (10)

where $r_{c}$ is the typical length scale of the radial profile of the boson
cloud. This implies

$\displaystyle v_{a}\sim\frac{GM\mu}{\hbar c}c\sim\alpha c\sim\alpha.$ (11)

The relaxation timescale for the boson can be approximated by

$\tau_{r}\sim\frac{R}{v_{a}}\sim\frac{R}{\alpha}.$ (12)

The time $\tau_{r}$ characterizes the timescale for the boson moving from one
black hole to the other. The period of the binary black holes is simply given
by

$T=\frac{2\pi}{\Omega}=2\pi\sqrt{\frac{R^{3}}{M(1+q)}}.$ (13)

Therefore the ratio between $\tau_{r}$ and $T$ is

$\displaystyle\frac{\tau_{r}}{T}\sim\frac{\sqrt{1+q}}{2\pi}\sqrt{\frac{M}{\alpha^{2}R}}\sim\frac{\sqrt{1+q}}{2\pi}\sqrt{\frac{r_{b}}{R}}.$
(14)

For $q=1$ and $R=32\,r_{b}$, $\tau_{r}/T\sim 1/8\pi\sim 0.04$. As we will see
in the following discussion, the mass transfer occurs when $R\gtrsim
50\,r_{b}$. This shows that for an intermediate orbital separation, the
relaxation time of the boson is much shorter than the period of the binary.
Therefore, the two black holes can be treated as quasi-static when considering
the transfer of bosons and the Born-Oppenheimer approximation applies.

The orbit of the binary black holes shrinks due to the emission of
gravitational waves Peters (1964),

$\displaystyle
R(t)=\bigg{[}\frac{M(1+q)}{\Omega_{0}^{2}}\bigg{]}^{1/3}\bigg{(}-\frac{t}{\tau_{0}}\bigg{)}^{1/4},$
(15)

where we set $t=0$ as the moment of merger, and $\tau_{0}$ is the time to
merger for an initial orbital frequency $\Omega_{0}$. The initial time
$\tau_{0}$ and initial orbital frequency $\Omega_{0}$ are related via

$\displaystyle\frac{\tau_{0}}{M(1+q)}=\frac{5}{256}\frac{(1+q)^{2}}{q}\bigg{[}\frac{1}{M(1+q)\Omega_{0}}\bigg{]}^{8/3}.$
(16)

The characteristic timescale for coalescence can be estimated via

$\tau_{p}=\bigg{|}\bigg{(}\frac{dR}{dt}\bigg{)}^{-1}\bigg{|}R\approx\frac{5}{64}\frac{R^{4}}{M^{3}}\frac{1}{q(1+q)}.$
(17)

This timescale is evidently much longer than the period of the binary $T$, and
it is thus also much longer than the boson relaxation timescale $\tau_{r}$.

$\displaystyle\frac{\tau_{r}}{\tau_{p}}\sim\frac{32}{5}\alpha^{5}\bigg{(}\frac{r_{b}}{R}\bigg{)}^{3}q(1+q).$
(18)

When the orbital separation of the two black holes is sufficiently large and
$q$ is not so large, these two ratios are both much smaller than one, so the
Born-Oppenheimer approximation is sound.

###  III.3 Eigenfunctions and energy eigenvalues

To find the exact eigenfunctions and eigenfrequencies, one in principle needs
to solve the Klein-Gordon equation of the boson in the background spacetime of
the binary black hole system, which is quite a challenging task. When the mass
of the boson is small, the boson cloud is far away from the black hole and
Newtonian approximation can be used to derive the eigenfunctions Baumann et
al. (2019). To the order of $1/r$, the wave function of the boson around a
single rotating black hole satisfies

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi(t,\boldsymbol{r})=\bigg{(}-\frac{1}{2\mu}\nabla^{2}-\frac{\alpha}{r}\bigg{)}\psi(t,\boldsymbol{r}),$
(19)

which is exactly in the same form as the Schrödinger equation for the electron
in a hydrogen atom. When the two black holes are not so close to each other
and the boson lingers around a regime far away from both black holes, then the
Newtonian approximation applies. To the order of $1/r$, the wave function of
the boson in the binary black hole system satisfies

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi(t,\boldsymbol{r})=\bigg{(}-\frac{1}{2\mu}\nabla^{2}-\frac{\alpha}{r_{1}}-\frac{\alpha}{r_{2}}\bigg{)}\psi(t,\boldsymbol{r}),$
(20)

where $r_{1}$ is the distance between the boson and the central black hole,
and $r_{2}$ is the distance between the boson and the companion black hole.
Equation (20) is exactly in the same form as the Schrödinger equation for the
electron in the hydrogen molecule ion, without including the interaction
between two protons.

In the Born-Oppenheimer approximation, the eigenfunctions of the boson at any
given time can be derived straightforwardly. They are exactly in the same form
as those given by Eqs. (III.1) and (III.1), with $\ket{\varphi_{n\ell m}}$ the
eigenfunctions of the boson in a single isolated rotating black hole.
Furthermore, the normalization constants $N_{i}$ depend on the configuration
of the BH-cloud-BH system. For the case that we consider here, the
normalization constants have no analytic expressions and are needed to be
evaluated numerically.

Once the eigenfunctions are known, the energy eigenvalues can be calculated
straightforwardly, which are simply the expectation values of the Hamiltonian.
From Eq. (20) it is evident that the Hamiltonian of the boson in the Newtonian
limit is

$\displaystyle\hat{H}=-\frac{1}{2\mu}\nabla^{2}-\frac{\alpha}{r_{1}}-\frac{\alpha}{r_{2}}.$
(21)

To simplify the calculation, we divide the Hamiltonian $\hat{H}$ into two
parts: the Hamiltonian of the boson in the central black hole and the
potential produced by the companion black hole, namely,
$\hat{H}=\hat{H}_{1}-\frac{\alpha}{r_{2}}$ with
$\hat{H}_{1}=-\frac{1}{2\mu}\nabla^{2}-\frac{\alpha}{r_{1}}$. Then the
eigenfrequencies are given by

$\displaystyle\omega_{1}$ $\displaystyle=$
$\displaystyle\bra{\sigma}\hat{H}\ket{\sigma}=\frac{1}{N_{1}^{2}}\bigg{[}-\frac{1}{8}\mu\alpha^{2}\big{(}1-\langle\varphi_{2p_{x}}^{1}\ket{\varphi_{2p_{x}}^{2}}\big{)}-\bra{\varphi_{2p_{x}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{x}}^{1}}+\bra{\varphi_{2p_{x}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{x}}^{2}}\bigg{]},$
$\displaystyle\omega_{2}$ $\displaystyle=$
$\displaystyle\bra{\sigma^{*}}\hat{H}\ket{\sigma^{*}}=\frac{1}{N_{2}^{2}}\bigg{[}-\frac{1}{8}\mu\alpha^{2}\big{(}1+\langle\varphi_{2p_{x}}^{1}\ket{\varphi_{2p_{x}}^{2}}\big{)}-\bra{\varphi_{2p_{x}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{x}}^{1}}-\bra{\varphi_{2p_{x}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{x}}^{2}}\bigg{]},$
$\displaystyle\omega_{3}$ $\displaystyle=$
$\displaystyle\bra{\pi}\hat{H}\ket{\pi}=\frac{1}{N_{3}^{2}}\bigg{[}-\frac{1}{8}\mu\alpha^{2}\big{(}1+\langle\varphi_{2p_{y}}^{1}\ket{\varphi_{2p_{y}}^{2}}\big{)}-\bra{\varphi_{2p_{y}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{y}}^{1}}-\bra{\varphi_{2p_{y}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{y}}^{2}}\bigg{]},$
$\displaystyle\omega_{4}$ $\displaystyle=$
$\displaystyle\bra{\pi^{*}}\hat{H}\ket{\pi^{*}}=\frac{1}{N_{4}^{2}}\bigg{[}-\frac{1}{8}\mu\alpha^{2}\big{(}1-\langle\varphi_{2p_{y}}^{1}\ket{\varphi_{2p_{y}}^{2}}\big{)}-\bra{\varphi_{2p_{y}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{y}}^{1}}+\bra{\varphi_{2p_{y}}^{1}}\frac{\alpha}{r_{2}}\ket{\varphi_{2p_{y}}^{2}}\bigg{]}.$
(22)

For the case that we consider here, there exists no analytic expressions for
$\omega_{i}$, so they have to be evaluated numerically. Note that in deriving
the eigenfunctions and eigenfrequencies, we have ignored the higher order
corrections to the Hamiltonian. This results in degeneracy of all four energy
levels when the two black holes are infinitely far away.

##  IV Cloud depletion

Assume that the central black hole has a boson cloud surrounding it while the
companion black hole does not. When the two black holes are very far away from
each other, the bosons move around the central black hole and cannot escape to
the companion. When the two black holes become closer to each other, the boson
orbits that belong to the central black hole and that belong to the companion
have more overlap, forming boson orbits that are analogous to the molecular
orbits of the hydrogen molecule ion. Therefore, the boson around the central
black hole may jump to the companion. This can have two consequences: first,
the bosons redistribute in the binary black hole system, changing its
quadruple and thus modifying the gravitational waveform; second, the bosons
may jump to the decaying modes of the companion, and therefore could be
absorbed by the companion black hole. To estimate the importance of the above
two consequences, one needs to model the process of boson transfer from the
central black hole to its companion.

The binary black hole system emits gravitational waves so its orbit shrinks,
namely, the orbital separation between the two black holes shrinks due to the
radiation of energy. The orbital separation, denoted as $R(t)$, is thus time
dependent. This implies that the eigenfunctions and eigenfrequencies are also
time dependent because they are solved at a given time assuming the orbital
separation is fixed. These solutions are meaningful only in the case where the
time scale of the orbital shrinking is much larger than that of the boson
relaxation. This is exactly the case as we have discussed in Sec. III.2. In
order to know the evolution of the boson cloud in the binary black hole
system, one needs to solve a time-dependent Schrödinger equation, in which the
Hamiltonian varies slowly. This sort of problem can be solved using the
adiabatic approximation if the energy levels are not degenerate and the energy
gap is sufficiently large.

###  IV.1 Adiabatic approximation

We numerically calculate the eigenfrequencies given by Eq. (III.3) and plot
them in Figs. 2 and 3. It can be seen from Fig. 2 that the four energy levels
split and the energy gap is larger when the two black holes are closer. The
gap decreases as the orbital separation increases, which is consistent with
the fact that the four energy levels are degenerate when $R\rightarrow\infty$.
Figure 3 shows that the energy gaps between orbits $\sigma$ and $\sigma^{*}$,
and orbits $\pi$ and $\pi^{*}$ decrease much faster than that between
$\sigma$’s and $\pi$’s orbits. Because the energy gaps between orbits $\sigma$
and $\sigma^{*}$, and orbits $\pi$ and $\pi^{*}$ are very close, a question
arises as to whether the adiabatic approximation is still valid.

Figure 2: Energy eigenvalues of the boson molecular orbits for small orbital
separation. Here we choose $q=1$ and $\alpha=0.1$ as an example. Figure 3:
Energy eigenvalues of the boson molecular orbits for large orbital separation.
Here we choose $q=1$ and $\alpha=0.1$ as an example.

In the adiabatic approximation, the eigenfunctions and eigenfrequencies vary
slowly with time, but the boson remains in the initial energy level during the
subsequent evolution, namely, no transition from one energy level to others.
In order to see whether the adiabatic approximation is valid for the problem
we consider here, we need to check whether the transition between orbits
$\sigma$ and $\sigma^{*}$, and orbits $\pi$ and $\pi^{*}$ is negligible.

Consider an arbitrary quantum state

$\displaystyle\ket{\psi(t)}=\sum_{n}c_{n}(t)\ket{\psi_{n}(t)},$ (23)

where $\ket{\psi_{n}(t)}$ and $c_{n}(t)$ are the time-dependent energy
eigenstates and their corresponding coefficients, respectively. In the
adiabatic approximation, we have $|c_{n}(t+dt)|=|c_{n}(t)|$, namely, there is
no transition between different energy levels. To estimate the accuracy of the
adiabatic approximation, one can check the amplitude of the time derivative of
the coefficient $c_{n}(t)$. It can be shown that (see Appendix C for details)

$\displaystyle\frac{dc_{k}}{dt}=-\sum_{n\neq
k}c_{n}\frac{1}{E_{n}-E_{k}}\bra{\psi_{k}}\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}-ic_{k}E_{k},$ (24)

where the summation characterizes the transition from the initial energy level
to other energy levels, and the last term represents the free evolution in the
initial energy level. Now consider the set of state $\\{\ket{\psi_{k}}\\}$ as
the four orbits $\ket{\sigma},\ket{\sigma^{*}},\ket{\pi}$ and $\ket{\pi^{*}}$.
Due to the symmetry of the system, it can be shown (see Appendix C for
details) that

$\displaystyle\bra{\psi_{k}}\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}=0,$ (25)

for $n\neq k$. This implies that the adiabatic approximation is a good
approximation to the time evolution of the boson.

Taking into account the emission of gravitational waves and the orbital
shrinking, the eigenfrequencies and eigenstates are all time dependent,
denoted as $\omega_{i}(t)$, $\ket{\sigma(t)},\ket{\sigma^{*}(t)},\ket{\pi(t)}$
and $\ket{\pi^{*}(t)}$, respectively. When the two black holes are infinitely
far away, namely, $t\rightarrow-\infty$ and $R\rightarrow\infty$, the
normalization constants defined in Eq. (III.1) are the same and equal to one.
This implies that the initial state of the boson can be written as

$\displaystyle\ket{\Phi(-\infty)}\equiv\ket{\varphi_{2,1,1}^{1}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}\ket{\sigma(-\infty)}+\ket{\sigma^{*}(-\infty)}+i\ket{\pi(-\infty)}+i\ket{\pi^{*}(-\infty)}\big{]}.$

According to the adiabatic approximation, the subsequent evolution of the
state is

$\displaystyle\ket{\Phi(t)}=\frac{1}{2}\bigg{[}e^{-i\int_{-\infty}^{t}\omega_{1}(\tau)\rm{d}\tau}\ket{\sigma(t)}+e^{-i\int_{-\infty}^{t}\omega_{2}(\tau)\rm{d}\tau}\ket{\sigma^{*}(t)}$
$\displaystyle+i\,e^{-i\int_{-\infty}^{t}\omega_{3}(\tau)\rm{d}\tau}\ket{\pi(t)}+i\,e^{-i\int_{-\infty}^{t}\omega_{4}(\tau)\rm{d}\tau}\ket{\pi^{*}(t)}\bigg{]},$
(27)

We now substitute Eqs. (III.1) and (III.1) into Eq. (IV.1) and obtain an
expansion of $\ket{\Phi(t)}$ in terms of the isolated orbits
$\ket{\varphi_{2,1,\pm 1}^{1}}$ and $\ket{\varphi_{2,1,\pm 1}^{2}}$. The
coefficient of the state $\ket{\varphi_{2,1,-1}^{2}}$ is given by

$\displaystyle
C_{-}(t)=-\frac{1}{4N_{1}}e^{-i\int_{-\infty}^{t}\omega_{1}(\tau)\rm{d}\tau}+\frac{1}{4N_{2}}e^{-i\int_{-\infty}^{t}\omega_{2}(\tau)\rm{d}\tau}$
$\displaystyle-\frac{1}{4N_{3}}e^{-i\int_{-\infty}^{t}\omega_{3}(\tau)\rm{d}\tau}+\frac{1}{4N_{4}}e^{-i\int_{-\infty}^{t}\omega_{4}(\tau)\rm{d}\tau}.$
(28)

The nonzero overlap between the isolated orbit $\ket{\varphi_{2,1,-1}^{2}}$
and $\ket{\Phi(t)}$ implies that bosons can jump to the orbits of the
companion black hole during their evolution. The modular square of $C_{-}(t)$,

$\displaystyle|C_{-}(t)|^{2}=\sum_{i}\frac{1}{16N_{i}^{2}}$
$\displaystyle+\sum_{i<j}(-1)^{i+j}\frac{1}{8N_{i}N_{j}}\cos\bigg{[}\int_{-\infty}^{t}\big{(}\omega_{i}(\tau)-\omega_{j}(\tau)\big{)}\rm{d}\tau\bigg{]},$

represents the probability of jumping to the isolated orbit
$\ket{\varphi_{2,1,-1}^{2}}$. If we only consider the effect of gravitational
wave emission and neglect other effects like back reaction of the boson cloud
onto the orbital evolution, then we can use the relation defined in Eq. (15)
to express the coefficient as a function of $R$, namely, $C_{-}(R)$; and the
time integration in Eq. (IV.1) can be replaced by an integration over the
orbital separation $R$.

Figure 4:  Probability for a boson jumping to the decaying mode
$\ket{\psi_{2,1,-1}^{2}}$ of the companion black hole. Here we choose $q=1$
and $\alpha=0.1$. The blue curve includes both the rapidly oscillating and
slowly oscillating terms, while the orange curve includes only the slowly
oscillating terms.

When $R$ is sufficiently large, the energy levels of orbits $\sigma$ and
$\sigma^{*}$ are almost degenerate ($\omega_{1}\approx\omega_{2}$), as well as
that of orbits $\pi$ and $\pi^{*}$ ($\omega_{3}\approx\omega_{4}$), as can be
seen from Fig. 3. In addition, the normalization factors satisfy
$N_{1}(R)\approx N_{2}(R)$ and $N_{3}(R)\approx N_{4}(R)$. This implies that
the coefficient $C_{-}(R)$ approaches to zero when the two black holes are
sufficiently far away from each other, which is consistent with the fact that
the boson cannot escape to a companion that is very far away. When the two
black holes are closer, the energy levels of orbits $\sigma$ and $\sigma^{*}$
become nondegenerate ($\omega_{1}\neq\omega_{2}$), as well as those of orbits
$\pi$ and $\pi^{*}$. The coefficient $C_{-}(R)$ gradually becomes nonzero,
indicating that the boson can transfer to the companion black hole. From Eq.
(IV.1) we can see that the probability $|C_{-}(R)|^{2}$ is the sum of six
oscillating terms having frequencies $|\omega_{i}-\omega_{j}|$, with
$i,j=\\{\sigma,\sigma^{*},\pi,\pi^{*}\\}$. The frequencies
$|\omega_{i}-\omega_{j}|$ with $i\in\\{\sigma,\sigma^{*}\\}$ and
$j\in\\{\pi,\pi^{*}\\}$, or $j\in\\{\sigma,\sigma^{*}\\}$ and
$i\in\\{\pi,\pi^{*}\\}$, are much higher than that when
$i,j\in\\{\sigma,\sigma^{*}\\}$ or $i,j\in\\{\pi,\pi^{*}\\}$.

Figure 4 shows the occupation probability $|C_{-}(R)|^{2}$ of the orbit
$\ket{\psi_{2,1,-1}^{2}}$, which is one of the decaying modes of the companion
black hole. In the numerical calculation, we use a sufficiently large distance
$R=80\,r_{b}$ as a replacement for the infinitely large distance ($t=-\infty$)
in the lower limit of the integration in Eq. (IV.1), which introduces a
negligible error. The blue curve takes into account the full expression of
$|C_{-}(R)|^{2}$ and includes highly oscillating terms. If we remove the
highly oscillating terms and only keep two terms with oscillating frequencies
$|\omega_{\pi}-\omega_{\pi^{*}}|$ and $|\omega_{\sigma}-\omega_{\sigma^{*}}|$,
we obtain the profile of the probability $|C_{-}(R)|^{2}$, which is the orange
curve shown in Fig. 4. The profile captures the change of the probability that
the boson occupies the decaying mode of the companion black hole during the
evolution of the binary black holes. When $R>63\,r_{b}$, the probability is
almost zero. It gradually increases as $R$ decreases and reaches its first
maximum around $R=55.5\,r_{b}$. The first maximum of the probability is about
$0.3$, showing that a significant amount of boson transfers to the companion
black hole and occupies the decaying mode $\ket{\psi_{2,1,-1}^{2}}$. As the
orbital separation $R$ further decreases, the profile of the probability
oscillates with a higher and higher frequency. Note that the orbital
separation corresponding to the first maximum is larger than the Loche limit,
which is about $10\,r_{b}$ for $q=1$.

The boson occupying the decaying mode may decay into the black hole, resulting
in the depletion of the boson cloud. The boson that transfers to the decaying
mode of the companion therefore may decay into the companion black hole. This
is another channel other than the hyperfine mixing and Bohr mixing that could
result in cloud depletion. Assuming that the back reaction is small and the
decay rate of the black hole is not affected by the presence of another black
hole, the time evolution of the cloud mass can be described as

$\displaystyle\frac{dM_{c}}{dt}=2\sum_{i=1}^{2}\Gamma_{2,1,-1}|C^{i}_{2,1,-1}|^{2}M_{c},$
(30)

where $M_{c}$ is the mass of the boson cloud and $\Gamma_{2,1,-1}$, defined by
Eq. (3), is the decay rate of the decaying mode $\ket{\psi_{2,1,-1}}$. Here
$|C^{1}_{2,1,-1}|^{2}$ is the occupation probability of the decaying mode
$\ket{\psi_{2,1,-1}^{1}}$ of the central black hole and $|C^{2}_{2,1,-1}|^{2}$
is the occupation probability of the decaying mode $\ket{\psi_{2,1,-1}^{2}}$
of the companion black hole, namely, $|C^{2}_{2,1,-1}|^{2}=|C_{-}(R)|^{2}$.
The summation takes into account the decay of the boson into both the central
and companion black holes.

To see the effect of mass transfer onto the cloud depletion, we calculate the
time evolution of the cloud mass under the assumption that the bosons only
decay into the companion black hole. Suppose the total mass of the cloud right
before the boson escapes to the companion, e.g., $R=80\,r_{b}$, is $M_{c,0}$.
We further assume that the spin of the two black holes are both
$\chi=\frac{4\alpha}{1+4\alpha^{2}}$, which is the critical spin for
$\ket{\varphi_{2,1,1}}$ to be saturated Baumann et al. (2019); and the
companion black hole has no bosons around it (Otherwise, we have to modify the
initial condition of the state, which in principle can also be handled in our
framework). Back reaction of the boson cloud onto the black holes are
neglected. The result is shown as the blue dashed curve in Fig. 5. We can see
that the boson cloud quickly decays around $R=60\,r_{b}$ and finally almost
all bosons are absorbed by the companion black hole. To be more specific, the
cloud remains unchanged when $R\gtrsim 67\,r_{b}$ and almost completely
disappears when $R\lesssim 57\,r_{b}$. This means all bosons decay into the
companion black hole when the probability $|C_{-}(R)|^{2}$ is approaching to
its first maximum. At first glance this seems impossible because the maximal
probability is about $0.3$, which is smaller than one. From the physical
perspective this is reasonable. Each boson moves between the central and
companion black holes. It may be absorbed by the companion black hole when it
occupies the orbit $\ket{\varphi_{2,1,-1}}$ of the companion. Though the
probability of being absorbed within one round trip is quite small, the boson
can travel many round trips when the orbital separation decreases from
$67\,r_{b}$ to $57\,r_{b}$. Therefore, with a very high probability, the boson
has been absorbed by the companion black hole when $R<57\,r_{b}$. This is also
the key difference from the depletion mechanism due to the hyperfine
resonance. It is true that at the resonance point, nearly all bosons jump to
the decaying mode, however, the bosons stay in the decaying mode only for a
short time, so the total mass that absorbed by the black hole may not be
large, as demonstrated in Ref. Baumann et al. (2019).

Figure 5: Evolution of the boson cloud mass for $\alpha=0.1$. The blue dashed
curve denotes the cloud depletion due to the mass transfer to the decaying
mode of the companion black hole, the orange curve denotes the cloud depletion
due to the hyperfine mixing of the boson in the central black hole and the
purple dashed-dotted curve denotes the total cloud depletion.

We can also include the cloud depletion due to the hyperfine mixing. Instead
of directly calculating the coefficients of the orbits $\ket{\varphi_{2,1,\pm
1}^{1}}$ from Eq. (IV.1), we use the method developed in Ref. Baumann et al.
(2019) to calculate the probability that the boson jumps to the decaying mode
of the central black hole due to the perturbation of the companion black hole.
Firstly, the effect of rotation is not included in Eq. (IV.1). Secondly,
according to the adiabatic theorem, the slow orbital shrinking cannot induce a
direct transition from the growing mode $\ket{\varphi_{2,1,1}^{1}}$ to the
decaying mode $\ket{\varphi_{2,1,-1}^{1}}$ of the central black hole. A boson
may jump to the decaying mode $\ket{\varphi_{2,1,-1}^{1}}$ by first moving to
the companion black and then back to the central black hole. However, this is
a second order effect and the probability is at the order of $|C_{-}(R)|^{4}$,
which therefore can be neglected.

For the parameters that we consider in this paper, the hyperfine resonance
occurs at the orbital separation $R\approx 660\,r_{b}$, which is much larger
than the orbital separation where the mass transfer occurs. We therefore
expect that the effect of hyperfine mixing is subdominant in the regime
$R<80\,r_{b}$, which is confirmed by Fig. 5. The orange curve represents the
cloud depletion only due to the hyperfine mixing, which decreases rather
slowly as compared to that due to the mass transfer to the companion black
hole. The reason is that the probability for a boson jumping to the decaying
mode of the central black hole is very small in this regime. Combining the
effects of the mass transfer and hyperfine mixing together, we obtain the
total cloud depletion, which is shown in Fig. 5. When $R\gtrsim 67\,r_{b}$,
the depletion is dominated by the hyperfine mixing; and when $R\lesssim
67\,r_{b}$, the depletion is dominated by the mass transfer to the companion
black hole and the decay into it. At about $R\sim 57\,r_{b}$, the cloud has
completely decayed into the two black holes, which is much earlier than the
prediction of hyperfine mixing alone.

In the above discussion we assume that a boson cloud exists when the orbital
separation is about $R=80\,r_{b}$. This could happen in several cases. The
boson cloud may form when the two black holes are close enough, in particular,
when the orbital separation is much smaller than the orbit separation when the
hyperfine resonance occurs. In this case the cloud depletes very slowly due to
the hyperfine mixing and a significant amount of bosons may still remain when
the mass transfer starts to dominate the depletion. The boson cloud may also
form when the two black holes are very far away, in particular, when the
orbital separation is larger than the orbit separation when the hyperfine
resonance occurs. There are two possibilities in this case. If the companion
black hole and the cloud are counter-rotating, then the hyperfine resonance
never occurs. The hyperfine and Bohr mixing are both weak so the cloud
depletes very slowly, as shown in Fig. 6, and therefore a substantial amount
of cloud may still remain. If the companion black hole and the cloud are co-
rotating, then the hyperfine resonance occurs. The hyperfine resonance may
cause a strong depletion of the cloud so that only a small amount of cloud is
left when the mass transfer starts to dominate the depletion. If this is the
case, then the effect of mass transfer does not play a dominant role in the
cloud depletion.

Figure 6: Evolution of the boson cloud mass for different $\alpha$ and initial
cloud mass $M_{c,0}=\alpha M$ Brito et al. (2017).

As an example, we consider the case where the companion black hole and the
cloud are co-rotating, and the cloud forms after the hyperfine resonance
occurs. The cloud has a life time $\tau_{c}$ that determined only by its own
gravitational radiation, which is given by Yoshino and Kodama (2014); Brito et
al. (2017); Baumann et al. (2019)

$\displaystyle\tau_{c}\sim
10^{7}\bigg{(}\frac{M}{3M_{\odot}}\bigg{)}\bigg{(}\frac{0.07}{\alpha}\bigg{)}^{15}~{}\text{years},$
(31)

where $M$ is the mass of the central black hole. We set the initial time of
the cloud evolution as the half lifetime of the cloud. The result is shown in
Fig. 6, where we have included the contribution from the hyperfine mixing,
Bohr mixing, mass transfer and gravitational radiation by the cloud itself.
For a fixed value of $\alpha$, the cloud depletion due to the hyperfine
mixing, Bohr mixing and gravitational radiation dominates at large orbital
separation, which proceeds very slowly. When the two black holes are close
enough and the mass transfer occurs, the cloud quickly depletes and almost all
remaining bosons are absorbed by the companion black hole. For a smaller value
of $\alpha$, the cloud depletes slower during the stage of hyperfine mixing
and the mass transfer occurs earlier, i.e., at a larger orbital separation.
This implies that a larger fraction of boson depletes due to the mass transfer
and is absorbed by the companion black hole.

##  V Conclusions

We develop a framework to study the transfer of bosons between two black holes
in a binary black hole system. The framework is formulated through the analogy
between the BH-cloud-BH system and the hydrogen molecule ion system in which
an electron moves in the potential generated by two protons. When two black
holes are sufficiently close, the bosons initially confined around the central
black hole can escape to the companion. In the language of quantum mechanics,
molecular orbits of the boson form and the boson moves back and forth between
two black holes. This results in cloud mass redistribution in the binary black
hole system. Furthermore, the boson escapes to the companion may occupy the
decaying mode and therefore may decay into the companion black hole. We find
that the boson cloud that exists right before the mass transfer occurs
completely disappears thgough this mechanism of cloud depletion. The cloud
depletion to the companion black hole may have important consequences to the
orbital evolution of the binary black hole system and the gravitational
waveform emitted from it.

Our framework is used to study the simplest model where the two black holes
have the same mass and spin, and their spin orientation is parallel. It can be
straightforwardly generalized to explore more realistic models, for example,
the two black holes may have unequal mass and spin, or their spin orientation
may be different. This would require a modification to the variational method
used to calculate the molecular orbits and the energy eigenvalues. Another
interesting case is that the companion may not be a black hole but a compact
star. Then there is only cloud mass redistribution but no cloud depletion due
to the decay into the companion. However, this could also have important
consequences to the evolution of the binary system and their gravitational
waveforms.

Acknowledgements: D. S. is supported by the Fundamental Research Funds for the
Central Universities, HUST (Grant No. 5003012068). Y. M. is supported by the
University start-up funding provided by Huazhong University of Science and
Technology.

##  Appendix A Choice of coordinates

To describe the orbits of the boson in the binary black hole system, we need
to set up an appropriate coordinate system, which is schematically shown in
Fig. 1. The origin of the coordinate system is located at the central black
hole, and the $z$-axis is parallel to the spin of the central black hole and
the $x$-axis is pointing towards the companion black hole. We use spherical
coordinates $(r_{1},\theta_{1},\varphi_{1})$ to represent the position of the
boson relative to the central black hole, and $(r_{2},\theta_{2},\varphi_{2})$
to denote the position of the boson relative to the companion black hole.
Whilst one set of coordinates is sufficient, we introduce the coordinates
$(r_{2},\theta_{2},\varphi_{2})$ only for convenience since the isolated wave
functions of the boson belonging to the companion black hole can be
conveniently expressed using the coordinates $(r_{2},\theta_{2},\varphi_{2})$.
The coordinates $(r_{2},\theta_{2},\varphi_{2})$ can be written in terms of
the coordinates $(r_{1},\theta_{1},\varphi_{1})$,

$\displaystyle r_{2}$ $\displaystyle=$
$\displaystyle\sqrt{r_{1}^{2}+R^{2}-2r_{1}R\sin\theta_{1}\cos\varphi_{1}}~{},$
$\displaystyle\cos\theta_{2}$ $\displaystyle=$
$\displaystyle\frac{r_{1}\cos\theta_{1}}{r_{2}},$
$\displaystyle\sin\theta_{2}$ $\displaystyle=$
$\displaystyle\sqrt{1-(\cos\theta_{2})^{2}}~{},$
$\displaystyle\cos\varphi_{2}$ $\displaystyle=$
$\displaystyle\frac{r_{1}\sin\theta_{1}\cos\varphi_{1}-R}{r_{2}\sin\theta_{2}},$
$\displaystyle\sin\varphi_{2}$ $\displaystyle=$
$\displaystyle\frac{r_{1}\sin\theta_{1}\sin\varphi_{1}}{r_{2}\sin\theta_{2}}.$
(32)

By using these relations we can carry out all the calculation using only one
set of coordinates, namely, the coordinates $(r_{1},\theta_{1},\varphi_{1})$.

##  Appendix B Wave functions and overlap integrals

The wave functions for the $n=2$ boson around an isolated black hole are given
by

$\displaystyle\varphi_{2,1,1}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{1}{\pi}}r_{b}^{-5/2}re^{-r/2r_{b}}\sin\theta
e^{i\varphi},$ $\displaystyle\varphi_{2,1,-1}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{1}{\pi}}r_{b}^{-5/2}re^{-r/2r_{b}}\sin\theta
e^{-i\varphi},$ $\displaystyle\varphi_{2p_{x}}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{2}{\pi}}r_{b}^{-5/2}re^{-r/2r_{b}}\sin\theta\cos\varphi,$
$\displaystyle\varphi_{2p_{y}}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{2}{\pi}}r_{b}^{-5/2}re^{-r/2r_{b}}\sin\theta\sin\varphi,$
(33)

where $r_{b}$ is the Bohr radius of the boson.

The wave functions for the molecular orbits of the boson in the binary black
hole system are given by

$\displaystyle\varphi_{\sigma}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8N_{1}}\sqrt{\frac{1}{\pi}}\,\bigg{(}r_{b1}^{-5/2}r_{1}e^{-r_{1}/2r_{b1}}\sin\theta_{1}\cos\varphi_{1}-r_{b2}^{-5/2}r_{2}e^{-r_{2}/2r_{b2}}\sin\theta_{2}\cos\varphi_{2}\bigg{)},$
$\displaystyle\varphi_{\sigma^{*}}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8N_{2}}\sqrt{\frac{1}{\pi}}\,\bigg{(}r_{b1}^{-5/2}r_{1}e^{-r_{1}/2r_{b1}}\sin\theta_{1}\cos\varphi_{1}+r_{b2}^{-5/2}r_{2}e^{-r_{2}/2r_{b2}}\sin\theta_{2}\cos\varphi_{2}\bigg{)},$
$\displaystyle\varphi_{\pi}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8N_{3}}\sqrt{\frac{1}{\pi}}\,\bigg{(}r_{b1}^{-5/2}r_{1}e^{-r_{1}/2r_{b1}}\sin\theta_{1}\sin\varphi_{1}+r_{b2}^{-5/2}r_{2}e^{-r_{2}/2r_{b2}}\sin\theta_{2}\sin\varphi_{2}\bigg{)},$
$\displaystyle\varphi_{\pi^{*}}({\bf r})$ $\displaystyle=$
$\displaystyle\frac{1}{8N_{4}}\sqrt{\frac{1}{\pi}}\,\bigg{(}r_{b1}^{-5/2}r_{1}e^{-r_{1}/2r_{b1}}\sin\theta_{1}\sin\varphi_{1}-r_{b2}^{-5/2}r_{2}e^{-r_{2}/2r_{b2}}\sin\theta_{2}\sin\varphi_{2}\bigg{)},$
(34)

where $r_{b1}$ and $r_{b2}$ are the Bohr radius for the boson in central and
companion black holes, respectively. The overlap integrals are given by

$\displaystyle\langle\varphi_{2p_{x}}^{1}\ket{\varphi_{2p_{x}}^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{32\pi}r_{b1}^{-\frac{5}{2}}r_{b2}^{-\frac{5}{2}}\int
r_{1}^{3}r_{2}e^{-\frac{r_{1}}{2r_{b1}}-\frac{r_{2}}{2r_{b2}}}\sin^{2}\theta_{1}\sin\theta_{2}\cos\varphi_{1}\cos\varphi_{2}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1},$
$\displaystyle\bigg{\langle}\varphi_{2p_{x}}^{1}\bigg{|}\frac{\alpha}{r_{2}}\bigg{|}\varphi_{2p_{x}}^{2}\bigg{\rangle}$
$\displaystyle=$
$\displaystyle\frac{\alpha}{32\pi}r_{b1}^{-\frac{5}{2}}r_{b2}^{-\frac{5}{2}}\int
r_{1}^{3}e^{-\frac{r_{1}}{2r_{b1}}-\frac{r_{2}}{2r_{b2}}}\sin^{2}\theta_{1}\sin\theta_{2}\cos\varphi_{1}\cos\varphi_{2}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1},$
$\displaystyle\bigg{\langle}\varphi_{2p_{x}}^{1}\bigg{|}\frac{\alpha}{r_{2}}\bigg{|}\varphi_{2p_{x}}^{1}\bigg{\rangle}$
$\displaystyle=$
$\displaystyle\frac{\alpha}{32\pi}r_{b1}^{-5}\int\frac{r_{1}^{4}}{r_{2}}e^{-\frac{r_{1}}{r_{b1}}}\sin^{3}\theta_{1}\cos^{2}\varphi_{1}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1},$
$\displaystyle\langle\varphi_{2p_{y}}^{1}\ket{\varphi_{2p_{y}}^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{32\pi}r_{b1}^{-\frac{5}{2}}r_{b2}^{-\frac{5}{2}}\int
r_{1}^{3}r_{2}e^{-\frac{r_{1}}{2r_{b1}}-\frac{r_{2}}{2r_{b2}}}\sin^{2}\theta_{1}\sin\theta_{2}\sin\varphi_{1}\sin\varphi_{2}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1},$
$\displaystyle\bigg{\langle}\varphi_{2p_{y}}^{1}\bigg{|}\frac{\alpha}{r_{2}}\bigg{|}\varphi_{2p_{y}}^{2}\bigg{\rangle}$
$\displaystyle=$
$\displaystyle\frac{\alpha}{32\pi}r_{b1}^{-\frac{5}{2}}r_{b2}^{-\frac{5}{2}}\int
r_{1}^{3}e^{-\frac{r_{1}}{2r_{b1}}-\frac{r_{2}}{2r_{b2}}}\sin^{2}\theta_{1}\sin\theta_{2}\sin\varphi_{1}\sin\varphi_{2}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1},$
$\displaystyle\bigg{\langle}\varphi_{2p_{y}}^{1}\bigg{|}\frac{\alpha}{r_{2}}\bigg{|}\varphi_{2p_{y}}^{1}\bigg{\rangle}$
$\displaystyle=$
$\displaystyle\frac{\alpha}{32\pi}r_{b1}^{-5}\int\frac{r_{1}^{4}}{r_{2}}e^{-\frac{r_{1}}{r_{b1}}}\sin^{3}\theta_{1}\sin^{2}\varphi_{1}\,{\rm
d}r_{1}{\rm d}\theta_{1}{\rm d}\varphi_{1}.$ (35)

##  Appendix C Quantum adiabatic theorem

In this appendix, we show that the adiabatic approximation can be applied in
our case. Since the potential generated by the two black holes changes very
slowly, we assume that the energy eigenvalues also evolve slowly and
continuously with time, and they satisfy the equation

$\displaystyle\hat{H}(t)\ket{\psi_{n}(t)}=E_{n}(t)\ket{\psi_{n}(t)}.$ (36)

We assume that the energy levels are not degenerate. An arbitrary state at a
given time can be written as

$\displaystyle\ket{\psi(t)}=\sum_{n}c_{n}(t)\ket{\psi_{n}(t)},$ (37)

and satisfies the Schrödinger equation

$\displaystyle i\frac{\partial}{\partial
t}\ket{\psi(t)}=\hat{H}(t)\ket{\psi(t)},$ (38)

By substituting Eq. (37) into Eq. (38), and multiplying both sides from the
left by $\bra{\psi_{k}}$, we have

$\displaystyle i\frac{\partial c_{k}}{\partial
t}+i\sum_{n}c_{n}\langle\psi_{k}\ket{\frac{\partial\psi_{n}}{\partial
t}}=c_{k}E_{k}.$ (39)

Now we are going to derive the expression for
$\langle\psi_{k}\ket{\frac{\partial\psi_{n}}{\partial t}}$. We start with a
time derivation of Eq. (36),

$\displaystyle\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}+\hat{H}\ket{\frac{\partial\psi_{n}}{\partial
t}}=\frac{\partial E_{n}}{\partial
t}\ket{\psi_{n}}+E_{n}\ket{\frac{\partial\psi_{n}}{\partial t}}.$ (40)

For $k\neq n$, we multiply both sides from the left by $\bra{\psi_{k}}$ and
obtain

$\displaystyle\bra{\psi_{k}}\frac{\partial\psi_{n}}{\partial
t}\rangle=\frac{1}{E_{n}-E_{k}}\bra{\psi_{k}}\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}.$ (41)

We also need to consider the case when $k=n$, namely, the expression for
$\bra{\psi_{k}}\frac{\partial\psi_{k}}{\partial t}\rangle$. Taking time
derivative of the normalization condition

$\displaystyle\bra{\psi_{k}}\psi_{k}\rangle=1,$ (42)

we have

$\displaystyle\bra{\psi_{k}}\frac{\partial\psi_{k}}{\partial
t}\rangle+\bra{\frac{\partial\psi_{k}}{\partial t}}\psi_{k}\rangle=0.$ (43)

It is evident that $\langle\psi_{k}\ket{\frac{\partial\psi_{k}}{\partial t}}$
is purely imaginary, which in principle can be canceled by appropriately
adding a phase. Therefore, the time derivative of the coefficient $c_{n}$ is
given by

$\displaystyle\frac{dc_{k}}{dt}=-\sum_{n\neq
k}c_{n}\frac{1}{E_{n}-E_{k}}\bra{\psi_{k}}\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}-ic_{k}E_{k}.$ (44)

We now show that $\bra{\psi_{k}}\frac{\partial\hat{H}}{\partial
t}\ket{\psi_{n}}=0$ when $n\neq k$, which is the equality given by Eq. (25).
Since $n,k\in\\{\sigma,\sigma^{*},\pi,\pi^{*}\\}$, so there are six off-
diagonal elements for $\frac{\partial\hat{H}}{\partial t}$. Let us first
consider $\bra{\sigma}\frac{\partial\hat{H}}{\partial t}\ket{\sigma^{*}}$ and
$\bra{\pi}\frac{\partial\hat{H}}{\partial t}\ket{\pi^{*}}$.

$\displaystyle\left\langle\sigma\left|\frac{\partial\hat{H}}{\partial
t}\right|\sigma^{\ast}\right\rangle$ (45) $\displaystyle=$
$\displaystyle\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{x}}^{1}\right\rangle-\left\langle\varphi_{2p_{x}}^{2}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{x}}^{2}\right\rangle$ $\displaystyle+$
$\displaystyle\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{x}}^{2}\right\rangle-\left\langle\varphi_{2p_{x}}^{2}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{x}}^{1}\right\rangle.$

In our simple model we assume that two black holes have exactly the same
parameters, so the isolated orbits of the boson are the same for the central
and companion black holes. Due to the symmetry of the configuration of the BH-
cloud-BH system, we have
$\bra{\varphi_{2p_{x}}^{1}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{x}}^{1}}=\bra{\varphi_{2p_{x}}^{2}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{x}}^{2}}$. The wave functions $\varphi_{2p_{x}}^{i}$ are
real and the time derivative of the Hamiltonian is Hermitian, so we have
$\bra{\varphi_{2p_{x}}^{1}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{x}}^{2}}=\bra{\varphi_{2p_{x}}^{2}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{x}}^{1}}$. Therefore, we find that
$\bra{\sigma}\frac{\partial\hat{H}}{\partial t}\ket{\sigma^{*}}=0$. By using
similar arguments we also have $\bra{\pi}\frac{\partial\hat{H}}{\partial
t}\ket{\pi^{*}}=0$.

To calculate other elements of $\frac{\partial\hat{H}}{\partial t}$, we need
to know its explicit expression, which is given by

$\displaystyle\frac{\partial\hat{H}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\partial\hat{H}}{\partial r_{2}}\frac{\partial
r_{2}}{\partial R}\frac{\partial R}{\partial
t}=\frac{\alpha}{r_{2}^{2}}\frac{\partial r_{2}}{\partial R}\frac{\partial
R}{\partial t}$ (46) $\displaystyle=$
$\displaystyle-\frac{64}{5}\alpha\,q(1+q)\frac{M^{3}}{R^{3}}\frac{R-r_{1}\sin\theta_{1}\cos\varphi_{1}}{r_{2}^{3}}.$

By using the explicit expressions of the orbits given by Eqs. (III.1) and
(III.1), we find that other elements of $\frac{\partial\hat{H}}{\partial t}$
can be expanded using
$\bra{\varphi_{2p_{x}}^{i}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{y}}^{j}}$, with $i,j\in\\{1,2\\}$. We now show that
$\bra{\varphi_{2p_{x}}^{i}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{y}}^{j}}=0$ for all choices of $i,j$.

Neglecting the constant factor that is independent of
$(r_{1},\theta_{1},\varphi_{1})$, we find

$\displaystyle\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{1}\right\rangle$ (47) $\displaystyle\sim$
$\displaystyle\int{\rm
d}V\,\frac{R-r_{1}\sin\theta_{1}\cos\varphi_{1}}{r_{2}^{3}}r_{1}^{2}e^{-r_{1}/r_{b}}\sin^{2}\theta_{1}\sin\varphi_{1}\cos\varphi_{1}$
$\displaystyle\sim$ $\displaystyle\int{\rm d}x{\rm d}y{\rm
d}z\,\frac{x(R-x)y}{(r_{1}^{2}+R^{2}-2Rx)^{3/2}}\,e^{-r_{1}/r_{b}},$

where we have used the Cartesian coordinates $(x,y,z)$,

$\displaystyle
x=r_{1}\sin\theta_{1}\cos\varphi_{1},~{}~{}y=r_{1}\sin\theta_{1}\sin\varphi_{1},~{}~{}z=r_{1}\cos\theta_{1}.$

It is evident that the integrand is an odd function of $y$, so the integral is
zero.

$\displaystyle\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{2}\right\rangle$ (48) $\displaystyle\sim$
$\displaystyle\int{\rm
d}V\frac{R-r_{1}\sin\theta_{1}\cos\varphi_{1}}{r_{2}^{3}}e^{-\left(r_{1}+r_{2}\right)/2r_{b}}r_{1}\sin\theta_{1}\cos\varphi_{1}$
$\displaystyle\times\,r_{2}\sin\theta_{2}\sin\varphi_{2}$ $\displaystyle\sim$
$\displaystyle\int{\rm d}x{\rm d}y{\rm
d}z\,\frac{x(R-x)y}{(r_{1}^{2}+R^{2}-2Rx)^{3/2}}\,e^{-(r_{1}+r2)/2r_{b}},$

where $r_{2}=\sqrt{r_{1}^{2}+R^{2}-2Rx}$, and we have used the equality
$y=r_{1}\sin\theta_{1}\cos\varphi_{1}=r_{2}\sin\theta_{2}\cos\phi_{2}$. It is
evident that the integrand is an odd function of $y$, so the integral is zero.

By using the symmetry of the BH-cloud-BH system, we have

$\displaystyle\left\langle\varphi_{2p_{x}}^{2}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{1}\right\rangle$ $\displaystyle=$
$\displaystyle-\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{2}\right\rangle=0,$
$\displaystyle\left\langle\varphi_{2p_{x}}^{2}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{2}\right\rangle$ $\displaystyle=$
$\displaystyle-\left\langle\varphi_{2p_{x}}^{1}\left|\frac{\partial\hat{H}}{\partial
t}\right|\varphi_{2p_{y}}^{1}\right\rangle=0.$

Therefore, We have $\bra{\varphi_{2p_{x}}^{i}}\frac{\partial\hat{H}}{\partial
t}\ket{\varphi_{2p_{y}}^{j}}=0$ for all choices of $i,j$. As a result, all
other elements of $\frac{\partial\hat{H}}{\partial t}$ are zero. In summary,
we have

$\displaystyle\left\langle\psi_{k}\left|\frac{\partial\hat{H}}{\partial
t}\right|\psi_{n}\right\rangle=0$ (50)

for $n\neq k$.

## References

  * Abbott et al. (2016) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
  * Luo et al. (2016) J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, et al., Classical and Quantum Gravity 33, 035010 (2016).
  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977), URL https://link.aps.org/doi/10.1103/PhysRevLett.38.1440.
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978), URL https://link.aps.org/doi/10.1103/PhysRevLett.40.223.
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978), URL https://link.aps.org/doi/10.1103/PhysRevLett.40.279.
  * Svrcek and Witten (2006) P. Svrcek and E. Witten, Journal of High Energy Physics 2006, 051 (2006), URL https://dx.doi.org/10.1088/1126-6708/2006/06/051.
  * Preskill et al. (1983) J. Preskill, M. B. Wise, and F. Wilczek, Physics Letters B 120, 127 (1983), ISSN 0370-2693, URL https://www.sciencedirect.com/science/article/pii/0370269383906378.
  * Abbott and Sikivie (1983) L. Abbott and P. Sikivie, Physics Letters B 120, 133 (1983), ISSN 0370-2693, URL https://www.sciencedirect.com/science/article/pii/037026938390638X.
  * Dine and Fischler (1983) M. Dine and W. Fischler, Physics Letters B 120, 137 (1983), ISSN 0370-2693, URL https://www.sciencedirect.com/science/article/pii/0370269383906391.
  * Li et al. (2014) B. Li, T. Rindler-Daller, and P. R. Shapiro, Phys. Rev. D 89, 083536 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.89.083536.
  * Hui et al. (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.043541.
  * Zel’Dovich (1971) Y. B. Zel’Dovich, Soviet Journal of Experimental and Theoretical Physics Letters 14, 180 (1971).
  * Zel’Dovich (1972) Y. B. Zel’Dovich, Soviet Journal of Experimental and Theoretical Physics 35, 1085 (1972).
  * Misner (1972) C. Misner, in _Bulletin of the American Physical Society_ (1972), vol. 17, p. 472.
  * Starobinskiǐ (1973) A. A. Starobinskiǐ, Soviet Journal of Experimental and Theoretical Physics 37, 28 (1973).
  * Starobinskiǐ and Churilov (1974) A. A. Starobinskiǐ and S. M. Churilov, Soviet Journal of Experimental and Theoretical Physics 38, 1 (1974).
  * East and Pretorius (2017) W. E. East and F. Pretorius, Phys. Rev. Lett. 119, 041101 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.041101.
  * East (2017) W. E. East, Phys. Rev. D 96, 024004 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.024004.
  * Brito et al. (2020) R. Brito, V. Cardoso, and P. Pani, _Superradiance_ , vol. 10 (Springer, 2020).
  * Cromb et al. (2020) M. Cromb, G. M. Gibson, E. Toninelli, M. J. Padgett, E. M. Wright, and D. Faccio, Nature Physics 16, 1069 (2020), URL https://doi.org/10.1038%2Fs41567-020-0944-3.
  * Arvanitaki and Dubovsky (2011) A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.044026.
  * Brito et al. (2017) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. D 96, 064050 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.064050.
  * Arvanitaki et al. (2015) A. Arvanitaki, M. Baryakhtar, and X. Huang, Physical Review D 91 (2015), URL https://doi.org/10.1103%2Fphysrevd.91.084011.
  * Arvanitaki et al. (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Physical Review D 81 (2010), URL https://doi.org/10.1103%2Fphysrevd.81.123530.
  * Baumann et al. (2019) D. Baumann, H. S. Chia, and R. A. Porto, Phys. Rev. D 99, 044001 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.99.044001.
  * Zhang and Yang (2019) J. Zhang and H. Yang, Phys. Rev. D 99, 064018 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.99.064018.
  * Berti et al. (2019) E. Berti, R. Brito, C. F. B. Macedo, G. Raposo, and J. a. L. Rosa, Phys. Rev. D 99, 104039 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.99.104039.
  * Takahashi et al. (2023) T. Takahashi, H. Omiya, and T. Tanaka, Phys. Rev. D 107, 103020 (2023), URL https://link.aps.org/doi/10.1103/PhysRevD.107.103020.
  * Baumann et al. (2022) D. Baumann, G. Bertone, J. Stout, and G. M. Tomaselli, Phys. Rev. Lett. 128, 221102 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.221102.
  * Siemonsen et al. (2023) N. Siemonsen, T. May, and W. E. East, Phys. Rev. D 107, 104003 (2023), URL https://link.aps.org/doi/10.1103/PhysRevD.107.104003.
  * Wong (2020) L. K. Wong, Phys. Rev. D 101, 124049 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.124049.
  * Ikeda et al. (2021) T. Ikeda, L. Bernard, V. Cardoso, and M. Zilhão, Phys. Rev. D 103, 024020 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.024020.
  * Detweiler (1980) S. Detweiler, Phys. Rev. D 22, 2323 (1980), URL https://link.aps.org/doi/10.1103/PhysRevD.22.2323.
  * Cohen-Tannoudji et al. (1986) C. Cohen-Tannoudji, B. Diu, and F. Laloe, Quantum Mechanics 2, 626 (1986).
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.136.B1224.
  * Yoshino and Kodama (2014) H. Yoshino and H. Kodama, Progress of Theoretical and Experimental Physics 2014, 043E02 (2014).


Nb | $(2.5\pm 0.2)\times 10^{4}$ | $1.09\pm 0.12$ | $23.0\pm 1.8$ | $0.15\pm 0.02$ | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $3.9\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$ | $3.1\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$
S1 | $(4.7\pm 0.3)\times 10^{4}$ | $1.06\pm 0.37$ | $31.1\pm 0.3$ | $0.14\pm 0.05$ | $1.4\begin{subarray}{c}+0.3\\\ -0.6\end{subarray}$ | $0.5\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
S2 | $(2.4\pm 0.2)\times 10^{4}$ | $0.74\pm 0.29$ | $31.4\pm 1.3$ | $0.07\pm 0.03$ | $0.7\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $1.8\begin{subarray}{c}+1.1\\\ -1.2\end{subarray}$ | $0.8\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$
N | $(3.4\pm 0.2)\times 10^{4}$ | $1.00\pm 0.31$ | $22.6\pm 0.4$ | $0.16\pm 0.05$ | $0.8\begin{subarray}{c}+0.2\\\ -0.3\end{subarray}$ | $0.8\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$ | $0.6\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$

Note. — Magnetic field strengths and mass-to-flux ratios derived from DCF
method. The uncertainties listed here are obtained from propagating the
observational uncertainties through the corresponding equations using a Monte
Carlo approach. Possible additional systematic uncertainties due to, e.g., the
unknown dust opacity, are not included.

Table 1: Magnetic field strengths and mass-to-flux ratios in NGC 2264, using
C18O.

Leaves | $n_{H_{2}}$ | $\sigma_{v,NT}$ | $\delta\phi$ | $B_{pos}$ | $\lambda$ | $T/W$ | $M/W$
---|---|---|---|---|---|---|---
| $\textrm{cm}^{-3}$ | (km s-1) | (deg) | (mG) |  |  |
S1a1 | $(3.7\pm 0.3)\times 10^{5}$ | $0.97\pm 0.25$ | $34.0\pm 0.4$ | $0.35\pm 0.09$ | $1.6\begin{subarray}{c}+0.4\\\ -0.5\end{subarray}$ | $0.4\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.1\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
S1a2 | $(6.2\pm 0.4)\times 10^{5}$ | $0.69\pm 0.18$ | $9.5\pm 1.4$ | $1.16\pm 0.35$ | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $0.8\begin{subarray}{c}+0.4\\\ -0.3\end{subarray}$ | $3.8\begin{subarray}{c}+2.1\\\ -2.1\end{subarray}$
S21 | $(5.6\pm 0.4)\times 10^{4}$ | $0.58\pm 0.09$ | $14.3\pm 2.1$ | $0.19\pm 0.04$ | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.5\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$ | $2.9\begin{subarray}{c}+1.2\\\ -1.2\end{subarray}$
N1 | $(1.5\pm 0.1)\times 10^{5}$ | $0.90\pm 0.07$ | $8.5\pm 0.4$ | $0.82\pm 0.08$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.2\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $7.0\begin{subarray}{c}+1.3\\\ -1.3\end{subarray}$
Na1 | $(2.2\pm 0.2)\times 10^{5}$ | $0.92\pm 0.11$ | $16.4\pm 2.0$ | $0.55\pm 0.09$ | $0.4\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.4\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $2.3\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$
Na2 | $(1.4\pm 0.1)\times 10^{5}$ | $1.22\pm 0.15$ | $11.4\pm 0.7$ | $0.80\pm 0.11$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $2.0\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$ | $6.4\begin{subarray}{c}+1.7\\\ -1.7\end{subarray}$
Na3 | $(1.6\pm 0.1)\times 10^{5}$ | $0.78\pm 0.07$ | $9.3\pm 0.9$ | $0.67\pm 0.09$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.7\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $8.0\begin{subarray}{c}+2.0\\\ -2.0\end{subarray}$
Nb1 | $(8.7\pm 0.6)\times 10^{4}$ | $0.70\pm 0.08$ | $28.0\pm 2.7$ | $0.15\pm 0.02$ | $0.6\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.6\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $0.8\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
S1a | $(1.5\pm 0.1)\times 10^{5}$ | $0.94\pm 0.24$ | $30.2\pm 0.3$ | $0.24\pm 0.06$ | $1.4\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$ | $0.3\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
Na | $(5.1\pm 0.4)\times 10^{4}$ | $1.04\pm 0.23$ | $21.5\pm 0.4$ | $0.22\pm 0.05$ | $0.6\begin{subarray}{c}+0.1\\\ -0.2\end{subarray}$ | $0.9\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$ | $0.8\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$
Nb | $(2.5\pm 0.2)\times 10^{4}$ | $0.72\pm 0.09$ | $23.0\pm 1.8$ | $0.10\pm 0.01$ | $0.5\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.8\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$ | $1.3\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$
S1 | $(4.7\pm 0.3)\times 10^{4}$ | $0.91\pm 0.25$ | $31.1\pm 0.3$ | $0.13\pm 0.04$ | $1.3\begin{subarray}{c}+0.5\\\ -0.3\end{subarray}$ | $0.4\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
S2 | $(2.4\pm 0.2)\times 10^{4}$ | $0.57\pm 0.09$ | $31.4\pm 1.3$ | $0.06\pm 0.01$ | $0.8\begin{subarray}{c}+0.1\\\ -0.2\end{subarray}$ | $1.0\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $0.4\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N | $(3.4\pm 0.2)\times 10^{4}$ | $1.01\pm 0.23$ | $22.6\pm 0.4$ | $0.16\pm 0.04$ | $0.7\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.7\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$ | $0.6\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$

Note. — Identical to Table 1.

Table 1: Magnetic field strengths and mass-to-flux ratios in NGC 2264, using
N2H+.

Leaves | $n_{H_{2}}$ | $\sigma_{v,NT}$ | $\delta\phi$ | $B_{pos}$ | $\lambda$ | $T/W$ | $M/W$
---|---|---|---|---|---|---|---
| $\textrm{cm}^{-3}$ | (km s-1) | (deg) | (mG) |  |  |
N1 | $(1.5\pm 0.2)\times 10^{5}$ | $0.76\pm 0.04$ | $8.5\pm 0.4$ | $0.68\pm 0.05$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $0.9\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $4.9\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$
Na1 | $(2.2\pm 0.2)\times 10^{5}$ | $0.61\pm 0.12$ | $16.4\pm 2.0$ | $0.34\pm 0.08$ | $0.6\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $0.7\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $1.0\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$
Na2 | $(1.4\pm 0.1)\times 10^{5}$ | $1.06\pm 0.18$ | $11.4\pm 0.7$ | $0.70\pm 0.13$ | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.5\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$ | $4.9\begin{subarray}{c}+1.7\\\ -1.7\end{subarray}$
Na3 | $(1.6\pm 0.1)\times 10^{5}$ | $0.59\pm 0.04$ | $9.3\pm 0.9$ | $0.48\pm 0.06$ | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1.0\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $4.2\begin{subarray}{c}+0.9\\\ -0.9\end{subarray}$
Nb1 | $(8.7\pm 0.6)\times 10^{4}$ | $0.41\pm 0.07$ | $28.0\pm 2.7$ | $0.09\pm 0.02$ | $1.0\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.7\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
Na | $(5.1\pm 0.4)\times 10^{4}$ | $0.86\pm 0.13$ | $21.5\pm 0.4$ | $0.18\pm 0.03$ | $0.7\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $0.6\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.6\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
Nb | $(2.5\pm 0.2)\times 10^{4}$ | $0.53\pm 0.15$ | $23.0\pm 1.8$ | $0.07\pm 0.02$ | $0.7\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $1.1\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$ | $0.8\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$
N | $(3.4\pm 0.2)\times 10^{4}$ | $0.87\pm 0.13$ | $22.6\pm 0.4$ | $0.14\pm 0.02$ | $0.8\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $0.5\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $0.4\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$

Note. — Identical to Table 1.

## References

  * Allen (1972) Allen, D. A. 1972, ApJ, 172, L55, doi: 10.1086/180890
  * Andersson et al. (2015) Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501, doi: 10.1146/annurev-astro-082214-122414
  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27, doi: 10.2458/azu_uapress_9780816531240-ch002
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102, doi: 10.1051/0004-6361/201014666
  * Arzoumanian et al. (2011) Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6, doi: 10.1051/0004-6361/201116596
  * Arzoumanian et al. (2021) Arzoumanian, D., Furuya, R. S., Hasegawa, T., et al. 2021, A&A, 647, A78, doi: 10.1051/0004-6361/202038624
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Berry et al. (2005) Berry, D. S., Gledhill, T. M., Greaves, J. S., & Jenness, T. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 343, Astronomical Polarimetry: Current Status and Future Directions, ed. A. Adamson, C. Aspin, C. Davis, & T. Fujiyoshi, 71
  * Buckle et al. (2012) Buckle, J. V., Richer, J. S., & Davis, C. J. 2012, MNRAS, 423, 1127, doi: 10.1111/j.1365-2966.2012.20941.x
  * Buckle et al. (2009) Buckle, J. V., Hills, R. E., Smith, H., et al. 2009, MNRAS, 399, 1026, doi: 10.1111/j.1365-2966.2009.15347.x
  * Buckner et al. (2020) Buckner, A. S. M., Khorrami, Z., González, M., et al. 2020, A&A, 636, A80, doi: 10.1051/0004-6361/201936935
  * Burkhart et al. (2015) Burkhart, B., Collins, D. C., & Lazarian, A. 2015, ApJ, 808, 48, doi: 10.1088/0004-637X/808/1/48
  * Busquet et al. (2013) Busquet, G., Zhang, Q., Palau, A., et al. 2013, ApJ, 764, L26, doi: 10.1088/2041-8205/764/2/L26
  * Cavanagh (2008) Cavanagh, B. 2008, Starlink User Note, 260
  * Chandrasekhar & Fermi (1953) Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113, doi: 10.1086/145731
  * Chapin et al. (2013) Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545, doi: 10.1093/mnras/stt052
  * Chung et al. (2022) Chung, E. J., Lee, C. W., Kwon, W., et al. 2022, AJ, 164, 175, doi: 10.3847/1538-3881/ac8a43
  * Chung et al. (2019) Chung, E. J., Lee, C. W., Kim, S., et al. 2019, ApJ, 877, 114, doi: 10.3847/1538-4357/ab12d1
  * Clarke et al. (2018) Clarke, S. D., Whitworth, A. P., Spowage, R. L., et al. 2018, MNRAS, 479, 1722, doi: 10.1093/mnras/sty1675
  * Crutcher et al. (2004) Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ, 600, 279, doi: 10.1086/379705
  * Cunningham et al. (2016) Cunningham, N., Lumsden, S. L., Cyganowski, C. J., Maud, L. T., & Purcell, C. 2016, MNRAS, 458, 1742, doi: 10.1093/mnras/stw359
  * Currie et al. (2014) Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay, 391
  * Davis (1951) Davis, L. 1951, Physical Review, 81, 890, doi: 10.1103/PhysRev.81.890.2
  * Dobashi et al. (2014) Dobashi, K., Matsumoto, T., Shimoikura, T., et al. 2014, ApJ, 797, 58, doi: 10.1088/0004-637X/797/1/58
  * Evans et al. (2009) Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321, doi: 10.1088/0067-0049/181/2/321
  * Federrath (2016) Federrath, C. 2016, MNRAS, 457, 375, doi: 10.1093/mnras/stv2880
  * Fiege & Pudritz (2000) Fiege, J. D., & Pudritz, R. E. 2000, MNRAS, 311, 105, doi: 10.1046/j.1365-8711.2000.03067.x
  * Frau et al. (2015) Frau, P., Girart, J. M., Alves, F. O., et al. 2015, A&A, 574, L6, doi: 10.1051/0004-6361/201425234
  * Friberg et al. (2016) Friberg, P., Bastien, P., Berry, D., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9914, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VIII, 991403, doi: 10.1117/12.2231943
  * Ganguly et al. (2023) Ganguly, S., Walch, S., Seifried, D., Clarke, S. D., & Weis, M. 2023, MNRAS, 525, 721, doi: 10.1093/mnras/stad2054
  * Gómez & Vázquez-Semadeni (2014) Gómez, G. C., & Vázquez-Semadeni, E. 2014, ApJ, 791, 124, doi: 10.1088/0004-637X/791/2/124
  * Gómez et al. (2018) Gómez, G. C., Vázquez-Semadeni, E., & Zamora-Avilés, M. 2018, MNRAS, 480, 2939, doi: 10.1093/mnras/sty2018
  * Goodman et al. (1993) Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528, doi: 10.1086/172465
  * Goodman et al. (2009) Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., et al. 2009, Nature, 457, 63, doi: 10.1038/nature07609
  * Henshaw et al. (2016a) Henshaw, J. D., Longmore, S. N., Kruijssen, J. M. D., et al. 2016a, MNRAS, 457, 2675, doi: 10.1093/mnras/stw121
  * Henshaw et al. (2016b) —. 2016b, SCOUSE: Semi-automated multi-COmponent Universal Spectral-line fitting Engine, Astrophysics Source Code Library, record ascl:1601.003.  http://ascl.net/1601.003
  * Henshaw et al. (2019) Henshaw, J. D., Ginsburg, A., Haworth, T. J., et al. 2019, MNRAS, 485, 2457, doi: 10.1093/mnras/stz471
  * Heyer et al. (2008) Heyer, M., Gong, H., Ostriker, E., & Brunt, C. 2008, ApJ, 680, 420, doi: 10.1086/587510
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Hildebrand et al. (2009) Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567, doi: 10.1088/0004-637X/696/1/567
  * Hwang et al. (2021) Hwang, J., Kim, J., Pattle, K., et al. 2021, ApJ, 913, 85, doi: 10.3847/1538-4357/abf3c4
  * Hwang et al. (2022) —. 2022, ApJ, 941, 51, doi: 10.3847/1538-4357/ac99e0
  * Ibáñez-Mejía et al. (2022) Ibáñez-Mejía, J. C., Mac Low, M.-M., & Klessen, R. S. 2022, ApJ, 925, 196, doi: 10.3847/1538-4357/ac3b58
  * Inoue et al. (2018) Inoue, T., Hennebelle, P., Fukui, Y., et al. 2018, PASJ, 70, S53, doi: 10.1093/pasj/psx089
  * Inoue & Inutsuka (2012) Inoue, T., & Inutsuka, S.-i. 2012, ApJ, 759, 35, doi: 10.1088/0004-637X/759/1/35
  * Inutsuka et al. (2015) Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49, doi: 10.1051/0004-6361/201425584
  * Inutsuka et al. (2016) Inutsuka, S.-I., Inoue, T., Iwasaki, K., Hosokawa, T., & Kobayashi, M. I. N. 2016, in From Interstellar Clouds to Star-Forming Galaxies: Universal Processes?, ed. P. Jablonka, P. André, & F. van der Tak, Vol. 315, 61–68, doi: 10.1017/S1743921316007262
  * Jones et al. (2015) Jones, T. J., Bagley, M., Krejny, M., Andersson, B. G., & Bastien, P. 2015, AJ, 149, 31, doi: 10.1088/0004-6256/149/1/31
  * Jones et al. (2016) Jones, T. J., Gordon, M., Shenoy, D., et al. 2016, AJ, 151, 156, doi: 10.3847/0004-6256/151/6/156
  * Jow et al. (2018) Jow, D. L., Hill, R., Scott, D., et al. 2018, MNRAS, 474, 1018, doi: 10.1093/mnras/stx2736
  * Koch et al. (2012a) Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2012a, ApJ, 747, 79, doi: 10.1088/0004-637X/747/1/79
  * Koch et al. (2012b) —. 2012b, ApJ, 747, 80, doi: 10.1088/0004-637X/747/1/80
  * Koch et al. (2013) —. 2013, ApJ, 775, 77, doi: 10.1088/0004-637X/775/1/77
  * Koch et al. (2018) Koch, P. M., Tang, Y.-W., Ho, P. T. P., et al. 2018, ApJ, 855, 39, doi: 10.3847/1538-4357/aaa4c1
  * Koch et al. (2014) —. 2014, ApJ, 797, 99, doi: 10.1088/0004-637X/797/2/99
  * Koch et al. (2022) —. 2022, ApJ, 940, 89, doi: 10.3847/1538-4357/ac96e3
  * Kumar et al. (2020) Kumar, M. S. N., Palmeirim, P., Arzoumanian, D., & Inutsuka, S. I. 2020, A&A, 642, A87, doi: 10.1051/0004-6361/202038232
  * Kwon et al. (2022) Kwon, W., Pattle, K., Sadavoy, S., et al. 2022, ApJ, 926, 163, doi: 10.3847/1538-4357/ac4bbe
  * Li et al. (2013) Li, H.-b., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707, doi: 10.1093/mnras/stt1849
  * Li et al. (2014) Li, H. B., Goodman, A., Sridharan, T. K., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 101, doi: 10.2458/azu_uapress_9780816531240-ch005
  * Li & Houde (2008) Li, H.-b., & Houde, M. 2008, ApJ, 677, 1151, doi: 10.1086/529581
  * Liu et al. (2023a) Liu, H.-L., Tej, A., Liu, T., et al. 2023a, MNRAS, 522, 3719, doi: 10.1093/mnras/stad047
  * Liu et al. (2019) Liu, J., Qiu, K., Berry, D., et al. 2019, ApJ, 877, 43, doi: 10.3847/1538-4357/ab0958
  * Liu et al. (2023b) Liu, J., Zhang, Q., Koch, P. M., et al. 2023b, ApJ, 945, 160, doi: 10.3847/1538-4357/acb540
  * Mairs et al. (2016) Mairs, S., Johnstone, D., Kirk, H., et al. 2016, MNRAS, 461, 4022, doi: 10.1093/mnras/stw1550
  * Mairs et al. (2021) Mairs, S., Dempsey, J. T., Bell, G. S., et al. 2021, AJ, 162, 191, doi: 10.3847/1538-3881/ac18bf
  * Margulis et al. (1988) Margulis, M., Lada, C. J., & Snell, R. L. 1988, ApJ, 333, 316, doi: 10.1086/166748
  * Maury et al. (2009) Maury, A. J., André, P., & Li, Z. Y. 2009, A&A, 499, 175, doi: 10.1051/0004-6361/200811442
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Mestel & Spitzer (1956) Mestel, L., & Spitzer, L., J. 1956, MNRAS, 116, 503, doi: 10.1093/mnras/116.5.503
  * Motte et al. (2018) Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41, doi: 10.1146/annurev-astro-091916-055235
  * Myers (2009) Myers, P. C. 2009, ApJ, 700, 1609, doi: 10.1088/0004-637X/700/2/1609
  * Nakamura & Li (2008) Nakamura, F., & Li, Z.-Y. 2008, ApJ, 687, 354, doi: 10.1086/591641
  * Nakamura et al. (2014) Nakamura, F., Sugitani, K., Tanaka, T., et al. 2014, ApJ, 791, L23, doi: 10.1088/2041-8205/791/2/L23
  * Nakano & Nakamura (1978) Nakano, T., & Nakamura, T. 1978, PASJ, 30, 671
  * Naranjo-Romero et al. (2022) Naranjo-Romero, R., Vázquez-Semadeni, E., & Loughnane, R. M. 2022, MNRAS, 512, 4715, doi: 10.1093/mnras/stac804
  * Nony et al. (2021) Nony, T., Robitaille, J. F., Motte, F., et al. 2021, A&A, 645, A94, doi: 10.1051/0004-6361/202039353
  * Oliphant (2006–) Oliphant, T. 2006–, NumPy: A guide to NumPy, USA: Trelgol Publishing.  http://www.numpy.org/
  * Ostriker et al. (2001) Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980, doi: 10.1086/318290
  * Padoan et al. (2001) Padoan, P., Juvela, M., Goodman, A. A., & Nordlund, Å. 2001, ApJ, 553, 227, doi: 10.1086/320636
  * Palmeirim et al. (2013) Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38, doi: 10.1051/0004-6361/201220500
  * Pattle et al. (2017) Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122, doi: 10.3847/1538-4357/aa80e5
  * Pattle et al. (2019) Pattle, K., Lai, S.-P., Hasegawa, T., et al. 2019, ApJ, 880, 27, doi: 10.3847/1538-4357/ab286f
  * Pattle et al. (2021) Pattle, K., Lai, S.-P., Di Francesco, J., et al. 2021, ApJ, 907, 88, doi: 10.3847/1538-4357/abcc6c
  * Peretto et al. (2006) Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979, doi: 10.1051/0004-6361:20053324
  * Peretto et al. (2013) Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112, doi: 10.1051/0004-6361/201321318
  * Pillai et al. (2020) Pillai, T. G. S., Clemens, D. P., Reissl, S., et al. 2020, Nature Astronomy, 4, 1195, doi: 10.1038/s41550-020-1172-6
  * Pineda et al. (2022) Pineda, J. E., Arzoumanian, D., André, P., et al. 2022, arXiv e-prints, arXiv:2205.03935, doi: 10.48550/arXiv.2205.03935
  * Planck Collaboration et al. (2015) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, A&A, 576, A105, doi: 10.1051/0004-6361/201424086
  * Planck Collaboration et al. (2016) —. 2016, A&A, 586, A138, doi: 10.1051/0004-6361/201525896
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A12, doi: 10.1051/0004-6361/201833885
  * Priestley & Whitworth (2022) Priestley, F. D., & Whitworth, A. P. 2022, MNRAS, 509, 1494, doi: 10.1093/mnras/stab2816
  * Rapson et al. (2014) Rapson, V. A., Pipher, J. L., Gutermuth, R. A., et al. 2014, ApJ, 794, 124, doi: 10.1088/0004-637X/794/2/124
  * Robitaille (2019) Robitaille, T. 2019, APLpy v2.0: The Astronomical Plotting Library in Python, doi: 10.5281/zenodo.2567476
  * Robitaille & Bressert (2012) Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library.  http://ascl.net/1208.017
  * Sadavoy et al. (2013) Sadavoy, S. I., Di Francesco, J., Johnstone, D., et al. 2013, ApJ, 767, 126, doi: 10.1088/0004-637X/767/2/126
  * Schneider et al. (2012) Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11, doi: 10.1051/0004-6361/201118566
  * Schneider & Elmegreen (1979) Schneider, S., & Elmegreen, B. G. 1979, ApJS, 41, 87, doi: 10.1086/190609
  * Seifried et al. (2020) Seifried, D., Walch, S., Weis, M., et al. 2020, MNRAS, 497, 4196, doi: 10.1093/mnras/staa2231
  * Soler (2019) Soler, J. D. 2019, A&A, 629, A96, doi: 10.1051/0004-6361/201935779
  * Soler & Hennebelle (2017) Soler, J. D., & Hennebelle, P. 2017, A&A, 607, A2, doi: 10.1051/0004-6361/201731049
  * Soler et al. (2013) Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128, doi: 10.1088/0004-637X/774/2/128
  * Sousbie (2011) Sousbie, T. 2011, MNRAS, 414, 350, doi: 10.1111/j.1365-2966.2011.18394.x
  * Storm et al. (2016) Storm, S., Mundy, L. G., Lee, K. I., et al. 2016, ApJ, 830, 127, doi: 10.3847/0004-637X/830/2/127
  * Sugitani et al. (2011) Sugitani, K., Nakamura, F., Watanabe, M., et al. 2011, ApJ, 734, 63, doi: 10.1088/0004-637X/734/1/63
  * Sung et al. (2009) Sung, H., Stauffer, J. R., & Bessell, M. S. 2009, AJ, 138, 1116, doi: 10.1088/0004-6256/138/4/1116
  * Suri et al. (2019) Suri, S., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 623, A142, doi: 10.1051/0004-6361/201834049
  * Tahani et al. (2019) Tahani, M., Plume, R., Brown, J. C., Soler, J. D., & Kainulainen, J. 2019, A&A, 632, A68, doi: 10.1051/0004-6361/201936280
  * Tahani et al. (2022a) Tahani, M., Glover, J., Lupypciw, W., et al. 2022a, A&A, 660, L7, doi: 10.1051/0004-6361/202243322
  * Tahani et al. (2022b) Tahani, M., Lupypciw, W., Glover, J., et al. 2022b, A&A, 660, A97, doi: 10.1051/0004-6361/202141170
  * Tobin et al. (2015) Tobin, J. J., Hartmann, L., Fűrész, G., Hsu, W.-H., & Mateo, M. 2015, AJ, 149, 119, doi: 10.1088/0004-6256/149/4/119
  * Van Loo et al. (2014) Van Loo, S., Keto, E., & Zhang, Q. 2014, ApJ, 789, 37, doi: 10.1088/0004-637X/789/1/37
  * Vázquez-Semadeni et al. (2017) Vázquez-Semadeni, E., González-Samaniego, A., & Colín, P. 2017, MNRAS, 467, 1313, doi: 10.1093/mnras/stw3229
  * Vázquez-Semadeni et al. (2019) Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061, doi: 10.1093/mnras/stz2736
  * Virtanen et al. (2019) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2019, arXiv e-prints, arXiv:1907.10121.  https://arxiv.org/abs/1907.10121
  * Wang et al. (2020a) Wang, J.-W., Koch, P. M., Galván-Madrid, R., et al. 2020a, ApJ, 905, 158, doi: 10.3847/1538-4357/abc74e
  * Wang et al. (2020b) Wang, J.-W., Lai, S.-P., Clemens, D. P., et al. 2020b, ApJ, 888, 13, doi: 10.3847/1538-4357/ab5c1c
  * Wang et al. (2019) Wang, J.-W., Lai, S.-P., Eswaraiah, C., et al. 2019, ApJ, 876, 42, doi: 10.3847/1538-4357/ab13a2
  * Wang et al. (2022) Wang, J.-W., Koch, P. M., Tang, Y.-W., et al. 2022, ApJ, 931, 115, doi: 10.3847/1538-4357/ac6872
  * Ward-Thompson et al. (2000) Ward-Thompson, D., Zylka, R., Mezger, P. G., & Sievers, A. W. 2000, A&A, 355, 1122
  * Ward-Thompson et al. (2017) Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66, doi: 10.3847/1538-4357/aa70a0
  * Wardle & Kronberg (1974) Wardle, J. F. C., & Kronberg, P. P. 1974, ApJ, 194, 249, doi: 10.1086/153240
  * Williams et al. (2018) Williams, G. M., Peretto, N., Avison, A., Duarte-Cabral, A., & Fuller, G. A. 2018, A&A, 613, A11, doi: 10.1051/0004-6361/201731587
  * Zhou et al. (2022) Zhou, J.-W., Liu, T., Evans, N. J., et al. 2022, MNRAS, 514, 6038, doi: 10.1093/mnras/stac1735
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51, doi: 10.1051/0004-6361/201936145