# Gravitational waves and kicks from the merger of
unequal mass, highly compact boson stars

Miguel Bezares1,2, Mateja Bošković1,2, Steven Liebling3, Carlos Palenzuela4,
Paolo Pani5, Enrico Barausse1,2 1SISSA, Via Bonomea 265, 34136 Trieste, Italy
and INFN Sezione di Trieste  2IFPU - Institute for Fundamental Physics of the
Universe, Via Beirut 2, 34014 Trieste, Italy  3Long Island University,
Brookville, New York 11548, USA  4Departament de Física $\&$ IAC3, Universitat
de les Illes Balears and Institut d’Estudis Espacials de Catalunya, Palma de
Mallorca, Baleares E-07122, Spain  5Dipartimento di Fisica, “Sapienza”
Università di Roma & Sezione INFN Roma1, P.A. Moro 5, 00185, Roma, Italy

###### Abstract

Boson stars have attracted much attention in recent decades as simple, self-
consistent models of compact objects and also as self-gravitating structures
formed in some dark-matter scenarios. Direct detection of these hypothetical
objects through electromagnetic signatures would be unlikely because their
bosonic constituents are not expected to interact significantly with ordinary
matter and radiation. However, binary boson stars might form and coalesce
emitting a detectable gravitational wave signal which might distinguish them
from ordinary compact object binaries containing black holes and neutron
stars. We study the merger of two boson stars by numerically evolving the
fully relativistic Einstein-Klein-Gordon equations for a complex scalar field
with a solitonic potential that generates very compact boson stars. Owing to
the steep mass-radius diagram, we can study the dynamics and gravitational
radiation from unequal-mass binary boson stars with mass ratios up to
$q\approx 23$ without the difficulties encountered when evolving binary black
holes with large mass ratios. Similar to the previously-studied equal-mass
case, our numerical evolutions of the merger produce either a nonspinning
boson star or a spinning black hole, depending on the initial masses and on
the binary angular momentum. We do not find any evidence of synchronized
scalar clouds forming around either the remnant spinning black hole or around
the remnant boson stars. Interestingly, in contrast to the equal-mass case,
one of the mechanisms to dissipate angular momentum is now asymmetric, and
leads to large kick velocities (up to a few $10^{4}\,{\rm km/s}$) which could
produce wandering remnant boson stars. We also compare the gravitational wave
signals predicted from boson star binaries with those from black hole
binaries, and comment on the detectability of the differences with ground
interferometers.

###### Contents

  1. I Introduction
  2. II Setup
    1. II.1 Einstein-Klein-Gordon equations
    2. II.2 Binary initial data
    3. II.3 Numerical setup and analysis
  3. III Dynamics for unequal-mass BS binaries
    1. III.1 Binary dynamics in the inspiral
    2. III.2 Final fate of the binary merger
    3. III.3 Remnants: Scalar clouds, blobs, and kicks
    4. III.4 Collision of boson and anti-boson stars
  4. IV Gravitational Wave Signal
    1. IV.1 Late inspiral and merger
    2. IV.2 Post-merger
    3. IV.3 Solitonic BSs in the LIGO/Virgo band
  5. V Conclusions

##  I Introduction

We are well into the era of gravitational wave (GW) astronomy with the rapidly
growing catalog of GW events detected by the LIGO-Virgo collaboration
LIGOScientific:2014pky ; TheVirgo:2014hva .

With the very recent release of the third GW transient catalog
LIGOScientific:2021djp , the total number of reported coalescences increased
to $90$. Some of the more remarkable events detected to date include:

  * •

GW190412 LIGOScientific:2020stg , a binary black hole (BBH) with asymmetric
component masses, showing evidence for higher harmonics in its GW signal;

  * •

GW190425 LIGOScientific:2020aai , identified with a binary neutron star (NS)
merger lacking evidence of an electromagnetic counterpart;

  * •

GW190521 LIGOScientific:2020iuh , a BBH with a total mass greater than 150
solar masses, which is the most massive binary yet detected, in which the
posterior distribution of the primary mass is nearly entirely in the pair-
instability supernova mass gap where BHs are not expected to form from the
collapse of massive stars;

  * •

GW190814 LIGOScientific:2020zkf , a highly asymmetric system consistent with
the merger of a 23 solar mass black hole (BH) with a 2.6 solar mass compact
object, making the latter either the lightest BH or the heaviest NS observed
in a compact binary;

  * •

GW200105 and GW200115 LIGOScientific:2021qlt , which are the first detections
consistent with a NS-BH merger.

The planned upgrades by the LIGO-Virgo collaboration and the addition of the
KAGRA detector KAGRA:2020agh  promise even more exciting observations in the
future.

A primary target of GW observations is the merger of very compact objects,
with BHs and NSs being the most natural candidates. However, a number of other
hypothetical compact objects have been proposed, called exotic compact objects
(ECOs) Giudice:2016zpa ; Cardoso:2019rvt . The motivations for various ECOs
arise both in beyond-Standard-Model theories and in modified-gravity
scenarios, and some of the most popular models include fuzzballs Mathur:2005zp
, gravastars Mazur:2001fv , wormholes Damour:2007ap , anisotropic stars
Raposo:2018rjn , and boson stars (BSs) Ruffini:1969qy . Phenomenological
studies of ECOs are required to perform actual searches for their signatures.
No evidence for such ECOs has yet been found, but, because they are expected
to be too dim electromagnetically, it is mostly through GW detections that we
can hope to observe them Cardoso:2019rvt .

In this work we study BSs, which are solutions of the Einstein equations
coupled to a complex scalar field with a harmonic time dependence describing a
macroscopic wave-function of a Bose-Einstein condensate (see Ref.
Schunck:2003kk ; liebpa ; Visinelli:2021uve  for reviews). BSs are
particularly promising as possible astrophysical objects because: (i) a
formation mechanism for BSs has been identified, known as gravitational
cooling Seidel:1993zk ; Guzman:2006yc , whereby BSs can be produced from
arbitrary scalar field configurations, (ii) their stability properties
resemble those of NSs so that static BSs below a critical mass are radially
stable Colpi:1986ye ; liebpa ; mace , and finally (iii) BSs have been invoked
in open problems in cosmological and particle physics, such as the nature of
the dark matter and the possibility of early Universe remnants. For instance,
the idea that dark matter is composed of ultra-light bosonic fields has
received significant attention recently Hui:2016ltb ; Marsh:2016rep ;
Hui:2021tkt . Although leading candidates for this kind of dark matter are
real scalars that are organized in time-dependent configurations Seidel:1991zh
, BSs can serve as a proxy for such configurations Brito:2015yfh . Some of
these scenarios can allow for compact BSs (or similar objects) to be produced
in the early Universe Troitsky:2015mda ; Krippendorf:2018tei ; Levkov:2018kau
; Cotner:2019ykd ; Amin:2019ums ; Widdicombe:2018oeo ; Arvanitaki:2019rax

Collisions of BSs have been studied extensively, including: head-on and
orbital mergers of mini-BSs pale1 ; pale2 , head-on mergers of oscillatons
Brito:2015yfh ; Helfer:2018vtq , orbital collisions of solitonic BSs bezpalen
; PhysRevD.96.104058 , and head-on and orbital mergers of Proca stars
Brito:2015pxa ; Sanchis-Gual:2017bhw  as a possible alternative explanation of
the GW190521 event PhysRevD.99.024017 ; Bustillo:2020syj . The merger of ECOs
can be studied within various dark matter scenarios as well, as for example:
mergers between a NSs and a star made of axions, one of the most popular dark
matter type candidates Dietrich:2018bvi ; Dietrich:2018jov ; Clough:2018exo ,
mergers of dark stars composed of bosonic fields Bezares:2018qwa , or mergers
of binary NSs containing a small fraction of dark matter Bezares:2019jcb
modeled using fermion-BSs suspalen .

Motivated by the recent GW detections of very unequal mass binary mergers, we
study here the coalescence of unequal mass BS binaries, focusing on their
dynamics and GW radiation. As in our previous works bezpalen ;
PhysRevD.96.104058 ; Bezares:2018qwa , we adopt the nontopological solitonic
BS potential frie  to construct our asymmetric binaries because: (i) it allows
for very compact configurations that reach a maximum compactness (see below
for its definition) in the stable branch of approximately $C\approx 0.35$
Boskovic:2021nfs ; Cardoso:2021ehg , and (ii) one can construct binaries with
a large mass ratio. Indeed, defining the mass ratio $q\equiv m_{1}/m_{2}$ such
that $m_{1}>m_{2},$ we can produce compact binaries with a mass ratio
ranging111Solitonic BSs in general admit two stable and two unstable branches
Tamaki:2011zza ; Boskovic:2021nfs . Here we focus on the more massive stable
branch, while the other stable branch corresponds to the weak-field regime of
mini BSs for our choice of the potential parameters Boskovic:2021nfs .
approximately from $1$ to $45$. Here, we focus on binaries within the range
$q\in[2,23]$. We note that in contrast to the difficulties encountered when
evolving BBH with large mass ratios Gonzalez:2008bi ; Lousto:2010qx ;
Lousto:2010ut ; BBHune , these evolutions require no change to the choice of
coordinates, namely gamma-driver shift condition, nor an exceptionally high
resolution. The reason for this difference is because the radii of solitonic
BSs even with vastly different masses are of the same order, whereas the
radius of the BH scales linearly with the mass, and therefore a large mass
ratio in a BH binary necessarily implies a large separation of length scales.

Our mergers of unequal mass solitonic BSs produce either a non-rotating BS or
a spinning BH, as in the equal-mass cases PhysRevD.96.104058 . In the former
cases, all the angular momentum is emitted to infinity through scalar field
and GW radiation, while in the latter case, after performing a very long-term
simulation, we find no indication of a scalar cloud synchronized with the
rotation of the remnant BH, as found in Ref. Sanchis-Gual:2020mzb . For one of
our simulations with large angular momentum, a blob of scalar field is ejected
after the merger, producing a significant kick velocity of the remnant. Note
that, this blob ejection has already been observed in solitonic BS binaries of
equal mass PhysRevD.96.104058 . Additionally, we study the dynamics and GW
radiation of a binary composed of a BS and an anti-boson (aBS) star, i.e. with
the opposite frequency, allowing some annihilation of the Noether charge
during the merger.

This work is organized as follows: in Sec. II, we review the evolution
equations describing BSs, followed by the construction of initial data for
binary BSs and numerical implementation. In Sec. III, the coalescence of
unequal-mass BS binaries is studied in detail. The GWs produced by these
systems are explored in Sec. IV, in particular, analyzing the imprint of
higher-order modes in the signal and the post-merger frequencies of the
remnant’s signal. In Sec. V, we summarize our results. We use geometric units
in which $G=1$ and $c=1,$ unless otherwise stated.

##  II Setup

In this section, we briefly summarize the evolution equations describing a
self-gravitating (complex) scalar field and the construction of binary BSs in
quasicircular orbits that constitute the initial data. We also outline the
numerical methods and grid setup employed to perform the simulations. Notice
that our setup is very similar to the one used in Ref. PhysRevD.96.104058
(hereafter Paper I) for studying equal-mass binary BSs.

###  II.1 Einstein-Klein-Gordon equations

Self-gravitating (complex) scalar-fields are described by the Einstein-Klein-
Gordon (EKG) equations

$\displaystyle R_{ab}-\frac{1}{2}g_{ab}R$ $\displaystyle=$ $\displaystyle
8\pi\,T_{ab}\,,$ (1) $\displaystyle g^{ab}\nabla_{a}\nabla_{b}\Phi$
$\displaystyle=$ $\displaystyle\frac{dV}{d|\Phi|^{2}}\Phi\,,$ (2)

where $R_{ab}$ is the Ricci tensor associated with the metric $g_{ab}$, $\Phi$
is a minimally coupled, complex scalar field, and $V\left(|\Phi|^{2}\right)$
is its associated self-interaction potential. The stress-energy tensor
$T_{ab}$ for the complex scalar field is given by

$T_{ab}=\nabla_{a}\Phi\nabla_{b}\Phi^{*}+\nabla_{a}\Phi^{*}\nabla_{b}\Phi-
g_{ab}\left[\nabla^{c}\Phi\nabla_{c}\Phi^{*}+V\left(|\Phi|^{2}\right)\right],$

where $\Phi^{*}$ is the complex conjugate of $\Phi$. Different BS models are
classified according to their scalar self-potential
$V\left(|\Phi|^{2}\right)$. Here we focus on the solitonic potential frie ,
which allows for highly compact BSs and is given by

$V\left(|\Phi|^{2}\right)=m_{b}^{2}|\Phi|^{2}\left[1-\frac{2|\Phi|^{2}}{\sigma_{0}^{2}}\right]^{2}\,,$
(3)

where $m_{b}$ and $\sigma_{0}$ are two free parameters. In our units (in which
the scalar field is dimensionless), $m_{b}$ has the dimensions of an inverse
length, $m_{b}\hbar$ is the bare mass of the scalar field, whereas
$\sigma_{0}$ is dimensionless. We define $\lambda\equiv\sigma_{0}\sqrt{8\pi}$
and set $m_{b}\lambda=1$ for the rest of the paper. However, in some
occurrences we shall re-insert the proper factors of $m_{b}\lambda$.

In the complex-$\Phi$ space the potential has the typical Mexican-hat shape,
with a maximum at $\Phi=0$ and a minimum (degenerate vacuum) at
$|\Phi|=\sigma_{0}/\sqrt{2}$. When $\sigma_{0}\ll 1$ the scalar profile is
roughly constant within the star and steeply vanishes over a lengthscale $\sim
1/m_{b}$ Kesden:2004qx ; Boskovic:2021nfs .

Due to the $U(1)$ invariance of the EKG action, BSs admit a conserved Noether
charge current

$j^{a}=ig^{ab}(\Phi^{*}\,\nabla_{b}\Phi-\Phi\,\nabla_{b}\Phi^{*})\,.$ (4)

