# Fixing the dynamical evolution of self-interacting vector fields

Marcelo Rubio1,2<EMAIL_ADDRESS>Guillermo Lara3<EMAIL_ADDRESS>Miguel
Bezares4,5<EMAIL_ADDRESS>Marco Crisostomi6,7
<EMAIL_ADDRESS>Enrico Barausse1,2<EMAIL_ADDRESS>1SISSA, Via Bonomea
265, 34136 Trieste, Italy and INFN Sezione di Trieste
2IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy
3Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, 14476 Potsdam, Germany
4Nottingham Centre of Gravity, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK
5School of Mathematical Sciences, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK
6TAPIR, Division of Physics, Mathematics, and Astronomy, California Institute
of Technology, Pasadena, CA 91125, USA
7Dipartimento di Fisica, Università di Pisa, Largo B. Pontecorvo 3, 56127
Pisa, Italy

###### Abstract

Numerical simulations of the Cauchy problem for self-interacting massive
vector fields often face instabilities and apparent pathologies. We explicitly
demonstrate that these issues, previously reported in the literature, are
actually due to the breakdown of the well-posedness of the initial-value
problem. This is akin to shortcomings observed in scalar-tensor theories when
derivative self-interactions are included. Building on previous work done for
$k$-essence, we characterize the well-posedness breakdowns, differentiating
between Tricomi and Keldysh-like behaviors. We show that these issues can be
avoided by “fixing the equations”, enabling stable numerical evolutions in
spherical symmetry. Additionally, we show that for a class of vector self-
interactions, no Tricomi-type breakdown takes place. Finally, we investigate
initial configurations for the massive vector field which lead to
gravitational collapse and the formation of black holes.

##  I Introduction

The detection of gravitational waves (GWs) from the merger of black holes
(BHs) and neutron stars (NSs) LV_Catalog2019 ; LIGO_GW170817 ; LIGO_GW190521 ;
LIGO_GW170817_2  has undoubtedly opened an exciting new window for testing
fundamental physics in extreme conditions, and for understanding the nature of
compact objects CardosoPani_LRR_2019 . Moreover, next-generation GW detectors
LISA_WP_Baker:2019nia ; LISA_Littenberg:2019grv ; ET_Branchesi:2023mws  will
be capable of probing modifications to General Relativity (GR) with higher
precision, and will provide new data that may prove crucial for unveiling
fundamental questions about the internal structure of NSs, the validity of the
Kerr hypothesis Israel:1967wq ; Hawking:1971vc ; Robinson:1975bv  and the
possible existence of “exotic” compact objects, many of which have been
proposed in the literature Bezares:2024btu ; Berti2006 ; ECO1 ; ECO2 ; ECO3 ;
ECO4 ; Kesden2005 .

An example of the latter are boson stars (BSs) JETZER1992163 ; BSs1 ;
Schunck:2003kk ; PhysRevD.80.084023 ; PhysRevLett.57.2485 , defined as regular
gravitating solitons PhysRev.172.1331 . While BSs often refer to spin-0 scalar
Bose-Einstein condensates PhysRevD.38.2376 ; PhysRevD.39.1257 ;
lee1992nontopological ; Lee:2008jp ; Lee:2008ab , spin-1 vector BSs, commonly
referred to as Proca stars (PSs), have also been put forward as models of
“black-hole mimickers” BRITO2016291 ; Herdeiro:2023lze ;
PhysRevLett.123.221101 ; Herdeiro:2023wqf . In fact, shadows and lensing
around PSs have been simulated and are compatible with current constraints
from the Event Horizon Telescope (EHT) Pheno_Sengo2022 ; Sengo:2024pwk . In
addition, dynamical mechanisms for the formation of PSs via gravitational
cooling have also been put forward in the literature Pheno2_Herdeiro2023 ;
PhysRevD.98.064044 ; Shen:2016acv .

The simplest models for compact exotic stars are given by a complex scalar or
vector field minimally coupled to gravity, and they lead to peculiar stability
properties. For instance, scalar BSs are generally stable under non-spherical
perturbations, but static, spherical PSs might display generic instabilities
Pheno2_Herdeiro2023 . Although some of these issues can be resolved for
stationary rotating stars PhysRevLett.123.221101 , the situation notably
changes when accounting for self-interactions. This fact has motivated deeper
investigations on more general Einstein-Proca systems (see for instance
Herdeiro:2023lze ). Additionally, explorations on self-interacting massive
vector fields reported superradiant instabilities due to exponential
amplification of bound state modes, and energy/angular momentum extraction
from spinning BHs EastPRL17 ; East17 ; Dolan18 .

A numerical study of self-interacting Proca fields on a Kerr background was
performed in Proca_CloughPRL2022 , also in the context of superradiance, but
reporting instabilities which were attributed to the development of ghost (or
tachyonic) modes. Similar issues were found afterwards, e.g. in
Proca_Coates1_PRL2022 ; Barvinsky:1989wy ; Zong-Gang22 . In Ref.
ShortPaper_Barausse2022 , some of us pointed out that these pathologies are
actually due not to ghost/tachyonic modes but to the breakdown of the well-
posedness of the corresponding Cauchy problem. Indeed, although the principal
part of the equations of motion governing vector theories might at first seem
the same as for Maxwell equations (and therefore innocuous), self-interactions
drastically modify it and potentially make the Cauchy problem ill-posed. This
was confirmed by further studies Proca_Coates2_PRL2023 ; Proca_Coates3_PRD2023
; Proca_Coates4_PRD2023 ; Proca_Silva2022 ; Aoki22 , which report
instabilities occurring when the dynamical fields reach a point where the
initial-value problem is mathematically no longer well-defined.

Another key point raised by Ref. ShortPaper_Barausse2022  is that the well-
posedness of the Cauchy Problem of self-interacting massive vector fields is
completely analogous to that of k-essence theories (i.e., scalar theories with
first-order derivative self-interactions), which was thoroughly studied and
characterized in the last few years Bernard19 ; Bezares_Kdynamics2020 ;
Lotte_PRL21 ; Bezares21KS ; Bezares:2021dma ; Figueras:2020dzx . Ref.
ShortPaper_Barausse2022  showed that this analogy becomes evident when
restoring the $U(1)$ gauge symmetry via a Stückelberg transformation. Indeed,
this transformation highlights that vector field self-interactions actually
hide derivative self-interactions of the Stückelberg scalar field (which
corresponds to the longitudinal mode of the vector field), radically modifying
the principal part of the evolution system. As a result, strategies similar to
those successfully adopted for k-essence could also be implemented for
simulating self-interacting vector fields, avoiding the aforementioned
breakdown of the initial-value problem.

In this work, and following ShortPaper_Barausse2022 , we carry out numerical
investigations of massive vector fields in the Stückelberg formulation, and
study the well-posedness of the associated initial-value problem. In
particular, we show an explicit example of a Tricomi-type evolution (i.e., one
leading to a change of character of the equations from hyperbolic to parabolic
or elliptic) and two ways of avoiding it. The first one is provided by the
“fixing-the-equations” technique Fix_CayusoPRD2017 , which consists of
modifying the principal part of the evolution system by introducing extra
dynamical fields satisfying ad hoc differential equations. This approach has
already been applied with success to several gravitational theories
Fix_Franchini2022 ; Bezares21KS ; CayusoNonIterative20 ; MemoUV2022 ;
Lara:2024rwa ; Corman:2024cdr . It was also considered in Ref.
PhysRevD.108.L101501 , in the context of massive fields with quadratic self-
interactions, performing evolutions in flat space and in one spatial
dimension. The second alternative explored in the present work is to account
for a cubic self-interacting term in the Lagrangian, without having to “fix”
the equations. In both cases, we perform full numerical relativity simulations
in spherical symmetry. We also explore configurations leading to gravitational
collapse and the formation of black holes in theories with cubic self-
interactions.

The outline of the paper is the following: in section II, we revisit the
theory of self-interacting massive vector fields in the Stückelberg language,
its corresponding Cauchy problem, and the “fixing-the-equations” technique. In
section III, we describe our numerical set-up, initial data and boundary
conditions in detail. The main results of the paper are reported in section
IV, which is divided into three parts. Firstly, we show numerically an example
of a Tricomi-type breakdown of the Cauchy problem and present how the “fixing”
approach restores the well-posedness during the evolution. Secondly, we allow
for a cubic self-interaction in the Lagrangian, and show that no Tricomi-type
breakdown is developed during the whole evolution. Finally, we study
gravitational collapse in the case of cubic self-interactions, reporting an
explicit example of an initial configuration whose evolution leads to the
formation of a black hole. An overall discussion of our results is left for
Section V. Throughout this work, we use the mostly-plus signature for the
spacetime metric, $(-,+,+,+)$, while the employed units are specified at the
beginning of Section IV.

