# Generic initial data for binary boson stars

Nils Siemonsen<EMAIL_ADDRESS>Perimeter Institute for
Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada  Arthur B. McDonald
Canadian Astroparticle Physics Research Institute, 64 Bader Lane, Queen’s
University, Kingston, ON K7L 3N6, Canada  Department of Physics & Astronomy,
University of Waterloo, Waterloo, ON N2L 3G1, Canada  William E. East
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

(August 28, 2024)

###### Abstract

Binary boson stars can be used to model the nonlinear dynamics and
gravitational wave signals of merging ultracompact, but horizonless, objects.
However, doing so requires initial data satisfying the Hamiltonian and
momentum constraints of the Einstein equations, something that has not yet
been addressed. In this work, we construct constraint-satisfying initial data
for a variety of binary boson star configurations. We do this using the
conformal thin-sandwich formulation of the constraint equations, together with
a specific choice for the matter terms appropriate for scalar fields. The free
data is chosen based upon a superposition of isolated boson star solutions,
but with several modifications designed to suppress the spurious oscillations
in the stars that such an approach can lead to. We show that the standard
approach to reducing orbital eccentricity can be applied to construct quasi-
circular binary boson star initial data, reducing the eccentricity of selected
binaries to the $\sim 10^{-3}$ level. Using these methods, we construct
initial data for quasi-circular binaries with different mass-ratios and spins,
including a configuration where the spin is misaligned with the orbital
angular momentum, and where the dimensionless spins of the boson stars exceeds
the Kerr bound. We evolve these to produce the first such inspiral-merger-
ringdown gravitational waveforms for constraint-satisfying binary boson stars.
Finally, we comment on how equilibrium equations for the scalar matter could
be used to improve the construction of binary initial data, analogous to the
approach used for quasi-equilibrium binary neutron stars.

##  I Introduction

In recent years, the detection of gravitational wave (GW) signals have
transformed the way we observe the universe and probe fundamental physics.
Predictions for the gravitational waveforms of inspiraling and merging binary
compact objects, such as black holes and neutron stars, play a crucial role in
the success of GW observation campaigns. While during the early inspiral one
can use perturbative methods, such an approach breaks down as nonlinear
effects take over at small separations of the binary. This necessitates the
use of numerical relativity to accurately predict the merger dynamics and
gravitational waveform through merger. Crucially, numerical evolutions of
binary compact objects rely on initial data that satisfies the constraints of
the Einstein equations as a starting point. For binary black holes, binary
neutron stars, and black hole-neutron star binaries, there has been extensive
work addressing the problem of constructing consistent initial data (see,
e.g., Refs. Cook (2000); Pfeiffer (2005); Gourgoulhon (2007); Baumgarte and
Shapiro (2010)), and a whole suite of formalisms, numerical methods, and
tricks have been developed (see Ref. Tichy (2017) for a recent review) to
obtain initial data suitable for binary evolutions that can directly be
compared to perturbative methods, both in the late inspiral, and the ringdown
phases of the coalescences (see, e.g., Refs. Pan  _et al._ (2014); Hannam  _et
al._ (2014); Bohé  _et al._ (2017)). Thus, binary initial data is a crucial
part of making the predictions that enable detecting GWs, characterizing GW
sources, and drawing conclusions for astrophysics and fundamental physics.

All GW observations of compact binary coalescence thus far are consistent with
arising from black hole and neutron stars. However, besides these two classes,
(ultra) compact and black hole mimicking objects motivated by extensions of
the Standard Model and models of quantum gravity have been proposed to solve
various problems in high energy and particle physics Cardoso and Pani (2019).
Boson stars (BSs) are particularly simple examples of such ultracompact
objects Kaup (1968); Ruffini and Bonazzola (1969); Yoshida and Eriguchi
(1997a); Seidel and Suen (1990, 1991); Brito  _et al._ (2016), and exhibit
many common features of black hole mimickers, such as ergoregions, stable
light rings, and the absence of horizons Friedberg  _et al._ (1987); Bošković
and Barausse (2022); Palenzuela  _et al._ (2017); Kleihaus  _et al._ (2005,
2008) (see Refs. Schunck and Mielke (2003); Visinelli (2021) for reviews).
Therefore, considerable effort has been put into the development of a program
to use binary BSs as a simple test bed to study the nonlinear and highly
dynamical regimes of the larger class of exotic compact objects (see Ref.
Liebling and Palenzuela (2017) for a review). Progress has been made in
understanding the early inspiral of binary BSs and the resulting GW emission
using perturbative methods Macedo  _et al._ (2013); Sennett  _et al._ (2017);
Krishnendu  _et al._ (2017); Johnson-Mcdaniel  _et al._ (2020); Vaglio  _et
al._ (2022); Adam  _et al._ (2022); Pacilio  _et al._ (2020); Vaglio  _et al._
(2023). There is also a significant body of work considering numerical
relativity evolutions of scalar and vector binary BSs, ranging from head-on
collisions of non-spinning and spinning binaries Lai (2004); Choptuik and
Pretorius (2010); Paredes and Michinel (2016); Bernal and Siddhartha Guzman
(2006); Schwabe  _et al._ (2016); Palenzuela  _et al._ (2007); Mundim (2010);
Bezares  _et al._ (2017); Sanchis-Gual  _et al._ (2019a); Helfer  _et al._
(2022); Sanchis-Gual  _et al._ (2022); Evstafyeva  _et al._ (2022), to orbital
inspirals Palenzuela  _et al._ (2008, 2017); Bezares  _et al._ (2022); Bezares
and Palenzuela (2018); Bezares  _et al._ (2017); Sanchis-Gual  _et al._
(2019a); Siemonsen and East (2023). However, despite the considerable amount
of work on the binary BS evolution problem, the problem of constructing
initial data for these compact binaries that is consistent with the Einstein
equations has yet to be addressed.

The common practice in the literature, when constructing binary BS initial
data, has been so far to simply superpose two boosted star solutions. More
recently, in Refs. Helfer  _et al._ (2022); Evstafyeva  _et al._ (2022), a
modified superposition trick was utilized to reduce (but not eliminate)
violations of the Hamiltonian and momentum constraints. While simple,
superposed initial data leads to large constraint violations even at moderate
separations, meaning that solutions to the evolution equations do not
accurately approximate solutions of the Einstein equations and effectively
excluding, for example, the quasi-circular binaries relevant to current GW
observations of inspiral-merger-ringdown111Note that this can not necessarily
be addressed by including constraint damping terms in the evolution equations.
See Appendix A. .

In this work, we develop and implement methods for constructing constraint-
satisfying binary scalar BS initial data for a wide variety of configurations.
(We also used these methods recently in Ref. Siemonsen and East (2023).) We
solve the constraint equations in the conformal thin-sandwich (CTS) formalism
York (1999), using free data based on superposing stationary BS solutions,
similar to what was done in Ref. East  _et al._ (2012a) for black hole and
fluid stars. However, considering scalar matter introduces several new
complications which we address here, including the choice of which matter
degrees of freedom to fix, as well as how to minimize spurious oscillations
which may be induced in the BSs. Our approach is very flexible, and we use it
to construct binary BSs with different mass-ratios, spin magnitudes, and spin
orientations. We evolve several such binaries, including several cases where
the BSs are super-spinning (i.e., have dimensionless spins exceed the Kerr
bound of unity), through inspiral and merger. We do this both for physical
interest, as well as to demonstrate that we can construct quasi-circular
binaries by adapting eccentricity reduction techniques.

The remainder of this paper is organized as follows. We briefly review the
relevant physics of isolated stationary BSs in Sec. II.1, and proceed in Sec.
II.2 to introduce our procedure to self-consistently solve the elliptic
constraint equations in the CTS formalism York (1999) numerically, given an
initial guess for the binary, utilizing the elliptic solver introduced in Ref.
East  _et al._ (2012a). To that end, we identify the most suitable
parameterization of the matter content of BSs in Sec. II.3. We comment on
possible equilibrium conditions for the scalar matter of these stars in Sec.
II.4, though we do not implement such an approach in this study. In Sec.
III.2, we analyze the quality of the constructed initial data and devise
methods to reduce spurious oscillations in each star, as well as in the
resulting gravitational radiation, and lastly, in Sec. III.3, we test
eccentricity reduction schemes in the context of binary BS inspirals and
comment on their possible limitations. We consider binary configurations in
two different scalar potentials, with equal and unequal masses, as well as
non-spinning and spinning constituent stars with aligned and misaligned spins.
Finally, in Sec. IV, we analyze the dynamics of selected eccentricity-reduced
binary configurations, and present inspiral-merger-ringdown gravitational
waveforms. We give details on the numerical evolution scheme, present
convergence results, and briefly compare the constraint satisfying initial
data we construct in this study to superposed initial data in Appendix A. We
provide some further details on spurious high frequency components to the GWs
and correcting for center-of-mass motion of the BS binaries in Appendices B
and C, respectively. We use $G=c=1$ units throughout.

##  II Methodology

###  II.1 Isolated boson star

Before turning to the construction of binary BS data, we briefly review
isolated BSs in their rest frames. Here and throughout, we focus entirely on
scalar BSs, and leave a full consideration of vector BSs to future work.
Scalar BSs are regular, asymptotically flat, stationary, and axisymmetric
solutions to the globally U(1)-invariant Einstein-Klein-Gordon theory

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{R}{16\pi}-g^{\alpha\beta}\nabla_{(\alpha}\bar{\Phi}\nabla_{\beta)}\Phi-V(|\Phi|)\right],$
(1)

where $R$ is the Ricci scalar of the spacetime $g_{\mu\nu}$, $\Phi$ is the
complex scalar field (with an overbar denoting complex conjugation), and
$V(|\Phi|)$ is the global U(1)-preserving scalar potential. We consider cases
where the lowest order term in $|\Phi|$ of $V$ is a mass term
$\mu^{2}|\Phi|^{2}$. The conserved Noether current associated with the $U(1)$
symmetry of the action (1) is

$\displaystyle\nabla_{\mu}j^{\mu}=0,$ $\displaystyle
j^{\mu}=-i(\bar{\Phi}\nabla^{\mu}\Phi-\Phi\nabla^{\mu}\bar{\Phi}).$ (2)

Intuitively, the charge $Q$ of a BS counts the number bosons in the solution.
This is in direct analogy to the conservation of baryon number in the case of
fluid stars.

The metric ansatz for isolated star solutions in Lewis-Papapetrou coordinates
takes the form

$\displaystyle\begin{aligned}
ds^{2}=-fdt^{2}+lf^{-1}\big{\\{}g(dr^{2}+r^{2}d\theta^{2})\\\
+r^{2}\sin^{2}\theta(d\varphi-\Omega r^{-1}dt)^{2}\big{\\}},\end{aligned}$ (3)

where the metric functions $f$, $l$, $g$, and $\Omega$ depend, in general, on
both $r$ and $\theta$. The boundary conditions on these functions for BS
solutions can be found in Ref. Kleihaus  _et al._ (2005); generally, the
solution is regular at the star’s center and asymptotically flat. The scalar
field ansatz for BSs is

$\displaystyle\Phi=\phi(r,\theta)e^{i(\omega t-m\varphi)},$ (4)

where $\phi$ is the magnitude, while $\omega$ is the star’s internal
frequency, and $m$ is the integer azimuthal index. At large distances, due to
the non-zero scalar mass, the field falls off exponentially as

$\displaystyle\lim_{r\rightarrow\infty}\phi\sim
e^{-\sqrt{\mu^{2}-\omega^{2}}r},$ (5)

ensuring that the solution is asymptotically flat. Spherically symmetric star
solutions attain a non-zero scalar field magnitude at the center,
$\phi|_{r=0}\neq 0$, while rotating solutions (i.e., those with $|m|\geq 1$)
exhibit a vortex line through their centers, and hence, are toroidal in shape
(see Ref. Siemonsen and East (2023) for a discussion of the vortex structure).
Details on the numerical construction of these solutions used in this work are
discussed in the appendices of Ref. Siemonsen and East (2021).

We are interested in testing our methods to construct binary BS initial data
in various physically interesting regimes. Therefore, we focus entirely on
scalar models with non-vanishing self-interactions, as these have shown to
lead to highly-compact Palenzuela  _et al._ (2017); Bošković and Barausse
(2022), as well as stably rotating BS solutions Siemonsen and East (2021). To
that end, we consider the solitonic scalar potential Friedberg  _et al._
(1987)

$\displaystyle
V(|\Phi|)=\mu^{2}|\Phi|^{2}\left(1-\frac{2|\Phi|^{2}}{\sigma^{2}}\right)^{2},$
(6)

controlled by the coupling $\sigma$, as well as the leading repulsive scalar
self-interaction

$\displaystyle V(|\Phi|)=\mu^{2}|\Phi|^{2}+\lambda|\Phi|^{4},$ (7)

in the following simply referred to as the repulsive potential (restricted to
$\lambda>0$). Imposing the appropriate asymptotically flat condition at
spatial infinity, and regularity conditions at the origin, in conjunction with
the ansätze (3) and (4), and the field equations resulting from (1), yields a
one-parameter family of solutions for each azimuthal index $m$ parameterized
by the internal frequency $\omega$222We focus solely on stars in their
respective radial ground states, i.e., the scalar field profiles have no
additional radial nodes.. For each member of the family of spacetimes, the
total mass $M$ is given by the Komar expression, coinciding with the ADM mass,
and the radius $R$ is defined as the areal (circular) radius $r_{*}$ at which
$99\%$ of the mass of the spherical (rotating) star solution lies within
$r_{*}$. The angular momentum $J$ is well-defined by means of the Komar
expression in this stationary and axisymmetric class of spacetimes. Finally,
the charge $Q$ of each star follows from the U(1)-Noether charge of the action
(1) (details can be found in Ref. Siemonsen and East (2021)). For these
isolated BS solutions, this Noether charge satisfies the “quantization”
relation $J=mQ$.

Figure 1: We plot the relationship between the mass $M$ and the radius $R$, in
units of the boson mass parameter $\mu$, for four families of BS solutions. We
show the solitonic scalar model (6) with coupling $\sigma=0.05$ (blue curves),
and the repulsive potential (7) with coupling $\lambda/\mu^{2}=10^{3}$ (orange
curves), including both spherically symmetric (dotted) and $m=1$ rotating
families of stars (solid). The regime with $M/R>1/2$ is labelled as “BH.”

In Fig. 1, we show the mass-radius relationships of the isolated BS solutions
considered in this work. As can be seen there, there is a regime where the
solitonic families have highly compact solutions; as a result, some of these
solutions exhibit regions of stable trapping of null geodesics, as well as
large ergoregions if $m>0$. Hence, this class of spacetimes exhibit highly
relativistic features potentially leading to interesting phenomenology in the
context of binary systems. In contrast to this, solutions of families in the
scalar model with repulsive self-interactions are mostly less compact, $C_{\rm
max}\approx 0.14$, enabling the study of BSs in the Newtonian regime.
Importantly, both scalar models contain spinning solutions with $m=1$ that
have been shown to be stable Siemonsen and East (2021) against the non-
axisymmetric instability discovered in Ref. Sanchis-Gual  _et al._ (2019b) on
long timescales.

###  II.2 Conformal thin-sandwich formulation

Moving now to the formalism used to construct constraint satisfying binary BS
initial data, we introduce the CTS formulation of the Hamiltonian and momentum
constraints of the Einstein equations. To that end, the spacetime is foliated
into a series of spacelike hypersurfaces $\Sigma_{t}$, parameterized by the
coordinate time $t$, with future-pointing unit-normal to the hypersurface
$n^{\mu}$. The tangent of lines of constant spatial coordinates is then

