# General Relativistic Polarized Proca Stars

Zipeng Wang<EMAIL_ADDRESS>Department of Physics and Astronomy, Johns
Hopkins University, Baltimore, MD 21218, USA  Thomas Helfer
<EMAIL_ADDRESS>Institute for Advanced Computational Science, Stony Brook
University, Stony Brook, NY 11794 USA  Mustafa A. Amin
<EMAIL_ADDRESS>Department of Physics and Astronomy, Rice University,
Houston, Texas 77005, USA

###### Abstract

Massive vector fields can form spatially localized, non-relativistic,
stationary field configurations supported by gravitational interactions. The
ground state configurations (p-solitons/vector solitons/dark photon
stars/polarized Proca stars) have a time-dependent vector field pointing in
the same spatial direction throughout the configuration at any instant of
time, can carry macroscopic amounts of spin angular momentum, and are
spherically symmetric and monotonic in the energy density. In this paper, we
include general relativistic effects, and numerically investigate the
stability of compact polarized Proca stars (linear and circularly polarized)
and compare them to hedgehog-like field configurations (with radially pointing
field directions). Starting with approximate field profiles of such stars, we
evolve the system numerically using 3+1 dimensional numerical simulations in
general relativity. We find that these initial conditions lead to stable
configurations. However, at sufficiently large initial compactness, they can
collapse to black holes. We find that the initial compactness that leads to
black hole formation is higher for circularly polarized stars (which carry
macroscopic spin angular momentum), compared to linearly polarized ones, which
in turn is higher than that for hedgehog configurations.

##  I Introduction

Massive vector fields (or dark photons) can constitute all or part of dark
matter. If their mass is $\lesssim 10\,\rm{\rm eV}$, the occupation numbers in
typical astrophysical/cosmological settings are large enough to admit a
classical field description. Their early universe production (for example,
Graham  _et al._ (2016); Agrawal  _et al._ (2020); Co  _et al._ (2019); Dror
_et al._ (2019); Bastero-Gil  _et al._ (2019); Long and Wang (2019); Adshead
_et al._ (2023); Kitajima and Takahashi (2023)), astrophysical/cosmological
phenomenology, as well as direct and indirect detection strategies are being
extensively explored (see Antypas  _et al._ (2022) for a recent review).
Numerical simulations investigating the nonlinear (and non-relativistic)
gravitational dynamics of such fields in an astrophysical setting have been
initiated recently Amin  _et al._ (2022); Gorghetto  _et al._ (2022); Jain
_et al._ (2023); Chen  _et al._ (2023).

Similar to the case of scalar fields, in the non-relativistic limit, one
expects massive vector fields to form spatially localized, non-relativistic,
stationary field configurations (solitons or Boson stars) supported by
gravitational interactions. At any instant of time, such polarized Proca stars
(also referred to as $p$-solitons, vector solitons, dark photon stars) have a
spatially constant orientation of the field polarization throughout the
configuration Adshead and Lozanov (2021); Jain and Amin (2022). Depending on
the polarization, they can carry macroscopic amounts of spin angular momentum
Jain and Amin (2022); Amin  _et al._ (2022). They are spherically symmetric in
energy density, but not in the field configuration (but are node-free). Non-
relativistic (fractionally) polarized solitons have been shown to from
generically from cosmological, as well as astrophysical initial conditions
Gorghetto  _et al._ (2022); Amin  _et al._ (2022); Jain  _et al._ (2023); Chen
_et al._ (2023).

General relativistic effects become necessary to consider if the vector field
configurations become sufficiently compact.111Current simulations show that
solitons forming from generic initial conditions are indeed fractionally
polarized (with macroscopic spin) Jain  _et al._ (2023), and are typically too
diffuse to warrant studying relativistic corrections. However, as such
solitons accrete fields from their surrounding, they can become increasingly
more compact. Such an analysis is relevant for understanding the detailed
nature of the compact configurations, including their maximal compactness,
intrinsic spin, stability and deformability. These properties can be critical
when considering (mergers of) such compact objects as gravitational wave
sources Sanchis-Gual  _et al._ (2022a); Calderón Bustillo  _et al._ (2021).
222The merger of scalar boson stars and their gravitational wave emission has
been studied extensively in the literature, see for example Palenzuela  _et
al._ (2007, 2008); Choptuik and Pretorius (2010); Bezares  _et al._ (2017);
Palenzuela  _et al._ (2017); Helfer  _et al._ (2019); Widdicombe  _et al._
(2020); Sanchis-Gual  _et al._ (2020); Diamond  _et al._ (2023); Evstafyeva
_et al._ (2023); Siemonsen and East (2023); Jaramillo  _et al._ (2022);
Bezares  _et al._ (2022). In this paper, we study such polarized Proca stars
within full numerical relativity.

We note that polarized Proca stars that we focus on here are unlike hedgehog-
like field configurations. Hedgehog-like configurations have spherically
symmetric field configurations with spatially varying field polarization, and
have a node in their field profiles at the origin. Hedgehog configurations
have been studied in detail including general relativistic corrections,
assisted by the spherical symmetry of the energy density and their field
configuration Brito  _et al._ (2016); Sanchis-Gual  _et al._ (2017). However,
such hedgehogs are higher energy states of the field for a given mass compared
to the polarized Proca stars/p-solitons mentioned above Adshead and Lozanov
(2021); Jain and Amin (2022). Unlike the polarized Proca stars, hedgegogs have
been shown to form only under a special set of spherically symmetric (field)
initial conditions and evolution Di Giovanni  _et al._ (2018).

As far as general relativistic polarized Proca stars are concerned, the $m=1$
case was studied in Brito  _et al._ (2016); Sanchis-Gual  _et al._ (2022a,
2019a) (for the complex vector field).333We thank William East and Nils
Siemonsen for guiding us to relevant existing work here. This case is similar
to circularly polarized solitons/Proc stars with macroscopic angular momentum
in Jain and Amin (2022). The analysis in Jain and Amin (2022) (which includes
circular, linear and fractionally polarized cases) is a non-relativistic
analysis, where the underlying relativistic vector field can be real valued.
To the best of our knowledge, general relativistic linearly polarized Proca
stars (with negligible angular momentum) have not been studied in the
literature before.

In this study, we explore the behavior of complex-valued polarized Proca stars
as their compactness is increased, which allows us to estimate a rough lower
bound of the maximum compactness these stars. We numerically investigate the
stability of both the linear and the circular polarization states, and compare
them to hedgehog-like field configurations when general relativistic effects
are included.444The maximal compactness for the hedgehog configurations, and
the circularly polarized ones was already provided in Brito  _et al._ (2016).
The maximal compactness for linearly polarized one, however, has not been
provided.