##  II Preliminaries

###  II.1 Massive vector fields

We study the dynamics of a self-interacting real vector field $A^{\mu}$ with
mass $m$, coupled to GR. The corresponding action is given by

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\left[\frac{R}{16\pi G}-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad-\frac{m^{2}}{2}A_{\mu}A^{\mu}+V(A_{\mu}A^{\mu})\right]\,,$

where $F_{\mu\nu}:=\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}$ is the vector
strength, and the self-interaction potential is

$V(A_{\mu}A^{\mu})=\frac{\beta}{4}(A_{\mu}A^{\mu})^{2}-\frac{\gamma}{8\Lambda^{2}}(A_{\mu}A^{\mu})^{3}.$
(2)

Here, $\Lambda$ is the energy scale suppressing the higher-order interactions,
whereas $\beta$ and $\gamma$ are $\mathcal{O}(1)$ dimensionless coupling
constants for the quadratic and cubic interactions, respectively. The field
$A_{\mu}$ propagates two transverse and one longitudinal modes. Although
action (II.1) lacks the $U(1)$ gauge symmetry $A_{\mu}\to
A_{\mu}+\nabla_{\mu}f$ of Maxwell electrodynamics, it is possible to restore
it at the cost of introducing an extra dynamical field. In fact, by performing
the Stückelberg transformation

$A_{\mu}\to A_{\mu}-\frac{1}{m}\nabla_{\mu}\phi,$ (3)

the product $A_{\mu}A^{\mu}$ transforms as

$A_{\mu}A^{\mu}\to
A_{\mu}A^{\mu}-\frac{2}{m}A^{\mu}\nabla_{\mu}\phi+\frac{1}{m^{2}}\nabla_{\mu}\phi\nabla^{\mu}\phi,$
(4)

and this particular mixing between $A_{\mu}$ and $\phi$ remains invariant
under the gauge transformation

$\begin{pmatrix}A_{\mu}\\\
\phi\end{pmatrix}\to\begin{pmatrix}A_{\mu}+\nabla_{\mu}f\\\
\phi+mf\end{pmatrix}.$ (5)

Thus, the theory obtained by replacing (3) into the action (II.1) results
invariant under (5). This symmetry induces a gauge freedom in the choice of
the dynamical fields, which can be fixed in different ways. Adopting for
instance the unitary gauge (i.e., $\phi=0$), one would recover the original
Lagrangian. The Lorenz gauge (i.e., $\nabla_{\mu}A^{\mu}=0$) is instead useful
for decoupling the vector and scalar degrees of freedom at high energies.

It is worthwhile to introduce the gauge invariant derivative of the scalar
field

$D_{\mu}\phi:=\nabla_{\mu}\phi-mA_{\mu}\,,$ (6)

so that, under the transformation (3), one has

$m^{2}A_{\mu}A^{\mu}\to D_{\mu}\phi D^{\mu}\phi=:X\,.$ (7)

Then, after the Stückelberg’s procedure, the action (II.1) can be re-written
as

$S=\int{d^{4}x\sqrt{-g}\left[\frac{R}{16\pi
G}-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+K(X)\right]}\,,$ (8)

where

$K(X):=-\frac{1}{2}X+\frac{\beta}{4m^{4}}X^{2}-\frac{\gamma}{8\Lambda^{2}m^{6}}X^{3}\,.$
(9)

Variations of (8) with respect to the dynamical fields
$\\{g_{\mu\nu},A^{\mu},\phi\\}$ lead to the equations of motion

$\displaystyle G_{\mu\nu}=8\pi
G\left(T^{\text{\tiny{(EM)}}}_{\mu\nu}+T^{(\phi)}_{\mu\nu}\right)\,,$
$\displaystyle\nabla_{\mu}F^{\mu\nu}=J_{(\phi)}^{\nu}\,,$ (10)
$\displaystyle\nabla_{\mu}\left[K^{\prime}(X)D^{\mu}\phi\right]=0\,,$

where $G_{\mu\nu}$ is the Einstein’s tensor,

$T^{\text{\tiny{(EM)}}}_{\mu\nu}=F_{\mu\rho}F_{\nu}{}^{\rho}-\frac{1}{4}F_{\rho\sigma}F^{\rho\sigma}g_{\mu\nu}$
(11)

is the energy-momentum tensor associated to $F_{\mu\nu}$, and

$T^{(\phi)}_{\mu\nu}=K(X)g_{\mu\nu}-2K^{\prime}(X)D_{\mu}\phi D_{\nu}\phi$
(12)

is the one associated with the Stückelberg field. The source of the equation
for $F^{\mu\nu}$ is

$J_{(\phi)}^{\mu}=-2mK^{\prime}(X)D^{\mu}\phi\,,$ (13)

where a prime denotes a derivative with respect to $X$.

Notice that, because $D_{\mu}\phi$ in Eq. (6) contains the vector field, the
action (8) still retains the full mixing between $\phi$ and $A^{\mu}$.
However, written in this form, it becomes easy to single out the respective
principal parts, which correspond to the high-energy/momentum decoupling
limit. Because this limit determines the well-posedness of the corresponding
initial-value problem, one can see from the action (8) that the scalar sector
exhibits first-order derivative self-interactions (as in $k$-essence
theories), which may lead to pathologies during dynamical evolutions Bernard19
; Bezares_Kdynamics2020 ; Bezares21KS ; MemoUV2022 ; ShortPaper_Barausse2022 ;
Lotte_PRL21 . This calls for a full well-posedness analysis of the equations
(II.1), particularly the scalar sector.

###  II.2 The Cauchy problem of $k$-essence theories

We will now briefly review the initial-value problem for $k$-essence theories,
following Refs. Bernard19 ; Bezares_Kdynamics2020 ; Bezares21KS ; Bezares21KS
; Figueras:2020dzx . Generally, we refer to the Cauchy problem of a system of
partial differential equations as well-posed Hadamard1908  if, for any given
initial-data set (e.g. at $t=t_{0}$), there exists a time interval
$I=[t_{0},T]$ in which: i) a solution exists; ii) it is unique; iii) it is a
continuous function of the initial data (with respect to the topologies where
the data and solutions are defined). The evolution is governed by the
principal part of the system of equations, as it contains all the information
about the propagation speeds of the different modes Kreiss70 ;
geroch1996partial . An algebraic strategy to assess these non-trivial
mathematical conditions is provided by the concept of hyperbolicity
Friedrichs54 ; Kreiss70 ; friedrichs1971systems ; geroch1996partial ; that is,
a set of algebraic conditions that the principal part must satisfy for the
system to admit a well-posed (sometimes said hyperbolic) initial-value
formulation. In the particular case of system (II.1), the scalar sector
governs the hyperbolicity of the system, and might be the cause of a possible
change of character during the evolution. In fact, its principal part can be
recast as a modified Klein-Gordon equation, with an effective metric
$\gamma^{\mu\nu}$ given by

$\gamma^{\mu\nu}=g^{\mu\nu}+\frac{2K^{\prime\prime}(X)}{K^{\prime}(X)}D^{\mu}\phi
D^{\nu}\phi\,.$ (14)

Thus, the system is strongly hyperbolic if and only if
$\det(\gamma^{\mu\nu})<0.$ Nevertheless, nothing prevents this condition from
potentially failing during the evolution, even when starting from regular
initial data sets. The highly non-linear evolution could cause the determinant
of the effective metric to become zero (implying that the equations become
parabolic), or it could give rise to very large (or even diverging)
characteristic speeds. The breakdown of the Cauchy problem due to these
shortcomings are usually referred to as Tricomi and Keldysh types,
respectively Bernard19 ; Bezares_Kdynamics2020 ; Ripley19 ; Ripley:2022cdh .

Let us analyze under which conditions these pathologies could appear, in the
particular case of a spherically symmetric configuration. The corresponding
line element can be expressed as

$ds^{2}=-\alpha(t,r)^{2}dt^{2}+g_{rr}(t,r)dr^{2}+r^{2}g_{\theta\theta}(t,r)d\Omega^{2},$
(15)

where $\alpha$ is the lapse function, $g_{rr}$ and $g_{\theta\theta}$ are
positive fields, and $d\Omega^{2}$ is the line element of the unitary sphere.
As pointed out in Ref. Bezares_Kdynamics2020 , the determinant of the
effective metric reads

$\det(\gamma^{\mu\nu})=-\frac{1}{\alpha^{2}g_{rr}}\left(1+\frac{2K^{\prime\prime}(X)}{K^{\prime}(X)}X\right),$
(16)

and its dynamical evolution can be assessed by examining the eigenvalues of
$\gamma^{\mu\nu}$, namely