$\displaystyle t^{\mu}=\alpha n^{\mu}+\beta^{\mu},$ (8)

with lapse function $\alpha$ and shift vector $\beta^{\mu}$, with
$n_{\mu}\beta^{\mu}=0$. Furthermore, let
$\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$ be the projector onto the
hypersurface $\Sigma_{t}$, such that $\gamma_{ij}$ is the spatial metric
induced on $\Sigma_{t}$. Lastly, the extrinsic curvature of $\Sigma_{t}$,

$\displaystyle K_{ij}=-\frac{1}{2}\mathcal{L}_{n}\gamma_{ij},$ (9)

is defined by means of the Lie derivative $\mathcal{L}_{n}$ along the
hypersurface normal. In this 3+1 language, the Hamiltonian and momentum
constraints, i.e., the projections of the Einstein equations along the
hypersurface normal, are with trace $K=\gamma^{ij}K_{ij}$, Ricci scalar
${}^{(3)}R$ and derivative $D_{i}$ defined with respect to the induced metric
$\gamma_{ij}$, and finally, the energy density $E=n_{\mu}n_{\nu}T^{\mu\nu}$
and momentum density $p^{i}=-\gamma^{i}_{\mu}n_{\nu}T^{\mu\nu}$ of the matter
content of the space in the Eulerian frame.

The CTS formulation York (1999) of the Hamiltonian and momentum constraints
(LABEL:eq:constraints) relies on relating the constraint satisfying metric
components $\gamma_{ij}$ to a freely specifiable conformal metric
$\tilde{\gamma}_{ij}$ as

$\displaystyle\gamma_{ij}=\Psi^{4}\tilde{\gamma}_{ij},$ (10)

with conformal factor $\Psi$. Furthermore, the traceless part of the extrinsic
curvature $A^{ij}$ is conformally decomposed as
$A^{ij}=\Psi^{-10}\hat{A}^{ij}$ with

$\displaystyle\hat{A}^{ij}=\frac{1}{2\tilde{\alpha}}\Big{[}(\tilde{\mathbb{L}}\beta)^{ij}+\partial_{t}\tilde{\gamma}^{ij}\Big{]},$
(11)

in terms of the conformal Killing form
$(\tilde{\mathbb{L}}\beta)^{ij}=\tilde{D}^{i}\beta^{j}+\tilde{D}^{j}\beta^{i}-2\tilde{\gamma}^{ij}\tilde{D}_{k}\beta^{k}/3$.
The conformal lapse and the time-derivative of the conformal metric are
$\tilde{\alpha}=\Psi^{-6}\alpha$ and
$\partial_{t}\tilde{\gamma}^{ij}=\Psi^{4}(\partial_{t}\gamma^{ij}-\gamma^{ij}\gamma_{kl}\partial_{t}\gamma^{kl}/3)$,
respectively. Utilizing this decomposition, the metric constraints
(LABEL:eq:constraints) are cast into the CTS equations

$\displaystyle\begin{aligned}
\tilde{D}_{i}\tilde{D}^{i}\Psi-\frac{\tilde{R}}{8}\Psi+\frac{\hat{A}_{ij}\hat{A}^{ij}}{8}\Psi^{-7}-\frac{K^{2}}{12}\Psi^{5}=&-2\pi\Psi^{5}E,\\\
\tilde{D}_{j}\hat{A}^{ij}-\frac{2}{3}\Psi^{6}\tilde{D}^{i}K=&\
8\pi\Psi^{10}p^{i}.\end{aligned}$ (12)

The geometric free data (i.e., those metric variables that need to be
specified) are comprised of the metric $\tilde{\gamma}_{ij}$ and its
coordinate time derivative $\partial_{t}\tilde{\gamma}^{ij}$, as well as the
trace $K$ of the extrinsic curvature and the lapse function $\tilde{\alpha}$.
This formulation is supplemented by a choice of energy and momentum densities,
$E$ and $p^{i}$, of the complex scalar matter, as well as the corresponding
scalar free data. This will be discussed in detail in the next section.

To specify the metric free data, we proceed as follows. First, we solve for
stationary isolated BSs in Lewis-Papapetrou coordinates, as outlined in Sec.
II.1. These solutions are subsequently transformed to Cartesian
coordinates333We transform from Lewis-Papapetrou to Cartesian coordinates by
applying the usual flat relations between spherical and Cartesian coordinates.
and boosted using initial coordinate velocities $v_{(A)}^{i}$, where
$A\in\\{1,2\\}$ labels each star in the binary, and placed at coordinate
positions $z_{(A)}^{i}$. Therefore, for each star we obtain the set of
variables $\gamma_{ij}^{(A)}$, $\partial_{t}\gamma_{ij}^{(A)}$,
$\alpha^{(A)}$, and $\beta^{i}_{(A)}$ (for both stars in Cartesian-type
coordinates). For all binary configurations presented in this work,
$z_{(A)}^{i}$ and $v_{(A)}^{i}$ are chosen such that the initial center-of-
mass location coincides with the origin of the numerical grid, and the initial
linear momentum of the center of mass vanishes (at the Newtonian level);
limitations of this approach are discussed in Appendix C. These two solutions
are then superposed as

$\displaystyle\begin{aligned} \gamma_{ij}^{\rm sup}=&\
\eta_{ij}+f_{(2)}\big{[}\gamma^{(1)}_{ij}-\eta_{ij}\big{]}+f_{(1)}\big{[}\gamma^{(2)}_{ij}-\eta_{ij}\big{]},\\\
\partial_{t}\gamma^{\rm sup}_{ij}=&\
f_{(2)}\partial_{t}\gamma_{ij}^{(1)}+f_{(1)}\partial_{t}\gamma_{ij}^{(2)},\\\
\alpha_{\rm sup}=&\
1+f_{(2)}\big{[}\alpha^{(1)}-1\big{]}+f_{(1)}\big{[}\alpha^{(2)}-1\big{]},\\\
\beta^{i}_{\rm sup}=&\
f_{(2)}\beta^{i}_{(1)}+f_{(1)}\beta^{i}_{(2)},\end{aligned}$ (13)

where $\eta^{i}{}_{j}=\delta^{i}_{j}$ is the flat 3-metric, and $f_{(A)}$ is
an attenuation function, which we introduce here for convenience and discuss
in detail in Sec. III.2. For now, we simply point out that the choice
$f_{(A)}\equiv 1$ corresponds to a plain superposition of the isolated stars.
As discussed in Ref. East  _et al._ (2012a), the metric free data is then
obtained from

$\displaystyle\begin{aligned} \tilde{\gamma}_{ij}=&\ \gamma_{ij}^{\rm sup},\\\
\partial_{t}\tilde{\gamma}^{ij}=&-\tilde{\gamma}^{ik}\tilde{\gamma}^{jl}\left(\partial_{t}\gamma_{kl}^{\rm
sup}-\frac{1}{3}\tilde{\gamma}_{kl}\tilde{\gamma}^{mn}\partial_{t}\gamma_{mn}^{\rm
sup}\right),\\\ \tilde{\alpha}=&\ \alpha_{\rm sup},\\\ K=&\
\frac{1}{2\tilde{\alpha}}\big{[}2\partial_{i}\beta^{i}_{\rm
sup}+\tilde{\gamma}_{ij}\partial_{t}\tilde{\gamma}^{ij}+\tilde{\gamma}^{ij}\beta^{k}_{\rm
sup}\partial_{k}\tilde{\gamma}_{ij}\big{]}.\end{aligned}$ (14)

The CTS equations admit solutions provided appropriate boundary conditions are
specified. In the context of binary BSs, i.e., asymptotically flat spacetimes,
we require that

$\displaystyle\lim_{|\textbf{x}|\rightarrow\infty}\Psi=1,$
$\displaystyle\lim_{|\textbf{x}|\rightarrow\infty}\beta^{i}=\beta^{i}_{\rm
sup}|_{\infty},$ (15)

where $\beta^{i}_{\rm sup}|_{\infty}$ is the shift of the free data at large
distances. With these boundary conditions, we solve the CTS equations
numerically using a multigrid scheme with fixed mesh refinement (further
details can be found in Ref. East  _et al._ (2012a)). In the context of
axisymmetry, we employ a generalized Cartoon method that provides derivatives
about the axis of symmetry by means of a the axisymmetric Killing field,
allowing also for harmonic azimuthal dependencies in the scalar sector.

###  II.3 Binary boson star sources

So far, we have left the precise parameterization of the scalar matter
sourcing the spacetime, $E$ and $p^{i}$, unspecified. In principle, various
choices of energy and momentum densities measured by an Eulerian observer, are
possible for time-dependent complex scalar field matter. However, we find the
precise choice to be crucial to achieve convergence of our numerical
implementation. Therefore, in the following we outline possible matter source
parameterizations, and especially, highlight the method we found to robustly
yield consistent binary BS initial data in any considered context.

We begin by introducing the necessary projections of the nonlinear complex
scalar energy-momentum tensor with respect to the foliation introduced in the
previous section. The latter is readily obtained from (1) in covariant form:

$\displaystyle T_{\mu\nu}=2\partial_{(\mu}\bar{\Phi}\partial_{\nu)}\Phi-
g_{\mu\nu}\Big{[}g^{\alpha\beta}\partial_{(\alpha}\bar{\Phi}\partial_{\beta)}\Phi+V(|\Phi|)\Big{]}.$
(16)

For convenience, we define the conjugate momentum of the complex field with
respect to a spatial slice as

$\displaystyle\eta:=\mathcal{L}_{n}\Phi=n^{\alpha}\partial_{\alpha}\Phi=\frac{1}{\alpha}(\partial_{t}\Phi-\beta^{i}D_{i}\Phi).$
(17)

With this, the scalar energy-momentum tensor can be written in the Eulerian
frame as

$\displaystyle\begin{aligned} E=&\
n_{\alpha}n_{\beta}T^{\alpha\beta}=\eta\bar{\eta}+D^{i}\Phi
D_{i}\bar{\Phi}+V(|\Phi|)\\\
p^{i}=&-\gamma^{i}{}_{\alpha}n_{\beta}T^{\alpha\beta}=-\eta
D^{i}\bar{\Phi}-\bar{\eta}D^{i}\Phi\\\ S^{ij}=&\
\gamma^{i}{}_{\alpha}\gamma^{j}{}_{\beta}T^{\alpha\beta}=2D^{(i}\Phi
D^{j)}\bar{\Phi}\\\ &\ +\gamma^{ij}\Big{[}\eta\bar{\eta}-D_{k}\Phi
D^{k}\bar{\Phi}-V(|\Phi|)\Big{]}.\end{aligned}$ (18)

Starting from these expressions, we discuss possible approaches to
parameterize the scalar matter, as well as the associated choices of scalar
free data accompanying the metric free data (14). To that end, and before
presenting the source parameterization that we found to work for any type of
binary BS configuration, it is instructive to consider two other approaches
that, while natural, exhibit fundamental issues.

####  II.3.1 Fixed energy and momentum densities

In the context of binary neutron star initial data, it is natural to choose
the conformal Eulerian energy $\tilde{E}$ and momentum densities
$\tilde{p}^{i}$ as free data for the CTS system of equations (12). The
corresponding physical momentum density $p^{i}$ is typically chosen to be
$p^{i}=\Psi^{-10}\tilde{p}^{i}$. One choice for the conformal scaling of the
energy is $E=\Psi^{-8}\tilde{E}$, which is motivated by uniqueness arguments
and the preservation of the dominant energy condition (i.e., if the free data
satisfies $\tilde{E}\geq\sqrt{\tilde{\gamma}^{ij}\tilde{p}_{i}\tilde{p}_{j}}$
so does the constraint satisfying initial data). In the case of fluid stars,
the initial physical pressure $P$ and density $\rho$ are recovered by means of
an algebraic relation between $E$, $p^{i}$, $P$, and $\rho$ derived from the
expression for the fluid energy-momentum tensor combined with the fluid
equation of state. This provides a means to reconstruct the constraint
satisfying fluid variable initial conditions directly from the free data and
constrained data.

A complex scalar field, on the other hand, has kinetic and gradient energy, in
addition to potential energy. The energy therefore depends on spatial
gradients $D^{i}\Phi$ and time-derivatives $\partial_{t}\Phi$ of the scalar
field, as can be seen in (18). Therefore, unlike in the binary neutron star
scenario, the relation between the physical energy and momentum densities, and
the matter field $\Phi$, is not purely algebraic, but rather of differential
form. This renders the reconstruction of the scalar field initial data $\Phi$
and $\partial_{t}\Phi$ (or equivalently $\Phi$ and $\eta$) from the constraint
satisfying energy and momentum densities non-trivial.

Irrespective of these shortcomings, we test this choice of matter source
variables with a set of single isolated non-spinning and spinning BSs. The
metric free data is constructed following the discussion in Sec. II.2 (setting
$\gamma_{ij}^{(2)}=0$ etc.), while the scalar sources $\tilde{E}$ and
$\tilde{p}^{i}$ are determined from (18). The CTS equations are then
numerically solved iteratively as outlined above. We succeeded in recovering
isolated, boosted, non-spinning and spinning BSs utilizing this approach,
i.e., the elliptic CTS solver removed truncation error of the isolated
solution to the precision allowed by the resolution of the discretization of
the CTS equations. This shows that, within our numerical setup, solutions to
isolated stars are in fact local attractors in the space of solutions using
this scalar matter parameterization.

####  II.3.2 Fixed scalar initial data

In order to circumvent the issue discussed in the previous section, i.e.,
instead of fixing the energy and momentum densities directly, one could
provide the scalar field initial data itself—$\Phi$ and $\partial_{t}\Phi$—as
free data to the CTS equations. This makes the reconstruction of the scalar
field trivial, ensuring that the metric and scalar initial data consistently
solve the Hamiltonian and momentum constraints. To that end, we rewrite (18)
in terms of the scalar fields $\Phi$ and $\partial_{t}\Phi$, leading to

$\displaystyle\begin{aligned} E=&\
\frac{\Psi^{-12}}{\tilde{\alpha}^{2}}|\partial_{t}\Phi-\beta^{i}\partial_{i}\Phi|^{2}+\Psi^{-4}\tilde{D}^{i}\bar{\Phi}\tilde{D}_{i}\Phi+V(|\Phi|),\\\
p^{i}=&-\frac{\Psi^{-10}}{\tilde{\alpha}^{2}}\Big{(}(\partial_{t}\Phi-\beta^{i}\partial_{i}\Phi)\tilde{D}^{i}\bar{\Phi}+(\partial_{t}\bar{\Phi}-\beta^{i}\partial_{i}\bar{\Phi})\tilde{D}^{i}\Phi\Big{)}\end{aligned}$
(19)

in terms of the conformal variables. The different scaling of the kinetic,
gradient, and potential energies with conformal factor $\Psi$, as well as the
dependencies on the shift vector $\beta^{i}$ indicate that this approaches
differs from providing $E$ and $p^{i}$ as free data by more than a simple
$\Psi$-rescaling.