We work in units where $\hbar=c=G_{N}=1$. In the captions we use the Planck
mass $M_{\rm pl}\equiv 1/\sqrt{G_{N}}$, and occasionally even $\hbar$, to make
the units explicit for clarity. We will use Greek letters ($\mu$,
$\nu,\ldots$) to represent four-dimensional indices, and Latin letters
($i$,$j,\ldots$) to represent three-dimensional spatial indices. We work with
the $-+++$ convention for the metric.

The rest of the paper is organized as follows. In section II, we provide the
underlying model for complex Proca fields in general relativity. In this
section we also describe the numerical relativity framework, construction of
initial data and provide some details of the numerical set-up. In section III,
we present and discuss the results of our simulations. We summarize our
results and provide a taste of their implications in section IV. Convergence
tests and level of constrain violations are discussed in an appendix.

Figure 1: Simulations of the three types of Proca stars. Top row: a hedgehog
configuration generated with initial compactness $\mathcal{C}\approx 0.04$
($\mu=0.04m$) shown on the $x$-$y$ plane at three different times. Middle row:
a linearly polarized Proca star with $\mathcal{C}\approx 0.06$ ($\mu=0.06m$)
shown on the $x$-$z$ plane. The Proca field vectors are polarized in the
$z$-direction. Bottom row: a circularly polarized star with initial
compactness $\mathcal{C}\approx 0.08$ ($\mu=0.10m$) shown on the $x$-$y$
plane. The energy density profiles, $\rho\times(M_{\rm pl}m)^{-2}$, are shown
as color plots. The black arrows show the direction and relative magnitude of
the real part of the spatial Proca vector field $\text{Re}(A_{i})$. The black
bars on top of each panel show the length scale of the plots in units of
$m^{-1}$. The time shown is in units of $m^{-1}$. The (real part of) the
vector field oscillates along the arrows in the top and middle panel, whereas
it rotates in the bottom one (with period $T=2\pi m^{-1}$). Note that the time
interval between snapshots is much longer than $T$; the changes in density
profiles due to perturbations happen on these longer timescales.

##  II Setup and numerical methods

###  II.1 Proca field in general relativity

We consider a vector field $X^{\alpha}$ in general relativity, with an action
Brito  _et al._ (2016)

$\mathcal{S}\\!=\\!\int
d^{4}x\sqrt{-g}\left(\frac{R}{16\pi}-\frac{1}{2}m^{2}X_{\alpha}\bar{X}^{\alpha}-\frac{1}{4}F^{\mu\nu}\bar{F}_{\mu\nu}\right).$
(1)

Here, $R$ is the Ricci scalar, $g=\text{det}(g_{\mu\nu})$, $m$ is the mass of
the vector field, $\bar{X}^{\alpha}$ is the complex conjugate of the vector
field ${X}^{\alpha}$, $F_{\mu\nu}=\partial_{\mu}X_{\nu}-\partial_{\mu}X_{\nu}$
is the field strength tensor, and $\bar{F}_{\mu\nu}$ is its complex conjugate.

Extremizing the action Eq. (1) with respect to variations in $g_{\mu\nu}$
yields the Einstein’s field equation

$G_{\mu\nu}=8\pi T_{\mu\nu}\,.$ (2)

Here, $G_{\mu\nu}$ is the Einstein tensor, and

$T_{\mu\nu}=\frac{1}{2}\left(F_{\mu\rho}\bar{F}_{\nu}{}^{\rho}+\bar{F}_{\mu\rho}F_{\nu}{}^{\rho}\right)-\frac{1}{4}F_{\rho\gamma}\bar{F}^{\rho\gamma}g_{\mu\nu}\\\
+\frac{m^{2}}{2}\left(X_{\mu}\bar{X}_{\nu}+\bar{X}_{\mu}X_{\nu}-g_{\mu\nu}\bar{X}^{\rho}X_{\rho}\right)\,$
(3)

is the stress-energy tensor associated with the Proca field. Similarly,
extremizing the action Eq. (1) with respect to $X_{\mu}$ leads to the Proca
field equation

$\nabla_{\mu}F^{\mu\nu}=m^{2}X^{\nu}\,.$ (4)

The above equation with $m\neq 0$, along with the antisymmetry of
$F^{\mu\nu}$, implies that the field $X^{\nu}$ must satisfy the Proca
constraint equation

$\nabla_{\nu}X^{\nu}=0\,.$ (5)

Eq. (2) and Eq. (4) govern the evolution of the Proca field and spacetime.

###  II.2 3+1 decomposition

Following Arnowitt  _et al._ (1959), we foliate spacetime with spatial slices
$\Sigma$ (with metric $\gamma_{ij}$), and connect the slices with each other
with a lapse function $\alpha$ and shift vector field $\beta^{i}$. The
spacetime metric can then be written as

$ds^{2}=-(\alpha^{2}-\beta^{i}\beta_{i})dt^{2}+2\beta_{i}dx^{i}dt+\gamma_{ij}dx^{i}dx^{j}\,,$
(6)

where $\beta_{i}\equiv\gamma_{ij}\beta^{j}$.

The unit future-directed normal vector of $\Sigma$:
$n_{\mu}=\left(-\alpha,0,0,0\right)$, and
$P_{\mu}{}^{\nu}\equiv\delta_{\mu}{}^{\nu}+n_{\mu}n^{\nu}$ is the projection
tensor that projects on to $\Sigma$ Croft (2022). The extrinsice curvature of
$\Sigma$ is $K_{ij}=-\mathcal{L}_{n}\gamma_{ij}/2$ with $\mathcal{L}_{n}$
being the Lie-derivative along the normal vector $n_{\mu}$. A decomposition of
the stress-energy tensor, Eq. (3), adapted to this foliation of spacetime is

$\rho\equiv
n_{\alpha}n_{\beta}T^{\alpha\beta},\quad\\!\\!\\!\\!S_{i}\equiv-\gamma_{i\alpha}n_{\beta}T^{\alpha\beta},\quad\\!\\!\\!\\!S_{ij}\equiv\\!\gamma_{i\alpha}\gamma_{j\beta}T^{\alpha\beta}.$
(7)

The evolution equation for $\gamma_{ij}$, $K_{ij}$, $\alpha$, and $\beta^{i}$
can be obtained from the Einstein equation (2); see Radia  _et al._ (2022) for
their explicit form.