The spatial integral of the time component of this current defines the
conserved Noether charge, $N$, which can be interpreted as the number of
bosonic particles in the star liebpa .

###  II.2 Binary initial data

Figure 1:  Isolated BS solutions. (Top) Compactness as a function of the
central value of the scalar field $\phi_{c}$ for isolated, nonrotating BSs
with $\sigma_{0}=0.05$. Circular markers refer to the equilibrium
configurations used in this work to construct initial data for BS binaries
[cf. Table 1]. The radius $R_{M}$ is defined as that containing $99\%$ of the
mass of the star. (Bottom) Profile of the scalar field as a function of the
isotropic radius for the different configurations.

The procedure to construct the initial data for a binary BS is the same as in
Paper I, that is, a superposition of two boosted, isolated, solitonic BSs.

The solution of a single solitonic BS is constructed as described in Ref. mace
, by adopting the usual harmonic ansatz for the scalar field
$\Phi=\phi(r)\,e^{-i\omega t}$ with a real frequency $\omega$. Assuming
stationarity and spherical symmetry, the EKG equations reduce to a set of
ordinary differential equations which can be solved numerically with a
shooting method. Integrating from the center with a given central value of the
scalar field $\phi_{c}$ and frequency $\omega$, one looks for solutions
satisfying regularity and boundary conditions. The resulting BS equilibrium
configurations can be characterized by their mass and radius. However, because
the scalar field only vanishes asymptotically as it decays exponentially, the
definition of its radius is necessarily somewhat ambiguous. Following previous
work, we can define the effective radius $R_{M}$ as the radius within which
$99\%$ of the total mass is contained, i.e. $m(R_{M})=0.99M$. Consequently, we
define the compactness as $C\equiv M/R_{M}$. As a reference, the compactness
for a Schwarzschild BH is $C=0.5$ and $C\approx 0.1-0.2$ for NSs. In numerical
simulations, it is however more convenient to estimate the radius of the final
remnant through the radius that contains $99\%$ of the Noether charge,
$R_{N}$, so we will use this definition when required. The radius of the
remnant is calculated with respect to its center of mass.

The maximum mass of static configurations in this model is

$M_{\rm max}\approx
5M_{\odot}\left[\frac{10^{-12}}{\sigma_{0}}\right]^{2}\left[\frac{500\,{\rm
GeV}}{m_{b}\hbar}\right]\,,$ (5)

where the scaling with $m_{b}^{-1}$ is exact, whereas the scaling with
$\sigma_{0}^{-2}$ is approximately valid only in the $\sigma_{0}\ll 1$ limit.
Thus, depending on $(m_{b},\sigma_{0})$ the model supports self-gravitating
configurations across a wide mass range.

Paper I presented a sequence of isolated BS solutions characterized by the
central value of the scalar field $\phi_{c}$ that we use to construct our
unequal mass binaries here. In the top panel of Fig. 1, the compactness $C$ is
shown as a function of $\phi_{c}$. The circular markers denote the five
representative BSs employed in this paper. The bottom panel of Fig. 1 displays
the radial profile of the scalar field for these isolated solutions, while
Table 1 lists the key properties of these configurations.

Notice that these solutions can be rewritten in terms of the following
dimensionless quantities mace

$\displaystyle
M(m_{b}\lambda),~{}~{}N(m_{b}\lambda)^{2},~{}~{}r(m_{b}\lambda),~{}~{}\omega/(m_{b}\lambda)\,,$
(6)

recalling that $\lambda=\sigma_{0}\sqrt{8\pi}$. In terms of these parameters,
the equations become independent of $m_{b}$, and hence $m_{b}$ serves to set
the units of the physical solution. Again, the linear scaling in $m_{b}$ in
the above expressions is exact, whereas that with respect to $\lambda$ is
approximately valid only in the $\sigma_{0}\ll 1$ limit. For the chosen value,
$\sigma_{0}=0.05$, this scaling is already a good approximation, and so
smaller values of $\sigma_{0}$ can be studied simply by applying such a
rescaling. Here we restrict ourselves to $\sigma_{0}=0.05$, which sufficiently
fulfills the condition $\sigma_{0}\ll 1$ and also allows for very compact,
stable configurations.

$C$ | $\phi_{c}/(\sigma_{0}/\sqrt{2})$ | $Mm_{b}\lambda$ | $N(m_{b}\lambda)^{2}$ | $(R_{M},R_{N})m_{b}\lambda$ | $\omega/(m_{b}\lambda)$ |  | $I/M^{3}$ | $k_{\rm tidal}$
---|---|---|---|---|---|---|---|---
0.03 | 1.065 | 0.0463 | 0.01653 | (1.507, 1.380) | 2.129620346 |  | $245.3$ | $136494$
0.06 | 1.045 | 0.1238 | 0.0605 | (2.0334, 1.8288) | 1.545745909 |  | $84.9$ | $8420$
0.12 | 1.030 | 0.3650 | 0.2551 | (3.0831, 2.8360) | 1.066612350 |  | $27.8$ | $332$
0.18 | 1.025 | 0.7835 | 0.7193 | (4.2572, 3.9960) | 0.790449025 |  | $12.5$ | $41$
0.22 | 1.025 | 1.0736 | 1.1147 | (4.9647, 4.7068) | 0.685760351 |  | $8.34$ | $20$
Table 1: Characteristics of solitonic BS models with $\sigma_{0}=0.05$. The
table shows: $\bullet$ the compactness $C$, $\bullet$ the central value of the
scalar field $\phi_{c}/(\sigma_{0}/\sqrt{2})$, $\bullet$ the ADM mass
$Mm_{b}\lambda$, $\bullet$ the Noether charge
$N\left(m_{b}\lambda\right)^{2}$, $\bullet$ the radius of the star (the radius
containing $99\%$ of either the mass or of the Noether charge for $R_{M}$ or
$R_{N}$, respectively), and $\bullet$ the angular frequency of the field in
the complex plane, $\omega/\left(m_{b}\lambda\right)$, in dimensionless units.
In the last two columns, we give $\bullet$ the normalized, Newtonian, moment
of inertia (where $I=\int\rho^{2}\,dm$, where $\rho$ is the distance from the
axis of rotation), and $\bullet$ the dimensionless tidal Love number, $k_{\rm
tidal}$, as computed in Refs. Cardoso:2017cfl ; Sennett:2017etc . For a NS
with an ordinary equation of state and $C\sim 0.1$, $k_{\rm tidal}={\cal
O}(200)$ while $k_{\rm tidal}=0$ for a BH.

The initial data for the BS binary follows the procedure described in Ref.
bezpalen  and Paper I. Once the isolated BSs are constructed in spherical
coordinates, the solution is extended to Cartesian coordinates, with the
centers of the stars located at along the $y$-axis at $(0,y^{j}_{c},0)$, so
that the center of mass of the system is located at the origin.222Here we
define the center of mass using the masses of isolated configurations listed
in Table 1. Constraint violation transient will change these masses, see the
discussion of “effective” configurations below. A Lorentz transformation is
performed to boost each star along the $\pm x$-directions, and finally the
boosted solutions for both stars are superposed to obtain our binary initial
data. Obviously, this superposition is only an approximate solution that does
not satisfy exactly the constraints at the initial time (see Ref.
Helfer:2021brt  for a partial solution in case of equal mass binaries of BSs).
However, our evolution scheme enforces an exponential decay of this constraint
violation dynamically (e.g., see Fig. 10 in Ref. bezpalen ).

In contrast with Paper I where the positions and initial velocities of each
binary were anti-symmetric (i.e., velocities with the same magnitude but
opposite direction), for these unequal cases we have set those parameters as
follows: given an initial separation we have calculated the 2nd order post-
Newtonian orbital velocity Mirshekari:2013vb  such that the system would be in
quasicircular orbit and the velocity of the center of mass would be close to
zero. Then, we modify these velocities by adding a tiny amount of linear drift
velocity to account for the finite initial orbital distance and higher-order
relativistic effects, and fix this drift velocity such that the velocity of
the binary center of mass is close to zero. The positions and velocities of
each binary system considered in this work, together with other parameters of
our simulations, are presented in Table 2.

Binaries | $q$ | $\nu$ | $y_{c}^{(1)}$ | $y_{c}^{(2)}$ | $v_{x}^{(1)}$ | $v_{x}^{(2)}$ | $M_{0}m_{b}\lambda$ | $J_{0}(m_{b}\lambda)^{2}$ | $t_{\mathrm{c}}$ | $t^{\mathrm{ret}}_{\mathrm{m}}$ | remnant | $M_{\rm r}m_{b}\lambda$ | $R_{\rm N}m_{b}\lambda$ | $M_{r}\omega^{0}_{\rm r}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
C003 - C022A | 23.2 | 0.039 |  $-$9.58 | 0.42 |  $-$0.34 | 0.02 | 1.16 | 0.229 | 790 | 811 | BS | 1.07 | 4.50 | 0.218
C003 - C022 | 23.2 | 0.039 |  $-$9.58 | 0.42 |  $-$0.34 | 0.02 | 1.16 | 0.229 | 790 | 808 | BS | 1.13 | 4.76 | 0.228
C006 - C022 | 8.6 | 0.093 |  $-$8.96 | 1.04 |  $-$0.36 | 0.05 | 1.34 | 0.668 | 510 | 539 | BS | 1.24 | 5.0 | 0.239
C012 - C022 | 2.9 | 0.189 |  $-$8.95 | 3.05 |  $-$0.33 | 0.136 | 1.90 | 2.388 | 370 | 402 | BH | 1.89 | 3.48 | 0.467
C012 - C018 | 2.1 | 0.21 |  $-$8.18 | 3.81 |  $-$0.26 | 0.135 | 1.36 | 1.488 | 660 | 684 | BS | 1.17 | 4.34 | 0.250
Table 2:  Properties of unequal binary BS models and of the final remnant.
Each case is characterized by: $\bullet$ the compactness $C$ of the individual
BSs in the binary, $\bullet$ the mass-ratio $q$, $\bullet$ the symmetric mass
ratio $\nu$, $\bullet$ the two initial centers $y_{c}^{(i)}$, $\bullet$ the
initial velocities of the boost $v_{x}^{(i)}$, $\bullet$ the ADM mass $M_{0}$
of the system, and $\bullet$ the orbital ADM angular momentum $J_{0}$ of the
system, after the constraint-violating transient respectively. The merger and
remnant are characterized by: $\bullet$ the coordinate time of contact of the
two stars $t_{c}$ (defined as the time at which the individual Noether charge
densities make contact for the first time), $\bullet$ the merger retarded time
(defined as the time when the maximum of the modulus of the $\Psi^{2,2}_{4}$
is produced minus the travel time to the sphere where it is measured),
$\bullet$ the type of final remnant, $\bullet$ the remnant mass
$M_{r}m_{b}\lambda$, $\bullet$ the remnant radius $R_{N}m_{b}\lambda$ (i.e.,
containing $99\%$ of the Noether charge), and $\bullet$ the main GW frequency
$M_{r}\omega_{r}^{0}$ in the post-merger. When the final remnant is a BH, the
radius is computed from the expression for Kerr BHs,
$R_{H}=M_{r}(1+\sqrt{1-a^{2}})$, where $a=J_{r}/M_{r}^{2}\approx 0.5$ is the
dimensionless spin.

As mentioned, our binary initial data is only approximate, but constraint
violations quickly propagate off the grid by our evolution scheme. Hence, it
makes sense to evaluate the global characteristics of the initial data not at
the initial time but instead just after the constraint-violating transient. We
therefore extract numerically the ADM mass, $M_{0}$, of the spacetime after
the transient, and, assuming that the mass ratio remains constant through the
transient, we decompose this mass into the constituent “effective” masses as

$\displaystyle\tilde{M}_{1}=\left(\frac{q}{q+1}\right)M_{0}\,,\quad\tilde{M}_{2}=\left(\frac{1}{q+1}\right)M_{0}\,.$
(7)

Notice that this calculation tacitly assumes that, even after the constraint
violation transient (approximately) ends, stars are sufficiently separated so
that GR nonlinearities are sub-leading. During this transient regime, we note
that the masses of the constituent stars increase which results in a decrease
in the number of orbits.

Furthermore, we can construct fitting formulae for the compactness, $C(M)$,
and particle number, $N(M)$, as functions of BS mass from the equilibrium
configurations of isolated BSs. We obtain

$\displaystyle C(M)$ $\displaystyle\approx$ $\displaystyle
0.0157+0.376M-0.3M^{2}+0.136M^{3}$ (8) $\displaystyle-0.0195M^{4}\,,$
$\displaystyle N(M)$ $\displaystyle\approx$
$\displaystyle-0.0187+0.6221M+0.3872M^{2}\,.$ (9)

With the above functions, one can calculate the “effective” Noether charges
and compactnesses of the stars in our binaries as a function of their
“effective” masses, respectively. In Table 3 , we provide this data for all
configurations consider in this work and Paper I. We also provide the relative
differences between the properties of the isolated initial data and the
“effective” ones. Comparing the total Noether charge in the system, $N_{0}$,
with the sum of the individually calculated charges,
$N(\tilde{M}_{1})+N(\tilde{M}_{2})$, provides a test of the consistency of
this approach. As explained below in Sec. III.2, the “effective” initial data
presented here agrees roughly with our initial data after the constraint-
violating transient.