We tested the above choice of sources within our numerical setup, similarly to
our tests of the formulation presented in Sec. II.3.1. We found robust
convergence of the numerical schemes in the case of isolated boosted non-
spinning BSs. However, we were unable to recover an isolated stationary
rotating BS solution using (19) in the CTS equations. Despite the CTS equation
residuals converging to zero at the expected order before the first iteration,
the solution of the elliptic solver moves away from the true solution
exponentially quickly with each iteration. This indicates that rotating BSs
are not attractors in the space of solutions when using (19) and our numerical
framework, or suggests a break-down in the uniqueness of this solution for the
given free and boundary data. Note, tests of uniqueness based on the maximum
principle (see, e.g., Ref. Gourgoulhon (2007)) are not applicable in this
case, since the momentum constraint is not trivially satisfied by stationary
rotating BS solutions444Even if the maximum principle could be applied in this
case, perturbations $\epsilon$ away from a solution $\Psi_{0}$ to the
Hamiltonian constraint follow to linear order the equation
$\tilde{D}_{i}\tilde{D}^{i}\epsilon=\Gamma\epsilon$, with
$\Gamma=\tilde{R}/2-14\pi\Psi_{0}^{-8}\tilde{\alpha}^{-2}|\partial_{t}\Phi-\beta^{i}\partial_{i}\Phi|^{2}+\dots$,
where we ignored all positive-definite terms. Hence, non-vanishing kinetic
energy and momentum densities $p^{i}$, present even in isolated rotating BSs,
may result in violations of the maximum principle, i.e., in $\Gamma<0$..

####  II.3.3 Fixed scalar kinetic energy

We turn now to the choice of scalar matter free data that we found to robustly
lead to constraint satisfying binary BS metric and scalar field initial data.
Similar choices of free data were considered recently in cosmological contexts
in Refs. Aurrekoetxea  _et al._ (2022); Corman and East (2022). Here, instead
of setting the scalar initial data $\\{\Phi,\partial_{t}\Phi\\}$ as free data
for the CTS equations, we replace $\partial_{t}\Phi$ by the scalar field’s
conformal conjugate momentum $\tilde{\eta}=\Psi^{6}\eta$ as free data. With
this, and in terms of the above conformal decomposition, the energy and
momentum densities turn into

$\displaystyle\begin{aligned} E=&\
\Psi^{-12}\tilde{\eta}\bar{\tilde{\eta}}+\Psi^{-4}\tilde{D}^{i}\Phi\tilde{D}_{i}\bar{\Phi}+V(|\Phi|),\\\
p^{i}=&\
-\Psi^{-10}(\tilde{\eta}\tilde{D}^{i}\bar{\Phi}+\bar{\tilde{\eta}}\tilde{D}^{i}\Phi).\end{aligned}$
(20)

The scalar data satisfying the constraint equations can then be recovered via
the algebraic relation

$\displaystyle\partial_{t}\Phi=\tilde{\alpha}\tilde{\eta}+\beta^{i}\partial_{i}\Phi,$
(21)

where $\beta$ is the solution to the vector CTS equation provided the free
data $\\{\Phi,\tilde{\eta}\\}$; notice
$\alpha\eta=\tilde{\alpha}\tilde{\eta}$. For completeness, we include here
also the expressions (20) in terms of $\Phi_{R}=(\Phi+\bar{\Phi})/2$, and
$\Phi_{I}=(\Phi-\bar{\Phi})/(2i)$, as these are the variables used in our
numerical implementation of the elliptic CTS solver, as well as the hyperbolic
evolution scheme:

$\displaystyle\begin{aligned} E=&\
\big{[}\Psi^{-12}\tilde{\eta}_{R}^{2}+\Psi^{-4}\tilde{D}_{i}\Phi_{R}\tilde{D}^{i}\Phi_{R}+(R\leftrightarrow
I)\big{]}\\\ &\qquad\qquad+V(\Phi_{R}^{2}+\Phi_{I}^{2}),\\\
p^{i}=&-2\Psi^{-10}\big{[}\tilde{\eta}_{R}\tilde{D}^{i}\Phi_{R}+(R\leftrightarrow
I)\big{]}.\end{aligned}$ (22)

We found this matter source parameterization to robustly recover any kind of
single BS solution in those tests outlined in Sec. II.3.1. Given this
parameterization passes these tests, we are now able to move to binary BS. To
that end, we construct the scalar free data in a similar fashion to the
superposed metric free data presented in (13). First, the scalar fields
$\Phi^{(A)}$ of each star are boosted with the same boost as the metric, then
the conjugate momenta $\tilde{\eta}^{(A)}$ are determined555Recall, the scalar
field has a non-trivial time-dependence [see also (4)], so that we first boost
the vector $\partial_{\mu}\Phi$, and then determine its conjugate momentum in
the boosted frame., and finally, the variables are superposed to obtain the
scalar free data as follows:

$\displaystyle\begin{aligned} \Phi^{\rm sup}=&\
\hat{f}_{(2)}\Phi^{(1)}+\hat{f}_{(1)}\Phi^{(2)},\\\ \tilde{\eta}^{\rm sup}=&\
\hat{f}_{(2)}\tilde{\eta}^{(1)}+\hat{f}_{(1)}\tilde{\eta}^{(2)}.\end{aligned}$
(23)

Here $\hat{f}_{(A)}$ are attenuation functions (directly analogous to
$f_{(A)}$, defined in Sec. II.2), which we discuss in detail below and simply
note here that $\hat{f}_{(A)}=1$ corresponds to a simple superposition of the
two star’s scalar field variables. Notice, while the constraint equations are
invariant under a global phase shift $\Phi\rightarrow\Phi e^{i\alpha}$, the
source functions (20) are not invariant under the phase of a single
constituent of a binary, e.g., $\Phi^{(1)}\rightarrow\Phi^{(1)}e^{i\alpha}$.

###  II.4 Scalar matter equilibrium

We focused so far on finding a scalar matter parameterization that robustly
yields constraint satisfying binary BS initial data. Since no assumptions on
the stars’ trajectories, spin orientations, or mass-ratio were built into the
formalism, this approach is very flexible. However, as we show in Sec. III,
the simple construction of the free data described above results in stars with
large internal oscillations and ejected scalar matter. In the case of binary
neutron stars, various methods have been introduced to alleviate these issues
by explicitly equilibriating the fluid and metric degrees of freedom Bonazzola
_et al._ (1997); Shibata (1998); Teukolsky (1998). These approaches are based
on assuming the existence of a helical Killing field $\ell^{\mu}$, which
provides a notion of equilibrium not just for the metric, i.e.,
$\mathcal{L}_{\ell}g_{\mu\nu}=0$, but also for the matter variables. In the
case of binary neutron stars, combining the conservation of rest mass density,
the conservation of the fluid’s energy-momentum, and matter equilibrium with
respect to $\ell^{\mu}$, results in an elliptic equation for the equilibriated
initial velocity of the fluid. In the following, we apply these arguments
qualitatively to the case of scalar field matter and outline approaches to
equilibriating the scalar matter. However, as we are not testing these
formalisms here explicitly, this is to be understood as a first step guiding
more thorough future analyses.

To that end, we assume the existence of a helical Killing field $\ell^{\mu}$,
such that $\ell^{\mu}=\alpha n^{\mu}+\tilde{V}^{\mu}$. In the asymptotically
inertial frame, the spatial velocity takes the form
$\tilde{V}^{\mu}=\beta^{\mu}+\tilde{\Omega}m^{\mu}$, where $m^{\mu}$ is
spacelike generating the azimuthal direction and $\tilde{\Omega}$ is the
orbital period of the helical field $\ell^{\mu}$. On the other hand, in the
co-rotating frame, the Killing field $\ell^{\mu}\rightarrow t^{\mu}$, such
that $\tilde{V}^{\mu}\rightarrow\beta^{\mu}$ (for a discussion on the
subtleties associated with this choice, see, e.g., Ref. Gourgoulhon  _et al._
(2002)). From the perspective of the observer associated with $n^{\mu}$, the
Noether-current, defined in (2), decomposes as

$\displaystyle j^{\mu}=$ $\displaystyle\ \rho n^{\mu}+J^{\mu},$ (24)
$\displaystyle\rho=$ $\displaystyle-
n_{\mu}j^{\mu}=i(\bar{\Phi}\eta-\Phi\bar{\eta}),$ (25)

with local boson number density $\rho$, and spatial current

$\displaystyle J^{\mu}=-i(\bar{\Phi}D^{\mu}\Phi-\Phi D^{\mu}\bar{\Phi}).$ (26)

In direct analogy to the rest mass conservation equation for fluids, the
global U(1) symmetry of the scalar theory implies the boson number
conservation [see (2)]. With this above decomposition of the current, the
conservation law (2) reduces to

$\displaystyle\mathcal{L}_{n}\rho=\rho K-\frac{1}{\alpha}D_{i}(\alpha J^{i}).$
(27)

Correspondingly, the evolution equation for the scalar field—the Klein-Gordon
equation—is readily obtained from (1):

$\displaystyle\left[\nabla_{\mu}\nabla^{\mu}-\partial_{|\Phi|^{2}}V(|\Phi|)\right]\Phi=0.$
(28)

Using the foliation defined by $n^{\mu}$, the Klein-Gordon equation takes the
form

$\displaystyle\mathcal{L}_{n}\eta=\frac{1}{\alpha}D_{i}(\alpha
D^{i}\Phi)+K\eta-\Phi\partial_{|\Phi|^{2}}V(|\Phi|),$ (29)

and similarly for the conjugate equation. To proceed, a series of equilibrium
conditions, utilizing the helical Killing vector, must be imposed on the
matter variables.

Contrary to the matter variables relevant for fluid stars, the scalar matter
making up BSs is not time-independent. To understand possible equilibrium
conditions, recall that the scalar field ansatz (4), with $\eta\sim
i(\omega-m\Omega/r)\Phi$, contains a harmonic time-dependence of the scalar
variables due to the linear time-dependence of the scalar phase. Therefore,
here we explore how the ansatz underlying the isolated BS solution, i.e.,

$\displaystyle\mathcal{L}_{\ell}\Phi=i\tilde{\omega}\Phi,$ (30)

may generalize to a binary system. To that end, we specialize to the co-
rotating frame for the remainder of this section. This implies that
$\mathcal{L}_{\ell}\alpha=0$ and $\mathcal{L}_{\ell}\beta^{\mu}=0$, which
together with (30) imply that $\mathcal{L}_{\ell}\eta=i\tilde{\omega}\eta$.
With these, the Klein-Gordon equation (29) reduces to the complex equation

$\displaystyle-\beta^{i}D_{i}\eta+i\tilde{\omega}\eta-\alpha
K\eta=D_{i}(\alpha D^{i}\Phi)-\alpha\Phi\partial_{|\Phi|^{2}}V(|\Phi|),$ (31)

where $\eta=(i\tilde{\omega}\Phi-\beta^{i}D_{i}\Phi)/\alpha$ in the co-
rotation frame. This second-order elliptic equation for $\Phi$ is directly
analogous to the scalar equation obtained from plugging the metric and scalar
field ansatz of an isolated stationary BS [given in (3) and (4), respectively]
into the Klein-Gordon equation (28). Binary BS boundary conditions are
$\lim_{|\textbf{x}|\rightarrow\infty}\Phi=0$.

In this approach, the frequency $\tilde{\omega}$ must be chosen. Even in the
equal-frequency regime, we expect the binary’s frequency $\tilde{\omega}$ to
differ from the isolated star’s frequency
$\omega_{1}=\omega_{2}\neq\tilde{\omega}$ due to the increase in binding
energy; hence, the naive expectation is that $\tilde{\omega}<\omega_{1}$. The
frequency $\tilde{\omega}$ could be iteratively adjusted to achieve a desired
property for the binary initial data (e.g. a value for the mass or scalar
charge), or a new binary BS solution may be found at fixed $\tilde{\omega}$
(analogous to the case of an isolated BS). For unequal-mass binaries, the
binary’s frequency is spatially dependent since $\omega_{1}\neq\omega_{2}$. A
simple choice is to construct a differentiable $\tilde{\omega}$ transitioning
from a fraction of $\omega_{1}$ around the first star to the same (or
different) fraction of $\omega_{2}$ around the second star666Instead of making
an ad-hoc choice for $\tilde{\omega}$, the total charge of the binary may be
kept fixed following the approach introduced in Ref. Kleihaus  _et al._ (2005)
for isolated BSs. While this approach is convenient for isolated solutions
with spatially constant frequency, a generalization to unequal-frequency
binary BSs with a spatially dependent frequency may be challenging..

The previous ansatz, based on assumption (30), requires specifying the binary
frequency $\tilde{\omega}$. An alternative is to specify some profile for the
scalar field (or conjugate momentum) and then assume only the equilibrium
conditions777Notice, these approaches are manifestly slice-dependent. Due to
the scalar time-dependence, the covariant equilibrium assumptions are
$\mathcal{L}_{\ell}[\bar{\Phi}\Phi]=0$, and
$\mathcal{L}_{\ell}[\bar{\Phi}\partial_{\mu}\Phi]=0$. This, in general,
implies $\mathcal{L}_{\ell}\rho\neq 0$ in the co-rotating frame.
$\mathcal{L}_{\ell}|\eta|^{2}=\mathcal{L}_{\ell}\rho=0$ [which are implied by
(30)]. Assuming these equilibrium conditions, the two equations (27) and (29)
reduce to

$\displaystyle D_{i}\left[\alpha(\bar{\Phi}D^{i}\Phi-\Phi
D^{i}\bar{\Phi})\right]=-(V^{i}D_{i}+\alpha
K)(\bar{\Phi}\eta-\Phi\bar{\eta}),$ (32)

as well as

$\displaystyle\begin{aligned} &\bar{\eta}D_{i}(\alpha D^{i}\Phi)+\eta
D_{i}(\alpha D^{i}\bar{\Phi})\\\ &=-(V^{i}D_{i}+2\alpha
K)|\eta|^{2}+\alpha(\bar{\eta}\Phi+\eta\bar{\Phi})\partial_{|\Phi|^{2}}V(|\Phi|),\end{aligned}$
(33)

respectively. Notice, both equations are real. In the parameterization of Sec.
II.3.3, it is natural to interpret these equations as elliptic equations for
either $\eta,\bar{\eta}$ or $\Phi,\bar{\Phi}$ (keeping the respective other
fixed). Here, no choice of the binary’s frequency is required, since the
$\tilde{\omega}$-dependence in (31) is canceled when adding (31) and its
complex conjugate to arrive at (33). Instead, the implicit assumption is that
the free data provides a sufficiently accurate profile for either the
conjugate momenta or scalar field to equilibriate the binary by adjusting the
other using (32) and (33). In principle, either approach could be applied to
constructing spinning binary BS initial data, assuming one starts with free
data that is sufficiently close to the desired solution.

Ultimately, only the direct implementation of these approaches may test their
applicability in reducing spurious perturbations and equilibriating the matter
in binary configurations, which we leave to future work.

##  III Quality of binary initial data

With the formalism in place to compute constraint satisfying binary initial
data, we now assess the quality of the constructed data for some specific
examples. Note, we are not equilibriating the scalar matter using the ansatz
discussed in the previous section. To that end, we compare the physical
properties of the free data with those of the constraint-satisfying initial
data in Sec. III.1. We proceed in Sec. III.2 by addressing the spurious
excitation of oscillation modes insides the BSs of the initial data, which
shows up in the emitted GW signal, by using two prescriptions to
systematically remove such artifacts. Finally, in Sec. III.3, we utilize the
standard procedure for reducing orbital eccentricity and discuss its
limitations in the context of binary BSs.

###  III.1 Superposed free data