$\lambda_{\pm}=\frac{1}{2}\left[\gamma^{tt}+\gamma^{rr}\pm\sqrt{(\gamma^{tt}-\gamma^{rr})^{2}+4(\gamma^{tr})^{2}}\right]\,.$
(17)

The characteristic speeds of the scalar equation of system (II.1) are given by

$V_{\pm}=-\frac{\gamma^{tr}}{\gamma^{tt}}\pm\sqrt{\frac{-\det(\gamma^{\mu\nu})}{(\gamma^{tt})^{2}}}\,.$
(18)

The hyperbolicity condition $\det(\gamma^{\mu\nu})<0$ restricts the possible
values of the coupling constants for the function $K(X)$ given in Eq. (9). As
an example, looking at Eq. (16) one can see that the choice $\beta>0$ and
$\gamma=0$ could give rise to the Tricomi pathology for certain initial data.
On the other hand, choosing for instance $\beta=0$ and $\gamma>0$ avoids it
altogether, although a Keldysh type behavior (i.e., $\gamma^{tt}\to 0$ in Eq.
(18)) can still occur. Nonetheless, we stress that the diverging speeds
characterizing the Keldysh type behavior do not violate hyperbolicity. The
system remains hyperbolic, but in practice numerical evolutions become
impossible as the Courant-Friedrichs-Lewy (CFL) condition forces the time step
to vanish if the evolution is performed with an explicit method. We also
notice that that Keldysh type behaviors depend on gauge; i.e., they can be (in
principle) avoided with a judicious gauge choice Bezares_Kdynamics2020 ;
Bezares:2021dma . In this sense, Keldysh type breakdown of well-posedness
could be a practical problem but not a fundamental pathology. In this work, we
numerically explore how these issues appear in the context of self-interacting
massive vector fields and show possible ways to cure them.

###  II.3 The fixing-the-equations approach

An operational method to deal with the “ill-suited” evolutions reviewed before
(which generally appear in the context of modified theories of gravity) is
through the so-called fixing-the-equations approach Fix_CayusoPRD2017 ;
Fix_Allwright2019 ; Fix_Cayuso2_PRD2020 . This possibility has been proven to
be useful in the strong, non-linear, and highly dynamical regime of
$k-$essence theories Bezares21KS ; MemoUV2022 .

Inspired by relativistic theories for dissipative fluids Israel761 ; Israel762
; Israel79  and by numerical methods to deal with interfaces between touching
numerical grids LuisPC , the idea of the fixing-the-equations approach is to
modify ad hoc the evolution equations (with particular focus on their higher-
derivative contributions), so that the high-frequency modes are controlled. In
this way, the norm of the solution remains finite and bounded at all times,
rendering the system well-posed. To do so, one introduces new auxiliary
variables and a timescale on which they settle to their true/physical values
through a driver equation. This method allows for capturing the main features
of the non-linear behavior of the system and can be easily implemented
numerically.

In MemoUV2022 , this method has been utilized for simulating spherical
collapse in quadratic k-essence after a Tricomi-type breakdown occurs. In
PhysRevD.108.L101501 , it was also implemented for evolving massive vector
fields with quartic self-interactions in one spatial dimension in flat space.
With a similar spirit, here we introduce a variable $\Sigma$, a timescale
$\tau$ and a driver equation in order to replace the evolution for the
Stückelberg field in (II.1) by the system

$\displaystyle\nabla_{\mu}\left(\Sigma D^{\mu}\phi\right)$ $\displaystyle=$
$\displaystyle 0\,,$ (19)
$\displaystyle\partial_{t}\Sigma+\frac{\Sigma-K^{\prime}(X)}{\tau}$
$\displaystyle=$ $\displaystyle 0.$ (20)

This modification damps all the modes with frequencies larger than $1/\tau$
and forces $\Sigma\to K^{\prime}(X)$ as $\tau\to 0$, resulting in a strongly
hyperbolic evolution. In fact, for the spherically symmetric Ansatz (15), and
for $\Sigma\neq 0$, the system admits the characteristic speeds

$V^{o}_{\scriptsize{\mbox{fixed}}}=0\,,\quad
V^{\pm}_{\scriptsize{\mbox{fixed}}}=\pm\frac{\alpha}{\sqrt{g_{rr}}}\,,$ (21)

which agree with the ones of the (spherically-symmetric) Klein-Gordon equation
in GR.

##  III Numerical Implementation

In this section, we give details about the methodology used for the numerical
simulations, including the evolution equations in spherical symmetry, the
initial data and the implemented boundary conditions.

We performed evolutions of both the original (II.1) and the fixed (19) systems
introduced in the previous section. For doing so, we extended our previous
code used in Bezares_Kdynamics2020 ; Bezares21KS ; Bezares21KS ; MemoUV2022
by including the equations for the electromagnetic sector. This code
implements a high-resolution shock-capturing (HRSC) finite-difference scheme,
originally developed in Alic:2007ev . It was used for simulating black holes
and for studying the dynamics of boson stars, fermion-boson stars and
anisotropic stars Z3_Bernal09 ; Z3_Valdez-Alvarado2013 ; Raposo:2018rjn . The
numerical scheme can be considered as a fourth-order finite-difference one
plus a third-order adaptive dissipation, with the dissipation coefficient
given by the maximum propagation speed at each grid point. Time evolution is
performed using the method of lines, with a third-order Runge-Kutta scheme. In
the context of this paper, the choice of an HRSC method is motivated by the
fact that $k-$essence theories might develop caustics/shocks during evolution,
even from smooth initial data, as reported in Felder:2002sv ; Reall:2014sla ;
Babichev:2017lrx ; Tanahashi:2017kgn .

###  III.1 Evolution equations

In order to set up evolution equations, we perform a $(3+1)$-decomposition by
introducing a foliation $\\{\Sigma_{t}\\}_{t\in\mathbb{R}}$ of spatial
hypersurfaces with normal $n_{\mu}=(-\alpha,\,0)$. We define the extrinsic
curvature of each $\Sigma_{t}$ as $K_{ij}=-\mathcal{L}_{n}h_{ij}/2$, where
$h_{ij}$ is the induced metric.

By taking the spherically-symmetric Ansatz (15) for the spacetime metric, the
trace of $K_{ij}$ reads

$K=K_{r}{}^{r}+2K_{\theta}{}^{\theta}.$ (22)

We denote the parallel and orthogonal projections of $A_{\mu}$ (with respect
to $n^{\mu}$) as

$\displaystyle A_{\parallel}$ $\displaystyle=$
$\displaystyle-n^{\mu}A_{\mu}\,,$ $\displaystyle A_{\perp}$ $\displaystyle=$
$\displaystyle\delta_{r}{}^{\mu}A_{\mu}\,,$ (23)

respectively.

####  III.1.1 Metric evolution

All the simulations were performed using the Z3 formulation PhysRevD.66.084013
, which is a strongly-hyperbolic first-order reduction for Einstein’s
equations. It can be obtained from the Z4 formulation by a symmetry breaking
procedure. We made use of the spherically-symmetry version developed in
Z3_Bernal09 ; Z3_Valdez-Alvarado2013 , and written in flux-conservative form
PhysRevD.40.1022 . The dynamical variables for the metric evolution are
$\\{\alpha_{r},D_{rr}{}^{r},D_{r\theta}{}^{\theta},Z_{r},K_{r}{}^{r},K_{\theta}{}^{\theta}\\}$,
where $\alpha_{r}=\partial_{r}\alpha/\alpha$,

$D_{ri}{}^{i}=\frac{g^{ii}}{2}\partial_{r}g_{ii},$ (24)

and $Z_{r}$ is the radial component of $Z_{i}$111We stress that the Z3
formulation considers and additional evolution field, $Z_{i}$, which accounts
for the Momentum constraint, and thus it should remain close to zero
throughout the whole numerical evolution. In spherical symmetry, only the
radial component $Z_{r}$ is dynamical. We refer the reader to Refs.
Z3_Bernal09 ; Z3_Valdez-Alvarado2013 , where the explicit set of field
equations in this formulation is displayed.. Finally, we fix the gauge freedom
by setting a “$1+\log$” slicing condition PhysRevLett.75.600  in normal
coordinates (with zero shift).

####  III.1.2 Scalar sector

We do a first-order reduction of the scalar field equation, by introducing the
variables

$\Pi=-\dfrac{\partial_{t}\phi}{\alpha}\,,\quad\Phi=\partial_{r}\phi.$ (25)

It is also useful to define the auxiliary “shifted” variables

