# Nonlinear treatment of a black hole mimicker ringdown

Nils Siemonsen<EMAIL_ADDRESS>Princeton Gravity Initiative,
Princeton University, Princeton NJ 08544, USA  Department of Physics,
Princeton University, Princeton, NJ 08544, USA

(August 28, 2024)

###### Abstract

We perform the first nonlinear and self-consistent study of the merger and
ringdown of a black hole mimicking object with stable light rings. To that
end, we numerically solve the full Einstein-Klein-Gordon equations governing
the head-on collisions of a series of binary boson stars in the large-mass-
ratio regime resulting in spinning horizonless remnants with stable light
rings. We broadly confirm the appearance of features in the extracted
gravitational waveforms expected based on perturbative methods: the signal
from the prompt response of the remnants approaches that of a Kerr black hole
in the large-compactness limit, and the subsequent emissions contain
periodically appearing bursts akin to so-called gravitational wave echoes.
However, these bursts occur at high frequencies and are sourced by
perturbations of the remnant’s internal degrees of freedom. Furthermore, the
emitted waveforms also contain a large-amplitude and long-lived component
comparable in frequency to black hole quasi-normal modes. We further
characterize the emissions, obtain basic scaling relations of relevant
timescales, and compute the energy emitted in gravitational waves.

##  I Introduction

The black hole is a remarkably successful paradigm explaining astrophysical
observations of electromagnetic and gravitational wave signals driven by
highly compact and dark objects. Despite this success, and motivated by
extensions to the Standard Model of particle physics or models of quantum
gravity, a large class of alternative objects has been developed Cardoso and
Pani (2019), challenging this paradigm. These objects—black hole mimickers—are
horizonless and far more compact than ordinary neutron stars, imitating many
or all of the observable signatures of black holes. In stark contrast to
ordinary black holes, unstable light rings around these mimickers are
accompanied by stable counterparts. Cardoso et al. pointed out that this light
ring structure plays a central role in the ringdown of any black hole mimicker
Cardoso  _et al._ (2016a, b); Cardoso and Pani (2017). They observed that the
gravitational wave emissions promptly after an extreme-mass-ratio merger of a
binary black hole mimicker is universally identical to that of a black hole
ringdown, independently of the mimicker’s internal structure. The emitted
signal deviates from the black hole ringdown only after a light-crossing time
of the interior of the mimicker and is characterized by repeated burst-like
gravitational wave echoes. As argued in Refs. Cardoso  _et al._ (2016a, b);
Cardoso and Pani (2017), the prompt response can be understood as the partial
reflection of the gravitational waves sourced by the plunging companion off of
the unstable light ring of the primary black hole mimicker. The subsequent
echoes are due to repeated leakage of the transmitted gravitational
perturbations traversing across the stable light ring in the mimicker’s
interior.

As a mechanism to discover as of yet unknown new physics, the ringdown of
black hole mimickers has received much attention. Broadly, efforts have been
devoted to understanding the ringdown of various classes of black hole
mimickers Cardoso  _et al._ (2016a, b); Holdom and Ren (2017); Bueno  _et al._
(2018); Barceló  _et al._ (2017); Urbano and Veermäe (2019); Raposo  _et al._
(2019); Pani and Ferrari (2018); Cardoso  _et al._ (2019); Oshita and Afshordi
(2019); Wang  _et al._ (2020); Oshita  _et al._ (2020); Maggio  _et al._
(2020); Dey  _et al._ (2020); Ikeda  _et al._ (2021), computing gravitational
waveforms and templates Mark  _et al._ (2017); Nakano  _et al._ (2017); Wang
and Afshordi (2018); Burgess  _et al._ (2018); Correia and Cardoso (2018);
Longo Micchi and Chirenti (2020); Maggio  _et al._ (2019); Srivastava and Chen
(2021); Annulli  _et al._ (2022); Xin  _et al._ (2021); Ma  _et al._ (2022),
determining detection prospects with gravitational wave detectors Maselli  _et
al._ (2017); Testa and Pani (2018); Tsang  _et al._ (2018); Longo Micchi  _et
al._ (2021), and searching for burst-like emissions in the gravitational wave
data following known binary coalescence events Abedi  _et al._ (2017); Ashton
_et al._ (2016); Westerweck  _et al._ (2018); Conklin  _et al._ (2018);
Nielsen  _et al._ (2019); Lo  _et al._ (2019); Tsang  _et al._ (2020);
Uchikata  _et al._ (2019); Abbott  _et al._ (2021a, b); Uchikata  _et al._
(2023) (see Ref. Abedi  _et al._ (2020) for a review). Despite the immense
progress of this program, results were largely obtained modeling the internal
structure of the mimicker by boundary conditions a small distance away from
the would-be horizon and treating it’s dynamical response at the test-field
and extreme-mass-ratio level. Thereby neglecting, (i) nonlinear gravitational
effects, (ii) coupling of the mimicker’s internal structure to gravitational
degrees of freedom even at the linear level, (iii) self-interactions of the
(effective) matter making up the object, and (iv) any finite size effects
(including spin) of the objects. Notably, a few approaches have been developed
in Refs. Carballo-Rubio  _et al._ (2018); Danielsson  _et al._ (2021);
Vellucci  _et al._ (2023); Dailey  _et al._ (2023) to address some of these
shortcomings. However, a fully nonlinear and self-consistent treatment of the
merger and ringdown of a black hole mimicker is still lacking. This leaves
many important questions unanswered. In particular, when including all effects
(i)–(iv), are the prompt emissions indeed identical to those of a ringing
black hole? What are the amplitudes of the quasi-normal modes of the remnant
excited during the merger? Specifically, is the waveform following the prompt
ringdown solely characterized by gravitational wave echoes? How does the
remnant’s spin impact these conclusions?

In this work, we address these questions fully self-consistently in the large-
mass-ratio and head-on merger setting. To that end, we drop all restrictions
(i)-(iv) by performing a series of four fully nonlinear numerical time-domain
evolutions of coalescing spinning binary boson stars—a particular set of black
hole mimickers—resulting in spinning horizonless remnants with stable light
rings. We find that the prompt dynamical response of the remnants, encoded in
the emitted gravitational waves, approaches that of Kerr black holes in the
large-compactness limit. The subsequent emissions contain both a high-
frequency burst-like component with frequency set the binary’s size-ratio, and
a large-amplitude and long-lived component indicative of excited trapped modes
in the remnant’s interiors.

##  II Model & Methods

$C$ | $\tilde{\omega}/\mu$ | $J/M^{2}$ | $M\mu$ | $q$ | $\eta$
---|---|---|---|---|---
$0.38$ | $0.105$ | $0.95$ | $18.7$ | $78$ | $4.4$
$0.36$ | $0.111$ | $0.98$ | $17.5$ | $73$ | $4.3$
$0.32$ | $0.130$ | $1.08$ | $13.9$ | $58$ | $3.9$
$0.28$ | $0.150$ | $1.22$ | $11.0$ | $46$ | $3.5$
Table 1: The compactness $C$, frequency $\tilde{\omega}$, anuglar momentum
$J$, and mass $M$ of the primary boson star in the sequence of binaries. The
secondary is the same in all binaries considered and has frequency
$\tilde{\omega}_{*}/\mu=0.85$, mass $M_{*}\mu=0.24$, and angular momentum
$J_{*}/M_{*}^{2}=12.1$. The last two columns show the mass- and size-ratios
$q$ and $\eta$ of the binaries, respectively.

Boson stars Kaup (1968); Ruffini and Bonazzola (1969) are the most developed
highly compact and black hole mimicking objects that can be treated within
numerical relativity Schunck and Mielke (2003); Liebling and Palenzuela
(2023). Those stars relevant for this work are stationary solutions in the
theory with Lagrangian density Friedberg  _et al._ (1976)

$\displaystyle\mathcal{L}=\frac{R}{16\pi}$
$\displaystyle-g^{\mu\nu}\partial_{(\mu}\bar{\Phi}\partial_{\nu)}\Phi-\mu^{2}|\Phi|^{2}\left(1-\frac{2|\Phi|^{2}}{\sigma^{2}}\right)^{2},$
(1)

where $g_{\mu\nu}$ is the metric with Ricci curvature scalar $R$, and $\Phi$
is the complex scalar field making up the star with mass parameter $\mu$ and
self-interaction strength $\sigma$. Here and in the following, we employ
$G=c=1$ units. Boson stars are characterized by a harmonic dependence on both
coordinate time $t$ and azimuthal coordinate $\varphi$, i.e., $\Phi\propto
e^{i(\tilde{\omega}t+\tilde{m}\varphi)}$, set by their internal frequency
$\tilde{\omega}$ and index $\tilde{m}$. Therefore, spinning stars, i.e., those
with $|\tilde{m}|>0$, exhibit toroidal surfaces of constant scalar field
magnitude $|\Phi|$. So far, no binary boson star evolution resulted in
remnants with stable light rings Cardoso  _et al._ (2016b); Palenzuela  _et
al._ (2007, 2008, 2017); Bezares  _et al._ (2017); Helfer  _et al._ (2022);
Bezares and Palenzuela (2018); Bezares  _et al._ (2022); Evstafyeva  _et al._
(2023); Siemonsen and East (2023a, b). We proceed by focusing entirely on
spinning stars Schunck and Mielke (1996); Volkov and Wohnert (2002); Kleihaus
_et al._ (2005) with $\tilde{m}=3$ in models with self-interaction strength
$\sigma=0.035$. The sequence of four binary boson stars considered in this
work have mass-ratio $q=M/M_{*}$ and size-ratio $\eta=R/R_{*}$ between the
primary (heavier) star of mass $M$ and radius $R$ and the secondary (lighter)
binary constituent of mass $M_{*}$ and radius $R_{*}$. We define $M_{0}\approx
M$ to be the ADM mass of the binaries for later convenience. In all cases, the
primary exhibits large compactness $C=M/R$, whereas the secondary is the same
spinning star in all cases with $C_{*}=0.02$. The properties of the binaries
and their constituents are summarized in Table 1. See Appendix A for a
discussion of this choice. Isolated star solutions are constructed numerically
using methods outlined in Ref. Siemonsen and East (2021) and used to generate
constraint-satisfying binary initial data as described in Ref. Siemonsen and
East (2023b). All primary stars, and due to the large mass-ratio all merger
remnants, exhibit a pair of counter rotating stable and unstable light rings
(see Appendix A for further details). In what follows, we identify the
binaries by the compactness of the primary constituent. As a point of
comparison, we also evolve a single black hole-boson star collision, where we
replaced the $C=0.38$ primary star by a black hole of the same mass and spin.
In all cases, the initial coordinate separation is $D/M_{0}=160$ and the
objects are boosted to Newtonian freefall velocities from infinity. The
numerical evolutions are performed imposing axisymmetry on the metric and
azimuthal symmetry $\partial_{\varphi}\Phi=3i\Phi$ on the scalar field; see
Appendix B for further details, as well as a discussion on numerical
resolution and convergence.