Binaries | $\tilde{M}_{1}m_{b}\lambda$ | $\Delta M_{1}/\tilde{M}_{1}$ | $\tilde{C}_{1}$ | $\Delta C_{1}/\tilde{C}_{1}$ | $\tilde{N}_{1}(m_{b}\lambda)^{2}$ | $\Delta N_{1}/\tilde{N}_{1}$ | $\tilde{M}_{2}m_{b}\lambda$ | $\Delta M_{2}/\tilde{M}_{2}$ | $\tilde{C}_{2}$ | $\Delta C_{2}/\tilde{C}_{2}$ | $\tilde{N}_{2}(m_{b}\lambda)^{2}$ | $\Delta N_{2}/\tilde{N}_{2}$
---|---|---|---|---|---|---|---|---|---|---|---|---
C006 - C006 | 0.13 | 0.066 | 0.061 | 0.0099 | 0.071 | 0.14 |  |  |  |  |  |
C012 - C012 | 0.43 | 0.15 | 0.13 | 0.093 | 0.32 | 0.21 |  |  |  |  |  |
C018 - C018 | 1.0 | 0.22 | 0.21 | 0.14 | 1.0 | 0.29 |  |  |  |  |  |
C022 - C022 | 1.6 | 0.32 | 0.28 | 0.20 | 1.9 | 0.42 |  |  |  |  |  |
C003 - C022 | 0.048 | 0.034 | 0.033 | 0.093 | 0.012 | 0.38 | 1.1 | 0.035 | 0.22 | 0.0028 | 1.2 | 0.034
C006 - C022 | 0.14 | 0.11 | 0.063 | 0.044 | 0.076 | 0.20 | 1.2 | 0.11 | 0.23 | 0.044 | 1.3 | 0.13
C012 - C022 | 0.49 | 0.25 | 0.14 | 0.16 | 0.38 | 0.32 | 1.4 | 0.24 | 0.25 | 0.14 | 1.6 | 0.32
C012 - C018 | 0.44 | 0.17 | 0.13 | 0.10 | 0.33 | 0.23 | 0.92 | 0.15 | 0.20 | 0.10 | 0.88 | 0.19
Table 3:  Effective properties of the individual stars within the binary after
the constraint violating transient. Tildes represent “effective“ quantities of
the stars in the binary. For the equal mass binaries of Paper I and the
unequal mass binaries studied here: $\bullet$ the mass from Eq. (7), $\bullet$
the compactness $C$ from Eq. (8), $\bullet$ the Noether charge from Eq. (9),
$\bullet$ for each of these, their fractional differences, $\Delta X/X$, with
respect to the initial data for the isolated star used in the construction of
the binary.

###  II.3 Numerical setup and analysis

The computational code, generated by the Simflowny platform ARBONA20132321 ;
ARBONA2018170 ; PALENZUELA2021107675 ; simflownywebpage , runs under the
SAMRAI infrastructure Hornung:2002 ; Gunney:2016 ; samraiwebpage , which
provides parallelization and the adaptive mesh refinement (AMR) required to
resolve the different scales in the problem. We use fourth-order spatial,
finite difference operators to discretize the EKG equations, which are evolved
in time using a fourth-order Runge-Kutta integrator Palenzuela:2018sly .

Our computational domain ranges within $[-264,264]^{3}$ and contains 8 levels
of refinement. Each level has twice the resolution of its coarser parent
level, achieving a resolution of $\Delta x_{8}=0.03125$ on the finest grid. We
use a Courant factor $\lambda_{c}\equiv\Delta t_{l}/\Delta x_{l}=0.4$ on each
refinement level $l$ to ensure the stability of the numerical scheme.

We analyze some relevant global physical quantities from our simulations, such
as the Arnowitt-Deser-Misner (ADM) and the Komar mass, the ADM angular
momentum, and the Noether charge, computed as in Ref. bezpalen . We focus our
attention mainly on the gravitational radiation represented by the strain $h$,
which is the quantity directly observable by GW detectors. We consider first
the Newman-Penrose scalar $\Psi_{4}$, which can be expanded in terms of spin-
weighted $s=-2$ spherical harmonics rezbish ; brugman  as

$r\Psi_{4}(t,r,\theta,\phi)=\sum_{l,m}\Psi_{4}^{l,m}(t,r)\,{}^{-2}Y_{l,m}(\theta,\phi),$
(10)

where the coefficients $\Psi_{4}^{l,m}$ are extracted and calculated on
spherical surfaces at different extraction radii. The relation between this
scalar and the two polarizations of the strain is given by
$\Psi_{4}=\ddot{h}_{+}-i\,\ddot{h}_{\times}$. The components of the strain in
the time domain can be calculated by performing the inverse Fourier transform
of the strain in the frequency domain, $h^{l,m}(t)\equiv{\cal
F}^{-1}[{\tilde{h}}^{l,m}(f)],$ where a high-pass filter has been applied in
the frequency domain in order to attenuate the signal with frequencies lower
than the initial orbital frequency 2011CQGra..28s5015R ; bezpalen . The
instantaneous angular frequency of each GW mode can be calculated easily from
$\Psi_{4}$ as

$\displaystyle\omega^{l,m}_{\rm{GW}}=-\frac{1}{m}\mathfrak{Im}\left(\frac{\dot{\Psi}_{4}^{l,m}}{\Psi_{4}^{l,m}}\right)~{},~{}~{}~{}f^{l,m}_{\rm{GW}}=\frac{\omega^{l,m}_{\rm{GW}}}{2\pi}.$
(11)

We will refer to $\omega_{\rm{GW}}$ as the one given by the dominant mode
$l=m=2$.

The mass, the angular momentum, and $\Psi_{4}$ are calculated on spherical
surfaces at different extraction radii between $R_{\rm{ext}}=100$ and
$R_{\rm{ext}}=200$, which are located far away from the sources in the wave
zone.

##  III Dynamics for unequal-mass BS binaries

We have evolved four unequal mass binary BS cases, $\\{$C003-C022, C006-C022,
C012-C022, C012-C018$\\}$, covering mass ratios $q\equiv m_{1}/m_{2}$ roughly
between $2$ and $23$. Additionally, we have studied a variation of the most
extreme case, C003-C022A, in which the heavier BS has been transformed into an
anti-BS. In what follows, we describe first qualitatively the dynamics for all
the cases and then analyze the GWs produced by these mergers in the next
section.

###  III.1 Binary dynamics in the inspiral

We display some representative snapshots along the equatorial plane to
characterize the dynamics of these binary evolutions. In particular, the
Noether charge densities in Fig. 2 show the dynamics of the condensed bosons,
whereas the scalar field norm in Fig. 3 shows the dynamics of the scalar field
generally.

The binaries in C003-C022A and C003-C022 complete five full orbits before
colliding, C006-C022 and C012-C018 complete three orbits, and C012-C022
performs just two. While such a short inspiral limits their use for guiding
templates, the inspiral is long enough for constraint violations resulting
from the construction of the initial data to propagate away.

Figure 2:  Dynamics of the Noether charge. Noether charge densities in the
equatorial plane $(z=0)$ at representative times. Each row represents one of
the cases (from top to bottom): $\\{$C003-C022A, C003-C022, C006-C022,
C012-C022, C012-C018$\\}$. The first column illustrates a time roughly one
orbit before the contact time $t_{c}$ (defined as the time at which the
individual Noether charge densities make contact for the first time), the
second column occurs at contact time, the third is roughly an orbit after the
contact time (except for the C012-C018 case, in order to visualize the ejected
blob), and the fourth one occurs at the end of our simulations. Note that the
final remnant for C012-C022 is a rotating BH which quickly engulfs the
surrounding scalar field (i.e., not visible on this natural scale). Figure 3:
Dynamics of the scalar field. Logarithm of the modulus of the scalar field
$||\Phi||$ in the equatorial plane $(z=0)$, at representative times. Each row
represents the unequal cases considered. Notice that there is only some scalar
emission around the contact time $t_{c}$ (defined as the time at which the
individual Noether charge densities make contact for the first time),
suggesting that the final object is an almost stationary BS (except for the
simulation on the third row, in which the remnant is a spinning BH).

During the inspiral, the spacetime curvature is dominated mainly by the
heavier BS, which moves in a spiral trajectory very close to the origin (i.e.,
see the leftmost column of Fig. 2), while the lighter one induces a
perturbation orbiting around the most massive object. This effect is
especially pronounced in the four most unequal mass cases in which the heavier
BS accounts for at least $75\%$ of the binary mass. During the inspiral, the
scalar field constituting each star has no significant overlap (see the first
column of Fig. 3), and therefore nonlinear scalar interactions only play a
significant role inside the stars. Roughly speaking, the BSs behave then like
point particles with moderate deviations produced by the tidal deformations.
As the mass ratio approaches unity, the binary behaves similarly to the equal-
mass cases of Paper 1. In particular, C012-C018 with $q=2.1$ resembles those
equal-mass cases.

The aforementioned deviations due to tidal deformations can be estimated by
looking at the quadrupole-moment tensor $Q_{ab}^{(i)}$ of the $i$-th object
induced by the tidal-field tensor $G_{ab}^{(j)}$ produced by the $j$-th object
($i,j=1,2$) Flanagan:2007ix ; PoissonWill ,

$Q_{ab}^{(i)}=\lambda_{i}G_{ab}^{(j)}\sim\lambda_{i}\frac{m_{j}}{r^{3}}\,,$
(12)

where $r$ is the orbital distance and
$\lambda_{i}=\frac{2}{3}m_{i}^{5}k^{(i)}_{\rm tidal}$ is the tidal Love number
of the $i$-th object, with $k^{(i)}_{\rm tidal}$ being its dimensionless
counterpart. Hence, the dimensionless quadrupole moment,
$\bar{Q}_{i}=|Q_{ab}^{(i)}|/m_{i}^{3}$, reads

$\displaystyle\bar{Q}_{1}$ $\displaystyle\sim k^{(1)}_{\rm
tidal}\frac{q^{2}}{(1+q)^{3}}\frac{M_{0}^{3}}{r^{3}}\,,$ (13)
$\displaystyle\bar{Q}_{2}$ $\displaystyle\sim k^{(2)}_{\rm
tidal}\frac{q}{(1+q)^{3}}\frac{M_{0}^{3}}{r^{3}}\,,$ (14)

where $M_{0}=m_{1}+m_{2}$ is the binary total mass. In the large mass-ratio
limit, $q\gg 1$, the tidally-induced quadrupole moments of the primary and of
the secondary are suppressed by a factor $q^{-1}$ and $q^{-2}$, respectively.
For example, for a fixed value of $k^{(i)}_{\rm tidal},$ the tidally-induced
quadrupole moment of the primary for $q=23$ is suppressed by a factor $3$
relative to $q=1$, whereas that of the secondary is even a factor $\sim 100$
smaller. Overall, tidal effects on the secondary object are less relevant than
those on the primary.

###  III.2 Final fate of the binary merger

If the system is sufficiently massive such that the remaining mass after
merger exceeds the maximum stable BS mass (i.e., $M_{r}\geq M_{\rm max}\approx
1.85/(m_{b}\lambda)$), one expects the system to collapse to a remnant BH. If
instead the total mass is below this threshold, a remnant BS is expected. In
the latter case, the possibility of forming a rotating BS should be
considered. At least two conditions appear to be required for such formation:
(i) because rotating BSs have quantized angular momentum, binaries need to
have angular momentum at the point of contact at least slightly larger than or
equal to the first discrete level of the rotating star, 333This argument
excludes some exotic possibility in which, say, GWs with some opposite angular
momentum are radiated copiously until the remnant achieves the sufficient
amount of angular momentum. and (ii) the rotating solution to which the
remnant might settle must be stable.

Once the stars contact each other, one expects scalar field interactions to
produce additional attractive forces that accelerate the merger (see the
discussion of the effective force with just a massive potential in Appendix B
of pale1 ). The newly formed, rotating object is initially largely
nonaxisymmetric, and, even by the end of our simulations, the remnant is a
highly perturbed BS (see the rightmost column of Fig. 2). Some general
features of the dynamics can be found in certain global quantities (mass,
Noether charge, and angular momentum) which are displayed in Fig. 4.

Figure 4:  Global quantities. ADM mass (top) , angular momentum $J_{\rm{z}}$
(middle) and Noether charge (bottom) as functions of time. The time has been
shifted such that contact time happens at $t=0$ and rescaled with the initial
total mass $M_{0}$ of each binary. Horizontal dashed lines in the middle panel
indicate the angular momentum of the ground state rotating BS corresponding to
the Noether charge measured at the contact time. Notice that the angular
momentum decreases monotonically (although slowly) after merger for all
remnants except for that of C012-C022 which collapses to a BH. This decrease
in angular momentum along with the fact that these binaries (except C012-C022
which collapses to a BH) have less angular momentum than any rotating BS with
the same Noether charge support our claim that all non-collapsed cases settle
to a nonrotating BS. The non-monotonic, brief drops in the mass and angular
momentum plots for the C003-C022A and C012-C018 cases correspond to the
passage of transients across the extraction surface on which mass and angular
momentum are calculated (the retarded time is used). The Noether charge is
computed as a volume integral and therefore less subject to such errors.

The mass and the Noether charge are unambiguously defined global quantities,
in contrast to the radius of the star. In the case of a complex field, the
$U(1)$ symmetry, which ensures the conservation of the Noether charge,
significantly restricts the ways in which the remnants might relax. Fig. 5
shows the mass-Noether charge phase space for two representative cases
C006-C012 and C012-C018. Here, we present several estimates of the initial and
final data along with families of isolated BSs, to facilitate the
understanding of the relaxation of the remnant.

The orange squares indicate the simplest estimate of the initial data,
$(M_{1}+M_{2},N_{1}+N_{2})$, obtained by adding the properties of the isolated
BSs used to construct the binary. These two squares fall far from our two
other estimates of the initial data. In particular, the total mass and Noether
charge measured by the numerics after the transient is shown in black circles.
We then construct the “effective” initial data (red crosses) by decomposing
the numerically obtained total mass via Eq. (7) and computing the charge of
each BSs from these individual masses (with Eq. (9) in Sec. II.2).