Label | Coupling | $m$ | $\omega/\mu$ | $C_{i}$ | $S_{i}/M^{2}_{i}$ | Setup | $\alpha$
---|---|---|---|---|---|---|---
$\mathcal{B}_{1}$ | $\sigma=0.05$ | $1$ | $0.4$ | $0.12$ | $2.0$ | Axi. | $0$
$\mathcal{B}_{2}$ | $\lambda/\mu^{2}=10^{3}$ | $0$ | $0.9$ | $0.08$ | $0$ | Axi. | $0$
$\mathcal{B}_{3}$ | $\sigma=0.05$ | $1$ | $0.3,0.35$ | $0.17,0.14$ | $1.37,1.70$ | 3D | $\frac{\pi}{3}$
$\mathcal{B}_{4}$ | $\lambda/\mu^{2}=10^{3}$ | $0$ | $0.86,0.9$ | $0.12,0.08$ | $0$ | 3D | $0$
$\mathcal{B}_{5}$ | $\sigma=0.05$ | $0$ | $0.25$ | $0.13$ | $0$ | 3D | $\frac{\pi}{2}$
$\mathcal{B}_{6}$ | $\sigma=0.05$ | $1$ | $0.3$ | $0.17$ | $1.37$ | Axi. | -
Table 1: The properties of the isolated constituents of the binaries used
throughout the remainder of this work. The configurations with coupling
$\lambda$ are solutions in the repulsive scalar model (7), while those with
coupling $\sigma$ are stars in the scalar theory with the solitonic potential
(6). Here, $\omega$ is the star’s frequency, $m$ is the azimuthal index,
$C_{i}$ is the compactness, and $S_{i}$ and $M_{i}$ are its individual spin
angular momentum and mass, respectively. Binaries $\mathcal{B}_{1,2,5,6}$
consist of identical stars, whereas $\mathcal{B}_{3,4}$ are made of two stars
with different frequencies (and hence, masses, spins, etc.). The mass-ratio of
the last mentioned binaries are $q=1.43$ and $q=1.13$, respectively. We also
consider non-zero initial complex scalar phase offsets $\alpha$ (as defined at
the end of Sec. II.3.3) between the two stars. Note, in the case of
$\mathcal{B}_{6}$, we vary the phase offset $\alpha\in[0,\pi]$ in Sec. IV.1;
hence, we leave $\alpha$ unspecified here. In the axisymmetric setup, the
stars are boosted by the Newtonian free-fall velocity at the given coordinate
separation, whereas in the 3D context, the stars are initialized with quasi-
circular orbital frequency and spins aligned with the orbital angular
momentum; note, however, we also consider a binary with the parameters of
$\mathcal{B}_{3}$ with misaligned spins in detail in Sec. IV.3.

To assess the quality of the constructed initial data, we consider a series of
binary BS configurations. For now, the focus is entirely on configurations
obtained with $f_{(A)}=\hat{f}_{(A)}=1$, i.e., the case of superposed free
data, as defined in Secs. II.2 and II.3.3. The main properties of the
considered configurations are summarized in Table 1. For binaries
$\mathcal{B}_{1,2,3}$, we analyze the impact of solving the constraint
equations on the physical properties of the resulting binary BS initial data
as a function of coordinate separation. To that end, in the top panels of Fig.
2, we compare the ADM masses $M$, charges $Q$, and angular momenta $J$ of the
constraint-satisfying data to the corresponding quantities of the binary at
infinite separation labeled $M_{0}=M_{1}+M_{2}$ and $Q_{0}=Q_{1}+Q_{2}$. The
physical properties of the initial data differ from the superposed
configuration at infinite separation by up to $\sim 10\%$; these differences
roughly scale as $\sim 1/D_{0}$ with the coordinate separation $D_{0}$ of the
binary. Furthermore, in the bottom left panel of Fig. 2, we show how much the
conformal factor $\Psi$ differs from unity. Generally, this difference is at
most a few percent for reasonable separations: $\max|\Psi-1|\lesssim 0.02$.

Figure 2: The properties of the constraint satisfying binary BS initial data
as a function of coordinate separation $D_{0}$ of the stars. Here, $M$ is the
ADM mass, $J$ the angular momentum (defined with respect to the axisymmetric
Killing field), and $Q$ the initial charge of the binary configurations with
properties summarized in Table 1. These are compared with the corresponding
quantities at infinite separation of the binary (e.g., $M_{0}=M_{1}+M_{2}$).
The amplitude of the spurious oscillations in the stars emerging during the
evolution of these binary initial data is defined in (41). Dotted lines
indicate the $\propto 1/D_{0}$ fall-off matched to the point with the largest
separations.

Since the scalar matter is not equilibriated, solving the constraint equations
(LABEL:eq:constraints) in this form with superposed free data leads to
spurious oscillations in the constituents of the binary. In the case of binary
neutron stars, these artifacts are identified by monitoring the central
density of the stars during the subsequent evolution of the initial data.
Here, we proceed analogously by tracking the global maximum of the magnitude
of the scalar field $\max|\Phi|$ throughout the first few oscillation periods
$T_{0}$ of excited modes in the stars. Specifically, we focus on those
oscillation modes of the normalized maximum on each time slice
$\phi_{m}(t)=\max|\Phi|_{t}/\max|\Phi|_{t=0}$ and quantify the amplitude of
these perturbations with

$\displaystyle\Delta\phi_{m}=\frac{\max_{t\in[t_{0},t_{0}+T_{0}]}\phi_{m}(t)}{\min_{t\in[t_{0},t_{0}+T_{0}]}\phi_{m}(t)}-1.$
(41)

In the case of binaries in the repulsive scalar model, measuring
$\Delta\phi_{m}$ with $t_{0}=0$ suffices, while for binary configurations in
the solitonic scalar theory, we typically extract $\Delta\phi_{m}$ with
$t_{0}>5T_{0}$ (once the binary settles into the dominant oscillation mode).
We find the maximum of the gauge-dependent U(1)-charge density to track these
oscillations equally well. In the bottom right panel of Fig. 2, we show
$\Delta\phi_{m}$ for different binary configurations and initial separations.
The amplitude of these oscillations increases with decreasing initial binary
separation. At large separations, $\Delta\phi_{m}\rightarrow 0$, indicating
that the construction of the initial data (without assuming matter
equilibrium) excites these oscillation modes. Notably, the magnitude of
$\Delta\phi_{m}$ is much smaller in the case of binaries in the solitonic
scalar model, compared with the repulsive scalar model. In the following
section, we analyze spurious oscillations of this kind in more detail, and
propose and test methods to help mitigate these effects.

###  III.2 Spurious Oscillations

We have seen that the naive choice for the metric and scalar free data, i.e.,
the superposed free data, leads to potentially significant spurious
oscillations in the individual stars in the subsequent evolution of the
initial data. To address this issue, it is instructive to consider possible
physical mechanisms unique to scalar BSs causing these artifacts. The
fundamental feature rendering the fluid star and the BS cases distinct is that
the microphysical scales of the latter are macroscopic, leading to wave-like
phenomena on scales of the star itself. Specifically, the BS can be thought of
as composed of a collection of bosons with Compton wavelength $\lambda\sim
1/\mu$ satisfying $M/\lambda\sim CR/\lambda\sim\mathcal{O}(1)$, since
$C\sim\mathcal{O}(0.1)$ in the relativistic regimes relevant for this work
(see, e.g., Fig. 1). Therefore, there may be distinct processes active in the
context of binary BSs affecting the quality of the initial data.

Self-gravitating solitonic solutions such as BSs consist of a single coherent
gravitationally bound state of bosons with energy888Note, here and in the
following, we use boson “energy” and “frequency” interchangeably, implicitly
setting $\hbar=1$. $\omega<\mu$. The marginally bound scenario, $\omega=\mu$,
separates the bound states from unbound and asymptotically free states with
energies $\omega>\mu$. Stationary isolated BSs are solutions with bosons of
energy $\omega$ precisely in equilibrium with the gravitational field.
However, perturbations introduced by superposing two stationary BSs and
solving the Hamiltonian and momentum constraints based on such non-
equilibriated free data disrupts this balance. Perturbations may elevate some
fraction of the bound bosons of energy $\omega$ to (i) completely free states
(dispersing away from the binary), (ii) states that are gravitationally bound
to the binary (as opposed to one of the constituents); analogous to a wave
dark matter halo with solitonic core (see, e.g., Ref. Hui (2021)), or (iii)
states gravitationally bound to a single star, but with energy $\omega$ that
is not at equilibrium with the gravitational sector. All these processes are,
in principle, able to excite oscillation modes inside the star, as well as
contaminate the GW signal at early times in a numerical evolution. Note, non-
linear scalar self-interactions may also disturb equilibrium configurations.
However, since since we mainly focus on high compact stars with
$\omega\ll\mu$, these effects are suppressed (see, e.g., Ref. Siemonsen and
East (2023)).

With these possible sources in mind, in the following we explore several
prescriptions for reducing spurious oscillations in binary BSs. We first
outline methods commonly used in the context of binary black hole and neutron
star initial data and discuss their effectiveness in the BS context, before
introducing and validating several methods specific to BSs.

####  III.2.1 Modified superposition

Figure 3: We show the temporal evolution of the normalized maximum on each
time slice $\phi_{m}=\max|\Phi|/\max|\Phi|_{t=0}$ and amplitude of the
perturbations $\Delta\phi_{m}$, defined in (41), for selected binary BS
configurations. In all cases shown here, the exponent $\gamma$, defined in
(42), is set to $\gamma=2$. (left top) The behavior of the maximum $\phi_{m}$
during the evolution of the binary $\mathcal{B}_{4}$ (see Table 1) with
initial coordinate separation $D_{0}=20M_{0}$, constructed with
$f_{(A)}=\hat{f}_{(A)}=1$ (labeled “superposition”) contrasted with the case,
where $\sigma_{(A)}/D_{0}=0.52$ and $\hat{\sigma}_{(A)}=0$ (labeled “modified
sup.”). (left bottom) The amplitude $\Delta\phi_{m}$ of the spurious
oscillations as functions of the lengthscale $\hat{\sigma}$. Cases indicated
with “(g)” correspond to only metric attenuation, $\hat{\sigma}_{(A)}=0$ and
$\hat{\sigma}=\sigma_{(A)}$, while for those labeled “(g+s)” both the metric
and scalar degrees of freedom are attenuated
$\hat{\sigma}=\hat{\sigma}_{(A)}=\sigma_{(A)}$. The binary $\mathcal{B}_{2}$
has initial coordinate separation $D_{0}=40M_{0}$. (top right) The behavior of
the maximum $\phi_{m}$ during the evolution of the binary $\mathcal{B}_{3}$
(see Table 1) with initial coordinate separation $D_{0}=12M_{0}$ (constructed
with $f_{(A)}=\hat{f}_{(A)}=1$) and rescaling the conformal kinetic energy in
(43) with $p=0$ as well as $p=-4$. (bottom right) We show the amplitude
$\Delta\phi_{m}$ of spurious oscillations emerging during the evolution of
binaries $\mathcal{B}_{1,3}$ with initial coordinate separations
$D_{0}=40M_{0}$ and $D_{0}=12M_{0}$, respectively. For the latter, we were
unable to construct binary BS initial data with $p<-4$.

We begin by returning to the largest spurious oscillations in the bottom right
panel of Fig. 2, i.e., those in the head-on collisions of non-rotating stars
in the repulsive scalar model. As we demonstrate below, some of these spurious
oscillation artifacts can be removed by choices of the attenuation functions
$f_{(A)}$ and $\hat{f}_{(A)}$, introduced in Secs. II.2 and II.3.3,
respectively. In the case of binary black hole or neutron star initial data,
it is common practice to remove the metric variables of one star at the
location of the other (analogous to approaches introduced in Refs. Marronetti
and Matzner (2000); Marronetti  _et al._ (2000); Bonning  _et al._ (2003)),
achieved by non-trivial attenuation functions $f_{(A)}$. This can reduce the
effect of the superposition on the individual stars. In order to remove the
metric and scalar free data due to one of the stars at the coordinate location
of the other, we choose

$\displaystyle
f_{(1)}(\textbf{x})=1-\exp\left[-\frac{|\textbf{x}-\textbf{z}_{(2)}|^{\gamma}}{\sigma_{(1)}^{\gamma}}\right],$
(42)

and the corresponding $(1)\leftrightarrow(2)$, as well as scalar attenuation
functions $\hat{f}_{(A)}$ with associated length scales $\hat{\sigma}_{(A)}$.
Here, $\textbf{z}_{(2)}$ is the initial coordinate position of the second star
(as defined in Sec. II.2), whereas the length scales $\sigma_{(A)}$ and
$\hat{\sigma}_{(A)}$ set the size of the attenuation region around each of the
constituents of the binary. We consider $\gamma=2$ and 4.

We find this approach to be effective in reducing spurious oscillations, as
measured by $\Delta\phi_{m}$, only for BS solutions in the repulsive scalar
model (7). In fact, applying this approach to stars in the solitonic models
worsens the artificial oscillations, requiring a different prescription to
handle the latter, as described in detail in the next section. In the left
panels of Fig. 3, we show the time-dependence of the maximum of the scalar
field magnitude, as well as the dependence of the amplitude of the spurious
oscillations on the lengthscales associated with the attenuation function
introduced in Eq. (42). As can be seen from the figure, increasing the
attenuation length scales (at fixed separation), decreases the amplitude
$\Delta\phi_{m}$. Around $\sigma_{(A)}/D_{0}\approx 0.5$, the spurious
oscillations are minimal, and increase in amplitude for
$\sigma_{(A)}/D_{0}\gtrsim 0.5$. Therefore, we find that in all cases
considered, tuning the length scales relevant in (42) results in binary BS
initial data with significantly reduced spurious oscillations. Finally,
considering $\gamma=4$ (as opposed to $\gamma=2$) results in no qualitative
difference to the behavior shown in left panels of Fig. 3.

Besides improving the binary BS initial data by removing spurious
oscillations, this modified superposition approach turns out to be necessary
in the case of highly compact binary inspirals at moderate separations within
the repulsive scalar model (7). Specifically, we focus on a binary BS with
$\omega_{1}=\omega_{2}=0.86\mu$, at initial coordinate separation
$D_{0}/M_{0}=20$ and quasi-circular boost velocities. For this, we find that
superposed free data, i.e., with $f_{(A)}=\hat{f}_{(A)}=1$, results in
premature collapse of each individual star to a black hole after $\approx
50M_{0}$. In contrast, with $\sigma_{(A)}=D_{0}/2$, the individual stars
remain stable throughout the inspiral of length $\sim\mathcal{O}(10^{3})M_{0}$
up to the point of contact. On the other hand, in Ref. Siemonsen and East
(2023), we found the simple choice $f_{(A)}=\hat{f}_{(A)}=1$ to be sufficient
to successfully evolve the binary with $\omega_{1}=\omega_{2}=0.9\mu$ in the
same family of solutions; hence, the attenuation is necessary for more
relativistic BS solutions. Note, similar premature collapse was observed in
Ref. Helfer  _et al._ (2019).

####  III.2.2 Conformally rescaled kinetic energy