$\displaystyle\tilde{\Pi}$ $\displaystyle\equiv$ $\displaystyle\Pi-
mA_{\parallel}$ (26) $\displaystyle\tilde{\Phi}$ $\displaystyle\equiv$
$\displaystyle\Phi-mA_{\perp},$ (27)

as they allow to express the kinetic term as

$X=-\tilde{\Pi}^{2}+g^{rr}\tilde{\Phi}^{2}\,.$ (28)

In order to deal with shocks, we write the evolution equations in flux-
conservative form,

$\partial_{t}\bm{U}_{(\phi)}+\partial_{r}\bm{F}_{\mathrm{(\phi)}}(\bm{U}_{\mathrm{(\phi)}})+\bm{S}_{\mathrm{(\phi)}}(\bm{U}_{\mathrm{(\phi)}})=\bm{0},$
(29)

where

$\bm{U}_{(\phi)}=\begin{bmatrix}\phi\\\ \Phi\\\
\sqrt{g_{rr}}g_{\theta\theta}K^{\prime}(X)\tilde{\Pi}\end{bmatrix}\,,$ (30)
$\bm{F}_{\mathrm{(\phi)}}(\bm{U}_{\mathrm{(\phi)}})=\begin{bmatrix}0\\\
\alpha\Pi\\\ \displaystyle\frac{\alpha
g_{\theta\theta}}{\sqrt{g_{rr}}}K^{\prime}(X)\tilde{\Phi}\end{bmatrix}\,,$
(31)

and

$\bm{S}_{(\phi)}(\bm{U}_{(\phi)})=\begin{bmatrix}\alpha\Pi\\\ 0\\\
\displaystyle\frac{2}{r}\frac{\alpha
g_{\theta\theta}}{\sqrt{g_{rr}}}K^{\prime}(X)\tilde{\Phi}\end{bmatrix}.$ (32)

Finally, the projections of the energy-momentum tensor (12) associated to the
scalar field are

$\displaystyle\tau^{(\phi)}$
$\displaystyle=-K(X)-2K^{\prime}(X)\tilde{\Pi}^{2},$ (33)
$\displaystyle{S^{(\phi)}}^{r}$
$\displaystyle=-2g^{rr}K^{\prime}(X)\tilde{\Pi}\tilde{\Phi},$ (34)
$\displaystyle{{S^{(\phi)}}_{r}}^{r}$
$\displaystyle=K(X)-2g^{rr}K^{\prime}(X)\tilde{\Phi}^{2},$ (35)
$\displaystyle{{S^{(\phi)}}_{\theta}}^{\theta}$ $\displaystyle=K(X).$ (36)

####  III.1.3 U(1) vector sector

The equations for the electromagnetic sector are usually given in terms of the
electric and magnetic fields $E^{\mu}=-t^{\nu}F_{\nu}{}^{\mu}$ and
$B^{\mu}=-t^{\nu}{}^{*}F_{\nu}{}^{\mu}$, respectively. These are measured with
respect to an Eulerian observer determined by a timelike, unitary and future-
pointing frame $t^{\mu}$. It is however natural to consider the evolution of
$A_{\mu}$, such that $F_{\mu\nu}=2\nabla_{[\mu}A_{\nu]}$. In spherical
symmetry, the only non-trivial components are $A_{\parallel}$ and $A_{\perp}$,
introduced in (III.1). This choice automatically implies that
$B^{i}=\varepsilon^{ijk}D_{j}A_{k}=0$, and the only nontrivial component for
the electric field is the radial one, namely $E^{r}$. Moreover, since the
scalar equation in system (II.1) is coupled with the vector field, we also
need to provide evolution equations for it. Like we did for the scalar sector,
we express them as

$\partial_{t}\,\bm{U}_{\mathrm{E}}+\partial_{r}\,\bm{F}_{\mathrm{E}}(\bm{U}_{\mathrm{E}})+\bm{S}_{\mathrm{E}}(\bm{U}_{\mathrm{E}})=\bm{0},$
(37)

for the variables

$\bm{U}_{\mathrm{E}}=\begin{bmatrix}E^{r}\\\ A_{\parallel}\\\
A_{\perp}\end{bmatrix},$ (38)

with radial flux

$\bm{F}_{\mathrm{E}}(\bm{U}_{\mathrm{E}})=\begin{bmatrix}0\\\
\displaystyle\frac{\alpha}{g_{rr}}\\\ \alpha A_{\parallel}\end{bmatrix},$ (39)

and sources

$\bm{S}_{\mathrm{E}}(\bm{U}_{\mathrm{E}})=\begin{bmatrix}\alpha(4\pi
j_{e}-KE^{r})\\\ \dfrac{\alpha
A_{\perp}}{g_{rr}}\left(\dfrac{2}{r}+2{D_{r\theta}}^{\theta}+{D_{rr}}^{r}\right)-\alpha
KA_{\parallel}\\\ \alpha\,g_{rr}\,E^{r}\end{bmatrix}.$ (40)

The Gauss constraint $D_{i}E^{i}=4\pi\rho_{e}$ reads

$\partial_{r}E^{r}+\left(\dfrac{2}{r}+2{D_{r\theta}}^{\theta}+{D_{rr}}^{r}\right)E^{r}-4\pi\rho_{e}=0,$
(41)

and the electromagnetic sources are explicitly given by

$\displaystyle 4\pi\rho_{e}$ $\displaystyle=$
$\displaystyle-2mK^{\prime}(X)\tilde{\Pi},$ (42) $\displaystyle 4\pi j_{e}$
$\displaystyle=$ $\displaystyle-2mK^{\prime}(X)g^{rr}\tilde{\Phi}.$ (43)

###  III.2 Boundary conditions

We consider approximate outgoing Sommerfeld conditions for the outer boundary
of our numerical domain, for which a detailed analysis of the characteristic
structure of the system is required. The electromagnetic sector admits the
eigenfields $\\{E^{r},\,A_{\pm}\\}$, where

$A_{\pm}=\frac{1}{2}\left(A_{\perp}\pm\sqrt{g_{rr}}A_{\parallel}\right),$ (44)

and with respective eigenvalues $\\{0,\,\pm\alpha/\sqrt{g_{rr}}\\}$. Inverting
Eq. (44), we get

$\displaystyle A_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{A_{+}-A_{-}}{\sqrt{g_{rr}}},$ (45) $\displaystyle
A_{\perp}$ $\displaystyle=$ $\displaystyle A_{+}+A_{-}.$ (46)

Following a similar idea implemented in Torres:2014fga , we apply outgoing
boundary conditions by requiring

$\displaystyle\partial_{t}A_{-}+v_{-}\partial_{r}A_{-}+\dfrac{v_{-}}{r}A_{-}=0,$
(47)

where $v_{-}=-\alpha/\sqrt{g_{rr}}$, and rewrite it in flux-conservative form.
For the rest of the system, we impose maximally dissipative boundary
conditions, for which all incoming fields at the outer boundary are
suppressed, thus damping spurious reflections as much as possible.

###  III.3 Initial data

For the initial-data set, we need to prescribe values for all the dynamical
fields in such a way that the Hamiltonian, momentum and Gauss constraints are
satisfied at $t=0$. For the metric fields, we set $\alpha(t=0,r)=1$, we write
the spatial metric in isotropic coordinates, yielding the line element
$dh^{2}=\psi^{4}(r)(dr^{2}+r^{2}d\Omega^{2})$, and solve for the conformal
factor $\psi(r)$ from the Hamiltonian constraint.

We consider a stationary initial configuration for which
$K_{r}{}^{r}=K_{\theta}{}^{\theta}=\Pi=0$, and set up the scalar field in two
ways:

  * •

ID Type I: we consider a pulse in $\Phi=\partial_{r}\phi$ given by

$\Phi(t=0,r)=A\,\exp\left[-\frac{(r-r_{c})^{2}}{\sigma^{2}}\right]\cos\left(\frac{\pi}{10}r\right);$
(48)

  * •

ID Type II: we set the pulse in the scalar field itself, with the form

$\phi(t=0,r)=A\,\exp\left[-\frac{(r-r_{c})^{2}}{\sigma^{2}}\right]\sin\left(\frac{r-r_{c}}{\sqrt{2}\,\sigma}\right),$
(49)

and thus $\Phi(t=0,r)=\partial_{r}\phi|_{t=0}$.

While the first configuration was used in MemoUV2022  to induce a Tricomi-type
breakdown in quadratic k-essence, the second one was implemented in
Bezares_Kdynamics2020  for simulating gravitational collapse when cubic
derivative self-interactions are taken into account.

For the electromagnetic variables, we choose the simplest possible initial
configuration, given by $E^{r}=A_{\parallel}=A_{\perp}=0$, for which the Gauss
constraint is automatically satisfied, as $\rho_{e}=0$ (see Eq. (42)). With
this choice for the initial fields, the momentum constraint is trivially
satisfied, while the Hamiltonian constraint yields a second-order elliptic
equation for the conformal factor $\psi(r)$, given by