We further note that, due to the nonlinearity of the function $N(M)$, some
amount of scalar and/or GW emission is needed during the merger in order for
the remnant to settle into either a static or rotating configuration. If the
remnant is assumed to be a BS that relaxes only by the emission of GWs, namely
no emission of scalar field to infinity, the evolutionary path of the binary
would follow a horizontal line in the $N$-$M$ phase space (blue dashed line on
Fig. 5), ultimately settling into the remnant BS occurring at the intersection
with the family of nonrotating BSs given by Eq. (9) (red dotted line). Our
simulations indicate emission of scalar field, in addition to GWs, a process
known as “gravitational cooling” Seidel:1993zk ; Guzman:2006yc . Indeed, the
path of the numerical evolution (green dots) indicates that the dynamics are
driving each system toward a stationary BS (red-dotted line). Although most of
these BS mergers ended before the remnant fully relaxed to stationarity, we
have established for C003-C022 and C012-C012 that the point $(M_{r},N_{r})$
(where $N_{r}\equiv N(t_{\mathrm{end}})$) indeed lies on the isolated BS
curve. However, the near constancy of the Noether charge in the late
postmerger (Fig. 4) and the close approach of the final simulation to the
isolated BS curve (Fig. 5) both indicate that the mergers that do not collapse
are forming a stable, nonrotating, solitonic BS.

If the late stage evolution is dominated by GW emission (since most of the
ambient scalar field has already been radiated), then we would expect the
final object to be that represented by the black solid diamond,
$(M_{r},N(M_{r}))$, itself a stationary BS, because the Noether charge would
not be changing.

---
Figure 5: Mass-Noether charge phase space. Evolution of the total mass and
charge of the system for two representative cases: C006-C012 (left) and
C012-C018 (right) which ejects a blob of scalar field after merger. The orange
square for each case represents the value of $(M_{1}+M_{2},N_{1}+N_{2})$
obtained from the individual stars as listed in Table 1. The black, open
circles are the masses and charges, $(M_{0},N_{0})$, obtained from the
numerical evolution just after the transient. The red crosses are the
“effective” values
$(\tilde{M}_{1}+\tilde{M}_{2},N(\tilde{M_{1}})+N(\tilde{M_{2}}))$, as
explained in Sec. II.2. The fact that the red crosses and black circles are
close to each other supports the validity of this approach. The green dots
trace the numerical evolution by showing the extracted quantities
$(M(t),N(t))$. The values characterizing the final state, $(M_{r},N_{r})$, of
the simulation are represented by a green triangle. The black diamond is the
point, $(M_{r},N(M_{r}))$, with the same mass as the green triangle but with
the charge obtained from the fit in Eq. (9). If one assumes that the remnant
is a BS that relaxes only via emission of gravitational waves, one obtains a
horizontal phase space trajectory (blue dashed line) through the initial data
(namely the black circle here). The family of nonrotating BSs given by Eq. (9)
are plotted with a (red dotted curve). Only the case on the right has angular
momentum comparable to the first rotating solution, and so for this case we
also show the family of $k=1$ rotating BS configurations for our same
potential from Ref. 2021PhRvD.103d4022S  with a solid green curve. Because
$J=kN$ for such rotating BSs, we also show the value of $N$ corresponding to
the angular momentum of the binary at contact time with a horizontal, yellow
dot-dashed line.
---
Figure 6: Scenario classification. The parameter space of solutions
represented in two different ways: (left) mass ratio versus total mass, and
(right) individual mass of one star versus the other. The outcomes of our
simulations are consistent with $M_{\max}$ being the simple delineator for the
BS/BH nature of the remnant. The blue region encloses configurations that
satisfy the angular momentum requirement of Eq. (15) (informed by our
evolutions that produce blobs), and thus we expect either blob ejection or the
formation of a rotating boson star. The red hashed region is a subset of this
region which we have not sampled, but where rotating BSs may form.

An important unresolved question is whether a merger of two BSs can produce a
rotating BS. The stability of rotating, solitonic BSs has been studied
recently. First, rotating BSs without scalar self-interactions were found to
be unstable due to a non-axisymmetric instability 2019PhRvL.123v1101S .
However, a subsequent study showed that this instability was quenched for the
solitonic model of the potential 2021PhRvD.103d4022S  (see also Ref.
Dmitriev:2021utv ) if $M>0.13/(m_{b}\lambda)$, for the value $\sigma_{0}=0.05$
considered here. Without stability, one would not expect formation of such
configurations from a merger.

Rotating BSs have quantized angular momentum, $J=kN$ for some integer $k$, and
one can calculate the function $N(M)$ for the $k=1$ family of rotating BSs
following Ref. 2021PhRvD.103d4022S  (see also Kleihaus:2005me ). We display
this family of solutions as a green solid curve in the right panel of Fig. 5
(case C012-C018) because this binary has angular momentum close to this first
quantized level. Actually, only two cases among those studied in this work and
Paper I (i.e., C012-C018 and C012-C012) are close to satisfying the
quantization condition, namely that the angular momentum is greater than or
equal to the Noether charge at the time of contact. In neither of these two
cases do we find a rotating remnant, and the angular momentum is primarily
reduced through emission of scalar “blobs.”

The case C012-C018 is shown in the right panel of Fig. 5. We display the
Noether charge equal to the binary’s angular momentum at the time of contact
with the horizontal, yellow dot-dashed line. However, as shown in the figure,
the point of intersection of the dynamical path of the binary, $(M(t),N(t))$
shown in green dots, with the curve indicating the $k=1$ family of rotating
BSs (solid green curve) occurs above this yellow line. Because these rotating
solutions have angular momentum equal to their charge, the evolution lacks
sufficient angular momentum to form the rotating BS indicated by this point of
intersection.

In an effort to understand the configuration space of binaries in terms of
possible endstates, in particular including formation of a rotating remnant or
a blob, we parameterize the quantization condition. We first compute a
Keplerian estimate of the angular momentum either at the time of first contact
$R_{c}=C_{1}/M_{1}+C_{2}/M_{2}$ or when the binary reaches the innermost
stable circular orbit, $R_{\rm ISCO}=6M_{0}$, whichever occurs first. We then
correct this estimate by including the relativistic effects of strong gravity.
Due to the pre-contact scalar emission, the total Noether charge in the binary
at the point of contact will be slightly smaller than the initial one. In
addition, we have observed blob emission in both cases where the total charge
of the binary is slightly higher than $J_{c}$. We incorporate these two
effects in our quantization condition for rotating boson stars, $J\geq N$, by
introducing two new parameters $\\{e_{N}\,,e_{J}\\}$ in the following way

$\frac{J_{c,K}(1+e_{J})}{N(M_{1})+N(M_{2})}>1+e_{N},$ (15)

where $J_{c,K}$ is the Keplerian estimate of the angular momentum at the
contact time, $e_{N}$ estimates either the amount of Noether charge radiated
during the merger $(e_{N}>0)$ or the difference between the critical angular
momentum and the charge at the point of contact that allows for blob emission
$(e_{N}<0)$. Finally, $e_{J}$ accounts for general relativistic corrections to
the Keplerian angular momentum calculation.

We use the above cases to estimate the value for the parameters $e_{J},e_{N}$.
To estimate $e_{J}$, we compute the differences between the Keplerian estimate
of the angular momentum at contact time and the numerical value, obtaining
$\sim 25\%$ in scenarios where we observe blob formation: C012-C012 and
C012-C018. In the low-mass regime, where solitonic BSs are in the weak-field
regime, we expect that $e_{J}\to 0$. Thus, we linearly interpolate $e_{J}$
between 0 and 0.25 for $M_{0}/2\in[M_{\rm min},M_{\rm C012}]$ and take the
constant value $e_{J}=0.25$ up to $M_{0}=M_{\rm max}$. Due to the initial data
constraint violation, we cannot estimate reliably how much of the Noether
charge is emitted before contact. For the sake of argument, we take
$e_{N}=0.01$ in this case, delineating a subset of the parameter space where
the strict form of the quantization condition is satisfied and where rotating
remnants may form. In addition, we require that the remnant has surpassed the
threshold mass estimate from Ref. 2021PhRvD.103d4022S .

In the two cases where the blobs are observed, one finds
$e_{N}=\\{-0.05,-0.02\\}$ for C012-C018 and C012-C012, respectively. Thus,
taking $e_{N}=-0.05$ would encompass both scenarios where blobs are found and
indicate the part of the parameter space where one can expect blobs
generically and even possibly rotating remnants (more restrictive condition).
We sketch the configuration space for these mergers in Fig. 6 in two ways: the
left panel plots the mass ratio versus the total mass, $(q,M)$, whereas the
right panel shows the space spanned by individual masses
$(M_{1}m_{b}\lambda,M_{2}m_{b}\lambda)$. Solutions exist only for binaries
constructed with stable BSs, $M_{i}<M_{\mathrm{max}}$, with regions outside
this indicated in white. For binaries with $M_{1}+M_{2}<M_{\mathrm{max}},$ the
formation of a rotating BS appears possible for the binaries that do not
collapse to a BH and possess angular momentum satisfying Eq. (15), although we
have not observed such formation (red hashed region).444Because the maximal
mass of rotating BSs is larger than that for non-rotating BSs, a priori, even
binaries with total mass slightly higher than the maximum mass for static
stars, $M_{\mathrm{max}}$, could allow for the formation of a rotating
remnant. Note, however, that the effective mass of C012-C022 is slightly
larger than the static maximum mass $M_{\mathrm{max}}$ and the configuration
collapses to BH. Whether this also happens for $q\to 1$ is an open question.
The set where we expect blob emission based on the results of C012-C018 and
C012-C012 cases (blue region) has a red hashed region as its subset. Note that
lacking an understanding of the physics of the blob formation, the blue region
should serve only for illustrative purposes. Both of these regions are
determined approximately and require more simulations in order to understand
their precise extent.

One expects qualitatively similar behavior near $M_{0}\approx
M_{\mathrm{max}}$ in the small $\lambda$ regime ($\lambda\ll 1$). In contrast,
when $M_{0}\to M_{\rm min}$ BSs behave as thick-walled Q-balls (where
“Q-balls” Coleman:1985ki  refers to the flatspace limit of solitonic BSs)
Boskovic:2021nfs , we can study the quantization condition (15) in detail in
this regime. We consider an equal-mass ($q=1$) binary with $N\approx\lambda M$
in which the objects collide at $R_{c}$ (for $q=1$ this happens when
$C<0.17$). Taking $e_{J}\approx e_{N}\approx 0$ (in $m_{b}=\lambda^{-1}$
units) and setting the angular velocity to the Keplerian estimate, it can be
shown with some algebra that Eq. (15) becomes

$C<\frac{M^{2}}{4\lambda^{2}}\,.$ (16)

Thus, for sufficiently small $\lambda$ (approximately an order of magnitude
smaller than the value in this work $\lambda=0.25$), the quantization
condition will be satisfied. This simple expression does not change
parametrically when a more precise description of the Q-balls is used
Boskovic:2021nfs . Although rotating Q-ball solutions have been constructed
Kleihaus:2005me , the non-axisymmetric instability (NAI) probably prevents one
from dynamically forming, based on the results of Ref. 2021PhRvD.103d4022S .
Whether in those cases blobs form or the non-axisymmetric instability would
kick in is an open question.

To conclude, we cannot rule out the formation of a rotating BS with the
solitonic potential although none has been formed. In any case, our parameter
space analysis indicates that the initial conditions would need significant
tuning, which may require more accurate initial data. Even in those cases
where the formation of rotating BS might be feasible, as suggested in Paper I,
the organization of the bosonic field into a rotating star from the very
nonlinear merger may be too difficult, particularly because the rotating BS
necessarily has a toroidal energy density 555Rotating Proca stars instead have
a spheroidal energy density and yet none of these have been formed from a
merger either PhysRevD.99.024017 . .

###  III.3 Remnants: Scalar clouds, blobs, and kicks

Another interesting possibility is the formation of a stable scalar cloud
surrounding the remnant rotating BH. A necessary condition for the scalar
field to remain around a spinning BH is the saturation of the superradiant
condition Brito:2015oca . In particular, this condition is saturated when the
phase oscillation frequency of the scalar field, $\omega$, is synchronized
with the angular frequency of the BH, $\Omega_{H}$, such that
$\Omega_{H}=\omega/m$ for some integer $m$. Synchronized scalar clouds were
not found originally from the mergers of Proca stars PhysRevD.99.024017 , but
more recent and detailed equal-mass binary simulations of non-solitonic
bosonic stars showed long-lived, scalar hair around a rotating horizon for a
small range of the initial angular momentum Sanchis-Gual:2020mzb .

We examine the case C012-C022 to determine whether any scalar field remains
after the remnant has collapsed to a BH. Visualizing the scalar field and its
associated Noether charged density reveals no significant remaining scalar
field (see the fourth column of the the C012-C022 case of Figs. 2 and 3).
Furthermore, one sees the total Noether charge drop quickly to zero after
merger in Fig. 4 and the bottom panel of Fig. 7.

To evaluate the possibility of formation of a synchronized scalar cloud, we
calculate the oscillation frequency of the scalar field for a numerical
comparison of the synchronization condition $\Omega_{H}=\omega/m$. We note
first that the final BH has mass $M=1.89/(m_{b}\lambda)$ and angular momentum
$J_{z}=1.92/(m_{b}\lambda)^{2}$, leading to a dimensionless spin
$a=J_{z}/M^{2}=0.537$. The radius and angular frequency of the BH can be
computed from expressions for Kerr-Schild BHs as
$R_{H}=M(1+\sqrt{1-a^{2}})\approx 3.5/(m_{b}\lambda)$ and
$\Omega_{H}=a/(2R_{H})\approx 0.077m_{b}\lambda$. We Fourier transform the
scalar field at an arbitrary point outside the BH but where the scalar field
is well above any numerical noise (roughly a distance of $12/(m_{b}\lambda)$
from the BH, as in the top panel of Fig. 7). We find a frequency
$\omega\approx 0.6m_{b}\lambda$ which, with the synchronization condition and
the estimate of the BH rotation $\Omega_{H}$, implies an azimuthal quantum
number $m\geq 8$. Interestingly, as shown in the top panel of Fig. 7, the
real, and similarly the imaginary (not shown), components of the late-time
scalar field configuration outside the BH resemble the high $m$-number
structure of a stationary cloud. However, the amplitude of the scalar field is
decreasing fast, consistent with the decrease in both the Komar mass and
Noether charge, displayed in the bottom panel of the same figure.