Figure 4: The GW amplitude (top) and phase (bottom) emitted by the (aligned-
spin) binary $\mathcal{B}_{3}$ with initial separation $D_{0}=12M_{0}$. The
two curves correspond to the evolution of initial data constructed using $p=0$
or $p=-4$ in (43). Specifically, the GW phase $\phi_{\rm GW}$ is the complex
phase of the Newman-Penrose scalar $\Psi_{4}$, whereas the GW amplitude is the
magnitude of the projection of $\Psi_{4}$ onto the $(\ell,m)=(2,2)$ mode of
$s=-2$ spin-weighted spherical harmonics on a coordinate sphere at radius
$r=100M_{0}$. Note, the binary orbits are eccentric with eccentricity
$e\approx 0.1$, resulting in modulations with period $\approx 250M_{0}$ as is
most striking in the top panel (eccentricity reduction is discussed separately
in Sec. III.3).

As alluded to above, we find attenuating the metric and scalar free data to
only be beneficial in binary BSs in the repulsive scalar theory. For the
solitonic case, we return to the notion of energy levels determined by the
frequency $\omega$ of the scalar bound state. In the context of the scalar
variables, this frequency is set by $\partial_{t}\Phi$, which enters the
kinetic energy $\sim|\eta|^{2}$ according to (17). For instance, for isolated
BSs $\partial_{t}\Phi=i\omega\Phi$, and similarly
$\eta\sim(\omega-m\Omega/r)\Phi$ (in the coordinates introduced in Sec. II.1).
To reduce the impact of the superposition on the BSs, as discussed in Sec.
III.2, we modify the kinetic energy of the scalar field999Note, in principle,
one could modify $\partial_{t}\Phi$ iteratively, even when working with the
parameterization of Sec. II.3.3.. The kinetic energy combines changes in the
frequencies of the individual stars with changes in the local linear and
angular momentum due to the orbital motion and the star’s intrinsic spin. As
such, increasing or decreasing the kinetic energy locally self-consistently
within the CTS setup may help address spurious oscillations. To incorporate
this in our initial data construction scheme discussed in Sec. II.3.3, we
rescale the physical conjugate moment $\eta$ by an additional power $p$ of the
conformal factor:

$\displaystyle\eta=\Psi^{-6-p}\tilde{\eta}.$ (43)

In the right panel of Fig. 3, we illustrate the impact of the choice (43) for
different $p$ on the spurious oscillations in the individual stars of two
types of binaries in the solitonic scalar model (in all binary BSs in the
repulsive scalar model, we set $p=0$). In the axisymmetric binary labelled
$\mathcal{B}_{1}$, we find that the amplitude of the spurious oscillations can
be minimized by adjusting the exponent $p$. In this case, a rescaling (43)
with $-4<p<-3$ minimizes the spurious oscillations measured by
$\Delta\phi_{m}$. These oscillations can be addressed also in the case of the
inspiraling binary $\mathcal{B}_{3}$; however, our numerical implementation
robustly relaxes into a solution to the constraint equations only for
$p\geq-4$. The oscillation amplitude is still decreasing with decreasing $p$
for $p=-4$, and this leaves a residual oscillation amplitude roughly a factor
of $5$ smaller compared with the $p=0$ initial data. While the spurious
oscillation amplitude of all binaries in the solitonic scalar model considered
is small, i.e., $|\Delta\phi_{m}|\sim 10^{-3}$, there is a correlation between
reducing these small artifacts and removing a high-frequency contamination
from the emitted gravitational waveform.

In Fig. 4, we show the GW amplitude and phase extracted from the binary
evolution of $\mathcal{B}_{3}$ constructed using (43) with $p=0$. We contrast
this with the signal extracted from the same binary, but with initial data
constructed using $p=-4$. Similar to the amplitude of the spurious
oscillations described above, the high-frequency and large amplitude
contamination of GWs at early times changes with the exponent $p$. As the
period of the high-frequency contamination matches the period of the spurious
oscillations in Fig. 3, we identify the latter as a source of the GW
contamination. Therefore, the adjustment of the kinetic energy by the
conformal rescaling (43) aids in reducing GW contamination. However, as for
the spurious oscillations, the trend of the magnitude of this contamination
with decreasing exponents $p<0$ suggests that if initial data with $p<-4$
could be constructed, the contamination would be further inhibited. In this
case, we find that the amplitude of the contamination is reduced by roughly a
factor of $2$ compared with the canonical choice $p=0$. Additionally, the
evolution of the GW phase shown in Fig. 4 indicates that in the $p=0$ case the
amplitude is dominated by the high-frequency contamination, while in the
$p=-4$ case, the phase evolution is determined primarily by the binary orbit,
as opposed to spurious oscillations. Finally, keeping the superposed free data
fixed, but varying the parameter $p$, results in constraint satisfying initial
data with different physical kinetic energy and momentum densities [according
to (43)]. This may lead to different orbital parameters for the binary,
explaining the shorter time-to-merger comparing the $p=0$ with the $p=-4$ case
in Fig. 4.

###  III.3 Eccentricity reduction

The flexibility of our approach to constructing the initial data allows us to,
in principle, find constraint satisfying binary BS data resulting in any
orbital configuration. Of particular interest in the case of compact binary
inspirals are low-eccentricity orbits. Hence, in a last step to improve the
quality of our binary BS initial data, we turn to applying common techniques
to reduce the eccentricity of compact binary initial data to the binary BSs
constructed here. To that end, we first define a notion of the BS coordinate
location valid throughout the evolution, we then briefly review the
eccentricity reduction methods introduced in Refs. Pfeiffer  _et al._ (2007);
Boyle  _et al._ (2007); Mroue  _et al._ (2010), and finally, we apply these
methods to selected spinning and non-spinning binary BSs.

To define the coordinate positions of the two BSs during a binary evolution,
we employ a two-step procedure: first, a rough estimate of the star’s position
restricted to the initial equatorial plane is obtained by finding the
coordinate locations of the local maxima of $|\Phi|$ for spherically symmetric
stars, and local minima in the case of rotating BSs associated with each star
(i.e., the intersection of the star’s vortex lines and the equatorial plane).
In a second step, the center of scalar field magnitude within a coordinate
sphere $B_{A}$ centered on the previously determined locations of extrema
enclosing star $A$,

$\displaystyle z^{i}_{(A)}(t)=\int_{B_{A}}d^{3}x|\Phi(t,\textbf{x})|x^{i},$
(44)

is used as the coordinate position at the given coordinate time $t$101010Note,
in the initial timeslice $z^{i}_{(A)}(0)$ agrees with the coordinate positions
$z^{i}_{(A)}$ defined in Sec. II.2 to a large degree; hence, we neglect any
potential difference between the two in the following.. The coordinate
separation $d(t)$ of the binary is then simply
$d(t)=|\textbf{z}_{(1)}(t)-\textbf{z}_{(2)}(t)|$. Note, the first step is
sufficient for stars in the repulsive scalar potential (7), as the locations
of the extrema are less prone to contamination by spurious oscillations within
the stars. In the case of BSs in the solitonic models, however, we find the
second step to be crucial particularly for non-spinning stars, as $|\Phi|$ is
roughly constant inside each star. In the remainder of this work, we use (44)
to define the BS position.

The general procedure to reduce eccentricity, following Ref. Pfeiffer  _et
al._ (2007), is to begin with a set of binary initial data, evolve these for a
sufficiently long time to be able to confidently fit for the binary’s orbital
parameters, and then compute a correction to the initial radial velocity and
orbital frequency to construct new initial data with lower eccentricity. This
process is repeated until the desired eccentricity is reached. Specifically,
we fit for the binary BS coordinate separation $d(t)$ using

$\displaystyle\hat{d}(t)=A_{-1}+A_{0}t+\frac{A_{1}}{2}t^{2}+\frac{B}{\omega}\sin(\omega
t+\varphi),$ (45)

and correspondingly for the binary’s radial velocity using the fit

$\displaystyle\hat{v}_{r}(t)=A_{0}+A_{1}t+B\cos(\omega t+\varphi).$ (46)

Based on this parameterization, the initial orbital frequency $\Omega_{0}$ and
initial radial velocity $v_{r}$ are corrected by Pfeiffer  _et al._ (2007)

$\displaystyle\begin{aligned}
\Omega_{0}&\rightarrow\Omega_{0}+\frac{B\omega\sin(\varphi)}{2\Omega_{0}d_{0}},\\\
v_{r}&\rightarrow v_{r}-B\cos(\varphi),\end{aligned}$ (47)

at each eccentricity iteration step, where $d_{0}=A_{-1}$111111Note, in
general $d_{0}\neq D_{0}$, as we show below explicitly.. Working entirely in
flat space, in order to translate these orbital parameters to the initial
coordinate positions $z^{i}_{(1)}$, $z^{i}_{(2)}$ and velocities
$v^{i}_{(1)}$, $v^{i}_{(2)}$ of the binary constituents, as defined in Sec.
II.2, we utilize the Newtonian center-of-mass expressions

$\displaystyle z^{i}_{(1)}=r^{i}\frac{M_{2}}{M_{0}},\qquad
z^{i}_{(2)}=-r^{i}\frac{M_{1}}{M_{0}}.$ (48)

Here, $r^{i}$ is the binary separation with $v^{i}=\partial_{t}r^{i}$, and the
corresponding velocities $v^{i}_{(1)}$ and $v^{i}_{(2)}$ are given by the time
derivatives of the above expressions. We decompose the center of mass velocity
$v^{i}$ into radial $v_{r}$ and tangential $v_{t}$ components as
$v^{i}=v_{r}n^{i}+v_{t}\lambda^{i}$, where $n^{i}=r^{i}/r$ and
$\lambda^{i}=\partial_{t}n^{i}$. The initial tangential component of the
center-of-mass velocity is set by the initial orbital velocity as well as
initial coordinate separation $d_{0}$ using $v_{t}=\Omega_{0}d_{0}$. The
constituents velocities are then reconstructed from $v^{i}$ utilizing the time
derivative of the expressions (48). In this context, the orbital eccentricity
is defined to be

$\displaystyle e=\frac{B}{\omega d_{0}}.$ (49)

The formulas (47) are based only on Newtonian gravity, with radiation reaction
and other general relativistic corrections assumed to be absorbed into the
linear and quadratic time dependence in (45); we note that BSs will also have
scalar interactions that become important during the late inspiral Siemonsen
and East (2023), particularly for stars with small compactness and
$\omega\approx\mu$, which may introduce extra complications in performing
eccentricity in this way at small separations.

Figure 5: (top panel) The initial radial velocity component $v_{r}$ and
orbital angular frequency $\Omega_{0}$ of the binary BS initial data
corresponding to $\mathcal{B}_{3}$ and $\mathcal{B}_{5}$ (see Table 1),
throughout the eccentricity reduction procedure, starting from iteration step
0. Here, $\Omega_{\infty}=(D_{0}^{3}/M_{0})^{1/2}$. (bottom panel) The
associated eccentricity, defined in (49), as a function of iteration step
$N_{e}$. We only perform a single iteration step for the $\mathcal{B}_{3}$
binary with misaligned spins (this case is discussed further in Sec. IV.3).

Applying this machinery to binary BSs, we find that spurious oscillations of
the stars and non-equilibriated scalar matter result in high-frequency
oscillations of the coordinate separation $d(t)$, limiting the eccentricity
reduction. Below eccentricities of $e\sim 10^{-2}$, the fit $\hat{v}_{r}(t)$
is too uncertain to confidently extract estimates for the subsequent iteration
step. Instead, in these cases, we resort to using (45) to determine the
corrections (47). For $e\lesssim 10^{-3}$, the amplitude of the modulation of
$d$ due to residual eccentricity reaches the amplitude of the oscillations in
$d$ introduced by these spurious oscillations. Therefore, even the fit
$\hat{d}(t)$ to the coordinate separation $d$ becomes uncertain, and we
terminate eccentricity reduction at $e\gtrsim 10^{-3}$. Furthermore, it is, of
course, crucial to minimize oscillations of the stars using the methods
discussed in Sec. III.2, i.e., find the exponent $p$ of (43) and length scales
$\\{\sigma_{(A)},\hat{\sigma}_{(A)}\\}$ in (42), before attempting to reduce
the eccentricity. Especially the choice of (43) alters the matter’s kinetic
energy and linear momentum, and hence, the orbital frequency and radial
velocities.

In Fig. 5, we show the orbital parameters throughout the iterative reduction
of the orbital eccentricity for two example binaries. In the case of the
spinning and unequal-mass binary $\mathcal{B}_{3}$, the eccentricity decreases
exponentially with the iteration step $N_{e}$ down to $e\sim 10^{-3}$, where
the fitting approach of (45) becomes unreliable. In the case of the non-
spinning and equal-mass binary $\mathcal{B}_{5}$, however, the eccentricity
decreases only slowly with each iteration step. Consulting the top panel of
Fig. 5, the first iteration step for $\mathcal{B}_{5}$ resulted in a positive
radial velocity. While $v_{r}$ is a coordinate quantity, and hence, carries
limited physical meaning, in binary black hole and neutron star initial data
constructions, this is found to consistently satisfy $v_{r}<0$. Both may be
due to fitting to the $d(t)$ time series before all spurious perturbations of
the initial data have decayed away sufficiently (e.g., fits of the $N_{e}=1$
iteration contained only roughly $3/2$ orbits), beyond which we cannot isolate
a cause of the slow convergence of $e$ for $\mathcal{B}_{5}$. Finally, in the
case of the precessing $\mathcal{B}_{3}$ (details can be found in Sec. IV.3),
we perform only a single iteration step and find similar convergence behavior
as for the aligned-spin binary $\mathcal{B}_{3}$. Note, in precessing cases
the spin-interactions (particularly for super-spinning compact objects as
considered in $\mathcal{B}_{3}$) induce physical oscillations of the binary
separation beyond residual eccentricity (see, e.g., Ref. Buonanno  _et al._
(2011)), which we, however, have not observed here. Lastly, in Appendix C, we
discuss the linear motion of the center of mass away from the center of the
numerical grid, and how to iteratively reduce this artifact, while
simultaneously reducing eccentricity.

##  IV Binary evolutions

In this section, we illustrate our method for constructing binary BS initial
data in the context of head-on collisions and quasi-circular inspirals,
including a precessing system. In the case of head-on collisions, we
demonstrate explicitly that equal-mass rotating binary BSs exhibit solitonic
behavior (resembling the dynamics found in Ref. Palenzuela  _et al._ (2007)),
bouncing off each other when colliding along their spin axes, if the phase
offset between the stars is precisely $\pi$. Furthermore, we consider two
eccentricity-reduced, quasi-circular inspiraling binary BSs: one non-spinning
equal-mass binary, and one super-spinning unequal-mass configuration. We
analyze their inspiral dynamics, show that non-trivial scalar interactions
result in strong deviations from the dynamics of binary black holes or neutron
stars (analogous to what was found in Ref. Siemonsen and East (2023)), and
characterize the GW strain. Finally, we consider a super-spinning, precessing
binary BS inspiral at moderately low eccentricity, analyze the merger
dynamics, and show the resulting gravitational waveform.

To that end, we employ the methods developed above in Sec. II. In particular,
we use the source parameterization introduced in Sec. II.3.3 (the scalar
matter is not equilibriated with approaches outlined in Sec. II.4). Since we
focus on BSs in the solitonic scalar model, we set the attenuation
coefficients (introduced in Sec. III.2.1), as well as the conformal rescaling
parameter $p$ (introduced in Sec. III.2.2), to zero unless otherwise stated.

###  IV.1 Head-on collisions