We also decompose the vector field as Zilhão  _et al._ (2015):

$X_{\mu}=A_{\mu}+n_{\mu}\varphi\,,$ (8)

where $\varphi=-n^{\mu}X_{\mu}$ is the component of the vector field normal to
the spatial slice Zilhão  _et al._ (2015), and
$A_{\mu}=P_{\mu}{}^{\nu}X_{\nu}$ are the components along spatial slice. The
electric field associated with the Proca field is defined as:

$E_{i}\equiv P_{i}{}^{\mu}n^{\nu}F_{\mu\nu}\,.$ (9)

Under this decomposition, the Proca constraint, Eq. (5), becomes

$\mathcal{C_{E}}=D_{i}E^{i}-m^{2}\varphi=0\,.$ (10)

where $D_{i}$ is the covariant derivative corresponding to $\gamma_{ij}$. The
Proca evolution equation ,(4), yields Zilhão  _et al._ (2015)

$\displaystyle\partial_{t}\varphi$
$\displaystyle=-A^{i}D_{i}\alpha+\alpha(K\varphi-
D_{i}A^{i}-Z)+\mathcal{L}_{\beta}\varphi\,,$ $\displaystyle\partial_{t}A_{i}$
$\displaystyle=-\alpha(E_{i}+D_{i}\varphi)-\varphi
D_{i}\alpha+\mathcal{L}_{\beta}A_{i}\,,$ $\displaystyle\partial_{t}E^{i}$
$\displaystyle=\alpha(KE^{i}+D^{i}Z+m^{2}A^{i}+{D^{k}D^{i}A_{k}-D^{k}D_{k}A^{i}})$
$\displaystyle+D^{j}\alpha(D^{i}A_{j}-D_{j}A^{i})+\mathcal{L}_{\beta}E^{i}\,,$
(11) $\displaystyle\partial_{t}Z$
$\displaystyle=\alpha(D_{i}E^{i}+m^{2}\varphi-\kappa
Z)+\mathcal{L}_{\beta}Z\,.$

In the above equation, $\mathcal{L}_{\beta}$ is the Lie derivative with
respect to the shift vector, $K=\gamma^{ij}K_{ij}$ is the mean curvature, and
$Z$ is an auxiliary variable introduced to keep $\mathcal{C_{E}}$ minimized
during evolution Zilhão  _et al._ (2015); Hilditch (2013); Palenzuela  _et
al._ (2010).

###  II.3 Initial data

In this section, we describe the construction of the initial data, which
approximates the full-GR solutions for three types of Proca stars.

####  II.3.1 Proca field

As initial data, we use the field profiles for stationary Proca stars in the
non-relativistic regime, $|\partial_{i}/m|\ll 1$. For details of this
construction, see Jain and Amin (2022). The three types of Proca stars under
consideration have a spatial vector field

${\bm{A}}(\tilde{t},\tilde{\bm{r}})=e^{i\tilde{t}}\frac{\mu}{m}\begin{cases}f^{\rm
lin}(\tilde{r})\hat{\bm{z}}&\textrm{linearly polarized},\\\ f^{\rm
cir}(\tilde{r})\dfrac{\hat{\bm{x}}+i\hat{\bm{y}}}{\sqrt{2}}&\textrm{circularly
polarized},\\\ f^{\rm
hh}(\tilde{r})\hat{\bm{r}}&\textrm{hedgehog}.\end{cases}$ (12)

Here, $\mu$ is the effective chemical potential with
$\mu/m\sim|\partial_{i}^{2}/m^{2}|\ll 1$ in the non-relativistic limit, and
the rescaled coordinates $\tilde{r}\equiv\sqrt{m\mu}r$,
$\tilde{t}\equiv(1-\mu/m)mt$. The profiles $f^{\rm lin}$, $f^{\rm cir}$ and
$f^{\rm hh}$ are approximately given by

$\displaystyle f^{\rm lin}(\tilde{r})$ $\displaystyle=f^{\rm
circ}(\tilde{r})\approx\frac{1.94}{(1+0.073\tilde{r}^{2})^{4}},$ (13)
$\displaystyle f^{\rm hh}(\tilde{r})$
$\displaystyle\approx\frac{0.76\tilde{r}}{(1+0.0096\tilde{r}^{2})^{16}}.$ (14)

These are fitting formulae for the profiles. More accurate profiles can be
obtained by numerically solving the corresponding profile equations Jain and
Amin (2022). Our fits deviate from these numerically obtained profiles by
$\sim 5\%$.

To specify the electrical field (9)

$E_{i}=\gamma^{-\frac{1}{3}}\left(\partial_{i}\varphi-\partial_{0}A_{i}\right)\,,$
(15)

where $\gamma\equiv\text{det}(\gamma_{ij})$. On the initial slice, we use Eq.
(12) to obtain $\partial_{0}A_{i}$. We ignore the $\partial_{i}\varphi$ since
it is suppressed by $|\partial_{i}|/m$. This $E^{i}=\gamma^{ij}E_{j}$ can then
be used in Eq. (10) to obtain $\varphi$.555We need not have ignored $\varphi$,
and could have solved for it using Eq. (10), together with the Hamiltonian and
momentum constraint equations. We found that this procedure was not
numerically stable.

####  II.3.2 Spacetime

On the initial slice, the gauge functions are assumed to be trivial,
$\alpha=1$ and $\beta^{i}=0$. The spatial metric is assumed have conformal
flatness: $\gamma_{ij}=\psi^{4}\delta_{ij}$, where
$\psi=[\text{det}(\gamma_{ij})]^{\frac{1}{12}}$ is the conformal factor of the
metric.

The fields from the previous subsection and the metric, must satisfy the
Hamiltonian and momentum constraint

$\displaystyle\mathcal{H}\equiv R-K_{ij}K^{ij}+K^{2}-16\pi\rho=0\,,$ (16)
$\displaystyle\mathcal{M}^{i}\equiv D_{j}\left(K^{ij}-\gamma^{ij}K\right)-8\pi
S^{i}=0\,,$ (17)

where $R$ is the three-dimensional Ricci scalar of $\gamma_{ij}$.

To solve the constraint equations, we follow the conformal-transverse-
traceless (CTT) formalism (see Baumgarte and Shapiro (2010) chap. 3 and
appendix B). We decompose the extrinsic curvature as
$K_{ij}=\psi^{-2}\bar{A}_{ij}+\gamma_{ij}K/3$, where $\bar{A}_{ij}$ is the
trace-free part of the extrinsic curvature. We further assume zero mean
curvature, $K=0$, on the initial slice. We decompose Eqns. (16) and (17) as