$\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial\psi}{\partial r}\right)-2\pi GK(X)\psi^{5}=0.$
(50)

We solve this equation by shooting in $\psi(r=0)$, requiring regularity at the
origin and imposing a Robin boundary condition at infinity; i.e.,

$\lim_{r\to 0}{\frac{\partial\psi}{\partial
r}}=0\,,\quad\lim_{r\to\infty}{\left(\psi+r\frac{\partial\psi}{\partial
r}\right)}=1\,.$ (51)

Finally, for the “fixed” evolution, we set the initial value of the auxiliary
field $\Sigma$ to

$\Sigma(t=0,r)=K^{\prime}(X)|_{t=0}\,.$ (52)

Figure 1: Dynamics of the Stückelberg scalar and vector fields. Comparison
between the evolutions of the theory in the original (solid black curve) and
fixed versions, with parameters $\tau=0.5$ (dotted red) and $\tau=0.1$ (dashed
blue). Left column: Derivative of the scalar field (top); time and radial
components of the vector field (mid and bottom) at $t=30$. Right column: same
fields, but at $t=59.9$ in which the Tricomi-type breakdown holds. Although a
small discrepancy in the scalar profiles is observed between the
approximations given by the fixed solutions, the profiles of the components of
the vector remain almost unaltered. The parameters of the initial pulse for
the radial derivative of the scalar field are $A=10^{-4}$, $\sigma=0.4$ and
$r_{c}=60$.

##  IV Results

In this section, we show the results of our simulations. For computational
convenience, we employ code units such that $c=1$ and
$\ell_{m}=t_{m}=m^{-2}\kappa^{-1/2}$ for length and time, where $\kappa=8\pi
G$.

We take $\Lambda=100$ for all the simulations, $m=0.1$ in subsections IV.1 and
IV.2, and $m=1$ in subsection IV.3, always in the above units.

###  IV.1 Fixing a Tricomi-type breakdown

We present an example of the dynamical evolution of the Proca field in the
Stückelberg formulation. We take $\beta=1$ and $\gamma=0$ for the coupling
constants. For this choice, a Tricomi-type breakdown of the strong
hyperbolicity is expected (for a wide variety of initial data) from previous
studies in k-essence Bernard19 ; Bezares_Kdynamics2020 ; MemoUV2022 . In fact,
this is exactly the example case where pathologies were initially reported in
Proca_CloughPRL2022 ; Proca_Coates1_PRL2022 ; Zong-Gang22 . We adopt Type I
initial data, which corresponds to a localized Gaussian pulse for the radial
derivative of the scalar field, with amplitude $A=10^{-4}$, width $\sigma=0.4$
and centered at $r_{c}=60$. The time and radial components of the vector
potential, as well as the radial component of the electric field, are all set
initially to zero. In particular, the Gauss constraint (41) is initially
exactly satisfied. For the numerical simulations, we consider a radial domain
of length $L=300$, a spatial resolution $\Delta r=0.01$, and a CFL factor
$C=0.25$, so that $\Delta t=C\Delta r$.

Figure 1 shows a comparison of the dynamics yielded by the original theory
(solid black) and the fixed one, with two values of the timescale; $\tau=0.5$
(doted red) and $\tau=0.1$ (dashed blue). The profiles of the scalar and
vector fields are displayed at two different times. The first column
corresponds to $t=30$, and represents initial stages of the evolution; i.e.,
far from the hyperbolicity breakdown. As expected, the profiles corresponding
to the original and fixed theories perfectly agree, showing the robustness of
the fixing, at least during the initial evolution. The second column, on the
other hand, displays the dynamics of the fields at (approximately) the time in
which a Tricomi-type breakdown occurs, which is
$t_{\scriptsize{\mbox{tricomi}}}\sim 60$. A small discrepancy is observed in
the scalar profile, between the fixed and original evolutions, while no major
discrepancies are reported for the vector profiles. It is expected that the
fixing may give rise to discrepancies in the evolutions close to the time in
which the breakdown of hyperbolicity happens. Nevertheless, we can assess the
robustness of this approximation by noticing that, as $\tau$ decreases, the
fixed evolution becomes closer to the physical solution. This trend is shown
in the top panel of the second column.

Figure 2 exhibits the evolution of the minimum and maximum of the eigenvalues
$\lambda_{\pm}$ of the effective metric $\gamma^{\mu\nu}$ governing the
dynamics of the scalar field. In analogy with Figure 1, we compare the
evolution of the original theory with the fixed approximations, taking the
same values for the parameter $\tau$ as before. The dashed red curve
corresponds to the evolution of the original theory, which ceases at
approximately $t\sim t_{\scriptsize{\mbox{tricomi}}}$; i.e., where the
breakdown of hyperbolicity occurs. In particular, we observe that
$\lambda_{+}\to 0$, meaning a Tricomi-type behavior. In order to verify that
this is actually the case, the evolution is followed up with more time
resolution, by decreasing the CFL factor. On the other hand, the dotted blue
and dotted-dashed green curves correspond to the evolutions with the fixed
equations. Once again, an agreement on the behavior of the eigenvalues is
observed until a few time steps before the breakdown takes place, while a
better approximation is reached for the smaller value of the timescale (i.e.,
for $\tau=0.1$). The lapse evolution was also monitored, observing that its
value remains very close to unity until $t_{\scriptsize{\mbox{tricomi}}}$.
This shows that, even for a weak initial pulse, a Tricomi-like breakdown of
hyperbolicity can arise.

For times greater than $t_{\scriptsize{\mbox{tricomi}}}$, the evolution
continues if followed with the fixed equations (19). Now, the characteristic
structure of the (fixed) scalar evolution corresponds to the one for the wave
equation in GR, yielding the eigenvalues

$\lambda^{\tiny{\mbox{(GR)}}}_{+}=\frac{1}{g_{rr}}\,,\quad\lambda^{\tiny{\mbox{(GR)}}}_{-}=-\frac{1}{\alpha^{2}}\,.$
(53)

After a short transitional time once the Tricomi-type breakdown has taken
place in the original theory, the prediction for the minimum and maximum of
the eigenvalues provided by the fixed evolution approaches their expected
behaviours (dashed orange lines).

Figure 2: Eigenvalues of the effective metric. Minimum and maximum of the
eigenvalues of $\gamma^{\mu\nu}$ for the original (dashed red) and fixed
theories, with $\tau=0.1$ (dotted-dashed green) and $\tau=0.5$ (dotted blue).
A Tricomi-type breakdown is observed at approximately
$t_{\scriptsize{\mbox{tricomi}}}\sim 60$, for which $\lambda_{+}\to 0$ (upper
red). The dotted blue and dashed-dotted green curves approach the eigenvalues
of the fixed equations (dashed orange) after a short transient time.

Figure 3: Characteristic evolution with a cubic self-interaction. Left panel:
maximum and minimum of the eigenvalues of the effective metric as a function
of time. From top to bottom: $\max(\lambda_{+})$ (dashed blue);
$\min(\lambda_{+})$ (solid blue); $\max(\lambda_{-})$ (dashed green) and
$\min(\lambda_{-})$ (solid green). A comparison with the case with only a
quadratic self-interaction is shown in dashed brown for the minimum of the
corresponding eigenvalues $\lambda_{\pm}^{\text{q}}$. Right panel: maximum and
minimum of the characteristic velocities $V_{\pm}$ of the scalar evolution
system. From top to bottom: $\max(V_{+})$ (dashed orange); $\min(V_{+})$
(solid orange); $\max(V_{-})$ (solid blue) and $\min(V_{-})$ (dashed
turquoise).

###  IV.2 Evolution with a cubic self-interaction

We also study the evolution of the Proca field with a cubic self-interaction,
taking $\beta=\gamma=1$ for the couplings. We evolve an initial data set very
similar to the one considered in the previous section. Although we set the
same parameters for the Gaussian pulse, the elliptic equation (50) needed to
obtain the initial metric fields depends on $K(X)$, which now includes a cubic
term. The numerical domain, spatial resolution and CFL factor are taken as in
the simulation showed before. From Eq. (16), it is straightforward to see
that, for this choice of the coupling constants, there can not be a breakdown
of hyperbolicity of Tricomi-type, since $\det(\gamma^{\mu\nu})<0$ for all $X$,
in agreement with ShortPaper_Barausse2022 ; Bezares_Kdynamics2020 . We have
confirmed this fact numerically, i.e. a Tricomi-type failure, which can arise
from the action considered in Proca_CloughPRL2022 ; Proca_Coates1_PRL2022 ;
Zong-Gang22 , can be easily avoided by adding a cubic term. Moreover, we
numerically find that a Keldysh-type behavior does not occur either.