In this section, we explore the merger dynamics of two rotating BSs during a
head-on collision along their respective spin axes, focusing on the
$\sigma=0.05$ solitonic scalar model (6) and binaries composed of $m=1$
rotating BSs. It is important to note that, in this setting, we evolve the
binary BSs using a generalized Cartoon method, which explicitly assumes the
scalar fields azimuthal dependence follows $\Phi\sim e^{im\varphi}$, in
addition to an axisymmetric metric (see Appendix A for details). As a result,
any modes violating this condition will not appear in the evolution. In
particular, this implies (i) any non-axisymmetric instability of the form
found in Ref. Sanchis-Gual  _et al._ (2019b) is suppressed, and (ii) the
vortex structure of the solution on the symmetry axis is conserved throughout
the evolution.

We begin by considering the predictions of the remnant map introduced in Ref.
Siemonsen and East (2023) for these head-on collisions. We do not repeat the
details of the construction of the remnant map here (which are found in Ref.
Siemonsen and East (2023)), and simply summarize the key features in the
context of the head-on collision of two $m=1$ BSs. This map assumes
U(1)-charge conservation ($Q_{\rm rem}=Q_{1}+Q_{2}$) during the merger of two
BSs to predict the qualitative and quantitative features of the merger
remnant. In order to use the remnant map, one must provide a plausible
candidate family of remnants. Due to our evolution methods, any binary
composed of two $m=1$ stars results in a remnant with $m=1$ vortex along the
symmetry axis. Therefore, it is natural to consider that the remnant of the
two rotating BSs is another rotating BS of the same vortex index (unless the
combined charge of the binary surpasses the maximum charge of the family of
$m=1$ BSs in this scalar model, in which case the remnant is likely a black
hole). In our axisymmetric setup, any known linear instability of the rotating
BSs in the $\sigma=0.05$ solitonic scalar model is suppressed, implying that
this condition allows the merger remnant to be a $m=1$ BS. Hence, we can map
all properties of any given $m=1$ binary BS, parameterized by their
frequencies $\omega_{1}$ and $\omega_{2}$, into the properties of the
resulting $m=1$ BS assuming charge conservation. In particular, in order to
consider the kinematics of the merger—and whether this favors the formation of
a $m=1$ BS remnant—we define the relative mass difference Siemonsen and East
(2023)

$\displaystyle\mathcal{M}=\frac{M_{1}+M_{2}-M_{\rm rem}}{M_{1}+M_{2}}.$ (50)

Here, $M_{1}$ and $M_{2}$ are the masses of the individual stars, while
$M_{\rm rem}$ is the mass of the rotating BS with charge $Q_{\rm
rot}=Q_{1}+Q_{2}$ obtained form the remnant map. If a binary configuration has
$\mathcal{M}>0$, then the formation of that $m=1$ BS remnant of mass $M_{\rm
rem}$ is energetically favorable.

Figure 6: The remnant properties of a $m_{1}=m_{2}=1$ spinning binary BS head-
on collision (with frequencies $\omega_{1}$ and $\omega_{2}$) in the
$\sigma=0.05$ solitonic scalar model (6), assuming the remnant is a $m=1$
rotating BS and using the remnant map of Ref. Siemonsen and East (2023) (i.e.,
assuming U(1)-charge conservation to predict the remnants properties for each
given binary configuration). The dimensionless angular momentum $J_{\rm
rem}/M_{\rm rem}^{2}$ and the associated remnant compactness $C_{\rm rem}$ are
shown as contours. The mass ratio $\mathcal{M}$ is defined in (50). We
restrict to the $\omega_{2}>\omega_{1}$ portion without loss of generality,
and indicate the regions with $Q_{\rm rem}>Q_{\max}$, where $Q_{\max}$ is the
maximum charge of the $m=1$ family of BSs, in black. Finally, we classify the
merger remnants into black holes (“BH”) and $m=1$ rotating BSs (“BS”), and
mark the binary that gave rise to the respective remnant with dots. [Notice,
the central (left) “BS” corresponds to binary $\mathcal{B}_{1}$
($\mathcal{B}_{6}$), see Table 1.]

In Fig. 6, we show the relative mass difference $\mathcal{M}$ across the
relevant binary BS parameter space. For all binary configurations shown, it is
indeed energetically favorable to form a $m=1$ rotating BS remnant after
merger. Hence, the remnant map predicts the merger remnant to be a $m=1$
rotating BS remnant. Note, if the initial phase offset $\alpha$ between the
stars is exactly $\alpha=\pi$, then the remnant is not a (parity-even) $m=1$
rotating BS, but rather a double rotating $m=1$ star (i.e., a single, parity-
odd rotating star), as we discuss in detail below.

To test this prediction, we perform a series of numerical evolutions of binary
BSs in the axisymmetric setting described above. The initial data is
constructed as discussed in the previous sections, where for simplicity, the
conformal rescaling power in (43) is set to $p=0$, and no modification of the
form give in (42) is applied. The initial velocities are chosen based on the
Newtonian free-fall velocity from rest at infinity, and the binary separation
is chosen to be $D_{0}=10M_{0}$ initially. Finally, here the phase offset
$\alpha$ (defined at the end of Sec. II.3.3) between the phases of the
rotating stars is set to vanish, $\alpha=0$. (We consider scenarios varying
the value of $\alpha$ below.) For each of the evolutions, we classify the
remnants as either BSs ($m=1$ rotating BSs) or spinning black holes. In Fig.
6, we show the binary configurations we numerically evolve, and indicate the
remnant type. First, in the case of equal-mass binaries, i.e., those with
$\omega_{1}=\omega_{2}$, the remnant is consistent with the prediction of the
remnant map, except for the case with $\omega_{1}=\omega_{2}=0.2\mu$. In this
case, and all other cases indicated as “BH” in Fig. 6, the merger product
collapses to a black hole during the nonlinear merger process. The fact that
the threshold for black hole formation is slightly lower than predicted by the
remnant map is likely due to the extra compression and kinetic energy due to
the collision. This explicitly demonstrates that the final remnant of the
head-on collision of two rotating BSs along their mutual vortex line results
in another rotating BS of the same vortex number (or a black hole if the
individual stars are sufficiently compact). In Ref. Nikolaieva  _et al._
(2022), a similar analysis was performed in the Newtonian limit dropping all
symmetry assumptions. Since they find that the central vortex line persists
throughout the merger, their results are consistent with our findings here,
and suggest that the latter generalize to 3D if linear instabilities are
absent in both the merging BSs and the remnant BS.

Figure 7: We show scalar field quantities on axial slices for the head-on
collision of binary $\mathcal{B}_{6}$ with initial phase-offset of
$\alpha=\pi$. The $z$-axis is the spin-axis and $z=0$ corresponds to the
center of mass of the system. (left) Surfaces of constant scalar field
magnitude at various times during the head-on collision. (right) The solution
after it has relaxed at late times ($t/M_{0}=856$). The black contour line
indicates the surfaces of constant scalar field magnitude, while the color
indicates the phase $\psi$ of the scalar field. This end state resembles the
parity-odd stationary double-BS solutions found in Refs. Kleihaus  _et al._
(2008); Gervalle (2022).

Let us now return to considering head-on collisions of two $m=1$ rotating
binary BS configurations while varying the initial phase offset $\alpha$. We
perform a series of evolutions of the binary configuration $\mathcal{B}_{6}$
(with initial separation of $D=10M_{0}$ and Newtonian free-fall velocities, as
before), with initial phase offsets of $\alpha/\pi\in\\{1/4,1/2,3/4,8/9,1\\}$.
Notice, the $\alpha=0$ case was found to result in a rotating $m=1$ BS as
indicated in Fig. 6. In Fig. 7, we show snapshots of the head-on collision of
the binary with maximal phase offset, i.e., $\alpha=\pi$. As evident there,
contrary to the expectation from the $\alpha=0$ scenario, the two stars bounce
off each other upon contact. After several bounces, the system settles into a
state of two spatially separated scalar distributions with a persistent phase
offset of $\pi$, as shown in the right panel of Fig. 7. The end state of this
merger is plausibly a stationary solution analogous to those found in Refs.
Kleihaus  _et al._ (2008); Gervalle (2022). These are parity-odd solutions
resembling two rotating BSs, where gravitational attraction is balanced by
scalar interactions. Physically, this solitonic behavior resembles the
dynamics during the head-on collision of two non-spinning BSs with phase
offset $\alpha=\pi$ reported in Ref. Palenzuela  _et al._ (2007) (and
associated stationary solutions of Ref. Yoshida and Eriguchi (1997b)).

Moving to the cases with $\alpha/\pi\neq 1$, we find that the remnant of the
corresponding head-on collision is always a single $m=1$ rotating BS at late
times. In the case of $\alpha/\pi=8/9$, the system performs a single bounce
upon collision, but then merges to a perturbed rotating BS. This demonstrates
that, similar to the head-on collision scenario of two non-spinning stars,
only the $\alpha=\pi$ configuration exhibits a final state different from a
$m=1$ rotating BS. This, of course, limits the validity of the remnant map to
those cases with $\alpha/\pi\neq 1$; however, as this is an edge case (similar
to two non-spinning stars with phase-offset of $\alpha=\pi$), the impact on
the applicability of the remnant map even in these head-on scenarios is
minimal.

###  IV.2 Quasi-circular binaries

Figure 8:  Trajectories and GWs from two quasi-circular binary BS mergers: an
equal-mass, non-spinning binary labelled $\mathcal{B}_{5}$ (top row), and a
unequal-mass aligned spin case $\mathcal{B}_{3}$ (bottom row). (The binary
properties are summarized in Table 1.) We show the BS trajectories in the
orbital plane [defined according to (44)] throughout the evolution up to the
point of contact (left panels), the binary coordinate separation $d$ as a
function of time (center panels), and the emitted $(\ell,m)=(2,2)$
$s=-2$-weighted spherical harmonic component of the GW strain $h_{+}^{(2,2)}$
(right panels). In the left and center panels, we indicate the initial time
(round markers) and the time of contact (square markers). In the bottom left
and center panels, we also indicate the time $t/M_{0}=1475$ (x-marker) and
show the subsequent in-plane star trajectories in the inset in the bottom left
panel up to the point of contact, which is discussed in the main text. The
dimensionless spins, $\chi$, of the binary constituents can be found in the
bottom right corners of the left panels. The legends in the center panels
indicate the eccentricity reduction step $N_{e}$ and the rescaling exponent
$p$ used in (43). For comparison, in the bottom-center panel only, we also
show two cases with different values of $p$ before eccentricity reduction
($N_{e}=0$). After eccentricity reduction, $e\approx 4\times 10^{-3}$
initially for the top row, and $e\approx 1.8\times 10^{-3}$ for the bottom
row. In the center top panel, we also show the time derivative of the
separation $d$ in the inset. Notice, the GW strain (right panels) contains
residual high-frequency contamination as discussed in Sec. III.2; this
contamination is shown in detail in Appendix B. At the point of contact of the
two stars, the GW frequency is roughly $\omega_{c}M_{0}\approx 0.1$ both for
$\mathcal{B}_{5}$ and $\mathcal{B}_{3}$). For the post-merger phase of binary
$\mathcal{B}_{5}$, we estimate the dominant frequency and exponential decay
timescale of $h_{+}^{(2,2)}$ (see top right panel) to be $\omega_{\rm
post}M_{0}\approx 0.23$ and $\tau_{\rm post}/M_{0}\approx 3\times 10^{2}$,
respectively. There is a slight drift of the center of mass that is barely
noticeable in bottom left panel corresponding to $v^{x}_{\rm
com}\approx-1.2\times 10^{-4}$ (all other components are $\lesssim 10^{-5}$;
see Appendix C for a discussion).

We return now to the eccentricity-reduced binary configurations discussed in
Sec. III.3, and provide further details on their orbital evolution as well as
GW emission. In particular, we consider the equal-mass, non-spinning binary
$\mathcal{B}_{5}$, and the aligned-spin super-spinning binary
$\mathcal{B}_{3}$ of mass-ratio $q=1.43$ with details provided in Table 1. In
both cases, we focus mainly on the last eccentricity-reduction iteration step
shown in Fig. 5 (i.e., $N_{e}=5$ and $N_{e}=3$ for the non-spinning and
spinning binaries, respectively) with eccentricity $e=4\times 10^{-3}$ in the
case of $\mathcal{B}_{5}$ and $e=1.8\times 10^{-3}$ in case of the non-
precessing binary $\mathcal{B}_{3}$. In Fig. 8, we show the orbital
trajectories, binary separations and radial velocities, as well as the emitted
GW strain for both configurations.

First, let us focus on the eccentricity reduced binary $\mathcal{B}_{5}$ (top
row of Fig. 8). After initial gauge dynamics (we utilize damped harmonic gauge
and the generalized harmonic formulation, see Appendix A for details), the
binary settles into a state with roughly 17% large coordinate separation of
$d/M_{0}\approx 14$. As can be seen from the top center panel of Fig. 8, the
time derivative $\dot{d}$ of the coordinate separation exhibits large, high-
frequency features throughout the inspiral. These are likely a result of the
residual low-amplitude perturbations of each star; the relative amplitude of
the oscillations in this case is at the $\Delta\phi_{m}\sim 10^{-3}$ level.
Clearly, further eccentricity reduction cannot rely on a fit to $\dot{d}$, but
even a fit to $d$ proves to be challenging at these eccentricities.
Nonetheless, the inspiral dynamics follows a quasi-circular trajectory up to
the point of contact, at a coordinate separation of roughly twice the star’s
radii: $d\approx 2R$. As the compactness of each star is $C=0.13$, we expect
qualitative similarities to the inspiral of a binary neutron star. After the
point of contact at $t/M_{0}=855$, the two non-spinning stars merge into
another non-spinning star. As shown in detail in Ref. Siemonsen and East
(2023), the binary $\mathcal{B}_{5}$ favors the formation of a $m=1$ rotating
remnant BS on purely energetic grounds. However, it lacks sufficient total
angular momentum, and with a phase offset of $\alpha/\pi=1/2$, is not expected
to form a rotating BS remnant, but instead a non-rotating star, with the
residual angular momentum shed in the form of scalar and gravitational
radiation. Hence, after merger, the system rings down towards a single, non-
spinning BS. This is reflected in the GW strain (top right panel of Fig. 8),
which exhibits a near exponential decay in the post-merger phase with a rough
decay timescale of $\tau_{\rm post}/M_{0}=3\times 10^{2}$ and a dominant
ringdown frequency of $\omega_{\rm post}M_{0}=0.23$. It should be noted that,
though not obvious from the figure, the GW strain shown in Fig. 8 does contain
some residual high-frequency GW contamination of the kind discussed in detail
in Sec. III.2. This is made more apparent in Appendix B.