##  III Results

Figure 1: The gravitational waves emitted during the merger of the binaries
labelled by their primary star compactness compared to the black hole-boson
star binary labelled “Kerr” (top and center panels contain the $C=0.38$ and
“Kerr” binaries only). We show the $\ell$th spin-weighted spherical harmonic
component of the Newman-Penrose scalar $\tilde{\Psi}_{4}^{\ell}=r_{\rm
extr.}M_{0}\Psi_{4}^{\ell}$ extracted on spheres of coordinate radius $r_{\rm
extr.}/M_{0}=100$ and with coordinate time $t$ (recall, $m=0$ in axisymmetry).
All waveforms are time shifted by aligning at the largest peak of the $\ell=2$
mode, whereas $t_{0}$ corresponds to the transition time from prompt to
subsequent emissions of the $C=0.38$ binary. The colored band indicates
uncertainties of our methods. Figure 2: The state of the $C=0.38$ binary
system at selected coordinate times $\tilde{t}=(t-t_{0}+r_{\rm extr.})/M_{0}$
during the plunge, merger, and ringdown. The main plots show the Newman-
Penrose scalar $\Psi_{4}$ at coordinate radii $r>3M_{0}$, which qualitatively
shows the local gravitational waves. We indicate the initial production and
subsequent propagation of the high-frequency component of the prompt response
as “PR” together with a dotted white circle (roughly corresponding to the
propagating wavefronts). Similarly, the first two gravitational wave bursts
are labelled $1^{\rm st}$ and $2^{\rm nd}$. In the region $r<3M_{0}$, we show
two (black and gray) surfaces of constant scalar field magnitude $|\Phi|$
(with a close-up in the insets); recall, spinning boson stars are toroidally
shaped. Arrows in the insets indicate the motion of the largest perturbation
of $|\Phi|$ between snapshots. The symmetry axis is at $x=0$.

Our first main result is that the prompt dynamical response of the merger
remnant encoded in the emitted gravitational waves approaches the response of
a Kerr black hole towards large-compactnesses. In Figure 1, we show the
waveforms emitted during the plunge and merger of the sequence of binaries and
compare these with the aforementioned black hole-boson star binary. While
there are differences in amplitude in the $\ell=2,3$ modes, the frequencies of
all $\ell=2,3,4$ modes broadly match with the black hole response. Since
$C=0.38$ is still a relatively low remnant compactness (with relatively short
light-crossing times; more details below), the peaks of the $\ell=2,3,4$
“Kerr”-waveforms are burried in the subsequent emissions (described below).
However, we find higher-$\ell$ modes to peak earlier leading to a separation
of the prompt response and subsequent emissions. Fitting for these $\ell=8$
frequencies and decay rates of the $C=0.38$ and “Kerr” cases (shown in Figure
1), we find that the former differ by $\approx 15\%$, while the latter are
consistent to within the uncertainties of our methods (though, the peak times
differ by $\approx 5M_{0}$). We discuss these differences below. Lastly, in
the bottom panel of Figure 1, we present the $\ell=4$ waveforms for all four
binaries. From there we conclude that the more compact the remnant (i.e., the
longer the interior’s light crossing time) the more and longer it’s response
is black hole-like. All in all, our nonlinear and self-consistent treatment of
the large-mass-ratio problem broadly confirms expectations based on
perturbative calculations.

Figure 3:  (left panels) The $\ell=2$ (top) and $\ell=8$ (bottom) components
of the gravitational waves emitted after the merger of the sequence of four
binary boson stars compared with that of the black hole-boson star binary
(labelled “Kerr”). The waveforms are aligned in time as in Figure 1. (central
panels) Close-ups of the first bursts in each of the two polar modes. See
Appendix B for details on uncertainties of these waveforms. (right panels) The
Fourier transform of the full $C=0.38$ signal as function of angular frequency
$\omega$, focusing only on dominant radiation components in each mode: the
long-lived component in the $\ell=2$ mode (top) and the burst-like high-
frequency component in the $\ell=8$ mode (bottom); we indicate $f_{\rm low}$
and $f_{\rm burst}$.

To understand the system’s dynamics during the plunge and merger, we present
snapshots of the evolution of the $C=0.38$ case in Figure 2. Around
$\tilde{t}\approx-25$ (the high-$\ell$ part of) the prompt Kerr response is
produced as the secondary approaches the primary111Note, while the $C=0.38$
star exhibits no strictly bound polar null geodesics (i.e., a light sphere),
there are “quasi”-bound such geodesics as detailed in Appendix A. (first
panel). While the latter propagates outwards, the binary merges and well-
localized perturbations propagate along the symmetry axis through the
remnant’s interior (second panel). As these perturbations reach the opposing
side and begin propagating poloidially around the remnant’s outer edge, the
first gravitational (and scalar) wave burst is emitted (third panel). Before
the second burst is emitted, the now less well-localized perturbations
propagate through the interior along the symmetry axis a second time as
indicated by the gray arrows (fourth panel). The perturbations come around a
third time in the last panel. At this stage, however, the initially localized
perturbations dispersed into a collection of gravitationally bound states in-
and outside the remnant as can be seen in the inset of the last panel of
Figure 2. As a point of comparison, the light crossing time of null geodesics
in the $C=0.38$ isolated solution traveling along the axis from $z=M/C$ to
$z=-M/C$ as seen by distant observers is $\tau_{\rm axis}/M=19.6$ (see
Appendix A).

Turning now to the gravitational waveform emitted after the prompt response.
Our second main result is that these emissions exhibit bursts akin to
gravitational wave echoes reappearing after roughly a light-crossing time
through the remnants interior and with frequency set by the secondary’s size.
In Figure 3, we present the $\ell=2$ and $\ell=8$ gravitational waveforms
after the prompt emissions. Both modes contain a set of high-frequency bursts
of frequency $f_{\rm burst}$ and separated by a timescale $\tau_{\rm burst}$.
The burst’s frequency can be understood as follows: The spatial scale
$R_{*}\approx M_{0}/(\eta C)$ of the secondary object sets the lengthscale
$\lambda$ of the perturbations sourced within the remnants interior:
$\lambda\lesssim R_{*}$. In Figure 4, we identify these perturbations as the
source of the bursts, by finding good agreement between the underlying
lengthscales, i.e., $f_{\rm burst}^{-1}\lesssim R_{*}$, for all four binaries
considered. The width of the peak of the Fourier transform in the bottom right
panel of Figure 3 around $f_{\rm burst}$ is linked to the decay rate of the
burst, which roughly aligns with $R_{*}$ as well. From the central panels of
Figure 3, we see that in contrast to the black hole ringdown frequencies, the
burst frequency is independent of the polar mode number $\ell$. The burst
period $\tau_{\rm burst}$ is less straightfoward to understand. At the test-
field level, high-frequency gravitational perturbations sourced by the merger
are propagating roughly along null geodesics. Therefore, most naturally the
burst period is compared to the light crossing time (with respect to far away
observers) of null geodesics traversing the remnants interior along the
direction of the perturbations sourced during the merger, i.e., $\tau_{\rm
axis}$, as defined above. In Figure 4, we compare $\tau_{\rm burst}$ against
this light crossing time $\tau_{\rm axis}$, again finding good agreement
$\tau_{\rm burst}\approx\tau_{\rm axis}$ for all remnant compactnesses.
However, purely considering null geodesics to determine the burst period
neglects that during the merger massless (gravitational) and massive (scalar)
modes are coupled and the relevant dynamics are not confined only to the axis.
Hence, $\tau_{\rm axis}$ is to be understood as a rough measure of crossing
times of high-frequency modes through these remnants more generally (see
Appendix A for other light crossing times). Finally, from Figure 3 we also
conclude that the amplitude of bursts following the first decreases with
increasing compactness. Overall, the appearence of burst-like emission
components separated by a light crossing time of the remnants interior is
consistent with the expectations based on test-field calculations. On the
other hand, the burst’s frequency $f_{\rm
burst}\gtrsim\mathcal{O}(1)M_{0}>f_{\ell m=20}^{\rm Kerr}$, at the considered
size-ratios, is much larger than the corresponding quasi-normal mode frequency
of a black hole with the same mass and spin Berti  _et al._ (2009, 2006), and
hence, larger than test-field computations predict the gravitational wave echo
frequency, and decay rates, to be Mark  _et al._ (2017); Wang and Afshordi
(2018); Maggio  _et al._ (2020); this is further discussed below.