Figure 4: Constraint propagation. $L^{2}$ norm of the Gauss (dashed brown),
Hamiltonian (dashed-dotted salmon) and Z3 (solid orange) constraints as a
function of time, when considering a cubic self-interaction.

Figure 3 illustrates essential aspects of the characteristic structure of the
evolution system. The left panel represents the maximum and minimum values of
the eigenvalues of the effective metric $\gamma^{\mu\nu}$ as a function of
time (light-blue dashed and blue continuous lines for, respectively, the
maximum and minimum of $\lambda_{+}$; light-green dashed and green continuous
lines for, respectively, the maximum and minimum of $\lambda_{-}$). Also, and
for illustrative purposes only, a comparison with the case $\gamma=0$ is
included (brown dashed curves). The eigenvalues remain far from zero,
exhibiting a well-posed and stable evolution. The right panel displays the
maximum and minimum values of the characteristic velocities of the scalar
equation in (II.1), verifying that no Keldysh divergence of the characteristic
speeds was found during the whole evolution.

The time evolution of the norm of the Gauss (dashed brown), Hamiltonian
(dashed orange) and Z3 (dotted-dashed salmon) constraints is displayed in
Figure 4. In particular, they remain less than $10^{-14}$, up to $t\sim 60$,
where the pulse reaches $r=0$. For that value of time, the norm grows to
approximately $10^{-11}$, which is expected from the boundary/regularity
conditions considered. After that, although the Z3 constraint considerably
decays, the Hamiltonian and Gaussian constraints remain bounded around
$10^{-14}$.

###  IV.3 Gravitational collapse with a cubic self-interaction

We report here an example of the dynamical evolution of the Proca field with a
cubic self-interaction leading to gravitational collapse, particularly in the
formation of a black hole. For the coupling constants we take $\beta=0$ and
$\gamma=1$, afresh preventing the determinant (16) from ever becoming zero.

We considered the Type II initial data given by Eq. (49), with amplitude
$A=0.093$, width $\sigma=0.942$ and center $r_{c}=55$. This configuration
yielded an Arnowitt-Deser-Misner (ADM) mass $M_{\tiny{\mbox{ADM}}}\simeq
1.590$. At the first stages of the simulation, the scalar pulse splits into
two modes: one moving towards the origin with an amplitude that increases,
approximately, as $1/r$ (as expected) and the other towards the outer
boundary. The spatial resolution, the CFL factor and the size of the numerical
domain were all taken to be the same as in the cases presented earlier.

Figure 5: Gravitational collapse with a cubic self-interaction. Profiles of
the lapse function for different values of time. Although it is initially set
to one, the “$1+\log$” slicing condition forces the lapse to vanish close to
the origin. The collapse occurs around $t_{\tiny{\mbox{AH}}}\sim 68.5$, where
the first apparent horizon is formed, and the lapse reaches zero soon after
$t_{\tiny{\mbox{AH}}}$.

Figure 6: Dynamics of the collapse. Top: Snapshots of the scalar field
(multiplied by $r$) in the region close to the apparent horizon, for different
times after the collapse takes place. The location of the apparent horizon at
each time is shown with a vertical line. Bottom: minimum and maximum values of
the characteristic velocities $V_{\pm}$ of the scalar equation as a function
of time. From top to bottom: $\max(V_{+})$ (dashed blue); $\min(V_{+})$ (solid
orange); $\max(V_{-})$ (solid green) and $\min(V_{-})$ (dashed red). The
maximum of the determinant of the effective metric $\gamma^{\mu\nu}$ is
represented in dotted black, and the time at which the first apparent horizon
is formed is meant by the dashed grey vertical line.

By letting the system evolve until $t_{\tiny{\mbox{final}}}\sim 100$, the
scalar field collapses and a black hole forms. This is signalled by the
formation of an apparent horizon (i.e., the outermost trapped surface) at
$t\sim 68.5$, whose position continuously increases in our coordinates. To
keep track of the dynamics of this surface, we considered the largest value of
$r$ for which the expansion of the outgoing null-ray congruence vanishes,
physically indicating a confinement of the latter. In spherical symmetry, this
condition reads Thornburg:2006zb

$D_{r\theta}{}^{\theta}+\frac{1}{r}-\sqrt{g_{rr}}K_{\theta}{}^{\theta}=0,$
(54)

where the function $D_{r\theta}{}^{\theta}$ is defined in Eq. (24). The
results for the lapse evolution are presented in Figure 5, showing profiles at
different times. Starting from $\alpha=1$, the lapse vanishes at approximately
$t\sim 69$, and the front propagation speed grows in time due to our gauge
choice. Also, in this particular gauge the singularity is avoided, allowing
one to proceed with the simulations after the appearance of the first apparent
horizon. The areal radius $r_{\tiny{\mbox{A}}}=r\,\sqrt{g_{\theta\theta}}$ is
also computed during the evolution. We estimate a black-hole area of
$A_{\tiny{\mbox{BH}}}\simeq 13.854$, and a mass of $M_{\tiny{\mbox{BH}}}\simeq
0.525$. The latter represents approximately $33\%$ of the total ADM mass of
the system. This apparent mass loss is consistent with the initial
configuration: just one half of the total scalar energy of the initial pulse
takes part of the physical process leading to collapse, while the other half
gets rapidly radiated away.

The top panel of Figure 6 shows profiles of the Stückelberg field close to the
horizon for different time values after the black hole has formed. The dashed
vertical lines show the location of the apparent horizon at these times, which
are meant by different colors. In the bottom panel, we display the maximum and
minimum values of the characteristic speeds of the scalar equation, as well as
the maximum of the determinant of the effective metric (dotted black). The
latter remains always bounded and different from zero, assuring no breakdown
of the Cauchy problem during evolution. After the black hole has formed, our
gauge choice “freezes-out” the region within the apparent horizon, thus not
giving any physical insights on the dynamics at the interior. Alternative and
more sophisticated numerical approaches to horizon formation would also be
possible; e.g., excision of the interior region from the numerical domain (at
the cost of putting boundary conditions at the horizon, which moves during the
evolution), or another gauge choice.

##  V Conclusions

In this paper, we conducted numerical simulations of self-interacting massive
vector fields coupled to General Relativity, in spherical symmetry. We
investigated the hyperbolicity of the corresponding Cauchy problem, displaying
numerical evidences that the breakdowns faced by the theory are analogous to
those occurring in scalar-tensor theories with first-order derivative self-
interactions (usually known as k-essence theories).

In particular, we gave an explicit example leading to a Tricomi-type
breakdown, when only quadratic self-interactions are accounted for. We showed
that this pathology can be avoided by suitably deforming the evolution
equations by means of a “fixing-the-equations” approach. We also explored the
dynamics of the vector field when cubic self-interactions are taken into
account. As expected, no Tricomi-type breakdown of the Cauchy problem
occurred. For the latter case, we also explored gravitational collapse. We
found suitable initial data giving rise to well-posed and stable evolutions
towards a black hole final state.

To conclude, we stress the advantage of our approach for performing numerical
simulations of massive vector fields, as it singles out the Stückelberg mode,
taming derivative self-interactions. This procedure can also be applied to
numerical simulations in massive gravity (i.e., a massive spin-2 field
breaking diffeomorphism invariance, which can be restored by a Stückelberg
transformation), for example in the context of BH superradiance Brito:2013wya
; brito2020superradiance . We leave these explorations for future work.

## Acknowledgements

We thank Ramiro Cayuso, Luis Lehner and Marc Schneider for discussions
throughout this work. MR acknowledges hospitality from the Perimeter Institute
for Theoretical Physics, where part of this work was carried out. MR and EB
acknowledge support from the European Union’s H2020 ERC Consolidator Grant
“GRavity from Astrophysical to Microscopic Scales” (Grant No. GRAMS-815673),
the PRIN 2022 grant “GUVIRP - Gravity tests in the UltraViolet and InfraRed
with Pulsar timing”, and the EU Horizon 2020 Research and Innovation Programme
under the Marie Sklodowska-Curie Grant Agreement No. 101007855. MC is funded
by the European Union under the Horizon Europe’s Marie Sklodowska-Curie
project 101065440.