We now turn to the inspiral of the binary $\mathcal{B}_{3}$ shown in the
bottom row of Fig. 8, which exhibits several new features fundamentally
different from either binary black hole or neutrons star coalescences. After
the initial gauge dynamics, the system has a coordinate separation of roughly
$d/M_{0}\approx 14$. Before $t/M_{0}=1475$ (indicated in the figure with an
x-marker), the binary exhibits a smooth inspiral with decreasing separation.
In the bottom center panel of Fig. 8, we compare the coordinate distance as
functions of time to those cases with no eccentricity reduction and no
conformal kinetic energy rescaling power $p$. Between the $N_{e}=0$ cases, the
binary with $p=-4$ conformally scaling for the scalar kinetic energy merges
before the otherwise identical binary with $p=0$. After several iterations of
eccentricity reduction ($N_{e}=3$), the merger occurs later. However, for
$t/M_{0}\gtrsim 1475$, the coordinate separation $d$ exhibits oscillations,
which culminate in the two stars moving away from each other, temporarily
increasing the coordinate separation by approximately $\sim 10$% to
$d/M_{0}>10$. After this sudden repulsion, the stars begin to merge at
$t/M_{0}=1700$ (indicated by a square marker in Fig. 8). In the inset of the
left bottom panel in Fig. 8, we show the in-plane trajectories of both stars
from $t/M_{0}=1475$ to merger. Qualitatively, the sudden repulsion of the two
stars results in a sudden increase of orbital eccentricity for the last $3/2$
orbits before merger (e.g., the orbital trajectories self-intersect). This
repulsive behavior is likely due to strong scalar interactions between the two
stars in the late stages of the inspiral. This behavior was first observed in
Ref. Siemonsen and East (2023) during the late inspiral and merger of two BSs
in a scalar theory with repulsive scalar self-interactions. Despite the
terminology, systems may exhibit repulsive behavior in both attractive and
repulsive scalar models (see e.g. Ref. Palenzuela  _et al._ (2007)). These
scalar interactions are exponentially suppressed by the separation of the
binary BS and, hence, only become active in the late stages of the inspiral,
depending on the constituent stars’ compactnesses and frequencies. In contrast
to the less relativistic binaries considered in Ref. Siemonsen and East
(2023), the constituents of binary $\mathcal{B}_{3}$ are highly compact, with
scalar densities that rapidly decay away from the individual stars. As a
result, the scalar interactions increase in importance over the gravitational
interaction only roughly two orbits before the point of contact. At this
stage, however, the strength of the scalar interactions surpasses that of the
gravity, likely resulting in the rapid increase of the coordinate separation
shown in the bottom center and left panels of Fig. 8. This is a feature absent
in mergers of compact binaries composed of black holes and neutron stars, and
may serve as a smoking gun signature to distinguish BS binaries from such
cases.

To understand the nonlinear merger dynamics of binary $\mathcal{B}_{3}$,
recall that the initial phase offset of this binary is $\alpha/\pi=1/3$ and
that $\omega_{1}\neq\omega_{2}$ (see Table 1). This latter renders the vortex
structure of the binary time-dependent even at the linear level. Hence, a
precise prediction and understanding of the merger outcome using the remnant
map of Ref. Siemonsen and East (2023) is challenging. However, the latter can
still be utilized to qualitatively analyze the merger process. Due to the
time-dependent scalar phase, there is no fixed vortex at the center of mass.
Thus, the vortex structure does not prevent the formation of a single BS, and
hence, the final remnant may be a combination of non-spinning stars including
possibly only a single spherically symmetric BS. Both transitions, two $m=1$
rotating BSs merging to two $m=0$ BSs of the same charge, or a single $m=0$ of
the same charge, are energetically favorable (i.e., the corresponding quantity
analogous to (50) satisfies $\mathcal{M}>0$ in the relevant part of the
parameter space). However, since the spins of the inspiraling super-spinning
BSs are aligned with the orbital angular momentum, the system contains large
amounts of angular momentum. While it is plausible that a single non-spinning
BS may shed this angular momentum sufficiently rapidly during the nonlinear
merger (based on findings of e.g., Ref. Palenzuela  _et al._ (2017)), we find
that the final remnant is instead a binary of non-spinning stars, which are
flung out away from the center of mass at high velocities (with the residual
angular momentum likely being converted into orbital angular momentum as the
stars come into contact). The binary separation continues to increase up until
we terminate the evolution, at which point the separation increased to roughly
$d/M_{0}\approx 40$. The breakup of the spinning binary at the point of
contact occurs at the locations of the vortices of each spinning star. As
such, the nonlinear merger process of $\mathcal{B}_{3}$ is qualitatively
similar to what in shown in Fig. 6 of Ref. Siemonsen and East (2023).

###  IV.3 Precessing binary

Finally, we turn to the precessing binary configuration $\mathcal{B}_{3}$ (see
Table 1) mentioned throughout this work, and analyze its nonlinear dynamics in
detail. The initial data for binary $\mathcal{B}_{3}$ is solved using an
initial coordinate separation $D_{0}=12M_{0}$ and a conformal rescaling
exponent $p=-4$. The spin of each star is chosen to make a 45∘-angle with the
plane containing the initial positions and velocities of the BSs, such that
the component parallel to this plane is in the opposite direction to the
initial tangential velocity. Recall, the dimensionless spins of both
individual stars are above the Kerr-bound $S_{i}/M^{2}_{i}>1$. As indicated in
Fig. 5, we perform a single eccentricity and center-of mass velocity reduction
step, resulting in an eccentricity of $e=0.013$ and linear center of mass
drift of $v_{\rm com}\approx 9\times 10^{-5}$ (with $|v^{z}_{\rm com}|\approx
7\times 10^{-5}$).

In Fig. LABEL:fig:prec_inspirals, we show snapshots of the evolution of these
binary BS initial data, as well as the gravitational waveform. Focusing first
on the inspiral, recall that the spin of a rotating BS points along the vortex
line through the center of the star, i.e., perpendicular to the torus formed
by surfaces of constant scalar field magnitude. As can be seen from the top
panel, the binary exhibits rapid precession of the star’s spin vectors
throughout the early inspiral. In fact, as the individual stars are super-
spinning, the strength of the spin-orbit and spin-spin interactions surpasses
that of any binary black hole during the inspiral. As the eccentricity is
relatively large compared to the other quasi-circular cases, the initial
separation is relatively close, and there is residual high-frequency GW
contamination of the type discussed in Sec. III.2, the typical modulation of
the GW amplitude due to precession cannot be seen in the extracted
gravitational waveform. The merger of this binary is qualitatively similar to
the aligned-spin case discussed in Sec. IV.2: the two rotating BSs collide to
form two non-spinning, highly perturbed BSs, which move outward from the
center of mass. This is shown in the last two panels of the top row of Fig.
LABEL:fig:prec_inspirals. The coordinate separation surpasses $d/M_{0}>40$
before we terminate the evolution, and is not clear whether this new binary is
gravitationally bound.

##  V Conclusion

In this work, we tackled the problem of robustly constructing binary BS
initial data satisfying the Hamiltonian and momentum constraints of the
Einstein equations utilizing the CTS formulation. We analyzed and tested
various approaches to specifying the scalar free data entering these equations
based on superposing isolated boosted BS solutions. Among these approaches, we
found considering the complex scalar field and its conformally rescaled
kinetic energy as free data to robustly lead to constraint satisfying initial
data that could be readily evolved without further reconstruction procedures.
Beyond a simple superposition of BS solutions, we also reduce the spurious
oscillations induced by non-equilibrium initial data using several methods. As
suggested in previous studies, we attenuate the superposed free data of one
compact object in the vicinity of the second compact object. In addition, here
we introduce a new approach where we change the relation between the initial
scalar kinetic energy and the conformally scaled version of this quantity
which is specified as free data. This reduces the scalar kinetic energy, and
hence, results in less perturbed scalar matter. Finally, we successfully
reduce the orbital eccentricity of various mass ratio binary BSs down to the
$e\sim 10^{-3}$ level.

Our procedure for constructing binary BS is highly generic, and thus, is
ideally suited to exploring a vast space of possible binary configurations.
Here we test it for head-on and quasi-circular inspiral cases, including
different mass-ratios, spin magnitudes, and spin orientations.

Ideally, instead of basing the scalar field configuration on superposed BS
solutions, one would solve additional equations imposing a quasi-equilibrium
of the scalar matter with respect to an approximate helical Killing field,
analogous to approaches to construction of equilibriated binary neutron star
initial data. In this work, we briefly discuss several ways in which one might
approach this problem, and some of the complications that may arise, in
particular if one wishes to tackle generic spinning binaries as considered
here. However, we leave an implementation and testing of such an approach to
future work.

While we were able to efficiently reduce spurious oscillations of the BSs in
the binary, particularly for star solutions in the scalar theory with
repulsive self-interactions, the residual perturbations limit the eccentricity
reduction and contaminate the extracted gravitational waveform. In the cases
we consider with eccentricities $\sim 10^{-3}$, we find further reduction of
this quantity to be challenging as the spurious oscillations in each star
induce high-frequency oscillations of the coordinate separation with
comparable or larger amplitude as the eccentricity. Additionally, even small-
amplitude perturbations in the scalar field making up BSs in scalar theories
with solitonic potential induce large-amplitude high-frequency contamination
in the gravitational radiation at early times during the numerical evolutions
of the initial data. Both artifacts may be suppressed by solving the CTS
constraint equations together with quasi-equilibrium scalar matter equations.

With the methods to construct binary BS initial data presented in this work,
accurate waveforms of low-eccentricity non-spinning and (super-)spinning
binary BSs can be obtained. Until the late inspiral, the binary evolution is
largely model-independent, i.e., the inspiral-dynamics is driven by
gravitational effects such as spin-interactions, rather than any mechanism
specific to the scalar matter making up the stars. The resulting waveforms
could be used in to validate and tune inspiral waveform models that would
allow for binary parameters outside the ranges allowed by black holes and
neutron stars (e.g., Ref. LaHaye  _et al._ (2022)). Likewise, current tests to
distinguish binary black holes and neutron stars from exotic alternatives
based on their GW signals Abbott  _et al._ (2021) could be validated with
accurate binary BS inspiral-merger-ringdown waveforms relying on the initial
data constructed in this work. Another interesting avenue for future work is
to study the impact of relativistic features such as stable light rings and
ergoregions of exotic compact objects on the inspiral gravitational waveform
using highly compact BSs.

###### Acknowledgements.

We would like to thank Luis Lehner for valuable discussions. We acknowledge
support from an NSERC Discovery grant. Research at Perimeter Institute is
supported in part by the Government of Canada through the Department of
Innovation, Science and Economic Development Canada and by the Province of
Ontario through the Ministry of Colleges and Universities. This research was
undertaken thanks in part to funding from the Canada First Research Excellence
Fund through the Arthur B. McDonald Canadian Astroparticle Physics Research
Institute. This research was enabled in part by support provided by SciNet
(www.scinethpc.ca), Calcul Québec (www.calculquebec.ca), and the Digital
Research Alliance of Canada (www.alliancecan.ca). Simulations were performed
on the Symmetry cluster at Perimeter Institute, the Niagara cluster at the
University of Toronto, and the Narval cluster at the École de technologie
supérieure in Montreal.

##  Appendix A Numerical setup

In this appendix, we discuss the details of the numerical evolution of the
binary BS initial data, which we perform using the same methods as in Refs.
Siemonsen and East (2021, 2023). We evolve the Einstein-Klein-Gordon system of
equations, derived from (1), using the generalized harmonic formulation of the
Einstein equations Pretorius (2005), including constraint damping terms
Gundlach  _et al._ (2005). The spatial fourth-order accurate discretization is
achieved using finite-difference stencils over a compactified grid containing
spatial infinity. At these boundaries, we impose asymptotically flat boundary
conditions both on the metric and the scalar variables. In order to track the
individual stars of a given binary, we utilize adaptive mesh refinement of the
Cartesian grid with a 2:1 refinement ratio (see Ref. East  _et al._ (2012b)).
The time-stepping is achieved using a fourth-order accurate Runge-Kutta
integration. As briefly mentioned in Sec. IV.1, the axisymmetric evolutions
are performed using a generalized Cartoon method Alcubierre  _et al._ (2001);
Pretorius (2005), which explicitly assumes an azimuthal Killing field
$k^{\mu}$, such that $\mathcal{L}_{k}g_{\mu\nu}=0$ and
$\mathcal{L}_{k}\Phi=im\Phi$. Typically, we require seven mesh refinement
levels with grid spacing $\Delta x/M_{0}=0.08$ on the finest level when
evolving an inspiraling binary BSs; in the case of axisymmetric evolutions the
grid spacing on the finest level is typically $\Delta x/M_{0}=0.01$. For all
evolutions considered here, the gauge is specified by setting the source
functions $H^{\mu}=\square x^{\mu}$ according to the damped harmonic gauge
Lindblom and Szilagyi (2009); Choptuik and Pretorius (2010).

Figure 10:  Convergence study of the numerical evolution of the axisymmetric
binary $\mathcal{B}_{1}$ with properties summarized in Table 1 (and initial
coordinate separation $D=10M_{0}$), at three different numerical resolutions,
where $h$ is the grid spacing of the lowest resolution. The point of contact
of this head-on collision is roughly at $t/M_{0}\approx 16$. The top panel
shows $\phi_{m}=\max|\Phi|/\max|\Phi|_{t=0}$, the normalized maximum of the
scalar field magnitude, while the lower panels show two different measures of
the constraint violation, $\max\mathcal{C}$ and $\mathcal{I}_{\mathcal{C}}$
(defined in the text). The constraint violation is converging to zero at
roughly third and fourth order in the middle and bottom panels, respectively.
Figure 11: The convergence behavior of the two norms, $\max\mathcal{C}$ and
$\mathcal{I}_{\mathcal{C}}$ (defined in the main text), of the constraint
violation of the binary initial data associated with $\mathcal{B}_{5}$ (see
Table 1) with increasing resolution $N$, the number of grid points in each
linear dimension, compared with the lowest resolution considered $N_{\rm min}$
extracted from subsequent numerical evolution at $t/M_{0}=5$. The solid lines
correspond to the constraint satisfying initial data constructed with methods
developed in the main text (with trivial attenuation functions and conformal
exponent $p=-4$). The dashed lines indicate the constraint violation of the
associated free data (i.e., the plain superposition of the two stars). In the
case of the constraint satisfying initial data, both norms exhibit (as
expected) roughly fourth-order accurate convergence towards zero as indicated
by the dash-dotted line. The constraint violation of the superposed initial
data, however, is not converging towards zero, even at very low resolutions
(note, $N/N_{\min}=1$, or equivalently $\Delta x/M_{0}=0.16$, corresponds to
roughly 10 points across each star).  Figure 12: The evolution of the
integrated constraint violation around the binary $\mathcal{B}_{5}$ starting
from plain superposed initial data at three different resolutions $N$ (number
of grid points in each linear dimension) compared to a minimal resolution
$N_{\min}$. In Fig. 11, we showed the non-convergence of the constraint
violation at early times.

We perform convergence tests in order to validate our numerical methods and
quantify the truncation error. Specifically, we measure the constraint
violation given by $\mathcal{C}=M_{0}\sum_{\mu}|(H^{\mu}-\square x^{\mu})|/4$
throughout the evolutions by computing $\max\mathcal{C}$, as well as the
integrated norm $\mathcal{I}_{\mathcal{C}}=M_{0}^{-2}\int
d^{3}x\sqrt{\gamma}\mathcal{C}$ in a coordinate sphere of radius $20M_{0}$. In
the generalized harmonic formulation, any violation of the Hamiltonian and
momentum constraints on the initial time slice will lead $\mathcal{C}$ to
evolve to a non-zero value (as does truncation error). In Fig. 10, we present
the typical convergence behavior of $\mathcal{I}_{\mathcal{C}}$ and
$\max\mathcal{C}$ for an axisymmetric binary BS evolution starting from
constraint satisfying initial data using the above numerical evolution setup.
In Fig. 11, we present a convergence study of the constraint violation of the
quasi-circular binary initial data associated with $\mathcal{B}_{5}$ (see
Table 1). In the axisymmetric cases, the maximum of the constraint violation
$\max\mathcal{C}$ is mainly set by the third-order accurate time-interpolation
performed on the mesh refinement boundaries, and as a result, converges at
roughly third order or better in Fig. 10. In the three-dimensional evolutions,
$\max\mathcal{C}$ is typically reached in the interior of the stars. Hence, in
those cases, we expect a roughly fourth-order convergence behavior, as shown
in Fig. 11. In both settings, however, the integrated norm
$\mathcal{I}_{\mathcal{C}}$ converges at approximately fourth order, as it is
less affected by the lower-order interpolation on the mesh refinement
boundaries.