$\displaystyle\Delta\psi+\frac{1}{8}\psi^{-7}\bar{A}_{ij}\bar{A}^{ij}=-2\pi\psi^{5}\rho\,,$
(18) $\displaystyle(\Delta_{L}W)^{i}=8\pi\psi^{10}S^{i}\,,$ (19)

where $\Delta\psi=\partial^{k}\partial_{k}\psi$ is the flat Laplacian of the
conformal factor, $W^{i}$ is the vector potential of $\bar{A}_{ij}$, and
$(\Delta_{L}W)^{i}=\partial^{j}\partial_{j}W^{i}+\frac{1}{3}\partial^{i}\partial_{j}W^{j}$
is the flat vector Laplacian of $W^{i}$.

We utilize a GRChombo-based elliptical solver that has adaptive mesh
refinement support to solve Eqns. (18) and (19) with $\rho$ and $S_{i}$ given
by Eq. (7) based on the vector field profiles. Additionally, the solver
requires an initial guess for the conformal factor $\psi$, which was set to be
the conformal factor of the full-GR stationary hedgehog stars in Brito  _et
al._ (2016). This solver improves this initial guess interatively, updating
$\psi$ and $W^{i}$ each time to reduce the Hamiltonian and momentum
constraints. At each iterative step, we update the $E^{i}$ components of the
field according to Eq. (15) (without $\partial_{i}\varphi$), and then update
the $\varphi$ component of the Proca field according to Eq. (10), to ensure
that the Proca constraint is still satisfied under the updated conformal
factor $\psi$. Note that since the Proca field distribution is compact, we use
the boundary conditions $\psi=1$ and $W^{i}=0$ and put the boundaries far from
the Proca star.

For the densest simulated Proca star generated with $\mu=0.10m$, we solve the
equations with side length $L=300m^{-1}$ with the number of points $N=96$ on
the coarsest level. We add three additional refinement levels enclosing the
star, with the finest resolution as $\Delta=0.2m^{-1}$. Under these
conditions, the procedure detailed above provides good convergence rates, with
$\mathcal{H}$ sufficiently small. In the Appendix, we demonstrate convergence
for $\mathcal{H}$ with different resolutions as an example. The momentum
constraint $\mathcal{M}^{i}$ and the Proca constraint $\mathcal{C_{E}}$
converge in a similar fashion.666However, we observed that, for denser initial
vector profiles with parameter $\mu>0.10m$, this procedure fails to converge.

For convenience, we call the Proca profiles detailed in II.3.1, along with the
spacetime metric solved with Eq. (18) and Eq. (19), “constructed Proca stars”.

A few comments are in order regarding our constructed Proca stars. We expect
our constructed Proca stars to be non-stationary. This is because the non-
relativistic profiles we use will deviate from the true relativistic solutions
as the compactness is increased. This is currently unavoidable for us because
unlike hedgehog-like Proca stars Brito  _et al._ (2016), constructing a
stationary solution at high compactness for polarized stars is difficult due
of the lack of spherical symmetry in the field configuration (and a likely
deviation from spherical symmetry in the energy density). Furthermore, in our
procedure, we used $K=0$ and conformal flatness, ignored the transverse
traceless part of the $K_{ij}$ and chose trivial functions for $\alpha$ and
$\beta^{i}$ (chosen as a numerical convenience) which might make the initially
constructed Proca stars deviate even more from their Newtonian counterparts at
low compactness. These shortcomings can be thought of as adding initial
perturbations to the possible stationary solution for each of the three stars.

###  II.4 Extraction of mass and angular momentum

Following Clough (2021); Croft (2022), we define the conserved mass ($M$) and
the $z$ component of the angular momentum ($J_{3}$) as

$\displaystyle Q=Q_{0}+\int_{0}^{t}\mathcal{S}dt\,,$ (20)

where $Q=M,J_{3}$, and

$\displaystyle Q_{0}$
$\displaystyle\equiv\int_{\Sigma}d^{3}x\sqrt{\gamma}\,n_{\nu}\zeta^{\mu}{T_{\mu}}{}^{\nu},$
(21) $\displaystyle\mathcal{S}$
$\displaystyle\equiv\int_{\Sigma}d^{3}x\sqrt{\gamma}\,\alpha
T^{\mu}{}_{\nu}\nabla_{\mu}\zeta^{\nu}\,.$ (22)

The above quantities differ for $Q=M,J_{3}$ in the choice of $\zeta^{\mu}$:

$\zeta^{\mu}=\begin{cases}(1,0,0,0)&\textrm{for}\quad Q=M\,,\\\
(0,-y,x,0)&\textrm{for}\quad Q=J_{3}\,.\end{cases}$ (23)

The explicit expressions for $M_{0}$ and $(J_{3})_{0}$ are given by

$\displaystyle M_{0}$
$\displaystyle=\int_{\Sigma}d^{3}x\sqrt{\gamma}\left(\alpha\rho-\beta_{j}S^{j}\right)\,$
(24) $\displaystyle(J_{3})_{0}$
$\displaystyle=\int_{\Sigma}d^{3}x\sqrt{\gamma}\left(yS_{x}-xS_{y}\right)$
(25)

where we $\rho,S_{i}$ are defined in Eq. (7).

In Fig. 2 we summarize the mass-radius relationship of the constructed Proca
stars solved using the CTT procedure in Sec. II.3.2, using $Q$ in Eq. (20) as
the measure of mass. We show that the initial data we obtain for all three
kinds of stars agree approximately with the mass-radius curve in the non-
relativistic limit Jain and Amin (2022). The compactness

$\displaystyle\mathcal{C}\equiv\frac{M}{R_{\rm 95}},$ (26)

of these stars ranges from $\approx 0.01$ to $0.1$. Here, $R_{95}$ is defined
as the radius containing $95\%$ of the mass. We note that the above measure
for “mass”, is really a measure of the total energy, including the rest mass.
It will agree with the rest mass (defined in the Newtonian solutions) at low
compactness, but can show deviations at larger ones.