## References

  * (1) B. P. Abbott et al., “Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs,” Phys. Rev. X, vol. 9, p. 031040, Sep 2019.
  * (2) B. P. Abbott et al., “Gw170817: Observation of gravitational waves from a binary neutron star inspiral,” Phys. Rev. Lett., vol. 119, p. 161101, Oct 2017.
  * (3) B. P. Abbott et al., “Gw190521: A binary black hole merger with a total mass of $150\text{ }\text{ }{M}_{\bigodot}$,” Phys. Rev. Lett., vol. 125, p. 101102, Sep 2020.
  * (4) B. P. Abbott et al., “Gw170817: Measurements of neutron star radii and equation of state,” Phys. Rev. Lett., vol. 121, p. 161101, Oct 2018.
  * (5) V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel., vol. 22, no. 1, p. 4, 2019.
  * (6) J. Baker et al., “The Laser Interferometer Space Antenna: Unveiling the Millihertz Gravitational Wave Sky,” 7 2019.
  * (7) T. B. Littenberg, K. Breivik, W. R. Brown, M. Eracleous, J. J. Hermes, K. Holley-Bockelmann, K. Kremer, T. Kupfer, and S. L. Larson, “Astro2020 Decadal Science White Paper: Gravitational Wave Survey of Galactic Ultra Compact Binaries,” 3 2019.
  * (8) M. Branchesi et al., “Science with the Einstein Telescope: a comparison of different designs,” JCAP, vol. 07, p. 068, 2023.
  * (9) W. Israel, “Event horizons in static vacuum space-times,” Phys. Rev., vol. 164, pp. 1776–1779, 1967.
  * (10) S. W. Hawking, “Black holes in general relativity,” Commun. Math. Phys., vol. 25, pp. 152–166, 1972.
  * (11) D. C. Robinson, “Uniqueness of the Kerr black hole,” Phys. Rev. Lett., vol. 34, pp. 905–906, 1975.
  * (12) M. Bezares and N. Sanchis-Gual, “Exotic compact objects: a recent numerical-relativity perspective,” 6 2024.
  * (13) E. Berti and V. Cardoso, “Supermassive black holes or boson stars? hair counting with gravitational wave detectors,” International Journal of Modern Physics D, vol. 15, no. 12, pp. 2209–2216, 2006.
  * (14) V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, “Gravitational-wave signatures of exotic compact objects and of quantum corrections at the horizon scale,” Phys. Rev. D, vol. 94, p. 084031, Oct 2016.
  * (15) P. V. P. Cunha, E. Berti, and C. A. R. Herdeiro, “Light-ring stability for ultracompact objects,” Phys. Rev. Lett., vol. 119, p. 251102, Dec 2017.
  * (16) E. Maggio, P. Pani, and V. Ferrari, “Exotic compact objects and how to quench their ergoregion instability,” Phys. Rev. D, vol. 96, p. 104047, Nov 2017.
  * (17) G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, “Anisotropic stars as ultracompact objects in general relativity,” Phys. Rev. D, vol. 99, p. 104072, May 2019.
  * (18) M. Kesden, J. Gair, and M. Kamionkowski, “Gravitational-wave signature of an inspiral into a supermassive horizonless object,” Phys. Rev. D, vol. 71, p. 044015, Feb 2005.
  * (19) P. Jetzer, “Boson stars,” Physics Reports, vol. 220, no. 4, pp. 163–227, 1992.
  * (20) S. L. Liebling and C. Palenzuela, “Dynamical boson stars,” Living Rev. Rel., vol. 26, no. 1, p. 1, 2023.
  * (21) F. E. Schunck and E. W. Mielke, “General relativistic boson stars,” Class. Quant. Grav., vol. 20, pp. R301–R356, 2003.
  * (22) F. S. Guzmán and J. M. Rueda-Becerril, “Spherical boson stars as black hole mimickers,” Phys. Rev. D, vol. 80, p. 084023, Oct 2009.
  * (23) M. Colpi, S. L. Shapiro, and I. Wasserman, “Boson stars: Gravitational equilibria of self-interacting scalar fields,” Phys. Rev. Lett., vol. 57, pp. 2485–2488, Nov 1986.
  * (24) D. J. Kaup, “Klein-gordon geon,” Phys. Rev., vol. 172, pp. 1331–1342, Aug 1968.
  * (25) M. Gleiser, “Stability of boson stars,” Phys. Rev. D, vol. 38, pp. 2376–2385, Oct 1988.
  * (26) M. Gleiser, “Erratum: Stability of boson stars [phys. rev. d 38, 2376 (1988)],” Phys. Rev. D, vol. 39, pp. 1257–1257, Feb 1989.
  * (27) T.-D. Lee and Y. Pang, “Nontopological solitons,” Physics Reports, vol. 221, no. 5-6, pp. 251–350, 1992.
  * (28) J.-W. Lee and S. Lim, “Minimum mass of galaxies from BEC or scalar field dark matter,” JCAP, vol. 01, p. 007, 2010.
  * (29) J.-W. Lee, “Is dark matter a BEC or scalar field?,” J. Korean Phys. Soc., vol. 54, p. 2622, 2009.
  * (30) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, “Proca stars: Gravitating Bose–Einstein condensates of massive spin 1 particles,” Phys. Lett. B, vol. 752, pp. 291–295, 2016.
  * (31) C. Herdeiro, E. Radu, and E. dos Santos Costa Filho, “Proca-Higgs balls and stars in a UV completion for Proca self-interactions,” JCAP, vol. 05, p. 022, 2023.
  * (32) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, “Nonlinear dynamics of spinning bosonic stars: Formation and stability,” Phys. Rev. Lett., vol. 123, p. 221101, Nov 2019.
  * (33) C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual, N. M. Santos, and E. dos Santos Costa Filho, “The non-spherical ground state of Proca stars,” Phys. Lett. B, vol. 852, p. 138595, 2024.
  * (34) I. Sengo, P. V. P. Cunha, C. A. R. Herdeiro, and E. Radu, “Kerr black holes with synchronised Proca hair: lensing, shadows and EHT constraints,” JCAP, vol. 01, p. 047, 2023.
  * (35) I. Sengo, P. V. P. Cunha, C. A. R. Herdeiro, and E. Radu, “The imitation game reloaded: effective shadows of dynamically robust spinning Proca stars,” JCAP, vol. 05, p. 054, 2024.
  * (36) C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual, N. M. Santos, and E. dos Santos Costa Filho, “The non-spherical ground state of Proca stars,” Phys. Lett. B, vol. 852, p. 138595, 2024.
  * (37) F. Di Giovanni, N. Sanchis-Gual, C. A. R. Herdeiro, and J. A. Font, “Dynamical formation of proca stars and quasistationary solitonic objects,” Phys. Rev. D, vol. 98, p. 064044, Sep 2018.
  * (38) T. Shen, M. Zhou, C. Bambi, C. A. R. Herdeiro, and E. Radu, “Iron K$\alpha$ line of Proca stars,” JCAP, vol. 08, p. 014, 2017.
  * (39) W. E. East and F. Pretorius, “Superradiant instability and backreaction of massive vector fields around kerr black holes,” Phys. Rev. Lett., vol. 119, p. 041101, Jul 2017.
  * (40) W. E. East, “Superradiant instability of massive vector fields around spinning black holes in the relativistic regime,” Phys. Rev. D, vol. 96, p. 024004, Jul 2017.
  * (41) S. R. Dolan, “Instability of the proca field on kerr spacetime,” Phys. Rev. D, vol. 98, p. 104006, Nov 2018.
  * (42) K. Clough, T. Helfer, H. Witek, and E. Berti, “Ghost instabilities in self-interacting vector fields: The problem with proca fields,” Phys. Rev. Lett., vol. 129, p. 151102, Oct 2022.
  * (43) A. Coates and F. M. Ramazanoğlu, “Intrinsic pathology of self-interacting vector fields,” Phys. Rev. Lett., vol. 129, p. 151103, Oct 2022.
  * (44) A. O. Barvinsky and Y. V. Gryzov, “Generalized causality condition in quantum field theory with torsion,” Sov. Phys. J., vol. 32, pp. 383–386, 1989.
  * (45) Z.-G. Mou and H.-Y. Zhang, “Singularity problem for interacting massive vectors,” Phys. Rev. Lett., vol. 129, p. 151101, Oct 2022.
  * (46) E. Barausse, M. Bezares, M. Crisostomi, and G. Lara, “The well-posedness of the Cauchy problem for self-interacting vector fields,” JCAP, vol. 11, p. 050, 2022.
  * (47) A. Coates and F. M. Ramazanoğlu, “Coordinate singularities of self-interacting vector field theories,” Phys. Rev. Lett., vol. 130, p. 021401, Jan 2023.
  * (48) A. Coates and F. M. Ramazanoğlu, “Pervasiveness of the breakdown of self-interacting vector field theories,” Phys. Rev. D, vol. 107, p. 104036, May 2023.
  * (49) K. I. Ünlütürk, A. Coates, and F. M. Ramazanoğlu, “Loss of hyperbolicity and tachyons in generalized proca theories,” Phys. Rev. D, vol. 108, p. 044022, Aug 2023.
  * (50) H. O. Silva, A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou, “Ghost of vector fields in compact stars,” Phys. Rev. D, vol. 105, p. 024046, Jan 2022.
  * (51) K. Aoki and M. Minamitsuji, “Resolving the pathologies of self-interacting proca fields: A case study of proca stars,” Phys. Rev. D, vol. 106, p. 084022, Oct 2022.
  * (52) L. Bernard, L. Lehner, and R. Luna, “Challenges to global solutions in horndeski’s theory,” Phys. Rev. D, vol. 100, p. 024011, Jul 2019.
  * (53) M. Bezares, M. Crisostomi, C. Palenzuela, and E. Barausse, “K-dynamics: well-posed 1+1 evolutions in K-essence,” JCAP, vol. 03, p. 072, 2021.
  * (54) L. ter Haar, M. Bezares, M. Crisostomi, E. Barausse, and C. Palenzuela, “Dynamics of screening in modified gravity,” Phys. Rev. Lett., vol. 126, p. 091102, Mar 2021.
  * (55) M. Bezares, L. ter Haar, M. Crisostomi, E. Barausse, and C. Palenzuela, “Kinetic screening in nonlinear stellar oscillations and gravitational collapse,” Phys. Rev. D, vol. 104, p. 044022, Aug 2021.
  * (56) M. Bezares, R. Aguilera-Miret, L. ter Haar, M. Crisostomi, C. Palenzuela, and E. Barausse, “No Evidence of Kinetic Screening in Simulations of Merging Binary Neutron Stars beyond General Relativity,” Phys. Rev. Lett., vol. 128, no. 9, p. 091103, 2022.
  * (57) P. Figueras and T. França, “Gravitational Collapse in Cubic Horndeski Theories,” Class. Quant. Grav., vol. 37, no. 22, p. 225009, 2020.
  * (58) J. Cayuso, N. Ortiz, and L. Lehner, “Fixing extensions to general relativity in the nonlinear regime,” Phys. Rev. D, vol. 96, p. 084043, Oct 2017.
  * (59) N. Franchini, M. Bezares, E. Barausse, and L. Lehner, “Fixing the dynamical evolution in scalar-Gauss-Bonnet gravity,” Phys. Rev. D, vol. 106, no. 6, p. 064061, 2022.
  * (60) R. Cayuso and L. Lehner, “Nonlinear, noniterative treatment of eft-motivated gravity,” Phys. Rev. D, vol. 102, p. 084008, Oct 2020.
  * (61) G. Lara, M. Bezares, and E. Barausse, “Uv completions, fixing the equations, and nonlinearities in $k$-essence,” Phys. Rev. D, vol. 105, p. 064058, Mar 2022.
  * (62) G. Lara, H. P. Pfeiffer, N. A. Wittek, N. L. Vu, K. C. Nelli, A. Carpenter, G. Lovelace, M. A. Scheel, and W. Throwe, “Scalarization of isolated black holes in scalar Gauss-Bonnet theory in the fixing-the-equations approach,” 3 2024.
  * (63) M. Corman, L. Lehner, W. E. East, and G. Dideron, “Nonlinear studies of modifications to general relativity: Comparing different approaches,” 5 2024.
  * (64) A. Coates and F. M. Ramazanoğlu, “Treatments and placebos for the pathologies of effective field theories,” Phys. Rev. D, vol. 108, p. L101501, Nov 2023.
  * (65) J. Hadamard, “Théorie des équations aux dérivées partielles linéaires hyperboliques et du problème de cauchy,” Acta Math., vol. 31, pp. 333–380, 1908.
  * (66) K. Heinz‐Otto, “Initial boundary value problems for hyperbolic systems,” Communications on Pure and Applied Mathematics, vol. 23, no. 3, pp. 277–298, 1970.
  * (67) R. Geroch, “Partial differential equations of physics,” General Relativity, Aberdeen, Scotland, pp. 19–60, 1996.
  * (68) F. K. O., “Symmetric hyperbolic linear differential equations,” Communications on Pure and Applied Mathematics, vol. 7, no. 2, pp. 345–392, 1954.
  * (69) K. O. Friedrichs and P. D. Lax, “Systems of conservation equations with a convex extension,” Proceedings of the National Academy of Sciences, vol. 68, no. 8, pp. 1686–1688, 1971.
  * (70) J. L. Ripley and F. Pretorius, “Hyperbolicity in Spherical Gravitational Collapse in a Horndeski Theory,” Phys. Rev. D, vol. 99, no. 8, p. 084014, 2019.
  * (71) J. L. Ripley, “Numerical relativity for Horndeski gravity,” Int. J. Mod. Phys. D, vol. 31, no. 13, p. 2230017, 2022.
  * (72) G. Allwright and L. Lehner, “Towards the nonlinear regime in extensions to GR: assessing possible options,” Class. Quant. Grav., vol. 36, no. 8, p. 084001, 2019.
  * (73) R. Cayuso and L. Lehner, “Nonlinear, noniterative treatment of eft-motivated gravity,” Phys. Rev. D, vol. 102, p. 084008, Oct 2020.
  * (74) W. Israel and J. Stewart, “Thermodynamics of nonstationary and transient effects in a relativistic gas,” Physics Letters A, vol. 58, no. 4, pp. 213–215, 1976.
  * (75) W. Israel, “Nonstationary irreversible thermodynamics: A causal relativistic theory,” Annals of Physics, vol. 100, no. 1, pp. 310–331, 1976.
  * (76) W. Israel and J. Stewart, “Transient relativistic thermodynamics and kinetic theory,” Annals of Physics, vol. 118, no. 2, pp. 341–372, 1979.
  * (77) L. Lehner personal communication, May 2024.
  * (78) D. Alic, C. Bona, C. Bona-Casas, and J. Masso, “Efficient implementation of finite volume methods in Numerical Relativity,” Phys. Rev., vol. D76, p. 104007, 2007.
  * (79) A. Bernal, J. Barranco, D. Alic, and C. Palenzuela, “Multi-state Boson Stars,” Phys. Rev. D, vol. 81, p. 044031, 2010.
  * (80) S. Valdez-Alvarado, C. Palenzuela, D. Alic, and L. A. Ureña López, “Dynamical evolution of fermion-boson stars,” Phys. Rev. D, vol. 87, no. 8, p. 084040, 2013.
  * (81) G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, “Anisotropic stars as ultracompact objects in General Relativity,” Phys. Rev., vol. D99, no. 10, p. 104072, 2019.
  * (82) G. N. Felder, L. Kofman, and A. Starobinsky, “Caustics in tachyon matter and other Born-Infeld scalars,” JHEP, vol. 09, p. 026, 2002.
  * (83) H. S. Reall, N. Tanahashi, and B. Way, “Shock Formation in Lovelock Theories,” Phys. Rev., vol. D91, no. 4, p. 044013, 2015.
  * (84) E. Babichev and S. Ramazanov, “Caustic free completion of pressureless perfect fluid and k-essence,” JHEP, vol. 08, p. 040, 2017.
  * (85) N. Tanahashi and S. Ohashi, “Wave propagation and shock formation in the most general scalar–tensor theories,” Class. Quant. Grav., vol. 34, no. 21, p. 215003, 2017.
  * (86) C. Bona, T. Ledvinka, and C. Palenzuela, “3+1 covariant suite of numerical relativity evolution systems,” Phys. Rev. D, vol. 66, p. 084013, Oct 2002.
  * (87) C. Bona and J. Massó, “Einstein’s evolution equations as a system of balance laws,” Phys. Rev. D, vol. 40, pp. 1022–1026, Aug 1989.
  * (88) C. Bona, J. Massó, E. Seidel, and J. Stela, “New formalism for numerical relativity,” Phys. Rev. Lett., vol. 75, pp. 600–603, Jul 1995.
  * (89) J. M. Torres and M. Alcubierre, “Gravitational collapse of charged scalar fields,” Gen. Rel. Grav., vol. 46, p. 1773, 2014.
  * (90) J. Thornburg, “Event and apparent horizon finders for 3+1 numerical relativity,” Living Rev. Rel., vol. 10, p. 3, 2007.
  * (91) R. Brito, V. Cardoso, and P. Pani, “Massive spin-2 fields on black hole spacetimes: Instability of the Schwarzschild and Kerr solutions and bounds on the graviton mass,” Phys. Rev. D, vol. 88, no. 2, p. 023514, 2013.
  * (92) R. Brito, V. Cardoso, P. Pani, et al., Superradiance New Frontiers in Black Hole Physics.  Springer, 2020.