Our third main result is that the low-multipole gravitational wave emissions
contain a large-amplitude and long-lived component. As evident from Figure 3,
the $\ell=2$ mode exhibits low-frequency oscillations between the high-
frequency burst. Physically, this quadrupolar contribution originates from
oscillation modes of the remnant star excited during the merger. In Figure 4,
we compare the fundamental frequency $f_{\rm low}$ of this component to the
oscillation period $\tau_{\rm osc.}$ of the remnants (measured from post-
merger oscillations of $\max|\Phi|$), finding good agreement as a function of
remnant compactness: $\tau_{\rm osc.}\approx f^{-1}_{\rm low.}$. The signal,
$\tilde{\Psi}_{4}$, is dominated by the $n=1$ harmonic of $f_{\rm low}$, where
each harmonic has frequency $(n+1)f_{\rm low}$ (see also the Fourier transform
in Figure 3, with evidence of a second set of overtones). Notice also that
$f_{\rm low}\approx f_{\ell m=20}^{\rm Kerr}/2$ (to within $\approx 25\%$).
Roughly, the dominant frequency in the $\ell$th mode of this second
gravitational wave component scales as $f^{\ell>2}_{\rm low}\sim\ell$ (similar
to black hole quasi-normal modes). As can be seen in Figure 3, this component
is long-lived, with longer decay timescales than our simulations. The equally
spaced frequency spectrum, the long-lived nature, together with the
observation $f^{-1}_{\rm low}\approx 2\tau_{\rm burst}$, are indicative of
trapped modes, which objects with stable light rings generally exhibit
Kokkotas and Schmidt (1999) (see also e.g., Ref. Macedo  _et al._ (2013); Wang
_et al._ (2020); Heidmann  _et al._ (2023)). Lastly, we find no evidence that
the amplitude of this long-lived component decreases with increasing
compactness.

Figure 4: The relevant timescales associated with the gravitational waveforms
emitted after the prompt emissions for all four binaries (identified by the
primary’s compactness $C$). We show the secondary’s size $R_{*}\approx
M_{0}/(\eta C)$, the frequency $f_{\rm burst}$ of the first burst, the light
crossing time $\tau_{\rm axis}$, the burst period $\tau_{\rm burst}$, the
dominant large-scale remnant oscillation mode period $\tau_{\rm osc.}$, and
the fundamental frequency $f_{\rm low}$ of the component of the waveform
between bursts.

The amplitude of $\tilde{\Psi}_{4}$, however, should be interpreted with care.
The observationally relevant strain $h$ scales as
$h(\omega)\sim\tilde{\Psi}_{4}(\omega)/\omega^{2}$ in the Fourier domain.
Hence, high-frequency features are suppressed by a factor of $\sim\omega^{-2}$
compared to low-frequency components of $\tilde{\Psi}_{4}$. Therefore, and
this must be emphasized, at the level of the strain the amplitude of the long-
lived component is larger than the amplitude of the bursts for all
$\ell\lesssim 6$ modes (and larger than the prompt response in all considered
modes). In turn, the amplitude of the bursts is larger than the prompt
response for $\ell\gtrsim 4$ only.

Lastly, the total gravitational wave energy $E_{<}(t)$ radiated through a
sphere of radius of $r_{\rm extr.}$ up to time $t$ for the $C=0.38$ binary is
$E_{<}(t_{0}+\tau_{\rm burst})/M_{0}\approx 6\times 10^{-6}$ and
$E_{<}(t_{0}+5\tau_{\rm burst})/M_{0}=6\times 10^{-5}$. The energy $E_{<}(t)$
still increases roughly linearly with rate $\dot{E}_{<}\approx 8\times
10^{-7}$ at this time. In contrast, the total radiated energy of the head-on
collision of a binary black hole with the same mass-ratio is
$\lim_{t\rightarrow\infty}E_{<}(t)=2\times 10^{-6}M_{0}$ Sperhake  _et al._
(2011).

##  IV Discussion and conclusion

In this work, we performed the first nonlinear and self-consistent study of
the ringdown of a black hole mimicker. To that end, we analyzed the dynamics
of spinning boson star remnants with stable light rings formed from large-
mass-ratio head-on collisions of binary boson stars. We broadly confirm the
appearence of features in the emitted gravitational waveform expected based on
test-field approaches: the prompt response of the remnant is consistent with
that of a Kerr black hole of the same mass and spin and the subsequent
emissions contain burst-like features. However, these bursts are sourced by
perturbations of the mimicker’s internal degrees of freedom, rather than a
propagation effect of massless waves through the mimicker’s interior.
Furthermore, we also identify a large-amplitude and long-lived emission
component in frequency comparable to the quasi-normal mode frequencies of the
Kerr spacetime that is indicative of excited trapped modes.

Although, the sequence of remnants, considered in this work, becomes
increasingly more compact (and e.g., the orbital frequency of the unstable
light ring asymptotes towards the Kerr value; see Appendix A), there are no
polar light spheres, and the primary’s exterior is (by construction) not
parametrically close to the Kerr exterior (due to the lack of a Birkhoff-like
theorem Birkhoff and Langer (1923); Jebsen (1921)). Therefore, it is
remarkable that the waveforms in Figure 1 match at all, and in principle, any
difference can be attributed to the differing exteriors (scalar interactions
may also impact the prompt ringdown waveform Palenzuela  _et al._ (2007);
Evstafyeva  _et al._ (2023); Siemonsen and East (2023a)). In analogy to the
black hole case Davis  _et al._ (1971); Sperhake  _et al._ (2011); East and
Pretorius (2013); East (2014), the small compactness of the secondary,
$C_{*}=0.02$, likely has small impact on the prompt response of the remnant in
low-$\ell$ modes.

Furthermore, we refrained from labeling the gravitational wave bursts as
“echoes”. Echoes—features of the scattering of massless waves in black hole
mimicker spacetimes neglecting any kind of backreaction—are expected to
exhibit frequencies and decay rates of the order of the Kerr quasi-normal
modes and appear after a light crossing time of the object. In this work, we
identified two distinct frequency components of the waveform emitted after the
prompt response. One of these is consistent with the quasi-normal mode
frequencies of a Kerr black hole of the same mass and spin, but long-lived,
while the other component is burst-like in nature, but with frequencies much
higher than those quasi-normal frequencies. Therefore, either component could
ambiguously be labelled “the echoes”. Since the bursts originate in
perturbations of the mimicker’s internal degrees of freedom and scale with the
size of the secondary, test-field approaches are generally insufficient to
capture this component. In fact, naively applying these perturbative methods,
one would not expect the waveforms to exhibit echoes, because the boson star’s
light-crossing time is comparable to Kerr quasi-normal mode decay timescales.

There are several interesting future directions: Our results may be
extrapolated to large remnant compactnesses to understand the rindown of
mimicker’s that are “closer” to black holes. This could also be used to inform
existing waveform templates and planned gravitational wave searches. A
systematic analysis of the quasi-normal mode spectrum of spinning boson stars
(as performed in the non-spinning case in Refs. Yoshida  _et al._ (1994);
Macedo  _et al._ (2013, 2016); Vásquez Flores  _et al._ (2019)), could guide
disentangling the two gravitational wave components. Lastly, we leave
performing a three-dimensional evolution of a quasi-circular large-mass-ratio
binary systems to future work.

###### Acknowledgements.

We would like to thank Niayesh Afshordi, Alejandro Cárdenas-Avendaño, Will
East, Suvendu Giri, Luis Lehner, Elisa Maggio, and Frans Pretorius for many
interesting discussions about aspects of this work. We especially thank Will
East and Frans Pretorius for comments on an earlier version of this draft. The
authors are pleased to acknowledge that the work reported on in this paper was
substantially performed using the Princeton Research Computing resources at
Princeton University which is a consortium of groups led by the Princeton In-
stitute for Computational Science and Engineering (PICSciE) and Office of
Information Technology’s Research Computing. This work used anvil at Purdue
University through allocation PHY230198 from the Advanced Cyberinfrastructure
Coordination Ecosystem: Services & Support (ACCESS) program Boerner  _et al._
(2023), which is supported by National Science Foundation grants #2138259,
#2138286, #2138307, #2137603, and #2138296.

##  Appendix A Boson stars

###  A.1 Parameter space

Figure 5: We show the masses and radii (measured with respect to the scalar
mass parameter $\mu$) of several $\tilde{m}=3$ families of rotating boson
stars in the scalar model with self-interaction strength $\sigma$ (see, eq.
(1)). In the decoupling limit $\sigma\rightarrow\infty$, the scalar potential
reduces to the mass term $\mu^{2}|\Phi|^{2}$. Regions of the parameter space
with compactness $C=M/R>1/2$ are labelled with “BH”. Two relevant radius
scalings are shown by dashed black lines in the respective regimes. Lastly, we
indicate the locations of those rotating stars considered in this work by
black squares and label each with their compactness (further properties of
these solutions can be found in Table 1).

The rotating boson stars relevant for this work are asymptotically flat
horizonless regular axisymmetric and stationary solutions to the Einstein-
Klein-Gordon equations associated with (1) Schunck and Mielke (1996); Volkov
and Wohnert (2002); Kleihaus  _et al._ (2005). As noted in the main text, the
full scalar field profile making up an isolated star takes the form
$\Phi=\phi(r,\theta)e^{i(\tilde{\omega}t+\tilde{m}\varphi)}$, with azimuthal
index $\tilde{m}$ and internal frequency $\tilde{\omega}$. The metric ansatz
for rotating solutions, $|\tilde{m}|\geq 1$, in Lewis-Papapetrou form is given
by