Figure 2: Initial mass and radius of Proca stars generated from non-
relativistic profiles. The dots show the Proca stars with Proca field
generated from Sec. II.3.1 and the spacetime metric solved from the CTT
procedure. The orange and blue dashed lines show the mass-radius curve as
predicted in the non-relativistic limit and under Newtonian gravity. Figure 3:
Total angular momentum against mass of the simulated Proca stars. For both the
hedgehog stars and the linearly polarized stars, the extracted angular
momentum from their respective simulations is consistent with zero. For the
circularly polarized stars, the extracted angular momentum satisfies the
relationship $J=\hbar M/m$, consistent with it being the spin angular momentum
discussed in Jain and Amin (2022).

##  III Results & Discussion

Figure 4: The mass-radius relationship of the simulated Proca stars. In all
panels the bars show ranges of the radial changes observed in simulations
(after an initial “settling in” period of $400m^{-1}$). For some stars, the
bars are replaced by arrows, indicating that the Proca star collapses into a
black hole. The grey dashed lines show the lines of constant compactness, and
the dark grey line show the compactness of a black hole, with its photosphere
as its radius, in isotropic coordinates. Left, Middle and Right panels show
the results for four hedgehog Proca stars, six linearly polarized Proca stars,
and five circularly polarized stars respectively. For $\mathcal{C}\lesssim
0.1$, the middle and left panels demonstrate the stability of compact,
gravitationally supported polarized stars. Near the upper bound of this range,
hedgehogs collapse at the lowest initial compactness ($\mathcal{C}\approx
0.06$), followed by linearly polarized ($\mathcal{C}\approx 0.08$), and then
(likely) circularly polarized stars ($\mathcal{C}>0.08$, although we were
unable to simulate collapse in circularly polarized stars). For non-collapsing
polarized stars, the mean of the radial variations provides insight into the
mass-radius relationship at these compactness.

We simulated the evolution of Proca stars (both polarized stars, as well as
hedgehog-like configurations for comparison), using constraint-fulfilling
initial Proca field profiles that are stationary under Newtonian gravity. The
evolution times were approximately $140$ cycles of the Proca field ($t_{\rm
sim}\approx 900m^{-1}$).

A sample evolution of the three different stars is shown in Fig. 1. During
their time evolution, all three stars exhibit radial oscillations in their
density, but do not disperse away for the duration of the simulation. The
period of this radial density oscillation is roughly $30-100$ times longer
than period of vector field cycle ($T=2\pi m^{-1}$). The radial oscillations
are likely excited due to the choice of initial data, see the last paragraph
of section II.3.

A summary of the mass radius-relationship of the three types of Proca stars
during the evolution (with varying initial compactess) is shown in Fig. 4. For
$\mathcal{C}\lesssim 0.1$, we see that all the stars show stable radial
oscillations with the amplitude of the oscillations being smallest at lowest
compactness. We take the survival of these stars for the duration of the
simulation, with perturbations introduced by imperfect initial data, as
evidence of the existence of long-lived, compact, polarized Proca stars within
full GR.

In Fig. 4, we can interpret the central values within the radial variations
(horizontal “error bars”) for each $M$ as defining the radius of a true
stationary solution at that mass. This should provide guidance in the
construction of the stationary solution at large compactness for the linearly
polarized stars, which is not known away from the Newtonian limit. However, we
urge caution here, especially as compactness gets high. For the hedgehog and
circularly polarized case, the mass radius relation is known at all
compactness Brito  _et al._ (2016). Using this known mass radius relationship
for hedgehog configurations, the radial variation does not include the
expected radius for $\mathcal{C}\approx 0.04$. We leave a more detailed
comparison with known stationary mass-radius relationships, as well as the
derivation of the stationary mass-radius curve for the linearly polarized case
in the high compactness regime for future work.