Previous studies suggest that initial data might need to be fine tuned in
order to form a stationary configuration, unless such a configuration is a
dynamical attractor as in the case of the superradiant instability
Brito:2015oca . It is worth noticing that Ref. Sanchis-Gual:2020mzb  found
Proca clouds with $m$ as high as $6$, but in the vector case the superradiant
instability develops much faster than in the scalar case at hand. In the small
$m_{b}M$ limit, the instability time scale of scalar fields is longer than
that of vectors by a factor $(m_{b}M)^{-2}\gg 1$ Brito:2015oca . Furthermore,
the instability is very suppressed for large azimuthal number $m$. The
imaginary part of the fundamental $(l,m)$ mode for a perturbation with spin
$s$ ($s=0,1$ for scalars and vectors, respectively) depends on the factor
$\beta_{ls}=\left[\frac{(l-s)!(l+s)!}{(2l)!(2l+1)!!}\right]^{2}$ Starobinskij2
; Brito:2015oca . This dependence is responsible for a suppression of the
superradiant instability time scale when $m\gg 1$. Indeed, for the most
relevant $l=m\gg 1$, and focusing only on the $m$-dependence, the instability
time scale reads

$\displaystyle\tau=\omega_{I}^{-1}\sim m^{2m}(m_{b}M)^{-4m-5+2s}\,,$ (17)

which quickly becomes extremely large as $m$ increases. Note that the fastest
growing superradiantly unstable mode in the scalar case ($l=m=1$) has an
instability time scale at least $\tau\sim 10^{5}M$ for a nearly extremal BH.
This is already much longer than the time scale of our simulations, and it
becomes much longer for $m\gg 1$ modes and moderately spinning BHs. This
discussion strongly suggests that it is unlikely that an $m=8$ superradiant
cloud could form dynamically over the limited time scale of our simulations.

It might be possible that mergers leading to smaller oscillation frequencies
of the remnant scalar field or with higher initial angular momentum (so that
the final BH is rapidly rotating) are more likely to produce clouds. Either of
these conditions would lead to a smaller required $m$, but may limit the
parameter space of cloud-generating solutions. Clearly, more work is needed to
answer this question.

Figure 7:  Details of the late-time behavior of case C012-C022. (Top) Snapshot
of the real part of the scalar field at the latest time. Some scalar field
remains around the BH (whose apparent horizon is represented by the black
sphere located at the center of the plot) with a large $m$ configuration. The
amplitude of the scalar field decays quickly during the timescale of our
simulations and does not reach the stationary state expected for a single-mode
synchronized scalar cloud. (Bottom) Noether charge and Komar mass (formulated
only in terms of the stress-energy tensor and hence measures only the scalar
energy and not that of the BH) as a function of time for the same case. After
the sudden drop at the merger, both quantities decay exponentially as the
scalar field falls into the BH.

We now consider the ejection of scalar blobs. As previously explained, the
case C012-C018 is the only one with contact angular momentum close to that of
the first quantized spinning BS configuration (namely, $J_{z}\gtrsim N$) that
does not collapse to a BH. Instead, whether the spheroidal energy density
formed in the merger somehow prevents the configuration from relaxing to the
toroidal shape of the rotating BS or not, the system relaxes instead to a
nonrotating BS. To do so, the system must shed its angular momentum.

In this case, the excess angular momentum is emitted in the form of a blob of
scalar field that is ejected from the remnant soon after the merger (see the
bottom row of Figs. 2 and 3). This blob travels outward on the grid, and its
passage across the spherical surface (i.e., around $(t-t_{c})/M_{0}\approx
100$) at which the system mass and angular momentum are computed disrupts the
assumptions of the calculation, seen as non-monotonicity in the global
quantities shown in Fig. 4.

Using the values before and after the drop in mass, we can estimate the blob’s
mass as $M_{\mathrm{blob}}\approx 0.12$. Despite the blob containing only a
small fraction of the total mass, it carries a significant fraction of the
total angular momentum due to its large velocity, $v_{\mathrm{blob}}\approx
0.5$ directed nearly tangentially away from the remnant and its distance from
the center of mass when ejected, $L\approx 7$. Indeed, using the same simple
estimate for the angular momentum as in Paper I, we obtain $J_{z}\approx
M_{\mathrm{blob}}v_{\mathrm{blob}}L\approx 0.4$, which is roughly equal to the
sharp decrease of angular momentum observed in the middle panel of Fig. 4. On
the time scale of our simulation, the blob appears bounded. In fact, the blob
satisfies the stability condition (in $m_{b}=\lambda^{-1}$ units) $\lambda
M_{\mathrm{blob}}<N_{\mathrm{blob}},$ with $N_{\mathrm{blob}}\approx
0.05/(m_{b}\lambda)^{2}$ (see Ref. Boskovic:2021nfs  and references therein
for a discussion of the stability regimes of solitonic BSs).

In addition to the unequal mass case C012-C018 presented here, the ejection of
condensed scalar field was observed in two equal-mass BS binary simulations,
one in Ref. bezpalen  and the other in Paper I. In those two cases, the
symmetry of the binary resulted in two, identical blobs propagating along
opposite directions. That three different studies found blob ejection suggests
that such ejection might be typical in solitonic BS binaries under certain
conditions.

The ejection of the blobs has important implications for the astrophysics of
BS mergers should such systems (or similar systems such as axion stars)
actually exist in nature. In contrast to the equal-mass case that ejects two
blobs in opposite directions, the ejection of a single blob generates a kick
on the remnant. For binaries with large enough mass ratios, the kick can be
large; large even compared to the superkicks of binary BHs (which are as large
as a ${\rm few}\times 10^{3}\,{\rm km/s}$ Boyle:2007sz ; Campanelli:2007cga ;
Gonzalez:2007hi ; Lousto:2011kp ) and larger than the typical escape
velocities of galaxies and of globular clusters (which are of ${\cal
O}(10^{2}-10^{3})\,{\rm km/s}$ and ${\cal O}(10)\,{\rm km/s}$, respectively).
For example, the linear momentum of the blob shown in the C012-C018 simulation
is roughly $M_{\rm blob}v_{\rm blob}\approx 0.06/(m_{b}\lambda)$ which, by
linear-momentum conservation, implies that the remnant with mass
$M_{r}=1.17/(m_{b}\lambda)$ recoils with a velocity $v_{r}\approx 0.05\sim
1.4\times 10^{4}\,{\rm km/s}$ . In practice, since the ejected scalar blobs
have a sizeable mass and relativistic speed, they induce remnant kicks much
more efficiently than GW emission in asymmetric binary BH systems Boyle:2007sz
; Campanelli:2007cga ; Gonzalez:2007hi ; Lousto:2011kp . These large kicks
would have important implications for the merger rate of BS binaries in the
universe, as they largely exceed the escape velocity from bound structures
(e.g. nuclear star clusters Antonini:2016gqe  and galaxies Bogdanovic:2007hp
). As a result, the rate of successive generations of mergers (which is
particularly important for supermassive objects, see e.g. Ref.
Barausse:2020gbp ) may be suppressed relative to the BH case. Moreover,
“‘stray” BSs moving at high speeds may be present in the intergalactic medium
as a result of ejections from the host galaxies.

Finally, one might be tempted to associate this disruption and blob ejection
to the nonaxisymmetric instability present in some rotating BSs
2019PhRvL.123v1101S . However, a recent study shows that the NAI should be
quenched for the solitonic potential for sufficiently compact BSs
2021PhRvD.103d4022S , which suggests that the NAI is not the cause of blob
ejection.

###  III.4 Collision of boson and anti-boson stars

The use of BSs as a model of compact objects allows easily for the study of
anti-stars Dupourque:2021poh . Here, the scalar field represents its own anti-
particle simply by changing the sign of its phase oscillation,
$\omega\rightarrow-\omega$ (this transformation also switches the sign of the
Noether charge associated with the BS). The case C003-C022A represents a BS
with a small “antimatter star” in a merging binary. The evolutions of
C003-C022A and C003-C022 differ only once the stars make contact, at which
point the antistar annihilates. In particular, the interaction of the
oppositely oscillating complex field annihilates such that the remaining field
lacks the harmonic oscillation in time and the Noether charge density
vanishes. Therefore, this scalar field interaction breaks the coherence of the
BS solution, the dispersive nature of the scalar field dominates over the
attraction of gravity, and the unbound scalar field is radiated to infinity.

The small difference between the final and initial mass of case C003-C022A,
$\Delta M\approx 0.09$, accounts roughly for twice the value of the lightest
BS, $M=0.0463$. This suggests that other energies, such as that of GW
radiation and binding energy, change little or are otherwise very small. This
simple calculation is also consistent with the lightest star being completely
annihilated during the merger. It is interesting, and a sign of the stability
of the solutions under a large perturbation, that the remnant still settles
into a stable BS, even after the annihilation of a significant fraction of its
Noether charge.

##  IV Gravitational Wave Signal

We now turn our attention to the analysis of the gravitational radiation
produced by unequal-mass BS binaries.

###  IV.1 Late inspiral and merger

Some of the most relevant $(l,m)$ modes of the gravitational radiation
represented by the strain, together with the angular frequency of the $(2,2)$
mode, are displayed in Fig. 8. A simple inspection of these profiles already
confirms that the dominant mode during the inspiral is always the $l=m=2$ for
our wide range of mass ratios. As expected, mass ratios closer to unity (i.e.,
such as the C012-C018 case), when the mass quadrupole moment is stronger,
displays a larger predominance of the $l=m=2$ mode. On the other hand, for
large mass ratios (i.e., such as the C003-C022 case), the importance of the
higher-order modes increases. It is interesting to note that after the merger
the amplitudes of the various modes are of the same order, without one clearly
dominating over the others.

Furthermore, as the mass ratio increases, the effects of tidal deformations on
the waveform become less relevant. This can be understood as follows. A
generic quadrupole-moment tensor, $Q_{ab}^{(i)}$, of the $i$-th object affects
the GW phase starting at second post-Newtonian order. The extra $1/r^{3}$
dependence of the tidally-induced quadrupole moment [see Eq. (12)] implies
that tidal effects enter the GW phase starting at the fifth post-Newtonian
order, with a phase correction Flanagan:2007ix ; PoissonWill

$\delta\phi_{\rm
tidal}=-\frac{117}{8}\frac{(1+q)^{2}}{q}\frac{\Lambda}{M_{0}^{5}}v^{5}~{},$
(18)

where $v=(\pi M_{0}f)^{1/3}$ is the orbital velocity, $f$ is the GW frequency,
and $\Lambda=\frac{1}{26}((1+12/q)\lambda_{1}+(1+12q)\lambda_{2})$ is the
weighted tidal deformability. When $q=1$,
$\Lambda=(\lambda_{1}+\lambda_{2})/2$ is simply the average of the two tidal
deformability parameters. However, in the large mass-ratio limit Pani:2019cyc
, we can write the correction as

$\delta\phi_{\rm
tidal}=-k_{1}\left(\frac{3}{8}v^{5}q+\ldots\right)-k_{2}\left(\frac{9}{2}v^{5}\frac{1}{q^{3}}+\ldots\right)\,,$
(19)

where we include for each of the tidal terms $k_{1}$ and $k_{2}$ only the
leading-order term in the $q\gg 1$ expansion. The above equation shows that
the tidal deformability of the primary is much more important than that of the
secondary, which is suppressed by a relative factor $\sim q^{-4}$. Thus, the
net contribution of the tidal deformability in the GW phase, compared to the
point-particle phase, depends on two competing effects: on the one hand, less
compact BSs have a large tidal Love number (see Table 1) but, on the other
hand, for binaries with very disparate mass stars the tidal Love number of the
secondary is negligible. The quantity $\Lambda/M_{0}^{5}$, which provides a
measure of the relevance of the tidal contribution compared to the leading-
order point particle phase, is presented in Table 4. For the binary systems
under consideration, the suppressing effect of large mass ratio more than
compensate for the large tidal Love number of the secondary, and hence the
quantity $\Lambda/M_{0}^{5}$ is larger for the smallest mass-ratio system in
the catalog.

Binaries | $q$ | $k_{1}$ | $k_{2}$ | $\Lambda/M_{0}^{5}$
---|---|---|---|---
C003 - C022A | 23.2 | 20 | 136494 | 0.75
C003 - C022 | 23.2 | 20 | 136494 | 0.75
C006 - C022 | 8.6 | 20 | 8420 | 0.99
C012 - C022 | 2.9 | 20 | 332 | 0.94
C012 - C018 | 2.1 | 41 | 332 | 1.79
Table 4:  Tidal properties of unequal binary BS models considered in our
simulations. For each binary studied here, the weighted tidal deformability
$\Lambda$ is shown. Despite the large tidal Love numbers in the most unequal
mass binaries, the binary deformability increases as $q\rightarrow 1$.