$\displaystyle\begin{aligned}
ds^{2}=-fdt^{2}&+\frac{l}{f}\bigg{[}g(dr^{2}+r^{2}d\theta^{2})\\\
&+r^{2}\sin^{2}\theta\Big{(}d\varphi-\frac{\Omega}{r}dt\Big{)}^{2}\bigg{]},\end{aligned}$
(2)

where $f,l,g$, and $\Omega$ depend on $r$ and $\theta$ only. Regularity on the
symmetry axis implies that the scalar field magnitude of rotating stars
vanishes there. At large distances, the magntiude $\phi$ decays exponentially
$\sim e^{-br}$, for some $b>0$. The mass $M$ and angular moment $J$ of these
stationary and axisymmetric spacetimes are given by their corresponding Komar
expressions. We define the radius of a given solution to be the circular
coordinate radius $R$ at which $99\%$ of the mass of the solutions lies at
$r<R$ Siemonsen and East (2021).

In Figure 5, we show the masses and radii of $\tilde{m}=3$ families of
spinning boson star solutions in scalar models with varying coupling $\sigma$.
In all cases, the scalar self-interactions are irrelevant in the Newtonian
limit, where solution’s masses scale as $M\sim R^{-1}$ and the compactness is
small. Towards the relativistic regime, the properties of families in models
with different coupling strength $\sigma$ diverge. Importantly, if self-
interactions are neglected, i.e., $\sigma\rightarrow\infty$, then the
solution’s mass scales roughly as in the Newtonian limit, even in the
relativistic regime. In contrast, the energy density of those rotating
solutions in strongly coupled scalar models is dominated by surface tensions
in this thin-wall limit Lee (1987); Friedberg  _et al._ (1987); Lee and Pang
(1992), scaling as $M\sim R^{2}$. Understanding these scalings is important
for the choice of binary boson star to model the ringdown of a black hole
mimicker in the large-mass-ratio regime, as we point out in the next section.
Families of solutions with higher $\tilde{m}$ are more rapidly spinning,
supporting more compact stars. This is the primary reason for considering
$\tilde{m}=3$, rather than lower-$\tilde{m}$ solutions. One drawback of this
choice is that $|\tilde{m}|>1$ solutions are likely linearly unstable to a
non-axisymmetric instability Sanchis-Gual  _et al._ (2019) even in the
small-$\sigma$ regime Siemonsen and East (2021). In our evolution setup
(enforcing the metric’s axisymmetry and the scalar azimuthal symmetry),
however, this instability is removed. Therefore, as there are no ergoregion
present in the considered boson stars, all star solutions relevant in this
work are stable against any known linear instabilities (in our evolution
setup).

###  A.2 Choice of binary boson stars

We are interested in studying the merger and ringdown of black hole mimicking
objects using boson stars. Utilizing these solutions, however, to model
mimickers comes with two major limitations: (i) as these are entirely
classical solutions, the merger remnant of two of these objects collapses to a
black hole roughly when its mass surpasses the hoop bound Misner  _et al._
(1973) or is above the maximum mass of the family of boson star solutions, and
(ii), within a given scalar theory (at fixed $\mu$ and $\sigma$), sequences of
solutions follow the scaling $M\sim R^{p}$ with $p\neq 1$ (i.e., not the black
hole scaling $M\sim R$). The former implies that one likely cannot consider
equal-mass mergers of ultra compact boson stars, since their merger products
are always black holes or less compact stars (see e.g., Ref. Cardoso  _et al._
(2016b) or Fig. 6 in Ref. Siemonsen and East (2023b)). Hence, in order to
circumvent this issue, we consider the large-mass-ratio regime, in which a
primary ultra compact star is perturbed by a lighter secondary star. We find
that a star of compactness $C\approx 0.39$, in the same sequence of solutions
detailed in Table 1, collapses to a black hole shortly after merger with the
$C_{*}=0.02$ secondary considered throughout this work.

Issue (ii) implies that, at fixed mass-ratio, only one of the binary’s
constituents may be ultra compact. In this work, we are interested in
constructing binaries, which model the merger of two ultra compact objects in
the large-mass-ratio regime the “closest”. Hence, ideally one would like to
use two objects of large compactness each; this, however, is not achievable in
the context of boson stars due to the restrictions shown in Figure 5 (i.e.,
nowhere does the boson star mass scale as $\sim R$). Let us consider binaries
built from stars with the Newtonian scaling, $M\sim R^{-1}$, in the large-
mass-ratio regime $M_{1}/M_{2}\gg 1$. Then the size of the lighter companion
star would be larger than the primary highly-compact star, since
$R_{1}/R_{2}\sim M_{2}/M_{1}\ll 1$. In this regime, even if the primary were a
black hole, the emitted gravitational radiation is not expected to exhibit a
black hole ringdown East (2014). Therefore, we focus on those solutions, which
are in the thin-wall regime of the parameter space instead. As pointed out
above, there the masses along a family of solutions scale as $M\sim R^{2}$,
which is much closer to the black hole $M\sim R$ scaling behavior. Ultimately,
a scalar potential yielding boson star solutions with $M\sim R$ in the
relativistic regime would be most relevant (see Ref. Pitz and Schaffner-
Bielich (2023) for a systematic study in the non-spinning case).

###  A.3 Null geodesics

####  A.3.1 Equatorial light rings

In order to understand the propagation of null geodesics through the boson
star spacetimes considered in this work, we begin by analyzing the relevant
potentials for radial motion of the latter confined to the equatorial plane.
To that end, let $\lambda$ be an affine parameter of the trajectory
$x^{\mu}(\lambda)$ of the null geodesic with tangent
$\dot{x}^{\mu}(\lambda)=dx^{\mu}/d\lambda=u^{\mu}(\lambda)$, such that
$u_{\mu}u^{\mu}=0$. The conserved energy, $E=-t_{\mu}u^{\mu}$, and angular
momentum, $L_{z}=\varphi_{\mu}u^{\mu}$, of the geodesic are measured with
respect to the Killing fields associated with stationarity and axisymmetry of
the isolated boson star spacetimes, $t^{\mu}$ and $\varphi^{\mu}$,
respectively. The normalization condition $u_{\mu}u^{\mu}=0$ yields

$\displaystyle\dot{r}^{2}=-\frac{f}{gl}V,$ $\displaystyle V=-\frac{\kappa
L_{z}^{2}}{f}(\kappa-H_{+})(\kappa-H_{-}),$ (3)

where the impact parameter is $\kappa=E/L_{z}$, and the potentials $H_{\pm}$
are Cunha  _et al._ (2016, 2017)

$\displaystyle H_{\pm}=\frac{\Omega}{r}\pm\frac{f}{r\sin\theta\sqrt{l}},$ (4)

in the Lewis-Papapetrou coordinates chosen above. Light rings inside the
equatorial plane ($\theta=\pi/2$ and $\dot{\theta}=0$) are defined by $V=0$
and $\partial_{r}V=0$ (the latter is equivalent to $\partial_{r}H_{\pm}=0$).
Furthermore, roots of $\partial_{r}H_{+}$ ($\partial_{r}H_{-}$) correspond to
co-(counter-)rotating light rings, and these are stable if
$\pm\partial_{r}^{2}H_{\pm}>0$ at the location of the light ring (or unstable
if $\pm\partial_{r}^{2}H_{\pm}<0$).

Figure 6: The potentials $H_{\pm}$, defined in (4), characterizing the
structure of light rings in the equatorial plane of a Kerr black hole with
spin $a/M=0.95$ (labelled “Kerr”), and the sequence of spinning boson stars of
increasing compactness considered throughout this work (labelled by their
compactness). The radial coordiante $\tilde{r}$ is the coordinate $r$ in the
case of boson stars and $r-r_{H}$ in the case of Kerr [in the coordinates of
(2)], where $r_{H}$ is the radial location of the horizon.

In Figure 6, we compare the potentials $H_{\pm}$ of a Kerr black hole with
spin $a/M=0.95$ with those of the sequence of stars of increasing compactness
detailed in Table 1. There exists a co- and counter rotating (unstable
equatorial) light ring around Kerr black holes. In this large-spin regime,
these light rings are located at $r/M\approx 1$ and $r/M\approx 4$
(corresponding to $\tilde{r}/M\approx 0$ and $3$ in Figure 6), respectively.
Turning now to the sequence of spinnning boson stars, all solutions in the
sequence exhibit a stable and unstable counter-rotating light ring in the
equatorial plane (and no co-rotating light ring). The orbital frequencies
$d\varphi/dt=\tilde{\Omega}$ of the stable light rings in the boson star
spacetimes are $M\tilde{\Omega}_{S}\times 10^{2}=2,5,9,11$, as well as
$M\tilde{\Omega}_{U}=0.1437,0.142,0.138,0.132$ for the unstable light rings,
in decreasing star compactness $C=0.38,0.36,0.32,0.28$, respectively.
Interestingly, the orbital frequency of the unstable light ring in a Kerr
spacetime of spin $a/M=0.95$ is $M\tilde{\Omega}_{\rm Kerr}=0.1446$.
Therefore, the relative difference of the unstable light ring’s frequency
between the $C=0.38$ rotating boson star and the Kerr spacetime with the same
angular momentum is $<1\%$. The absence of a co-rotating light ring may be
understood intuitively as follows: Along the considered sequence of rotating
boson stars, and with increasing compactness, the outer light ring (i.e., the
counter-rotating ring) appears first, whereas the inner, co-rotating, light
ring may appear only for more compact solutions.

####  A.3.2 Polar null orbits