As the compactness approaches $\mathcal{C}=0.1$, we start seeing qualitatively
different behavior for the three types of stars. We start seeing collapse to
black holes in some types of stars. At $\mathcal{C}\approx 0.08$, the
circularly polarized star shows large radial oscillations, but does not
collapse to a black hole (left most point in the right panel of Fig. 4. The
linearly polarized one (middle panel), however does collapse at this same
compactness. The hedghog star (left panel) collapses at an even smaller
compactness of $\approx 0.06$ (where both linear and circularly polarized
stars are stable). That is, the initial compactness range $(0,\mathcal{C}$)
where hedgehogs are stable is smaller compared with that of the polarized
stars. Between the two polarized stars, the stable region of the linearly
polarized stars is smaller than that of the circularly polarized ones. We were
unable simulate any type of Proca stars with $\mathcal{C}>0.08$ because the
constraint equations solver Eqns. (18) and (19) failed to converge for these
configurations.

Limitations: We cannot control the magnitude of the perturbations around the
stationary solution induced by the imperfect initial conditions. Therefore, it
is possible that the perturbation is larger in the case of hedgehogs, which
might be a confounding factor leading to collapse to black holes at smaller
compactness. For this reason, we cannot “prove” that that polarized stars are
more stable than the hedgehog configuration. A controlled quantitative
analysis will be possible, after the stationary solutions (like the $m=1$ case
in Brito  _et al._ (2016)) of polarized Proca stars are found at high
compactness.

At large field amplitudes corresponding to highest compactness explored here
($|A_{i}|\sim 0.1M_{\rm pl}$, the self-interactions of the vector field might
not be ignorable. While polarized Proca stars with self-interactions for
relatively small compactness have been explored in the literature Zhang  _et
al._ (2022); Jain (2022); Jain and Amin (2023); Zhang and Ling (2023), the
large amplitude here might bring additional complications, see Mou and Zhang
(2022); Clough  _et al._ (2022); Coates and Ramazanoğlu (2022, 2023).

##  IV Summary & Implications

We simulated two types of polarized Proca stars (linear and circularly
polarized), along with hedgehog-like Proca stars for comparison, using general
relativistic field equations. The initial conditions were based on field
profiles of related Proca star solutions in Newtonian gravity Jain and Amin
(2022)(see our Fig. 2), scaled to a higher compactness.

Our key results are as follows (see Fig. 4):

  * •

We provided evidence that high-compactness polarized stars can be stable for
$\mathcal{C}\lesssim 0.1$.

  * •

As we increase the initial compactness from approximately $0.01$ to $0.1$, the
linearly polarized, circularly polarized, and hedgehog stars evolve away from
their initial configurations and towards new, and slightly different fixed
points.

  * •

At sufficiently high compactness, some types of stars collapse to black holes.
We found that circularly polarized stars avoid collapse to black holes up to
higher initial compactness than linearly polarized ones, which in turn avoid
collapse up to a higher initial compactness than hedgehog-like stars. The
large intrinsic spin angular momentum of circularly polarized stars (see Fig.
3) might be playing a role in their relative robustness to collapse.

For circularly polarized stars, we did not observe collapse to a black hole up
to $\mathcal{C}=0.08$. We were unable to simulate stars with initial
compactness $\gtrsim 0.08$ due to numerical limitations. An improved procedure
for constructing the initial data which allows for control of perturbations
away from the stationary solution is needed. This can be done using an
improved initial data formulation such as the one in Brito  _et al._ (2016).

We hope our findings provide new phenomenology that can be incorporated in the
search for “exotic” compact objects Cardoso and Pani (2019) through
gravitational and electromagnetic radiation. Polarized Proca stars can form in
dark photon/ vector dark matter fields Gorghetto  _et al._ (2022); Amin  _et
al._ (2022); Jain  _et al._ (2023); Chen  _et al._ (2023), potentially
providing access the nature of the dark sector.

For the purpose of gravitational wave physics, both the increased compactness,
and the polarization of the stars, can have important implications. The
increased maximal compactness of polarized stars in this paper (compared to,
for example, hedgehog stars), suggests that they will get closer before they
merge, resulting in the emitted gravitational radiation being different from
hedgehog stars. The polarization of the star can also impact the dynamics of
the binary merger of such stars through finite size effects such as tidal
deformability $\Lambda\propto\mathcal{C}^{-5}$ Poisson and Will (2014), before
and during merger.777See Herdeiro  _et al._ (2020) for the tidal deformability
of hedgehog Proca stars. More generally, the relevance of large tidal
deformabilty of boson stars for gravitational wave emission is discussed in
detail in Chia  _et al._ (2023) (also see Cardoso  _et al._ (2017); Sennett
_et al._ (2017)). In addition, circularly polarized stars with maximal
intrinsic spin can lead to spin-orbit and spin-spin effects before they merge.
During the final phase of the merger, the generated gravitational waves can
also be directly impacted by the polarization state of the star. Analysis of
mergers of compact scalar boson with angular momentum leads to rich dynamics
(see, for example Siemonsen and East (2023); Sanchis-Gual  _et al._ (2019b)).
A similar analysis is warranted for polarized Proca stars; for related recent
work see Sanchis-Gual  _et al._ (2022a); Calderón Bustillo  _et al._ (2021).

We have focused on stars constructed out of complex valued Proca fields for
convenience. Similar constructions can be carried out for real valued fields
(which might have a different lifetime). As in the case of axion stars
Hertzberg and Schiappacasse (2018); Levkov  _et al._ (2020); Amin and Mou
(2021); Amin  _et al._ (2021); Sanchis-Gual  _et al._ (2022b); Chung-Jukko
_et al._ (2023), such polarized Proca stars can also emit electromagnetic
radiation, with the novelty that the properties of the radiation now depends
on the polarization state of the Proca star Amin  _et al._ (2023)(for effects
on gravitational radiation, see Nakayama  _et al._ (2023)). In particular, the
polarization patterns in the outgoing radiation could provide a new handle on
the nature of the underlying dark fields. It could be interesting to construct
multimessenger signals (gravitational and electromagnetic waves) from merging
polarized Proca stars.

###### Acknowledgements.

The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing HPC, visualization, database, or
grid resources that have contributed to the research results reported within
this paper Stanzione  _et al._ (2020). URL: http://www.tacc.utexas.edu. MA
acknowledges early discussions of related work with Peter Adshead, Mudit Jain,
Kaloian Lozanov, Helvi Witek and Hong-Yi Zhang. We also acknowledge many
fruitful conversations with Emanuele Berti, Robin Croft, Liina M. Chung-Jukko,
Tamara Evstafyeva, Eugene Lim and Ulrich Sperhake. Z.W and T.H. are supported
by NSF Grants No. AST-2006538, PHY-2207502, PHY-090003 and PHY20043, and NASA
Grants No. 19-ATP19-0051, 20-LPS20-0011 and 21-ATP21-0010. MA is supported by
a DOE grant DE-SC0021619. This research project was conducted using
computational resources at the Maryland Advanced Research Computing Center
(MARCC).

##  Appendix A Convergence tests

On the initial spatial slice, the Hamiltonian constraint Eq. (16), the
momentum constraint Eq. (17), and the Proca constraint Eq. (10) must be
satisfied. We use the conformal-transverse-traceless (CTT) formalism to reduce
Eq. (16) and Eq. (17) into elliptical equations of the conformal factor $\psi$
and the extrinsic curvature $K_{ij}$ (see Baumgarte and Shapiro (2010) chap. 3
and appendix B).

Figure 5: Hamiltonian constraint violation for the densest
($\mathcal{C}=0.08$) circularly polarized Proca star with two different
resolutions $\Delta_{1}=0.40m^{-1}$ and $\Delta_{2}=0.20m^{-1}$. Here $r$ is
the radial distance (in code units) from the center of the star. The dashed
lines show the predicted Hamiltonian constraint of the high-resolution run for
second-order and fourth-order convergence. Figure 6: Convergence test for the
densest circularly polarized Proca star with three different resolutions. The
top panel shows $\rho$ at the center of the star (averaged over a sphere with
radius $0.2m^{-1}$). The bottom panel shows the difference between simulations
of medium and low (M-L) resolution (black). The dashed line shows the
predicted difference, assuming fourth-order convergence, based on the
difference between the high-resolution run and the medium-resolution run.

To show the validity of the solutions obtained from the CTT equations, we
performed convergence tests on the initial data of the densest circularly
polarized Proca star (generated with $\mu=0.10m$). In Fig. 5, we show the
Hamiltonian constraint $\mathcal{H}$ with two different resolutions. We see
that $\mathcal{H}$ from the high-resolution run is smaller than that of the
low-resolution run, and is consistent with second-order convergence towards
zero. The momentum constraint $\mathcal{M}$ and the Proca constraint
$\mathcal{C_{\epsilon}}$ both behave similarly to $\mathcal{H}$, and are
consistent with second-order convergence.

Aside from the initial data test, We also performed a convergence test for the
Proca field evolution scheme. Using the same circularly polarized Proca star,
We perform three runs with resolutions $\Delta_{1}=0.146m^{-1}$,
$\Delta_{2}=0.117m^{-1}$, and $\Delta_{3}=0.098m^{-1}$ at the position of the
star. We then plot the energy density $\rho$ at the center of the star
(averaged over several cells) in Fig. 6. On the top panel, the radial
oscillations of the star are visible in the form of changes in its central
density. Here, the differences between simulations with different resolutions
are negligible compared to the physical density variations.

##  Appendix B Numerical Details

We performed a series of runs, as shown in Figure 4, where we selected
profiles with $\mu/m=0.01,0.02,0.04,\text{and }0.06$ for the hedgehog Proca
stars. For the generated linearly polarized
$\mu/m=0.01,0.02,0.04,0.06,0.087,\text{ and }0.10$ and circularly polarized
Proca stars, we used $\mu/m=0.01,0.02,0.04,0.06,\text{and }0.10$.

To determine the $R_{95}$ radius, we integrated the energy density $\rho$ over
the volume. We extracted the energy density along the $x$ axes and assumed
that these profiles are symmetric for each axis. The mass $M(R)$ is given by
the integral:

$M(R)=4\pi\int_{0}^{R}\rho_{(x)}(r)r^{2}dr.$ (27)

Next, we determined the radius $R_{95\%}$ that contains 95% of the total mass:
$M(R_{95\%})=0.95M(\infty)$. We calculated the mass based on all three
directions $x$, $y$ and $z$ and in most cases, the three radii agreed and did
not vary significantly, only when approaching highly relativistic regimes
before collapse did we observe noticeable deviations from sphericity. Lastly,
this measure is dependent on gauge, so it should be interpreted with care.

## References

  * Graham  _et al._ (2016) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] .
  * Agrawal  _et al._ (2020) P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Phys. Lett. B 801, 135136 (2020), arXiv:1810.07188 [hep-ph] .
  * Co  _et al._ (2019) R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, Phys. Rev. D 99, 075002 (2019), arXiv:1810.07196 [hep-ph] .
  * Dror  _et al._ (2019) J. A. Dror, K. Harigaya, and V. Narayan, Phys. Rev. D 99, 035036 (2019), arXiv:1810.07195 [hep-ph] .
  * Bastero-Gil  _et al._ (2019) M. Bastero-Gil, J. Santiago, L. Ubaldi, and R. Vega-Morales, JCAP 04, 015 (2019), arXiv:1810.07208 [hep-ph] .
  * Long and Wang (2019) A. J. Long and L.-T. Wang, Phys. Rev. D 99, 063529 (2019), arXiv:1901.03312 [hep-ph] .
  * Adshead  _et al._ (2023) P. Adshead, K. D. Lozanov, and Z. J. Weiner, Phys. Rev. D 107, 083519 (2023), arXiv:2301.07718 [hep-ph] .
  * Kitajima and Takahashi (2023) N. Kitajima and F. Takahashi, “Resonant production of dark photons from axion without a large coupling,” (2023), arXiv:2303.05492 [hep-ph] .
  * Antypas  _et al._ (2022) D. Antypas _et al._ , “New Horizons: Scalar and Vector Ultralight Dark Matter,” (2022), arXiv:2203.14915 [hep-ex] .
  * Amin  _et al._ (2022) M. A. Amin, M. Jain, R. Karur, and P. Mocz, JCAP 08, 014 (2022), arXiv:2203.11935 [astro-ph.CO] .
  * Gorghetto  _et al._ (2022) M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, JCAP 08, 018 (2022), arXiv:2203.10100 [hep-ph] .
  * Jain  _et al._ (2023) M. Jain, M. A. Amin, J. Thomas, and W. Wanichwecharungruang, Phys. Rev. D 108, 043535 (2023), arXiv:2304.01985 [astro-ph.CO] .
  * Chen  _et al._ (2023) J. Chen, X. Du, M. Zhou, A. Benson, and D. J. E. Marsh, “Gravitational Bose-Einstein Condensation of Vector/Hidden Photon Dark Matter,” (2023), arXiv:2304.01965 [astro-ph.HE] .
  * Adshead and Lozanov (2021) P. Adshead and K. D. Lozanov, Phys. Rev. D 103, 103501 (2021), arXiv:2101.07265 [gr-qc] .
  * Jain and Amin (2022) M. Jain and M. A. Amin, Phys. Rev. D 105, 056019 (2022), arXiv:2109.04892 [hep-th] .
  * Sanchis-Gual  _et al._ (2022a) N. Sanchis-Gual, J. Calderón Bustillo, C. Herdeiro, E. Radu, J. A. Font, S. H. W. Leong, and A. Torres-Forné, Phys. Rev. D 106, 124011 (2022a), arXiv:2208.11717 [gr-qc] .
  * Calderón Bustillo  _et al._ (2021) J. Calderón Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, Phys. Rev. Lett. 126, 081101 (2021), arXiv:2009.05376 [gr-qc] .
  * Palenzuela  _et al._ (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D 75, 064005 (2007), arXiv:gr-qc/0612067 .
  * Palenzuela  _et al._ (2008) C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. D 77, 044036 (2008), arXiv:0706.2435 [gr-qc] .
  * Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [gr-qc] .
  * Bezares  _et al._ (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [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] .
  * 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] .
  * Widdicombe  _et al._ (2020) J. Y. Widdicombe, T. Helfer, and E. A. Lim, JCAP 01, 027 (2020), arXiv:1910.01950 [astro-ph.CO] .
  * Sanchis-Gual  _et al._ (2020) N. Sanchis-Gual, M. Zilhão, C. Herdeiro, F. Di Giovanni, J. A. Font, and E. Radu, Phys. Rev. D 102, 101504 (2020), arXiv:2007.11584 [gr-qc] .
  * Diamond  _et al._ (2023) M. D. Diamond, D. E. Kaplan, and S. Rajendran, JHEP 01, 136 (2023), arXiv:2112.09147 [hep-ph] .
  * Evstafyeva  _et al._ (2023) T. Evstafyeva, U. Sperhake, T. Helfer, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Class. Quant. Grav. 40, 085009 (2023), arXiv:2212.08023 [gr-qc] .
  * Siemonsen and East (2023) N. Siemonsen and W. E. East, “Binary boson stars: Merger dynamics and formation of rotating remnant stars,” (2023), arXiv:2302.06627 [gr-qc] .
  * Jaramillo  _et al._ (2022) V. Jaramillo, N. Sanchis-Gual, J. Barranco, A. Bernal, J. C. Degollado, C. Herdeiro, M. Megevand, and D. Núñez, Phys. Rev. D 105, 104057 (2022), arXiv:2202.00696 [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] .
  * 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] .
  * Sanchis-Gual  _et al._ (2017) N. Sanchis-Gual, C. Herdeiro, E. Radu, J. C. Degollado, and J. A. Font, Phys. Rev. D 95, 104028 (2017), arXiv:1702.04532 [gr-qc] .
  * Di Giovanni  _et al._ (2018) F. Di Giovanni, N. Sanchis-Gual, C. A. R. Herdeiro, and J. A. Font, Phys. Rev. D 98, 064044 (2018), arXiv:1803.04802 [gr-qc] .
  * Sanchis-Gual  _et al._ (2019a) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, Phys. Rev. Lett. 123, 221101 (2019a).
  * Arnowitt  _et al._ (1959) R. L. Arnowitt, S. Deser, and C. W. Misner, Phys. Rev. 116, 1322 (1959).
  * Croft (2022) R. Croft, “Local Continuity of Angular Momentum and Noether Charge for Matter in General Relativity,” (2022), arXiv:2203.13845 [gr-qc] .
  * Radia  _et al._ (2022) M. Radia, U. Sperhake, A. Drew, K. Clough, P. Figueras, E. A. Lim, J. L. Ripley, J. C. Aurrekoetxea, T. França, and T. Helfer, Class. Quant. Grav. 39, 135006 (2022), arXiv:2112.10567 [gr-qc] .
  * Zilhão  _et al._ (2015) M. Zilhão, H. Witek, and V. Cardoso, Class. Quant. Grav. 32, 234003 (2015), arXiv:1505.00797 [gr-qc] .
  * Hilditch (2013) D. Hilditch, Int. J. Mod. Phys. A 28, 1340015 (2013), arXiv:1309.2012 [gr-qc] .
  * Palenzuela  _et al._ (2010) C. Palenzuela, L. Lehner, and S. Yoshida, Phys. Rev. D 81, 084007 (2010), arXiv:0911.3889 [gr-qc] .
  * Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro,  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge University Press, 2010).
  * Clough (2021) K. Clough, Class. Quant. Grav. 38, 167001 (2021), arXiv:2104.13420 [gr-qc] .
  * Zhang  _et al._ (2022) H.-Y. Zhang, M. Jain, and M. A. Amin, Phys. Rev. D 105, 096037 (2022), arXiv:2111.08700 [astro-ph.CO] .
  * Jain (2022) M. Jain, Phys. Rev. D 106, 085011 (2022), arXiv:2205.03418 [hep-ph] .
  * Jain and Amin (2023) M. Jain and M. A. Amin, JCAP 04, 053 (2023), arXiv:2211.08433 [astro-ph.CO] .
  * Zhang and Ling (2023) H.-Y. Zhang and S. Ling, JCAP 07, 055 (2023), arXiv:2305.03841 [astro-ph.CO] .
  * Mou and Zhang (2022) Z.-G. Mou and H.-Y. Zhang, Phys. Rev. Lett. 129, 151101 (2022), arXiv:2204.11324 [hep-th] .
  * Clough  _et al._ (2022) K. Clough, T. Helfer, H. Witek, and E. Berti, Phys. Rev. Lett. 129, 151102 (2022), arXiv:2204.10868 [gr-qc] .
  * Coates and Ramazanoğlu (2022) A. Coates and F. M. Ramazanoğlu, Phys. Rev. Lett. 129, 151103 (2022), arXiv:2205.07784 [gr-qc] .
  * Coates and Ramazanoğlu (2023) A. Coates and F. M. Ramazanoğlu, Phys. Rev. Lett. 130, 021401 (2023), arXiv:2211.08027 [gr-qc] .
  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Poisson and Will (2014) E. Poisson and C. M. Will,  _Gravity: Newtonian, Post-Newtonian, Relativistic_ (Cambridge University Press, 2014).
  * Herdeiro  _et al._ (2020) C. A. R. Herdeiro, G. Panotopoulos, and E. Radu, JCAP 08, 029 (2020), arXiv:2006.11083 [gr-qc] .
  * Chia  _et al._ (2023) H. S. Chia, T. D. P. Edwards, D. Wadekar, A. Zimmerman, S. Olsen, J. Roulet, T. Venumadhav, B. Zackay, and M. Zaldarriaga, (2023), arXiv:2306.00050 [gr-qc] .
  * Cardoso  _et al._ (2017) V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, Phys. Rev. D 95, 084014 (2017), [Addendum: Phys.Rev.D 95, 089901 (2017)], arXiv:1701.01116 [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] .
  * Sanchis-Gual  _et al._ (2019b) N. Sanchis-Gual, C. Herdeiro, J. A. Font, E. Radu, and F. Di Giovanni, Phys. Rev. D 99, 024017 (2019b), arXiv:1806.07779 [gr-qc] .
  * Hertzberg and Schiappacasse (2018) M. P. Hertzberg and E. D. Schiappacasse, JCAP 11, 004 (2018), arXiv:1805.00430 [hep-ph] .
  * Levkov  _et al._ (2020) D. Levkov, A. Panin, and I. Tkachev, Phys. Rev. D 102, 023501 (2020), arXiv:2004.05179 [astro-ph.CO] .
  * Amin and Mou (2021) M. A. Amin and Z.-G. Mou, JCAP 02, 024 (2021), arXiv:2009.11337 [astro-ph.CO] .
  * Amin  _et al._ (2021) M. A. Amin, A. J. Long, Z.-G. Mou, and P. Saffin, JHEP 06, 182 (2021), arXiv:2103.12082 [hep-ph] .
  * Sanchis-Gual  _et al._ (2022b) N. Sanchis-Gual, M. Zilhão, and V. Cardoso, Phys. Rev. D 106, 064034 (2022b), arXiv:2207.05494 [gr-qc] .
  * Chung-Jukko  _et al._ (2023) L. M. Chung-Jukko, E. A. Lim, D. J. E. Marsh, J. C. Aurrekoetxea, E. de Jong, and B.-X. Ge, (2023), arXiv:2302.10100 [hep-ph] .
  * Amin  _et al._ (2023) M. A. Amin, A. J. Long, and E. D. Schiappacasse, JCAP 05, 015 (2023), arXiv:2301.11470 [hep-ph] .
  * Nakayama  _et al._ (2023) K. Nakayama, F. Takahashi, and M. Yamada, JCAP 08, 058 (2023), arXiv:2306.12961 [hep-ph] .
  * Stanzione  _et al._ (2020) D. Stanzione, J. West, R. T. Evans, T. Minyard, O. Ghattas, and D. K. Panda, in  _Practice and Experience in Advanced Research Computing_, PEARC ’20 (Association for Computing Machinery, New York, NY, USA, 2020) p. 106–111.