Let us now consider the waveform’s higher modes. We can characterize the
effect of the mass ratio on the higher modes by examining the ratio between
some relevant modes $h^{l,m}$ and the dominant $h^{2,2}$ mode. The top panel
of Fig. 9 displays this ratio for the two extremely unequal cases, C012-C018
with $q=2.1$ (dashed lines) and C003-C022 with $q=23.2$ (solid lines). The
panel shows clearly that $h^{3,3}$ is much larger for the more unequal case
while the $h^{4,4}$ case is less clear. A more quantitative comparison can be
performed by averaging the ratios $h^{l,m}/h^{2,2}$ over the last few orbits,
corresponding roughly to the range of orbital frequencies
$M_{0}\omega_{\rm{GW}}\in(0.05,0.1)$. In this way, we are then excluding both
the early inspiral, contaminated with significant constraint violations, and
the post-merger phase. The bottom panel of Fig. 9 displays these ratios as a
function of the mass ratio. The mode $h^{3,3}$ increases relative to the
dominant one by almost a factor $1.7$ when passing from $q=2$ to $q=23$, while
$h^{4,4}$ barely changes. Even for our largest mass ratio, the amplitude of
the $h^{3,3}$ mode is at most $25\%$ of the dominant one $h^{2,2}$.

Figure 8:  GWs in the coalescence. Different modes $l=m$ of the strain as a
function of time, together with the frequency of the dominant mode $l=m=2$.
Clearly, the $l=m=2$ mode is always much larger than the others, even for the
largest mass ratio. The vertical, gray lines show the merger time.

Figure 9:  GWs in the inspiral stage. (Top) Ratio between the $l=m\geq 2$
modes of the strain and the dominant $l=m=2$ for the two most extreme cases
C012-C018 ($q=2.1$, dashed lines) and C003-C022 ($q=23.2$, solid lines).
(Bottom) The ratios averaged over the last few orbits for all the mass ratios
studied. Notice that the $l=m=3$ mode increases by a factor 1.7 in this range
of mass ratios, while the $l=m=4$ barely increases with mass ratio.

Figure 10:  GWs in the post-merger stage. (Top) Power spectral density of the
main mode (i.e., $l=m=2$) of the strain for the post-merger. (Bottom)
Frequencies of the first peak of the $l=m=2$ mode for all the cases
represented as a function of the compactness of the remnant $C_{N}\equiv
M_{r}/R_{N}$ (filled circles). Reliable results from the equal-mass mergers of
Paper I are shown with filled squares. The dashed line corresponds to the
lowest quasinormal mode (QNM) frequencies $\omega_{\mathrm{QNM}}$ of isolated
BSs. The fact that the frequencies of the merger remnants are in good
agreement with those of isolated BSs provides further evidence that the
remnants have settled to stable BSs. Notice that we have also included the
cases that collapse to rotating BHs (those with $C_{N}\geq 0.5$).

###  IV.2 Post-merger

We analyze the post-merger frequencies of the gravitational signal of the
remnant, showing the power spectral density of the dominant $l=m=2$ mode in
the top panel of Fig. 10. In the bottom panel of Fig. 10, we display the
frequency of the dominant mode for all the cases studied in this paper and
Paper I, together with the fundamental mode for isolated BS stars as a
function of the compactness of the remnant.

An analysis from Paper I indicates a correspondence between the frequency of
the first peak with the quasi-normal mode (QNM) of isolated BSs. We scrutinize
this hypothesis further by considering the post-merger behavior of all
configurations from both papers as well as the QNM of isolated solitonic BSs
with $C=\\{0.06,0.12,0.18,0.22\\}$ calculated in Paper I. We fit the spectral
lines with a Lorentzian function, i.e.

$\tilde{h}^{2,2}(\omega)\sim\frac{\omega_{I}}{\left(\omega^{2}_{I}+(\omega-\omega_{\rm
R})^{2}\right)}\,$ (20)

to determine the peak frequency of the main mode $\omega_{\rm R}$ and the
inverse decay time $\omega_{I}.$ In line with the discussion on the relaxation
of the remnant from Sec. III.2, one can construct quadratic fits for
$\omega_{\rm R}(M_{r})\,,\omega_{\rm I}(M_{r})$, where $M_{r}=M(N_{r})$ [Eq.
(9)],

$\displaystyle M_{r}\omega_{\rm R}$ $\displaystyle\approx$ $\displaystyle
0.05+0.3M_{r}-0.13(M_{r})^{2}\,,$ (21) $\displaystyle M_{r}\omega_{\rm I}$
$\displaystyle\approx$ $\displaystyle 0.013+0.007M_{r}\,,$ (22)

for isolated scenarios from Paper I. As shown in Fig. 11 (left panel),
excluding the C018-C018 case from Paper I where the post-merger behavior is
not reliable, real parts of the post-merger main mode frequencies agree well
with the isolated solitonic BS QNM fit.

However, in the case of the imaginary frequency (see Fig. 11, right panel) all
remnants produced in the binary coalescence have an offset with respect to the
isolated QNMs. We notice that the three configurations in which blobs do not
form have lower imaginary frequencies (longer decay times) compared to the
isolated configurations. In contrast, in the case of blob formation,
frequencies are higher (shorter decay times) than expected from the isolated
QNMs. Note the three cases with $m_{b}\lambda M_{r}\approx 1.07-1.1$ that have
almost identical real frequencies but vastly different imaginary components.
Understanding this peculiar behavior lies beyond the scope of this paper. We
speculate that the excess angular momentum requires longer decay times in
contrast to the isolated configurations, except in the case of blob formation
that removes the excess (rotational) energy more efficiently than in the
isolated case, thus shortening the decay time.

We have also compared the fits with the tabulated BH QNMs rdweb . The BH
remnant from Paper I, i.e. the C022-C022 case $(a=0.698)$, has tabulated value
$(M_{r}\omega_{R},M_{r}\omega_{I})=$ $(0.532,-0.081)$, while we find
$(0.469,-0.083)$, which is close to the time-domain fit from Paper I where
$(0.5,-0.07)$. For the C022-C012 case $(a=0.5)$, we find the tabulated value
$(0.464,-0.086)$, while the fit gives $(0.475,-0.061)$. This mild discrepancy
between the fit and the predicted ones for BHs may originate from the
numerical precision of the ADM mass/angular momentum extraction and the fit,
the presence of some remnant scalar surrounding the BH, or the fact that the
frequency estimate depend on the choice of the post-merger time. Nonetheless,
the overall agreement corroborates the conclusion that the remnant is a BH.

---
Figure 11: QNMs and post-merger spectrum. (Left) Real and (right) imaginary
parts of the frequencies of the first peak of the $l=m=2$ mode with respect to
remnant mass $M_{r}$ for binaries that form a remnant BS, with the equal mass
binaries of Paper I (filled) and the unequal mass binaries studied here
(open). Binaries that eject a blob are denoted with squares while circles
denote those that do not. The frequencies of the QNMs of the four isolated BSs
used for initial data are marked with stars, and the curves (solid) represent
the fits to the real and imaginary components of these QNM frequencies from
Eqs. (21) and (22).

###  IV.3 Solitonic BSs in the LIGO/Virgo band

In this subsection, we quantify the difference between the GW signal expected
from BS binaries and from binary BHs, focusing on the LIGO/Virgo band. In
particular, we assess whether analyzing LIGO/Virgo data with binary BH
templates can lead to missed detections or to biases on the estimate of the
parameters of the source, under the assumption that the latter consists of a
BS binary.

As a preliminary test of this, we consider the BS binary waveforms extracted
from the unequal-mass simulations of this paper, focusing on the $l=m=2$ mode
alone. Actually, each of these simulations can be taken to represent a binary
of any total mass, as long as frequencies and strain amplitudes are properly
rescaled, i.e. each simulation actually corresponds to a one-parameter family
of systems with varying binary mass $M$, but with fixed dimensionless product
$m_{b}M/\hbar$. We choose therefore to vary $M$ in a range likely to yield
observable effects in the LIGO/Virgo frequency band, i.e. we choose $M$ in the
interval $[M_{\rm min},100]M_{\odot}$, where $M_{\rm min}$ is such that the
smaller progenitor is always heavier than $\sim 3M_{\odot}$. For each BS
waveform obtained in this way, we rescale the (2,2)-mode strain amplitude to
correspond to a fiducial luminosity distance of 400 Mpc. (We recall that
choosing a slightly different distance will simply rescale strains and signal-
to-noise ratios by a linear factor, at leading order.) We then compare the BS
signal obtained to SEOBNRv4 BH binary waveforms Bohe:2016gbl , as implemented
in the Pycbc python package alex_nitz_2021_5347736 . The component masses and
luminosity distance of the BH binary waveform are chosen to match those of the
BS binary, the component BH spins are set to zero, and the initial phase and
merger time are chosen so as to minimize the “difference” of the two signals.
In particular, we minimize the signal-to-noise ratio of the difference of the
two signals,

$\rho(\Delta)=\left[4\int\frac{|\tilde{\Delta}(f)|^{2}}{S_{n}(f)}df\right]^{1/2},$
(23)

with $\Delta(t)\equiv h_{\rm BS}(t)-h_{\rm BH}(t)$ the residual, i.e. the
difference between BS and BH signals (computed for optimal detector
orientation and sky position), and with a tilde denoting a Fourier transform.
The (single-sided) power spectral density of the noise, $S_{n}$, is chosen to
be that of a single LIGO detector. More precisely, we consider both the case
in which $S_{n}$ corresponds to the Livingston detector in O3b
LIGOScientific:2021djp , or to the zero-detuning, high laser power design
sensitivity curve aligo . Accounting for the second LIGO interferometer and
for Virgo will further increase the signal-to-noise ratio, roughly by a factor
$\lesssim\sqrt{3}$ (with the $\lesssim$ due to the fact that the source can
only be optimally placed relative to one detector at a time, and that Virgo is
less sensitive than LIGO in O3b).

Two examples of BS binary waveforms, qualitatively representative of the two
possible post-merger scenarios (i.e. BH or BS remnant), are shown in Fig. 12,
where they are compared to the “most similar” BH binary waveforms identified
with this procedure.

Figure 12: Comparison between $(2,2)$-mode BS waveforms produced from our
simulations and BH binary SEOBNRv4 templates, for a system leading to a BS
($q=8.6$, $M=28M_{\odot}$; top) or a BH ($q=2.9$, $M=30M_{\odot}$; bottom).
The template’s initial phase and merger time are chosen to minimize the
residual signal-to-noise ratio. Both systems are optimally oriented at a
luminosity distance of 400 Mpc, and the BH component masses and spins are set
equal to the BS component masses and to zero, respectively.
---
Figure 13: Inspiral-merger-ringdown BH binary SEOBNRv4 signal-to-noise ratio
$\rho(h_{\rm BH})$ (top) and residual signal-to-noise ratio $\rho(\Delta)$
between BS and BH waveforms (bottom), as functions of the total binary mass.
The residual signal-to-noise ratio is minimized over initial phase and merger
time and then normalized by $\rho(h_{\rm BH})$, for the BS waveforms extracted
from our unequal mass simulations. Both the BS and BH binaries are optimally
oriented at a luminosity distance of 400 Mpc, and the BH component masses and
spins are set respectively to the BS component masses and to zero. For both
signals and templates only the $(2,2)$ mode is included. The signal-to-noise
ratios are computed with the O3b single-detector sensitivity (solid lines),
and with the single-detector design LIGO sensitivity in the zero-detuning,
high laser power configuration (dashed lines).

The signal-to-noise ratio $\rho(h_{\rm BH})$ of the BH binary waveform best
matching each BS signal is shown in the top panel of Fig. 13, as a function of
$M$ and for both the O3b and design LIGO configurations. The signal-to-noise
ratio is computed by using the aforementioned SEOBNRv4 waveforms, which
(unlike our short BS signals) include inspiral, merger and ringdown. In the
bottom panel, we show instead the residual signal-to-noise ratio
$\rho(\Delta)$, minimized over initial phase and merger time for all the
simulations that we have at our disposal.

This residual signal-to-noise ratio is computed by comparing BH and BS
waveforms that are both cut below the minimum frequency at which constraint
violations are significant in our simulations, in order to avoid biasing the
comparison. In this way, the residual signal-to-noise ratio includes
differences between BH and BS waveforms that occur in the post-merger phase
and also in the late inspiral, thus including at least some contribution from
tidal effects, while keeping the impact of initial constraint violations
subleading 666 Improvements in the initial data would change the residual
signal-to-noise ratio only marginally at high masses, while at low masses they
would allow for simulating a longer portion of the inspiral phase. Our
residual signal-to-noise ratios should then be regarded as lower bounds. . We
then normalize $\rho(\Delta)$ by the full inspiral-merger-ringdown BH signal-
to-noise ratio $\rho(h_{\rm BH})$, since we have no access to the full
inspiral-merger-ringdown BS waveforms. As can be seen,
$\rho(\Delta)/\rho(h_{\rm{BH}})$ is always very large and grows with $M$, as
expected because for massive binaries only the merger signal is in the band of
terrestrial interferometers (i.e. for those binaries the BS-BH differences in
the post-merger have a larger relative impact). Also note that
$\rho(\Delta)/\rho(h_{\rm{BH}})$ is smallest (although still quite
significant) in the case where the remnant is a BH $(q=2.9)$. Again, this is
expected: The collision of two non-rotating BHs with $q=2.9$ leads to a
rotating remnant with $a\approx 0.52$ Hofmann:2016yih , close to the value of
the BH spin produced by the BS binary in our simulation ($a\approx 0.5$). As
such, the differences in the merger-ringdown, where most of the in-band power
resides (at least for moderately high masses), are small, c.f. e.g. Fig. 12.

As a rough rule of thumb, residual signal-to-noise ratios
$\rho(\Delta)\lesssim 8/\sqrt{3}\sim 5$ may allow for claiming a BS binary
detection (as opposed to a BH binary one), provided that an accurate
determination of the component masses and spins is available (e.g. thanks to a
long inspiral). In the absence of a sufficiently long detected inspiral, large
residual signal-to-noise ratios may merely lead to biases in the estimation of
the parameters of the source (i.e. one could mistake a BS post-merger signal
for a BH ringdown with remnant mass different from the actual one, and/or non-
zero spin), or even missed detections. As can be seen from Fig. 13, this
second possibility seems the most likely at high masses, for which most of the
inspiral is out of band and BH templates miss most of the signal’s power for
binaries producing a BS remnant. Whether this leads to a bias on the recovered
parameters or just a missed detection should be ascertained by considering
templates with varying BH progenitor spins. However, given the long duration
of the BS post-merger signal (c.f. e.g. Fig. 12), it seems unlikely that it
can be detected by any one BH template, i.e. we expect mainly missed
detections at high $M$, at least for second-generation detectors and for
systems that lead to a BS remnant. For systems that instead lead to BH
formation (e.g. the $q=2.9$ case in Fig. 13), using BH templates may simply
produce a bias on the parameter estimation.