Also in Fig. 11, we contrast this convergence with the behavior of the
constraint violation with increasing resolution for initial data composed of
the plain superposition of the two stars, i.e. when not solving the constraint
equations. At the lowest resolution, the two are comparable in magnitude.
However, the constraint violation of the free data is not converging to zero
with increasing resolution. Resolutions of $N/N_{\min}<1$ are not suitable for
evolutions, as numerical instabilities associated with unresolved spatial
scales (such as the boson mass $\mu$) appear. Furthermore, we find no evidence
that the constraint damping terms we include in our evolution equations act to
damp away the constraint violations from the superposed initial data down to a
level where they set by truncation error. This is illustrated in Fig. 12.
Note, the damping timescale associated with the constraint damping is
$\tau_{\rm constr.}/M_{0}=0.5$.

##  Appendix B GW contamination

Figure 13: We show the $(\ell,m)=(2,2)$ spin-weighted spherical harmonic
components of the Newman-Penrose scalar $\Psi_{4}$ extracted on a coordinate
sphere of radius $100M_{0}$. The top panel corresponds to the $N_{e}=5$ and
$p=-4$ binary $\mathcal{B}_{5}$ shown in the top row of Fig. 8, while the
bottom panel shows the $N_{e}=3$ and $p=-4$ binary $\mathcal{B}_{3}$ shown in
the bottom row of Fig. 8. This shows the high-frequency contamination of the
gravitational waveform from the binaries at early times due to residual
spurious oscillations and unbound scalar matter in and around the constituents
of the binaries.

In this appendix, we briefly return to the high-frequency contamination in the
GW emission from the binaries presented in Sec. IV.2. This contamination
emerges from the residual perturbations present in the eccentricity reduced
binaries $\mathcal{B}_{3}$ and $\mathcal{B}_{5}$ constructed with conformally
rescaled kinetic energy using (43) with $p=-4$. In Fig. 13, we present the
$(\ell,m)=(2,2)$ spherical harmonic component of the Newman-Penrose scalar
$\Psi_{4}$, which can be compared to the corresponding GW strain shown in the
right panels of Fig. 8. Since $\Psi_{4}$ is related to the strain by two time
derivatives, it accentuates the high-frequency component from the perturbed
BSs. As shown in Sec. III.2, this contamination is reduced by means of the
rescaling (43) of the conformal kinetic energy. However, with the lowest
exponent with which the elliptic solver was able to find a solution, $p=-4$,
some oscillations remain in the stars, and lead to the high frequency
component to $\Psi_{4}$ evident at early times in Fig. 13.

##  Appendix C Center-of-mass motion

Within our approach, the initial linear momentum is set to zero using
Newtonian expressions for the initial boost velocities (as discussed in Sec.
III.3). We find this to be sufficient for non-spinning binaries, i.e, the
center-of-mass velocity throughout the evolution of the initial data remains
below $v_{\rm com}<10^{-7}$. For highly-spinning quasi-circular binary initial
data, the center of mass of the system exhibits larger drifts with constant
velocity away from origin of the numerical grid. For the aligned-spin binary
$\mathcal{B}_{3}$, the magnitude of the in-orbital-plane coordinate velocity
of the center of mass is $v_{\rm com}=6\times 10^{-3}$ (the out-of-plane
component is $<10^{-10}$). Leading post-Newtonian corrections to the center-
of-mass and center-of-momentum velocities, used to initialize the binaries in
this work, are roughly an order of magnitude too small to account for $v_{\rm
com}$. This strong drift may be the result of spurious gravitational and
scalar radiation emitted during the first few light crossing times of the
binary, as well as the large spins of the super-spinning binary
$\mathcal{B}_{3}$. We address this by measuring the in-plane components
$v_{\rm com}^{i}$ of the center-of-mass coordinate velocity and subtracting
this from the initial binary velocities $v_{(A)}^{i}$ of the free data as
defined in Sec. II.2. This is done in tandem with the eccentricity reduction.
Hence, the linear momentum of the binary can be iteratively reduced in this
way. After two iteration steps, the velocity is reduced by more than an order
of magnitude to $v_{\rm com}=1\times 10^{-4}$ for the binary $\mathcal{B}_{3}$
with aligned spins. In the case of the precessing binary $\mathcal{B}_{3}$,
discussed in Sec. IV.3, the center-of-mass motion in the direction of the
orbital angular momentum dominates; this we treat iteratively in precisely the
same manner as the in-plane center-of-mass velocity.

## References

  * Cook (2000) G. B. Cook, Living Rev. Rel. 3, 5 (2000), arXiv:gr-qc/0007085 .
  * Pfeiffer (2005) H. P. Pfeiffer, J. Hyperbol. Diff. Equat. 2, 497 (2005), arXiv:gr-qc/0412002 .
  * Gourgoulhon (2007) E. Gourgoulhon, (2007), arXiv:gr-qc/0703035 .
  * Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro,  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge University Press, 2010).
  * Tichy (2017) W. Tichy, Rept. Prog. Phys. 80, 026901 (2017), arXiv:1610.03805 [gr-qc] .
  * Pan  _et al._ (2014) Y. Pan, A. Buonanno, A. Taracchini, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 89, 084006 (2014), arXiv:1307.6232 [gr-qc] .
  * Hannam  _et al._ (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [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).
  * Yoshida and Eriguchi (1997a) S. Yoshida and Y. Eriguchi, Phys. Rev. D 56, 762 (1997a).
  * Seidel and Suen (1990) E. Seidel and W.-M. Suen, Phys. Rev. D 42, 384 (1990).
  * Seidel and Suen (1991) E. Seidel and W. M. Suen, Phys. Rev. Lett. 66, 1659 (1991).
  * Brito  _et al._ (2016) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Phys. Lett. B 752, 291 (2016), arXiv:1508.05395 [gr-qc] .
  * Friedberg  _et al._ (1987) R. Friedberg, T. Lee, and Y. Pang, Phys. Rev. D 35, 3658 (1987).
  * Bošković and Barausse (2022) M. Bošković and E. Barausse, JCAP 02, 032 (2022), arXiv:2111.03870 [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] .
  * Kleihaus  _et al._ (2005) B. Kleihaus, J. Kunz, and M. List, Phys. Rev. D 72, 064002 (2005), arXiv:gr-qc/0505143 .
  * Kleihaus  _et al._ (2008) B. Kleihaus, J. Kunz, M. List, and I. Schaffer, Phys. Rev. D 77, 064025 (2008), arXiv:0712.3742 [gr-qc] .
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Visinelli (2021) L. Visinelli, Int. J. Mod. Phys. D 30, 2130006 (2021), arXiv:2109.05481 [gr-qc] .
  * Liebling and Palenzuela (2017) S. L. Liebling and C. Palenzuela, Living Rev. Rel. 20, 5 (2017), arXiv:1202.5809 [gr-qc] .
  * 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] .
  * Sennett  _et al._ (2017) N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno, and S. Ossokine, Phys. Rev. D 96, 024002 (2017), arXiv:1704.08651 [gr-qc] .
  * Krishnendu  _et al._ (2017) N. V. Krishnendu, K. G. Arun, and C. K. Mishra, Phys. Rev. Lett. 119, 091101 (2017), arXiv:1701.06318 [gr-qc] .
  * Johnson-Mcdaniel  _et al._ (2020) N. K. Johnson-Mcdaniel, A. Mukherjee, R. Kashyap, P. Ajith, W. Del Pozzo, and S. Vitale, Phys. Rev. D 102, 123010 (2020), arXiv:1804.08026 [gr-qc] .
  * Vaglio  _et al._ (2022) M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, Phys. Rev. D 105, 124020 (2022), arXiv:2203.07442 [gr-qc] .
  * Adam  _et al._ (2022) C. Adam, J. Castelo, A. García Martín-Caro, M. Huidobro, R. Vázquez, and A. Wereszczynski, Phys. Rev. D 106, 123022 (2022), arXiv:2203.16558 [gr-qc] .
  * Pacilio  _et al._ (2020) C. Pacilio, M. Vaglio, A. Maselli, and P. Pani, Phys. Rev. D 102, 083002 (2020), arXiv:2007.05264 [gr-qc] .
  * Vaglio  _et al._ (2023) M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, (2023), arXiv:2302.13954 [gr-qc] .
  * Lai (2004) C.-W. Lai,  _A Numerical study of boson stars_ , Other thesis (2004), arXiv:gr-qc/0410040 .
  * Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [gr-qc] .
  * Paredes and Michinel (2016) A. Paredes and H. Michinel, Phys. Dark Univ. 12, 50 (2016), arXiv:1512.05121 [astro-ph.CO] .
  * Bernal and Siddhartha Guzman (2006) A. Bernal and F. Siddhartha Guzman, Phys. Rev. D 74, 103002 (2006), arXiv:astro-ph/0610682 .
  * Schwabe  _et al._ (2016) B. Schwabe, J. C. Niemeyer, and J. F. Engels, Phys. Rev. D 94, 043513 (2016), arXiv:1606.05151 [astro-ph.CO] .
  * Palenzuela  _et al._ (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D 75, 064005 (2007), arXiv:gr-qc/0612067 .
  * Mundim (2010) B. C. Mundim,  _A Numerical Study of Boson Star Binaries_ , Ph.D. thesis, British Columbia U. (2010), arXiv:1003.0239 [gr-qc] .
  * Bezares  _et al._ (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [gr-qc] .
  * Sanchis-Gual  _et al._ (2019a) N. Sanchis-Gual, C. Herdeiro, J. A. Font, E. Radu, and F. Di Giovanni, Phys. Rev. D 99, 024017 (2019a), arXiv:1806.07779 [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] .
  * Sanchis-Gual  _et al._ (2022) N. Sanchis-Gual, J. Calderón Bustillo, C. Herdeiro, E. Radu, J. A. Font, S. H. W. Leong, and A. Torres-Forné, (2022), arXiv:2208.11717 [gr-qc] .
  * Evstafyeva  _et al._ (2022) T. Evstafyeva, U. Sperhake, T. Helfer, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, (2022), arXiv:2212.08023 [gr-qc] .
  * Palenzuela  _et al._ (2008) C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. D 77, 044036 (2008), arXiv:0706.2435 [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] .
  * Bezares and Palenzuela (2018) M. Bezares and C. Palenzuela, Class. Quant. Grav. 35, 234002 (2018), arXiv:1808.10732 [gr-qc] .
  * Siemonsen and East (2023) N. Siemonsen and W. E. East, (2023), arXiv:2302.06627 [gr-qc] .
  * York (1999) J. W. York, Jr., Phys. Rev. Lett. 82, 1350 (1999), arXiv:gr-qc/9810051 .
  * East  _et al._ (2012a) W. E. East, F. M. Ramazanoglu, and F. Pretorius, Phys. Rev. D 86, 104053 (2012a), arXiv:1208.3473 [gr-qc] .
  * Siemonsen and East (2021) N. Siemonsen and W. E. East, Phys. Rev. D 103, 044022 (2021), arXiv:2011.08247 [gr-qc] .
  * Sanchis-Gual  _et al._ (2019b) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. Font, and E. Radu, Phys. Rev. Lett. 123, 221101 (2019b), arXiv:1907.12565 [gr-qc] .
  * Aurrekoetxea  _et al._ (2022) J. C. Aurrekoetxea, K. Clough, and E. A. Lim, (2022), arXiv:2207.03125 [gr-qc] .
  * Corman and East (2022) M. Corman and W. E. East, (2022), arXiv:2212.04479 [gr-qc] .
  * Bonazzola  _et al._ (1997) S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. Rev. D 56, 7740 (1997), arXiv:gr-qc/9710031 .
  * Shibata (1998) M. Shibata, Phys. Rev. D 58, 024012 (1998), arXiv:gr-qc/9803085 .
  * Teukolsky (1998) S. A. Teukolsky, Astrophys. J. 504, 442 (1998), arXiv:gr-qc/9803082 .
  * Gourgoulhon  _et al._ (2002) E. Gourgoulhon, P. Grandclement, and S. Bonazzola, Phys. Rev. D 65, 044020 (2002), arXiv:gr-qc/0106015 .
  * Hui (2021) L. Hui, Ann. Rev. Astron. Astrophys. 59, 247 (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Marronetti and Matzner (2000) P. Marronetti and R. A. Matzner, Phys. Rev. Lett. 85, 5500 (2000), arXiv:gr-qc/0009044 .
  * Marronetti  _et al._ (2000) P. Marronetti, M. Huq, P. Laguna, L. Lehner, R. A. Matzner, and D. Shoemaker, Phys. Rev. D 62, 024017 (2000), arXiv:gr-qc/0001077 .
  * Bonning  _et al._ (2003) E. Bonning, P. Marronetti, D. Neilsen, and R. Matzner, Phys. Rev. D 68, 044019 (2003), arXiv:gr-qc/0305071 .
  * Helfer  _et al._ (2019) T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, Phys. Rev. D 99, 044046 (2019), arXiv:1802.06733 [gr-qc] .
  * Pfeiffer  _et al._ (2007) H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom, G. Lovelace, and M. A. Scheel, Class. Quant. Grav. 24, S59 (2007), arXiv:gr-qc/0702106 .
  * Boyle  _et al._ (2007) M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, Phys. Rev. D 76, 124038 (2007), arXiv:0710.0158 [gr-qc] .
  * Mroue  _et al._ (2010) A. H. Mroue, H. P. Pfeiffer, L. E. Kidder, and S. A. Teukolsky, Phys. Rev. D 82, 124016 (2010), arXiv:1004.4697 [gr-qc] .
  * Buonanno  _et al._ (2011) A. Buonanno, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, and A. Taracchini, Phys. Rev. D 83, 104034 (2011), arXiv:1012.1549 [gr-qc] .
  * Nikolaieva  _et al._ (2022) Y. O. Nikolaieva, Y. M. Bidasyuk, K. Korshynska, E. V. Gorbar, J. Jia, and A. I. Yakimenko, (2022), arXiv:2212.06116 [nlin.PS] .
  * Gervalle (2022) R. Gervalle, Phys. Rev. D 105, 124052 (2022), arXiv:2206.03982 [gr-qc] .
  * Yoshida and Eriguchi (1997b) S. Yoshida and Y. Eriguchi, Phys. Rev. D 55, 1994 (1997b).
  * LaHaye  _et al._ (2022) M. LaHaye, H. Yang, B. Bonga, and Z. Lyu, (2022), arXiv:2212.04657 [gr-qc] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2112.06861 [gr-qc] .
  * Pretorius (2005) F. Pretorius, Class. Quant. Grav. 22, 425 (2005), arXiv:gr-qc/0407110 [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 .
  * East  _et al._ (2012b) W. E. East, F. Pretorius, and B. C. Stephens, Phys. Rev. D 85, 124010 (2012b), arXiv:1112.3094 [gr-qc] .
  * 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 .
  * Lindblom and Szilagyi (2009) L. Lindblom and B. Szilagyi, Phys. Rev. D80, 084019 (2009), arXiv:0904.4873 [gr-qc] .