In the case of ultra compact and spherically symmetric boson stars Palenzuela
_et al._ (2017); Bošković and Barausse (2022); Collodel and Doneva (2022) the
presence of equatorial light rings implies the existence of complete photon
spheres. As rotating boson stars are not connected to their non-rotating
counterparts by continuously increasing the angular momentum Kobayashi  _et
al._ (1994); Kleihaus  _et al._ (2005), the existence of polar null orbits,
i.e., those bound null geodesics with $L_{z}=0$, is not implied by the
presence of equatorial light rings. In the following, we briefly investigate
the existence and properties of the bound polar null geodesics in the sequence
of rotating boson stars (and beyond) presented in Table 1.

Figure 7: The polar and radial coordinates $\theta$ and $r$, respectively, of
the polar null geodesic, which maximizes $\tau_{\rm polar}$. The boson stars
are identified by their compactness $C$. Those with $C\leq 0.38$ are listed in
Table 1, and have been evolved in the context of binaries in the bulk of this
work. Those rotating stars with $C\geq 0.39$ are slightly more compact
representatives of the same family of solutions, which have frequencies
$\tilde{\omega}/\mu=0.101$ and $\tilde{\omega}/\mu=0.095$. Furthermore, dashed
lines indicate the coordinate radius of each star $r/M=1/C$ and the black
dotted (solid) line, corresponds to the outer (inner) closed null polar
geodesic of the $C=0.41$ spacetime.

Due to the lack of a Killing tensor, and associated constant of motion, we
proceed by numerically solving the geodesic equation to find closed polar null
geodesics (with further details on the numerical implementation deferred to
Appendix B). We integrate starting on the spin axis (with vanishing initial
radial velocity) of each of the stars listed in Table 1. This is performed for
various initial positions along the symmetry axis. We define $\tau_{\rm
polar}$ as the coordinate time such a geodesic spends inside the star,
$r<M/C$. More generally, and for later convenience, this coordinate time is
defined as

$\displaystyle\tau=\int_{\gamma}\frac{dt}{d\lambda}d\lambda.$ (5)

Then for those initial positions along the symmetry axis approaching a bound
null geodesic the time $\tau_{\rm polar}$ diverges. If a bound orbit exists,
then iteratively increasing this time results in the convergence of the
initial position along the symmetry axis towards that orbit. If no bound polar
null geodesics exist, this proceedure yields the geodesic with largest
$\tau_{\rm polar}$.

In Figure 7, we show those polar null geodesics, which maximize $\tau_{\rm
polar}$. The time $\tau_{\rm polar}$ increases along the sequence of boson
stars with increasing compactness. In the case of the $C=0.38$ star, the polar
null geodesic maximizing $\tau_{\rm polar}$ passes through the spin axis twice
before escaping to infinity. The sequence of stars listed in Table 1 are
elements of a family of solutions with more compact stars (see also Figure 5).
Continuing further along this family of solutions, the $C=0.39$ boson star
exhibits at least one polar light ring222With our numerical methods we can
confidently identify only the existence of a single polar orbit. It is
plausible that two distinct orbits exit, which we are unable to distinguish
due to their relative proximity., whereas the one with $C=0.41$ has two
distinct bound polar null orbits. This demonstrates that the family of
solutions described in Table 1 develops bound polar null orbits in the high-
compactness limit. On the other hand, those rotating boson stars evolved
within binaries in this work, exhibit only “quasi”-bound null orbits, with
large, but finite, $\tau_{\rm polar}$. Explicitly, the coordinate times are
$\tau_{\rm polar}/M=28.8,17.5,11.7,9.8$ for the stars of compactness
$C=0.38,0.36,0.32,0.28$, respectively.

####  A.3.3 Light-crossing times

As we saw in Figure 2, the high-frequency perturbations sourced during the
merger of the binary propagate through the remnant star’s interior along the
symmetry axis, partially reflect on the opposite side and propagate backwards
both along the symmetry axis and along the surface of the remnant. In the
eikonal limit, massless perturbations propagate along null geodesics through
the spacetime Goebel (1972); Isaacson (1968). Within our numerical setup, we
enforce axisymmetry of the metric and a $\tilde{m}=3$ azimuthal symmetry for
the scalar field333Recall, this corresponds to
$\mathcal{L}_{k}\Phi=i\tilde{m}\Phi$ with $\tilde{m}=3$ and the Lie derivative
$\mathcal{L}_{k}$ along the Killing vector $k^{\mu}$ generating the
spacetime’s axisymmetry.. Therefore, geodesics with vanishing angular
momentum, $L_{z}=0$, are particularly important to our discussion, as these
correspond to the propagation of high-frequency and $m=0$ massless modes (the
equatorial light rings correspond to modes with $m\gg 1$).

Accurately describing the wave dynamics in a merger setting as shown in Figure
2 is non-trivial. Therefore, in the following we determine the light crossing
times of the interior of the sequence of rotating boson stars listed in Table
1 along $L_{z}=0$ null geodesics simply as a rough measure of the time delay
of the gravitational wave bursts emerging during the ringdown shown in Figure
3. To that end, we focus on three classes of such geodesics: (I) “quasi”-bound
polar null geodesics (as defined above), (II) equatorial null geodesics
propagating from the unstable light ring through the spin axis to the opposite
side, and (III) massless mode propagation along the symmetry axis. In all
cases, we use (5) to determine the light crossing times with respect to
observers at infinity.

Figure 8: The light crossing time $\tau_{\rm axis}$ as a function of
start/endpoint $|z_{s}|$ along the symmetry axis for all stationary
ultracompact solutions considered in this work. The coordinate radii
$z_{s}/M=1/C$ are indicated by vertical dashed lines.

The crossing times of geodesics of type I were determined in the previous
section to be $\tau_{\rm polar}\sim\mathcal{O}(10)M$. The crossing time of
type-II or those of geodesics within the equatorial plane emerging from the
unstable null orbit (as found in Figure 6), passing through the symmetry axis,
and ending on the unstable null orbit (on the opposite side of the object).
These times are $\tau_{\rm equa.}/M=12.2,11.8,10.1$, and $9.2$, corresponding
to the stars of compactness $C=0.38,0.36,0.32$, and $0.28$, respectively.
Note, the light-crossing time through the interior increases with decreasing
radius $R=M/C$. Lastly, the light crossing time $\tau_{\rm axis}$ associated
with geodesics of category III are less naturally bounded. Therefore, we
determine $\tau_{\rm axis}$ for geodesics along the spin axis starting at
various coordinates $z_{s}>0$ above the equatorial plane and terminating at
$-z_{s}$. In Figure 8, we show the dependence of $\tau_{\rm axis}$ on this
choice of $z_{s}$. Those values shown in Figure 4 correspond to the choice
$\tau_{\rm axis}(z_{s}=M/C)$.

##  Appendix B Numerical implementation & uncertainties

###  B.1 Numerical evolution methods & convergence

The isolated rotating boson stars making up the binaries were obtained using
methods developed in Refs. Siemonsen and East (2021). Specifically, we use
Newton-Raphson relaxation techniques in conjunction with fifth-order finite
difference methods to solve the set of two-dimensional elliptic partial
differential equations emerging, when plugging the scalar field ansatz and (2)
into the Einstein-Klein-Gordon equations derived from (1). The equations are
discretized uniformly in the polar coordinate $\theta$ and compactified radial
coordinate $\bar{r}=r/(1+r)\in(0,1)$. Obtaining boson star solutions in the
thin-wall regime is challenging due to large spatial gradients developing at
the star’s surface. We overcome some of these issues by rescaling the radial
coordinate $\bar{r}$ with lengthscale $L$, i.e., $\bar{r}\rightarrow
L\bar{r}$, as well as all other dimensionful quantities, in order for the
star’s surface to reside (roughly) at $\bar{r}\approx 0.5$. This ensures that
the gradients at the star’s surface are best-resolved at fixed resolution. The
resolutions used for all boson stars in the radial and polar directions are
$N_{\bar{r}}\times N_{\theta}=1500\times 100$. For all highly compact binary
constituents listed in Table 1, we use $L=0.025$. Note, in the case of the
$C_{*}=0.02$ secondary companion, $L=1$ is entirely sufficient.

These isolated solutions are combined into binaries using techniques developed
in Refs. Siemonsen and East (2023b); East  _et al._ (2012a). In particular, we
solve for constraint satisfying binary boson star initial data using the
conformal thin-sandwich formulation of the Hamiltonian and momentum
constraints of the Einstein equations. The necessary free data is obtained by
superposing two displaced and boosted star solutions. The scalar kinetic
energy of the binary (entering the constraint equations) is conformally
rescaled with power $p=-3$, as defined in detail in Ref. Siemonsen and East
(2023b), to reduce spurious oscillations of the stars in the subsequent
evolution of the initial data. Note, for the black hole-boson star binary, we
use plain superposed initial data without solving the constraints explicitly.
In this latter case, we verified that the total charge $Q$ [associated with
the global U(1) symmetry of (1)] of the initial data agrees with the total
charge of the isolated secondary to within $<1\%$ relative difference (and is
stable throughout the entire evolution until merger). The black hole’s
properties differ from those parameters initialized in the initial data only
by $<0.1\%$ throughout the evolution up until the merger.