The situation will be more favorable for third-generation interferometers
Kalogera:2021bya  such as the Einstein Telescope or Cosmic Explorer, which
will observe many more inspiral cycles. Not only will this allow for a better
measurement of progenitor masses and spins (which will reduce degeneracies
when comparing the post-merger signal to BH templates), but it may also allow
for measuring the tidal Love number in the late inspiral Cardoso:2017cfl ;
Pacilio:2020jza ; Pacilio:2021jmq . This will provide additional hints on the
BS versus BH nature of the system. We will explore the discovery space of
these detectors, and at the same time refine our analysis, in future work.

##  V Conclusions

The coalescence of BSs allows us to study not only the binary dynamics of one
of the most viable and better motivated models of ECOs, but also the two-body
problem in General Relativity for large mass ratios. The soft dependence of
the BS radius with its mass, at least for the solitonic potential used here,
facilitates the numerical simulations of binaries with very different
compactness, as compared to the more challenging case of asymmetric BH
binaries. Taking advantage of this feature of solitonic BSs, we have studied
numerically the coalescence of unequal-mass binaries with mass ratios ranging
between 2 and 23. The analysis of our simulations, which extends the equal-
mass binaries considered in Paper I (i.e., Ref. PhysRevD.96.104058 ), confirms
many of the findings obtained in that previous study.

The fate of these binary mergers is either a nonrotating BS or a Kerr BH, as
confirmed not only by global quantities and by the structure of the solution,
but also by the gravitational QNMs of the remnant. As in Paper I, we once
again find no evidence that any of these binaries form a rotating BS or a
scalar cloud synchronized in its rotation about a spinning BH. The asymmetry
introduced by the unequal mass of the constituent stars perhaps makes the
formation of either of these remnants less likely.An analysis of the parameter
space indicates the need to refine the initial configurations to assess
whether a rotating remnant can be formed.

For a certain range of the initial angular momentum, the remnant undergoes a
process similar to a tidal disruption in NSs, and a blob of scalar field is
ejected. This process has already been observed in the equal-mass binaries of
Paper I, although the symmetry in that case induced the ejection of two blobs
in opposite directions instead of a single blob observed here. The ejection of
a single blob produces a large recoil of the remnant. In our C012-C018 case,
the estimate of the recoil velocity is more than $10^{4}$ km/s, larger than
the superkicks of binary BHs and large enough to have significant implications
for the expected dynamics of BSs in the universe. Because recent studies
suggest that rotating solitonic BSs should be stable against the
nonaxisymmetric instability 2021PhRvD.103d4022S , the ejection of the scalar
blob is not likely a result of such an instability.

We also evolve an unequal mass binary with one of the stars transformed to an
anti-star. This anti-star completely annihilates upon contact, dispersing
scalar field to infinity. The remaining scalar field settles to a lower mass,
static BS in a clear demonstration of the stability of these solutions under a
strong perturbation.

Regarding the GWs emitted during the coalescence, we have found results
comparable to those of binary BHs: the $l=m=2$ mode of the strain is always
dominant, although higher-order modes become more relevant as the mass ratio
increases. We have estimated how the ratio of the modes depends on the mass
ratio during the last few orbits of the coalescence.

We have also analyzed the prospect of detecting differences between binary BS
and binary BH gravitational signals with ground interferometers. We have found
that while the merger portion of the signal is significantly different between
the two classes of sources (at least if the final merger remnant is a BS),
distinguishing between the two might be difficult with second-generation
detectors due to degeneracies between merger and inspiral parameters. However,
this task will ease considerably with third-generation interferometers, such
as Cosmic Explorer or the Einstein Telescope.

Many interesting questions remain to be addressed, especially regarding the
final state of the remnant. Evolutions of solitonic BSs have yet to produce
either a spinning BS or a synchronized scalar cloud. More accurate and longer
simulations together with improved initial data may shed light on such
questions, or perhaps some a priori analysis will indicate whether and under
what conditions such end states will result.

### Acknowledgments

It is a pleasure to thank Nico Sanchis-Gual for helpful comments and
discussions on the manuscript, as well as Guillermo Lara for useful comments.
M.B, M.B., and E.B. acknowledge support from the European Union’s H2020 ERC
Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” (Grant
No. GRAMS815673). This work was supported by the EU Horizon 2020 Research and
Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No.
101007855. This work was supported by the NSF under grants PHY-1912769 and
PHY-2011383 (SLL). P.P. acknowledges financial support provided under the
European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480. We
also acknowledge support under the MIUR PRIN and FARE programmes (GW-NEXT,
CUP: B84I20000100001), and from the Amaldi Research Center funded by the MIUR
program “Dipartimento di Eccellenza” (CUP: B81I18001170001). This work was
supported by European Union FEDER funds, the Ministry of Science, Innovation
and Universities and the Spanish Agencia Estatal de Investigación grant
PID2019-110301GB-I00 (C.P.). Numerical calculations have been made possible
through a CINECA-INFN agreement, providing access to resources on MARCONI at
CINECA, as well as resources provided by XSEDE. We acknowledge the use of
CINECA HPC resources thanks to the agreement between SISSA and CINECA.