We proceed to evolve these initial data with the same methods as used in Refs.
Siemonsen and East (2021, 2023a). We numerically solve the full Einstein-
Klein-Gordon equations of motion, derived from (1), using the generalized
harmonic formulation Pretorius (2005). The equations are discretized using
fourth-order accurate finite differences in conjunction with fourth-order
Runge-Kutta time stepping. All evolutions are performed assuming axisymmetry
for the metric, and azimuthal symmetry for the scalar field (as defined
above), which is achieved numerically by means of a generalized Cartoon method
Pretorius (2005); Alcubierre  _et al._ (2001). Essential is the role of the
adaptive mesh refinement (AMR) with refinement ratio 2:1 East  _et al._
(2012b). Due to the large separation of scales of the problem, we employ nine
mesh refinement levels with a medium spatial resolution of $\Delta x_{\rm
finest}/M_{0}\approx 3\times 10^{-3}$ on the finest level (covering both the
primary and secondary of each binary). Crucially, due to the high-frequency
nature of the emitted burst-like gravitational waves after the merger of the
binary, we ensure the wave extraction zone (out to radial coordinate distances
$r/M_{0}\approx 133$) is well-resolved with grid spacing of $\Delta x_{\rm
wavezone}/M_{0}\leq 0.025$ at medium resolution. For all evolutions, we
utilize constraint damping terms with damping rate $\kappa_{D}M_{0}=3$ and
damping constant $\rho_{D}=0.75$ Gundlach  _et al._ (2005). For all boson star
binaries, we use stationary gauge throughout the evolution to minimize
spurious gauge dynamics. For the black hole-boson star we employ the damped
harmonic gauge condition Choptuik and Pretorius (2010); Lindblom and Szilagyi
(2009).

Figure 9: The convergence behavior of constraint violations, measured by
$\max\mathcal{C}$ and $\mathcal{I}$ (defined in the main text), throughout the
evolution of the binary $C=0.38$ with decreasing grid spacing from $h$ at low
resolutions to $h/2$ at high resolution. (top) The global maximum of the
scalar field magnitude throughout the plunge and merger of the binary (all
curves are aligned at the time of merger $t^{\prime}_{0}$). The evolution of
the global maximum $\max\mathcal{C}$ (center row) and the integrated norm
$\mathcal{I}$ (bottom row) of the constraint violations during the initial
phase before the plunge (left column), as well as during and after the merger
(right column). Note, the merger occurs at $t/M_{0}\approx 10^{3}$.

To validate our findings, we perform a convergence study on the boson star
binary with the most compact constituent (i.e., the $C=0.38$ binary), which
exhibits the largest separation of scales and most extreme remnant. To that
end, we use the supremiums norm $\max\mathcal{C}$ and integrated norm
$\mathcal{I}=\int d^{3}x\sqrt{\gamma}\mathcal{C}/M_{0}^{3}$ (integrated over a
ball of coordinate radius $64M_{0}$ centered around the center of mass), with
$\mathcal{C}=\sum_{\mu}|H_{\mu}-\square x_{\mu}|M_{0}/4$, to track the
constraint violations throughout the evolution. In Figure 9, these quantities
are shown as a function of time at low, medium (default) and high resolutions.
Early on in the evolution, the constraint violations converge towards zero at
convergence order between third and the expected fourth order. Note, the
convergence of $\max\mathcal{C}$ slows towards higher resolutions. This is
primarily due to the resolution limitations to solve for the isolated solution
of the ultra compact constituent as pointed out above. At merger, $t\approx
t^{\prime}_{0}$, short-wavelength modes are excited in the remnant resulting
in worse convergence behavior (i.e., the wavelength of some of these modes
approaches the grid scale) in both $\mathcal{I}$ and $\max\mathcal{C}$. At
this time, the integrated norm converges roughly at first order between the
low and medium resolution runs (zeroth order for the supremums norm), and at
second order between medium and high resolution cases (first order for the
supremums norm). This suggests that, likely driven by high-frequency modes,
the low resolution is not yet in the converging regime. We use the medium
resolution as the default resolution for all binaries considered in this work.
In the following section, we discuss the impact of resolution on the extract
gravitational waveforms.

###  B.2 Gravitational waveform uncertainties

To understand the numerical uncertainties associated with observables
presented in the bulk of this work, we estimate the impact of different
sources of error on the extracted gravitational waveform. In particular, due
to the largest mass- and size-ratio the $C=0.38$ binary configuration is the
most challenging to treat within our numerical setup. Hence, we focus on the
uncertainties of this binary merger as a representative example for all
considered configurations.

A major source of contamination of the emitted gravitational waves originates
from spurious perturbations introduced by the initial data. As shown in detail
in Ref. Siemonsen and East (2023b), even small amplitude long-lived
oscillations induced in each star can lead to the emission of spurious
radiation. The prompt response shown in Figure 1, is particularly prone to
these contaminations, since its amplitude scales linearly in the binary’s
mass-ratio (and $q=78$ in the most relativistic considered binary; see Table
1). The spurious gravitational wave emissions themselves scale as $\sim
D^{-1}$ with the initial coordinate separation. Therefore, increasing the
initial binary separation (in addition to the conformal rescaling as outlined
above, and with details in Ref. Siemonsen and East (2023b)) reduces spurious
gravitational wave contaminations. This is the primary reason for considering
such a large initial separation of $D/M_{0}=160$. Note, due to the scalar
field’s time dependence, comparing the gravitational radiation from the merger
of identical binary boson stars with varying initial separation is non-
trivial. Scalar interactions, active during the merger and ringdown, are
sensitive to the phase-offset of the boson stars at the point of contact,
which in general depends on the initial separation. We leave a study of these
effects to future work, and hence, focus entirely on the separation
$D/M_{0}=160$.

Figure 10: The convergence behavior of the spin-weighted spherical harmonic
components $\tilde{\Psi}_{4}^{\ell}=r_{\rm extr.}M_{0}\Psi_{4}^{\ell}$,
extracted at a coordinate sphere of radius $r_{\rm extr.}/M_{0}=100$, with
grid spacing $h$, corresponding to the prompt response of the $C=0.38$ binary
boson star (shown in Figure 1). All polar components are aligned, such that
the first peaks of the $\ell=2$ modes align in time (as before, $t_{0}$
indicates the end time of the prompt response).

The gravitational waveforms are obtained by evaluating the Newman-Penrose
scalar $\Psi_{4}$ on coordinate spheres with radii $r_{\rm
extr.}/M_{0}\in\\{50,76,100\\}$. We check explicitly that finite extraction
radius effects are smaller than the uncertainty due to the discretization of
the evolution grid. To quantify the latter, we compare the waveforms of the
$C=0.38$ binary boson star merger at three different resolutions at $r_{\rm
extr.}/M_{0}=100$. To that end, in Figure 10 we show the convergence behavior
of the prompt response of the remnant produced in this binary merger with
increasing resolution. While the lowest resolution differs significantly, both
the medium and high resolution waveforms align relatively well for all polar
modes. In the case of $\ell=8$, there is a relative timeshift of the waveform
of $\approx 0.6M_{0}$ between the medium and high resolutions. Since the
alignment of these waveforms is of importance for the interpretation of the
results, we present the medium-resolution waveforms in both Figure 1 and
Figure 3, with the exception of the central panel in Figure 1, where we
present the high-resolution $\ell=8$ waveform; note, non of the conclusions
drawn in the main text are impacted by this choice, as it implies merely a
timeshift of the waveform by the aforementioned $\approx 0.6M_{0}$ (i.e., no
change in frequency or decay rates). The uncertainties shown in Figure 1 are
obtained from a time-averaged difference between the medium and high
resolution waveforms shown in Figure 10.

Figure 11: As in Figure 10, however, now focusing on the convergence of the
waveform emitted after the prompt ringdown of the $C=0.38$ remnant (these
correspond to those waveforms presented in Figure 3 in the main text).

In Figure 11, we show the convergence behavior of the gravitational waveforms
after the prompt response. Analogously to the discussion above, the lowest
resolution waveform deviates significantly, while both the medium and high
resolution gravitational waveforms align reasonably well. Among others, the
burst’s frequency, the burst period, burst amplitude, the amplitude of the
long-lived gravitational wave component, as well as the frequencies of this
long-lived component are reasonably converged in both the $\ell=2$ and
$\ell=8$ polar modes. Naively, one may worry that the bursts in the waveform
originate (at least partially) from the reflection of outgoing waves off AMR
boundaries. However, for one in Figure 2, we explicitly present the physical
origin of these bursts; to double-check, we compared the burst waveforms at
different extraction radii finding good agreement of $\tau_{\rm burst}$ (i.e.,
if bursts originated from such reflections, then $\tau_{\rm burst}$ would
strongly depend on the extraction radius). The first burst’s frequency $f_{\rm
burst}$ is obtained from the centers of the peak of the Fourier transform of
the $\ell=8$ waveform shown for $C=0.38$ in Figure 3 (though the same result
is obtained from a fit to $\ell=2$). The burst period $\tau_{\rm burst}$ is
the time $t$ between the peak of each burst (if there are more than three
bursts, we use the average of the times between bursts). The star’s
oscillation timescale $\tau_{\rm osc.}$ is obtained from fitting
$\sin(at+b)+c$ to post-merger oscillations of $\max|\Phi|$ (as shown in Figure
9) for all four binaries considered in this work. Finally, the frequency
$f_{\rm low}$ is obtained directly from Fourier transform of the type shown in
Figure 3 (for $C=0.28$ and $C=0.32$, the uncertainties of this method increase
compared with the more compact binaries). Overall, we estimate the
uncertainties of these timescales, shown in Figure 4, to be at the $\lesssim
10\%$ level. The uncertainties of the computed total radiated energy in
gravitational waves we estimate to be $\approx 20\%$.

###  B.3 Geodesic integrator

To determine the light crossing times along various classes of null geodesics
in the stationary boson star spacetimes, we solve the geodesic equations
numerically. The latter equation is given by

$\displaystyle\frac{dx^{\mu}}{d\lambda}=u^{\mu},$
$\displaystyle\frac{du^{\mu}}{d\lambda}=-\Gamma^{\mu}_{\nu\gamma}u^{\nu}u^{\gamma},$
(6)

with Christoffel symbol $\Gamma^{\mu}_{\nu\gamma}$ associated with the boson
star spacetime $g_{\mu\nu}$. Following from the stationarity and axisymmetry
of $g_{\mu\nu}$, we have defined the conserved energy $E$ and angular momentum
$L_{z}$ of the geodesics above. Together with the normalization condition
$u_{\mu}u^{\mu}=0$, these are conserved along the geodesic. Initial conditions
are provided by choosing an initial location $x^{i}(0)$ and spatial velocity
$u^{i}(0)$. The temporal component $u^{0}(0)$ is then obtained from the
normalization condition, which fixes the parameterization. The geodesic
equations (6) are solved numerically using a fourth-order accurate Runge-Kutte
method. We find it convenient to perform the integration on a Cartesian
spatial grid, applying the Euclidean spherical to Cartesian coordinate
transformations to the Lewis-Papapetrou coordinates of (2). To determine
$\Gamma^{\mu}_{\nu\gamma}$, we use fifth-order accurate finite differences for
the metric derivatives in the $r$ and $\theta$ directions, and employ linear
order interpolations to obtain values of $f,l,g$ and $\Omega$ at an arbitrary
location in $(r,\theta)\in(0,\infty)\times(0,\pi)$. During the integration of
(6), we monitor the constancy of $E$ and $L_{z}$ and assess the convergence
based on the violations of the constraint $u^{\mu}u_{\mu}=0$ along the
geodesic. The convergence of the latter is entirely limited by the resolution
of the background spacetime (with grid spacing $\Delta r$ and $\Delta\theta$),
as the typical stepsize of the Runge-Kutta integrator satisfies
$\Delta\lambda\ll\Delta r,\Delta\theta$.