## References

  * (1) LIGO Scientific Collaboration, J. Aasi et al., “Advanced LIGO,” Class. Quant. Grav. 32 (2015) 074001, arXiv:1411.4547 [gr-qc].
  * (2) VIRGO Collaboration, F. Acernese et al., “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32 no. 2, (2015) 024001, arXiv:1408.3978 [gr-qc].
  * (3) LIGO Scientific, VIRGO, KAGRA Collaboration, R. Abbott et al., “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,” arXiv:2111.03606 [gr-qc].
  * (4) LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses,” Phys. Rev. D 102 no. 4, (2020) 043015, arXiv:2004.08342 [astro-ph.HE].
  * (5) LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$,” Astrophys. J. Lett. 892 no. 1, (2020) L3, arXiv:2001.01761 [astro-ph.HE].
  * (6) LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$,” Phys. Rev. Lett. 125 no. 10, (2020) 101102, arXiv:2009.01075 [gr-qc].
  * (7) LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object,” Astrophys. J. Lett. 896 no. 2, (2020) L44, arXiv:2006.12611 [astro-ph.HE].
  * (8) LIGO Scientific, KAGRA, VIRGO Collaboration, R. Abbott et al., “Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences,” Astrophys. J. Lett. 915 no. 1, (2021) L5, arXiv:2106.15163 [astro-ph.HE].
  * (9) KAGRA Collaboration, T. Akutsu et al., “Overview of KAGRA: Calibration, detector characterization, physical environmental monitors, and the geophysics interferometer,” arXiv:2009.09305 [gr-qc].
  * (10) G. F. Giudice, M. McCullough, and A. Urbano, “Hunting for Dark Particles with Gravitational Waves,” JCAP 10 (2016) 001, arXiv:1605.01209 [hep-ph].
  * (11) V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel. 22 no. 1, (2019) 4, arXiv:1904.05363 [gr-qc].
  * (12) S. D. Mathur, “The Fuzzball proposal for black holes: An Elementary review,” Fortsch. Phys. 53 (2005) 793–827, arXiv:hep-th/0502050 [hep-th].
  * (13) P. O. Mazur and E. Mottola, “Gravitational condensate stars: An alternative to black holes,” arXiv:gr-qc/0109035 [gr-qc].
  * (14) T. Damour and S. N. Solodukhin, “Wormholes as black hole foils,” Phys. Rev. D76 (2007) 024016, arXiv:0704.2667 [gr-qc].
  * (15) G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, “Anisotropic stars as ultracompact objects in General Relativity,” Phys. Rev. D 99 no. 10, (2019) 104072, arXiv:1811.07917 [gr-qc].
  * (16) R. Ruffini and S. Bonazzola, “Systems of selfgravitating particles in general relativity and the concept of an equation of state,” Phys. Rev. 187 (1969) 1767–1783.
  * (17) F. E. Schunck and E. W. Mielke, “General relativistic boson stars,” Class. Quant. Grav. 20 (2003) R301–R356, arXiv:0801.0307 [astro-ph].
  * (18) S. L. Liebling and C. Palenzuela, “Dynamical Boson Stars,” Living Reviews in Relativity 15 (Dec., 2012) 6, arXiv:1202.5809 [gr-qc].
  * (19) L. Visinelli, “Boson Stars and Oscillatons: A Review,” arXiv:2109.05481 [gr-qc].
  * (20) E. Seidel and W.-M. Suen, “Formation of solitonic stars through gravitational cooling,” Phys. Rev. Lett. 72 (1994) 2516–2519, arXiv:gr-qc/9309015.
  * (21) F. S. Guzman and L. A. Urena-Lopez, “Gravitational cooling of self-gravitating Bose-Condensates,” Astrophys. J. 645 (2006) 814–819, arXiv:astro-ph/0603613.
  * (22) M. Colpi, S. L. Shapiro, and I. Wasserman, “Boson Stars: Gravitational Equilibria of Selfinteracting Scalar Fields,” Phys. Rev. Lett. 57 (1986) 2485–2488.
  * (23) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, “Astrophysical signatures of boson stars: Quasinormal modes and inspiral resonances,” Phys. Rev. D 88 no. 6, (Sept., 2013) 064046, arXiv:1307.4812 [gr-qc].
  * (24) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D95 no. 4, (2017) 043541, arXiv:1610.08297 [astro-ph.CO].
  * (25) D. J. E. Marsh, “Axion Cosmology,” Phys. Rept. 643 (2016) 1–79, arXiv:1510.07633 [astro-ph.CO].
  * (26) L. Hui, “Wave Dark Matter,” arXiv:2101.11735 [astro-ph.CO].
  * (27) E. Seidel and W. M. Suen, “Oscillating soliton stars,” Phys. Rev. Lett. 66 (1991) 1659–1662.
  * (28) R. Brito, V. Cardoso, C. F. B. Macedo, H. Okawa, and C. Palenzuela, “Interaction between bosonic dark matter and stars,” Phys. Rev. D93 no. 4, (2016) 044045, arXiv:1512.00466 [astro-ph.SR].
  * (29) S. Troitsky, “Supermassive dark-matter Q-balls in galactic centers?,” JCAP 11 (2016) 027, arXiv:1510.07132 [hep-ph].
  * (30) S. Krippendorf, F. Muia, and F. Quevedo, “Moduli Stars,” JHEP 08 (2018) 070, arXiv:1806.04690 [hep-th].
  * (31) D. G. Levkov, A. G. Panin, and I. I. Tkachev, “Gravitational Bose-Einstein condensation in the kinetic regime,” Phys. Rev. Lett. 121 no. 15, (2018) 151301, arXiv:1804.05857 [astro-ph.CO].
  * (32) E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, “Analytic Description of Primordial Black Hole Formation from Scalar Field Fragmentation,” JCAP 1910 (2019) 077, arXiv:1907.10613 [astro-ph.CO].
  * (33) M. A. Amin and P. Mocz, “Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe,” Phys. Rev. D100 no. 6, (2019) 063507, arXiv:1902.07261 [astro-ph.CO].
  * (34) J. Y. Widdicombe, T. Helfer, D. J. E. Marsh, and E. A. Lim, “Formation of Relativistic Axion Stars,” JCAP 1810 no. 10, (2018) 005, arXiv:1806.09367 [astro-ph.CO].
  * (35) A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, “Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter,” Phys. Rev. D101 no. 8, (2020) 083014, arXiv:1909.11665 [astro-ph.CO].
  * (36) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, “Head-on collisions of boson stars,” Phys. Rev. D 75 no. 6, (Mar., 2007) 064005, gr-qc/0612067.
  * (37) C. Palenzuela, L. Lehner, and S. L. Liebling, “Orbital dynamics of binary boson star systems,” Phys. Rev. D 77 no. 4, (Feb., 2008) 044036, arXiv:0706.2435 [gr-qc].
  * (38) T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, “Gravitational Wave Emission from Collisions of Compact Scalar Solitons,” Phys. Rev. D 99 no. 4, (2019) 044046, arXiv:1802.06733 [gr-qc].
  * (39) M. Bezares, C. Palenzuela, and C. Bona, “Final fate of compact boson star mergers,” Phys. Rev. D 95 (Jun, 2017) 124005. https://link.aps.org/doi/10.1103/PhysRevD.95.124005.
  * (40) C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, “Gravitational wave signatures of highly compact boson star binaries,” Phys. Rev. D 96 (Nov, 2017) 104058. https://link.aps.org/doi/10.1103/PhysRevD.96.104058.
  * (41) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, “Proca stars: Gravitating Bose–Einstein condensates of massive spin 1 particles,” Phys. Lett. B 752 (2016) 291–295, arXiv:1508.05395 [gr-qc].
  * (42) N. Sanchis-Gual, C. Herdeiro, E. Radu, J. C. Degollado, and J. A. Font, “Numerical evolutions of spherical Proca stars,” Phys. Rev. D 95 no. 10, (2017) 104028, arXiv:1702.04532 [gr-qc].
  * (43) N. Sanchis-Gual, C. Herdeiro, J. A. Font, E. Radu, and F. Di Giovanni, “Head-on collisions and orbital mergers of proca stars,” Phys. Rev. D 99 (Jan, 2019) 024017. https://link.aps.org/doi/10.1103/PhysRevD.99.024017.
  * (44) J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, “GW190521 as a Merger of Proca Stars: A Potential New Vector Boson of $8.7\times 10^{-13}$ eV,” Phys. Rev. Lett. 126 no. 8, (2021) 081101, arXiv:2009.05376 [gr-qc].
  * (45) T. Dietrich, S. Ossokine, and K. Clough, “Full 3D numerical relativity simulations of neutron star–boson star collisions with BAM,” Class. Quant. Grav. 36 no. 2, (2019) 025002, arXiv:1807.06959 [gr-qc].
  * (46) T. Dietrich, F. Day, K. Clough, M. Coughlin, and J. Niemeyer, “Neutron star–axion star collisions in the light of multimessenger astronomy,” Mon. Not. Roy. Astron. Soc. 483 no. 1, (2019) 908–914, arXiv:1808.04746 [astro-ph.HE].
  * (47) K. Clough, T. Dietrich, and J. C. Niemeyer, “Axion star collisions with black holes and neutron stars in full 3D numerical relativity,” Phys. Rev. D 98 no. 8, (2018) 083020, arXiv:1808.04668 [gr-qc].
  * (48) M. Bezares and C. Palenzuela, “Gravitational Waves from Dark Boson Star binary mergers,” Class. Quant. Grav. 35 no. 23, (2018) 234002, arXiv:1808.10732 [gr-qc].
  * (49) M. Bezares, D. Viganò, and C. Palenzuela, “Gravitational wave signatures of dark matter cores in binary neutron star mergers by using numerical simulations,” Phys. Rev. D 100 no. 4, (2019) 044049, arXiv:1905.08551 [gr-qc].
  * (50) S. Valdez-Alvarado, C. Palenzuela, D. Alic, and L. A. Ureña-López, “Dynamical evolution of fermion-boson stars,” Phys. Rev. D 87 no. 8, (Apr., 2013) 084040, arXiv:1210.2299 [gr-qc].
  * (51) R. Friedberg, T. D. Lee, and Y. Pang, “Scalar soliton stars and black holes,” Phys. Rev. D 35 (June, 1987) 3658–3677.
  * (52) M. Bošković and E. Barausse, “Soliton boson stars, Q-balls and the causal Buchdahl bound,” arXiv:2111.03870 [gr-qc].
  * (53) V. Cardoso, C. F. B. Macedo, K.-i. Maeda, and H. Okawa, “ECO-spotting: looking for extremely compact objects with bosonic fields,” arXiv:2112.05750 [gr-qc].
  * (54) T. Tamaki and N. Sakai, “How does gravity save or kill Q-balls?,” Phys. Rev. D 83 (2011) 044027, arXiv:1105.2932 [gr-qc].
  * (55) J. A. Gonzalez, U. Sperhake, and B. Bruegmann, “Black-hole binary simulations: The Mass ratio 10:1,” Phys. Rev. D 79 (2009) 124006, arXiv:0811.3952 [gr-qc].
  * (56) C. O. Lousto, H. Nakano, Y. Zlochower, and M. Campanelli, “Intermediate-mass-ratio black hole binaries: Intertwining numerical and perturbative techniques,” Phys. Rev. D 82 (2010) 104057, arXiv:1008.4360 [gr-qc].
  * (57) C. O. Lousto and Y. Zlochower, “Orbital Evolution of Extreme-Mass-Ratio Black-Hole Binaries with Numerical Relativity,” Phys. Rev. Lett. 106 (2011) 041101, arXiv:1009.0292 [gr-qc].
  * (58) D. Müller, Numerical Simulations of Black Hole Binaries with Unequal Masses.  PhD thesis, PhD Thesis, 2011, 2011.
  * (59) N. Sanchis-Gual, M. Zilhão, C. Herdeiro, F. Di Giovanni, J. A. Font, and E. Radu, “Synchronized gravitational atoms from mergers of bosonic stars,” Phys. Rev. D 102 no. 10, (2020) 101504, arXiv:2007.11584 [gr-qc].
  * (60) M. Kesden, J. Gair, and M. Kamionkowski, “Gravitational-wave signature of an inspiral into a supermassive horizonless object,” Phys. Rev. D 71 (2005) 044015, arXiv:astro-ph/0411478.
  * (61) V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, “Testing strong-field gravity with tidal Love numbers,” Phys. Rev. D 95 no. 8, (2017) 084014, arXiv:1701.01116 [gr-qc]. [Addendum: Phys.Rev.D 95, 089901 (2017)].
  * (62) N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno, and S. Ossokine, “Distinguishing Boson Stars from Black Holes and Neutron Stars from Tidal Interactions in Inspiraling Binary Systems,” Phys. Rev. D96 no. 2, (2017) 024002, arXiv:1704.08651 [gr-qc].
  * (63) T. Helfer, U. Sperhake, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, “Malaise and remedy of binary boson-star initial data,” arXiv:2108.11995 [gr-qc].
  * (64) S. Mirshekari and C. M. Will, “Compact binary systems in scalar-tensor gravity: Equations of motion to 2.5 post-Newtonian order,” Phys. Rev. D 87 no. 8, (2013) 084070, arXiv:1301.4680 [gr-qc].
  * (65) A. Arbona, A. Artigues, C. Bona-Casas, J. Massó, B. Miñano, A. Rigo, M. Trias, and C. Bona, “Simflowny: A general-purpose platform for the management of physical models and simulation problems,” Computer Physics Communications 184 no. 10, (2013) 2321 – 2331. http://www.sciencedirect.com/science/article/pii/S0010465513001471.
  * (66) A. Arbona, B. Miñano, A. Rigo, C. Bona, C. Palenzuela, A. Artigues, C. Bona-Casas, and J. Massó, “Simflowny 2: An upgraded platform for scientific modelling and simulation,” Computer Physics Communications 229 (2018) 170–181. https://www.sciencedirect.com/science/article/pii/S0010465518300870.
  * (67) C. Palenzuela, B. Miñano, A. Arbona, C. Bona-Casas, C. Bona, and J. Massó, “Simflowny 3: An upgraded platform for scientific modeling and simulation,” Computer Physics Communications 259 (2021) 107675. https://www.sciencedirect.com/science/article/pii/S0010465520303271.
  * (68) “Simflowny project website.,” 2021.  https://bitbucket.org/iac3/simflowny/wiki/Home.
  * (69) R. D. Hornung and S. R. Kohn, “Managing application complexity in the samrai object-oriented framework,” Concurrency and Computation: Practice and Experience 14 no. 5, (2002) 347–368. http://dx.doi.org/10.1002/cpe.652.
  * (70) B. T. Gunney and R. W. Anderson, “Advances in patch-based adaptive mesh refinement scalability,” Journal of Parallel and Distributed Computing 89 (2016) 65 – 84. http://www.sciencedirect.com/science/article/pii/S0743731515002129.
  * (71) “Samrai project website.,” 2021.  https://computation.llnl.gov/project/SAMRAI/.
  * (72) C. Palenzuela, B. Miñano, D. Viganò, A. Arbona, C. Bona-Casas, A. Rigo, M. Bezares, C. Bona, and J. Massó, “A Simflowny-based finite-difference code for high-performance computing in numerical relativity,” Class. Quant. Grav. 35 no. 18, (2018) 185007, arXiv:1806.04182 [physics.comp-ph].
  * (73) N. T. Bishop and L. Rezzolla, “Extraction of gravitational waves in numerical relativity,” Living Reviews in Relativity 19 (Oct., 2016) 2, arXiv:1606.02532 [gr-qc].
  * (74) B. Brügmann, J. A. González, M. Hannam, S. Husa, U. Sperhake, and W. Tichy, “Calibration of moving puncture simulations,” Phys. Rev. D 77 no. 2, (Jan., 2008) 024027, gr-qc/0610128.
  * (75) C. Reisswig and D. Pollney, “Notes on the integration of numerical relativity waveforms,” Classical and Quantum Gravity 28 no. 19, (Oct., 2011) 195015, arXiv:1006.1632 [gr-qc].
  * (76) E. E. Flanagan and T. Hinderer, “Constraining neutron star tidal Love numbers with gravitational wave detectors,” Phys. Rev. D 77 (2008) 021502, arXiv:0709.1915 [astro-ph].
  * (77) E. Poisson and C. M. Will, Gravity: Newtonian, Post-Newtonian, Relativistic.  Cambridge University Press, 2014.
  * (78) N. Siemonsen and W. E. East, “Stability of rotating scalar boson stars with nonlinear interactions,” Phys. Rev. D 103 no. 4, (Feb., 2021) 044022, arXiv:2011.08247 [gr-qc].
  * (79) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, “Nonlinear Dynamics of Spinning Bosonic Stars: Formation and Stability,” Phys. Rev. Lett.  123 no. 22, (Nov., 2019) 221101, arXiv:1907.12565 [gr-qc].
  * (80) A. S. Dmitriev, D. G. Levkov, A. G. Panin, E. K. Pushnaya, and I. I. Tkachev, “Instability of rotating Bose stars,” Phys. Rev. D 104 no. 2, (2021) 023504, arXiv:2104.00962 [gr-qc].
  * (81) B. Kleihaus, J. Kunz, and M. List, “Rotating boson stars and Q-balls,” Phys. Rev. D 72 (2005) 064002, arXiv:gr-qc/0505143.
  * (82) S. R. Coleman, “Q-balls,” Nucl. Phys. B 262 no. 2, (1985) 263. [Addendum: Nucl.Phys.B 269, 744 (1986)].
  * (83) R. Brito, V. Cardoso, and P. Pani, “Superradiance: New Frontiers in Black Hole Physics,” Lect. Notes Phys. 906 (2015) pp.1–237, arXiv:1501.06570 [gr-qc].
  * (84) A. A. Starobinskij and S. M. Churilov, “Amplification of electromagnetic and gravitational waves scattered by a rotating black hole.,” Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki 65 (1973) 3–11.
  * (85) L. Boyle, M. Kesden, and S. Nissanke, “Binary black hole merger: Symmetry and the spin expansion,” Phys. Rev. Lett. 100 (2008) 151101, arXiv:0709.0299 [gr-qc].
  * (86) M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, “Maximum gravitational recoil,” Phys. Rev. Lett. 98 (2007) 231102, arXiv:gr-qc/0702133.
  * (87) J. A. Gonzalez, M. D. Hannam, U. Sperhake, B. Bruegmann, and S. Husa, “Supermassive recoil velocities for binary black-hole mergers with antialigned spins,” Phys. Rev. Lett. 98 (2007) 231101, arXiv:gr-qc/0702052.
  * (88) C. O. Lousto and Y. Zlochower, “Hangup Kicks: Still Larger Recoils by Partial Spin/Orbit Alignment of Black-Hole Binaries,” Phys. Rev. Lett. 107 (2011) 231102, arXiv:1108.2009 [gr-qc].
  * (89) F. Antonini and F. A. Rasio, “Merging black hole binaries in galactic nuclei: implications for advanced-LIGO detections,” Astrophys. J. 831 no. 2, (2016) 187, arXiv:1606.04889 [astro-ph.HE].
  * (90) T. Bogdanovic, C. S. Reynolds, and M. C. Miller, “Alignment of the spins of supermassive black holes prior to merger,” Astrophys. J. Lett. 661 (2007) L147, arXiv:astro-ph/0703054.
  * (91) E. Barausse and A. Lapi, “Massive Black Hole Mergers,” arXiv:2011.01994 [astro-ph.GA].
  * (92) S. Dupourqué, L. Tibaldo, and P. Von Ballmoos, “Constraints on the antistar fraction in the Solar System neighborhood from the 10-year Fermi Large Area Telescope gamma-ray source catalog,” Phys. Rev. D 103 no. 8, (2021) 083016, arXiv:2103.10073 [astro-ph.HE].
  * (93) P. Pani and A. Maselli, “Love in Extrema Ratio,” Int. J. Mod. Phys. D 28 no. 14, (2019) 1944001, arXiv:1905.03947 [gr-qc].
  * (94) Webpage with Mathematica notebooks and numerical quasinormal mode Tables:
http://centra.tecnico.ulisboa.pt/network/grit/files/
http://www.phy.olemiss.edu/~berti/ringdown/
http://www.roma1.infn.it/~pani/ .

  * (95) A. Bohé et al., “Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors,” Phys. Rev. D 95 no. 4, (2017) 044028, arXiv:1611.03703 [gr-qc].
  * (96) A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, T. Dent, L. Pekowsky, A. R. Williamson, G. S. C. Davies, S. De, M. Cabero, B. Machenschalk, P. Kumar, D. Macleod, S. Reyes, dfinstad, F. Pannarale, T. Massinger, S. Kumar, M. Tápai, L. Singer, S. Khan, S. Fairhurst, A. Nielsen, S. Singh, K. Chandra, shasvath, and B. U. V. Gadre, “gwastro/pycbc:,” Aug., 2021.  https://doi.org/10.5281/zenodo.5347736.
  * (97) “Advanced LIGO anticipated sensitivity curves : LIGO document T0900288-v3.” https://dcc.ligo.org/LIGO-T0900288/public.
  * (98) F. Hofmann, E. Barausse, and L. Rezzolla, “The final spin from binary black holes in quasi-circular orbits,” Astrophys. J. Lett. 825 no. 2, (2016) L19, arXiv:1605.01938 [gr-qc].
  * (99) V. Kalogera et al., “The Next Generation Global Gravitational Wave Observatory: The Science Book,” arXiv:2111.06990 [gr-qc].
  * (100) C. Pacilio, M. Vaglio, A. Maselli, and P. Pani, “Gravitational-wave detectors as particle-physics laboratories: Constraining scalar interactions with a coherent inspiral model of boson-star binaries,” Phys. Rev. D 102 no. 8, (2020) 083002, arXiv:2007.05264 [gr-qc].
  * (101) C. Pacilio, A. Maselli, M. Fasano, and P. Pani, “Ranking the Love for the neutron star equation of state: the need for third-generation detectors,” arXiv:2104.10035 [gr-qc].