## References

  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Cardoso  _et al._ (2016a) V. Cardoso, E. Franzin, and P. Pani, Phys. Rev. Lett. 116, 171101 (2016a), [Erratum: Phys.Rev.Lett. 117, 089902 (2016)], arXiv:1602.07309 [gr-qc] .
  * Cardoso  _et al._ (2016b) V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Phys. Rev. D 94, 084031 (2016b), arXiv:1608.08637 [gr-qc] .
  * Cardoso and Pani (2017) V. Cardoso and P. Pani, Nature Astron. 1, 586 (2017), arXiv:1709.01525 [gr-qc] .
  * Holdom and Ren (2017) B. Holdom and J. Ren, Phys. Rev. D 95, 084034 (2017), arXiv:1612.04889 [gr-qc] .
  * Bueno  _et al._ (2018) P. Bueno, P. A. Cano, F. Goelen, T. Hertog, and B. Vercnocke, Phys. Rev. D 97, 024040 (2018), arXiv:1711.00391 [gr-qc] .
  * Barceló  _et al._ (2017) C. Barceló, R. Carballo-Rubio, and L. J. Garay, JHEP 05, 054 (2017), arXiv:1701.09156 [gr-qc] .
  * Urbano and Veermäe (2019) A. Urbano and H. Veermäe, JCAP 04, 011 (2019), arXiv:1810.07137 [gr-qc] .
  * Raposo  _et al._ (2019) G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, Phys. Rev. D 99, 104072 (2019), arXiv:1811.07917 [gr-qc] .
  * Pani and Ferrari (2018) P. Pani and V. Ferrari, Class. Quant. Grav. 35, 15LT01 (2018), arXiv:1804.01444 [gr-qc] .
  * Cardoso  _et al._ (2019) V. Cardoso, V. F. Foit, and M. Kleban, JCAP 08, 006 (2019), arXiv:1902.10164 [hep-th] .
  * Oshita and Afshordi (2019) N. Oshita and N. Afshordi, Phys. Rev. D 99, 044002 (2019), arXiv:1807.10287 [gr-qc] .
  * Wang  _et al._ (2020) Q. Wang, N. Oshita, and N. Afshordi, Phys. Rev. D 101, 024031 (2020), arXiv:1905.00446 [gr-qc] .
  * Oshita  _et al._ (2020) N. Oshita, Q. Wang, and N. Afshordi, JCAP 04, 016 (2020), arXiv:1905.00464 [hep-th] .
  * Maggio  _et al._ (2020) E. Maggio, L. Buoninfante, A. Mazumdar, and P. Pani, Phys. Rev. D 102, 064053 (2020), arXiv:2006.14628 [gr-qc] .
  * Dey  _et al._ (2020) R. Dey, S. Chakraborty, and N. Afshordi, Phys. Rev. D 101, 104014 (2020), arXiv:2001.01301 [gr-qc] .
  * Ikeda  _et al._ (2021) T. Ikeda, M. Bianchi, D. Consoli, A. Grillo, J. F. Morales, P. Pani, and G. Raposo, Phys. Rev. D 104, 066021 (2021), arXiv:2103.10960 [gr-qc] .
  * Mark  _et al._ (2017) Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, Phys. Rev. D 96, 084002 (2017), arXiv:1706.06155 [gr-qc] .
  * Nakano  _et al._ (2017) H. Nakano, N. Sago, H. Tagoshi, and T. Tanaka, PTEP 2017, 071E01 (2017), arXiv:1704.07175 [gr-qc] .
  * Wang and Afshordi (2018) Q. Wang and N. Afshordi, Phys. Rev. D 97, 124044 (2018), arXiv:1803.02845 [gr-qc] .
  * Burgess  _et al._ (2018) C. P. Burgess, R. Plestid, and M. Rummel, JHEP 09, 113 (2018), arXiv:1808.00847 [gr-qc] .
  * Correia and Cardoso (2018) M. R. Correia and V. Cardoso, Phys. Rev. D 97, 084030 (2018), arXiv:1802.07735 [gr-qc] .
  * Longo Micchi and Chirenti (2020) L. F. Longo Micchi and C. Chirenti, Phys. Rev. D 101, 084010 (2020), arXiv:1912.05419 [gr-qc] .
  * Maggio  _et al._ (2019) E. Maggio, A. Testa, S. Bhagwat, and P. Pani, Phys. Rev. D 100, 064056 (2019), arXiv:1907.03091 [gr-qc] .
  * Srivastava and Chen (2021) M. Srivastava and Y. Chen, Phys. Rev. D 104, 104006 (2021), arXiv:2108.01329 [gr-qc] .
  * Annulli  _et al._ (2022) L. Annulli, V. Cardoso, and L. Gualtieri, Class. Quant. Grav. 39, 105005 (2022), arXiv:2104.11236 [gr-qc] .
  * Xin  _et al._ (2021) S. Xin, B. Chen, R. K. L. Lo, L. Sun, W.-B. Han, X. Zhong, M. Srivastava, S. Ma, Q. Wang, and Y. Chen, Phys. Rev. D 104, 104005 (2021), arXiv:2105.12313 [gr-qc] .
  * Ma  _et al._ (2022) S. Ma, Q. Wang, N. Deppe, F. Hébert, L. E. Kidder, J. Moxon, W. Throwe, N. L. Vu, M. A. Scheel, and Y. Chen, Phys. Rev. D 105, 104007 (2022), arXiv:2203.03174 [gr-qc] .
  * Maselli  _et al._ (2017) A. Maselli, S. H. Völkel, and K. D. Kokkotas, Phys. Rev. D 96, 064045 (2017), arXiv:1708.02217 [gr-qc] .
  * Testa and Pani (2018) A. Testa and P. Pani, Phys. Rev. D 98, 044018 (2018), arXiv:1806.04253 [gr-qc] .
  * Tsang  _et al._ (2018) K. W. Tsang, M. Rollier, A. Ghosh, A. Samajdar, M. Agathos, K. Chatziioannou, V. Cardoso, G. Khanna, and C. Van Den Broeck, Phys. Rev. D 98, 024023 (2018), arXiv:1804.04877 [gr-qc] .
  * Longo Micchi  _et al._ (2021) L. F. Longo Micchi, N. Afshordi, and C. Chirenti, Phys. Rev. D 103, 044028 (2021), arXiv:2010.14578 [gr-qc] .
  * Abedi  _et al._ (2017) J. Abedi, H. Dykaar, and N. Afshordi, Phys. Rev. D 96, 082004 (2017), arXiv:1612.00266 [gr-qc] .
  * Ashton  _et al._ (2016) G. Ashton, O. Birnholtz, M. Cabero, C. Capano, T. Dent, B. Krishnan, G. D. Meadors, A. B. Nielsen, A. Nitz, and J. Westerweck, (2016), arXiv:1612.05625 [gr-qc] .
  * Westerweck  _et al._ (2018) J. Westerweck, A. Nielsen, O. Fischer-Birnholtz, M. Cabero, C. Capano, T. Dent, B. Krishnan, G. Meadors, and A. H. Nitz, Phys. Rev. D 97, 124037 (2018), arXiv:1712.09966 [gr-qc] .
  * Conklin  _et al._ (2018) R. S. Conklin, B. Holdom, and J. Ren, Phys. Rev. D 98, 044021 (2018), arXiv:1712.06517 [gr-qc] .
  * Nielsen  _et al._ (2019) A. B. Nielsen, C. D. Capano, O. Birnholtz, and J. Westerweck, Phys. Rev. D 99, 104012 (2019), arXiv:1811.04904 [gr-qc] .
  * Lo  _et al._ (2019) R. K. L. Lo, T. G. F. Li, and A. J. Weinstein, Phys. Rev. D 99, 084052 (2019), arXiv:1811.07431 [gr-qc] .
  * Tsang  _et al._ (2020) K. W. Tsang, A. Ghosh, A. Samajdar, K. Chatziioannou, S. Mastrogiovanni, M. Agathos, and C. Van Den Broeck, Phys. Rev. D 101, 064012 (2020), arXiv:1906.11168 [gr-qc] .
  * Uchikata  _et al._ (2019) N. Uchikata, H. Nakano, T. Narikawa, N. Sago, H. Tagoshi, and T. Tanaka, Phys. Rev. D 100, 062006 (2019), arXiv:1906.00838 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2112.06861 [gr-qc] .
  * Uchikata  _et al._ (2023) N. Uchikata, T. Narikawa, H. Nakano, N. Sago, H. Tagoshi, and T. Tanaka, Phys. Rev. D 108, 104040 (2023), arXiv:2309.01894 [gr-qc] .
  * Abedi  _et al._ (2020) J. Abedi, N. Afshordi, N. Oshita, and Q. Wang, Universe 6, 43 (2020), arXiv:2001.09553 [gr-qc] .
  * Carballo-Rubio  _et al._ (2018) R. Carballo-Rubio, F. Di Filippo, S. Liberati, and M. Visser, Phys. Rev. D 98, 124009 (2018), arXiv:1809.08238 [gr-qc] .
  * Danielsson  _et al._ (2021) U. Danielsson, L. Lehner, and F. Pretorius, Phys. Rev. D 104, 124011 (2021), arXiv:2109.09814 [gr-qc] .
  * Vellucci  _et al._ (2023) V. Vellucci, E. Franzin, and S. Liberati, Phys. Rev. D 107, 044027 (2023), arXiv:2205.14170 [gr-qc] .
  * Dailey  _et al._ (2023) C. Dailey, N. Afshordi, and E. Schnetter, Class. Quant. Grav. 40, 195007 (2023), arXiv:2301.05778 [gr-qc] .
  * Kaup (1968) D. J. Kaup, Phys. Rev. 172, 1331 (1968).
  * Ruffini and Bonazzola (1969) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Liebling and Palenzuela (2023) S. L. Liebling and C. Palenzuela, Living Rev. Rel. 26, 1 (2023), arXiv:1202.5809 [gr-qc] .
  * Friedberg  _et al._ (1976) R. Friedberg, T. D. Lee, and A. Sirlin, Phys. Rev. D 13, 2739 (1976).
  * Palenzuela  _et al._ (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D 75, 064005 (2007), arXiv:gr-qc/0612067 .
  * Palenzuela  _et al._ (2008) C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. D 77, 044036 (2008), arXiv:0706.2435 [gr-qc] .
  * Palenzuela  _et al._ (2017) C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, Phys. Rev. D 96, 104058 (2017), arXiv:1710.09432 [gr-qc] .
  * Bezares  _et al._ (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [gr-qc] .
  * Helfer  _et al._ (2022) T. Helfer, U. Sperhake, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Class. Quant. Grav. 39, 074001 (2022), arXiv:2108.11995 [gr-qc] .
  * Bezares and Palenzuela (2018) M. Bezares and C. Palenzuela, Class. Quant. Grav. 35, 234002 (2018), arXiv:1808.10732 [gr-qc] .
  * Bezares  _et al._ (2022) M. Bezares, M. Bošković, S. Liebling, C. Palenzuela, P. Pani, and E. Barausse, Phys. Rev. D 105, 064067 (2022), arXiv:2201.06113 [gr-qc] .
  * Evstafyeva  _et al._ (2023) T. Evstafyeva, U. Sperhake, T. Helfer, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Class. Quant. Grav. 40, 085009 (2023), arXiv:2212.08023 [gr-qc] .
  * Siemonsen and East (2023a) N. Siemonsen and W. E. East, Phys. Rev. D 107, 124018 (2023a), arXiv:2302.06627 [gr-qc] .
  * Siemonsen and East (2023b) N. Siemonsen and W. E. East, Phys. Rev. D 108, 124015 (2023b), arXiv:2306.17265 [gr-qc] .
  * Schunck and Mielke (1996) F. E. Schunck and E. W. Mielke, “Rotating boson stars,” in  _Relativity and Scientific Computing: Computer Algebra, Numerics, Visualization_, edited by F. W. Hehl, R. A. Puntigam, and H. Ruder (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996) pp. 138–151.
  * Volkov and Wohnert (2002) M. S. Volkov and E. Wohnert, Phys. Rev. D 66, 085003 (2002), arXiv:hep-th/0205157 .
  * Kleihaus  _et al._ (2005) B. Kleihaus, J. Kunz, and M. List, Phys. Rev. D 72, 064002 (2005), arXiv:gr-qc/0505143 .
  * Siemonsen and East (2021) N. Siemonsen and W. E. East, Phys. Rev. D 103, 044022 (2021), arXiv:2011.08247 [gr-qc] .
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Berti  _et al._ (2006) E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Macedo  _et al._ (2013) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D 88, 064046 (2013), arXiv:1307.4812 [gr-qc] .
  * Heidmann  _et al._ (2023) P. Heidmann, N. Speeney, E. Berti, and I. Bah, Phys. Rev. D 108, 024021 (2023), arXiv:2305.14412 [gr-qc] .
  * Sperhake  _et al._ (2011) U. Sperhake, V. Cardoso, C. D. Ott, E. Schnetter, and H. Witek, Phys. Rev. D 84, 084038 (2011), arXiv:1105.5391 [gr-qc] .
  * Birkhoff and Langer (1923) G. D. Birkhoff and R. E. Langer,  _Relativity and modern physics_ (1923).
  * Jebsen (1921) J. T. Jebsen, Ark. Mat. Ast. Fys. 15 (1921).
  * Davis  _et al._ (1971) M. Davis, R. Ruffini, W. H. Press, and R. H. Price, Phys. Rev. Lett. 27, 1466 (1971).
  * East and Pretorius (2013) W. E. East and F. Pretorius, Phys. Rev. D 87, 101502 (2013), arXiv:1303.1540 [gr-qc] .
  * East (2014) W. E. East, Astrophys. J. 795, 135 (2014), arXiv:1408.1695 [gr-qc] .
  * Yoshida  _et al._ (1994) S. Yoshida, Y. Eriguchi, and T. Futamase, Phys. Rev. D 50, 6235 (1994).
  * Macedo  _et al._ (2016) C. F. B. Macedo, V. Cardoso, L. C. B. Crispino, and P. Pani, Phys. Rev. D 93, 064053 (2016), arXiv:1603.02095 [gr-qc] .
  * Vásquez Flores  _et al._ (2019) C. Vásquez Flores, A. Parisi, C.-S. Chen, and G. Lugones, JCAP 06, 051 (2019), arXiv:1901.07157 [hep-ph] .
  * Boerner  _et al._ (2023) T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth, and J. Towns, PEARC ’23: Practice and Experience in Advanced Research Computing (2023), 10.1145/3569951.3597559.
  * Lee (1987) T. D. Lee, Phys. Rev. D 35, 3637 (1987).
  * Friedberg  _et al._ (1987) R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D 35, 3658 (1987).
  * Lee and Pang (1992) T. D. Lee and Y. Pang, Phys. Rept. 221, 251 (1992).
  * Sanchis-Gual  _et al._ (2019) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, Phys. Rev. Lett. 123, 221101 (2019), arXiv:1907.12565 [gr-qc] .
  * Misner  _et al._ (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (1973).
  * Pitz and Schaffner-Bielich (2023) S. L. Pitz and J. Schaffner-Bielich, Phys. Rev. D 108, 103043 (2023), arXiv:2308.01254 [astro-ph.HE] .
  * Cunha  _et al._ (2016) P. V. P. Cunha, J. Grover, C. Herdeiro, E. Radu, H. Runarsson, and A. Wittig, Phys. Rev. D 94, 104023 (2016), arXiv:1609.01340 [gr-qc] .
  * Cunha  _et al._ (2017) P. V. P. Cunha, E. Berti, and C. A. R. Herdeiro, Phys. Rev. Lett. 119, 251102 (2017), arXiv:1708.04211 [gr-qc] .
  * Bošković and Barausse (2022) M. Bošković and E. Barausse, JCAP 02, 032 (2022), arXiv:2111.03870 [gr-qc] .
  * Collodel and Doneva (2022) L. G. Collodel and D. D. Doneva, Phys. Rev. D 106, 084057 (2022), arXiv:2203.08203 [gr-qc] .
  * Kobayashi  _et al._ (1994) Y. Kobayashi, M. Kasai, and T. Futamase, Phys. Rev. D 50, 7721 (1994).
  * Goebel (1972) C. J. Goebel, The Astrophysical Journal 172, L95 (1972).
  * Isaacson (1968) R. A. Isaacson, Phys. Rev. 166, 1263 (1968).
  * East  _et al._ (2012a) W. E. East, F. M. Ramazanoglu, and F. Pretorius, Phys. Rev. D 86, 104053 (2012a), arXiv:1208.3473 [gr-qc] .
  * Pretorius (2005) F. Pretorius, Class. Quant. Grav. 22, 425 (2005), arXiv:gr-qc/0407110 .
  * Alcubierre  _et al._ (2001) M. Alcubierre, S. Brandt, B. Bruegmann, D. Holz, E. Seidel, R. Takahashi, and J. Thornburg, Int. J. Mod. Phys. D 10, 273 (2001), arXiv:gr-qc/9908012 .
  * East  _et al._ (2012b) W. E. East, F. Pretorius, and B. C. Stephens, Phys. Rev. D 85, 124010 (2012b), arXiv:1112.3094 [gr-qc] .
  * Gundlach  _et al._ (2005) C. Gundlach, J. M. Martin-Garcia, G. Calabrese, and I. Hinder, Class. Quant. Grav. 22, 3767 (2005), arXiv:gr-qc/0504114 .
  * Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [gr-qc] .
  * Lindblom and Szilagyi (2009) L. Lindblom and B. Szilagyi, Phys. Rev. D 80, 084019 (2009), arXiv:0904.4873 [gr-qc] .