# Instability of vectorized stars

Ekrem S. Demirboğa<EMAIL_ADDRESS>Andrew Coates<EMAIL_ADDRESS>Fethi M. Ramazanoğlu<EMAIL_ADDRESS>Department of Physics, Koç
University, Rumelifeneri Yolu, 34450 Sariyer, Istanbul, Turkey

###### Abstract

In recent papers it has been shown that a large class of vectorization
mechanisms in gravity, which involve the vector fields becoming apparently
tachyonic in some regime, are actually dominated by ghosts and non-
perturbative behavior. Despite this, vectorized compact object solutions have
previously been found, which raises the question of how, and if, the newly
discovered ghosts are quenched in these cases. Here we develop the tools to
study the perturbations of vectorized compact objects, and demonstrate that
they suffer from ghosts and gradient instabilities as well. Thus, these
vectorized objects do not represent the stable end point of a quenched
instability unlike their scalarized counterparts in the spontaneous
scalarization literature.

##  I Introduction

General Relativity (GR) is incredibly successful from a phenomenological point
of view, having accounted for previously unexplained observations and made
many incredible predictions which have since been verified. Perhaps the most
spectacular of these is the confirmation of GR’s prediction of the existence
of gravitational waves. It matches observations of binary pulsars so well as
to provide an extremely stringent constraint on massless (and sufficiently
light) scalar fields [1]. It is also so far consistent with all of direct
gravitational wave observations performed by the LIGO-Virgo-KAGRA
collaboration (e.g [2]).

Given the success of GR, especially in the solar system [3], the most
interesting theories come with some mechanism to mimic GR in the solar system,
while allowing for large deviations elsewhere. Such mechanisms are termed
“screening mechanisms” and the prototypical examples are the Damour Esposito-
Farèse (DEF) model of spontaneous scalarization [4, 5] and the Vainshtein
mechanism [6] (see Ref. [7] for a review). We are interested in the former in
this study, which has an action of the form

$\displaystyle S$ $\displaystyle=\frac{1}{16\pi}\int{\rm
d}^{4}x\sqrt{-g}\,\left[R-2\left(\partial_{\mu}\phi\right)\left(\partial^{\mu}\phi\right)-2V(\phi)\right]$
$\displaystyle\quad+S_{\rm m}\left[f_{\rm
m};\,{\Omega}_{\phi}^{2}(\phi)\,g_{\mu\nu}\right]\ .$ (1)

Here $V(\phi)$ is a potential for the scalar field, $f_{\rm m}$ stands for any
matter fields and $\Omega_{\phi}^{2}(\phi)$ is a conformal factor through
which the scalar field and matter are coupled. Given some conditions on
$\Omega_{\phi}$ and $V$ one can arrange for the theory to be solved by GR
solutions with constant $\phi=\phi_{0}$. Perturbations of the scalar field
around such a solution have an effective mass squared of the form

$m_{\rm eff}^{2}=m_{\phi}^{2}-4\pi\beta T,$ (2)

where $m_{\phi}^{2}=\left.V^{\prime\prime}\right|_{\phi_{0}}$,
$\beta=\nicefrac{{{\rm d}^{2}\log\Omega_{\phi}^{2}}}{{{\rm
d}\phi^{2}}}|_{\phi_{0}}$ and $T$ is the trace of the stress energy tensor.
Choosing $m_{\phi}$ and $\beta$ correctly one can then arrange that the solar
system is stably given by GR while (say) neutron stars become unstable to
scalar perturbations. An instability due to a wrong sign of the effective
square mass is termed “tachyonic”. When this effective square mass is constant
throughout spacetime, it is sufficient that it be negative to trigger an
instability, whereas in this case it depends on the matter configuration.

If $m_{\rm eff}^{2}<0$ in some region of characteristic size $R$, a heuristic
condition for the instability is that $1/R^{2}+m_{\rm eff}^{2}\lesssim 0$ (so
that some mode which can probe this length scale is unstable within it) [8].
When $m_{\phi}^{2}$ is small, one can see that the controlling parameter for
the onset of the instability will be the compactness of the object. When the
mass is sufficiently small, binary pulsar tests heavily constrain the model
[9], and for very large masses the mechanism is never activated for
astrophysical objects [8]. However, there remains an interesting range of
parameters for light scalars in various astrophysical scenarios [10].

There is no explicit role of the scalar nature of the field in the spontaneous
scalarization mechanism we have presented, hence, there have been efforts to
generalize it to other fields, the most commonly studied case being vectors
[11, 12, 13, 14]. The common theme of these theories is that the vector field
spontaneously grows in the vicinity of compact objects similarly to the scalar
of the DEF model, and the vectorized objects generically feature large
deviations from GR.

Our analysis will focus on the direct generalization of the DEF model, Ref.
[11], though some of our results will cover any model of vectorization that
have perturbation equations of the form

$\nabla_{\mu}F^{\mu\nu}=\Xi^{\nu}{}_{\mu}A^{\mu}\ ,$ (3)

such as Refs. [12, 13, 14].

It was recently shown that theories of this type possess ghost instabilities
if they contain what one might term a “naive tachyon” [15, 16].111See Ref.
[17] for a similar earlier work on cosmology. If ${\Xi}$ were to be a good
mass squared tensor, it should behave, in some key ways, like
$\mu^{2}\delta^{\mu}{}_{\nu}$, the mass term of the Proca field. A “naive
tachyon” would then be when it does not meet these requirements, for example,
if ${\Xi}$ has negative eigenvalues as a matrix, one would expect some
tachyonic behaviour. As we will demonstrate later, ${\Xi^{\mu\nu}}$ plays an
even more prominent role than that of a mass squared tensor, being the inverse
of the effective metric for the scalar mode of the vector field, and so
requirements of causality will also play a role in deciding what forms of
${\Xi}$ allows healthy dynamics of the vector (see Ref. [18]).

As a concrete example, the spontaneous vectorization theory of Ref. [11] is
given by the action222Also see Ref. [19] for a similar cosmological theory.

$\displaystyle S$ $\displaystyle=\frac{1}{16\pi}\int{\rm
d}^{4}x\sqrt{-g}\,\left(R-F_{\mu\nu}F^{\mu\nu}-2m_{\rm
V}^{2}A_{\mu}A^{\mu}\right)$ $\displaystyle\quad+S_{\rm m}\left[f_{\rm
m};\,{\Omega_{\rm V}}^{2}(A_{\mu}A^{\mu})\,g_{\mu\nu}\right],$ (4)

where $A_{\mu}$ is a vector field with bare mass $m_{\rm V}$,
$F_{\mu\nu}=\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}$, and as before $f_{\rm
m}$ is a generic label for any matter field. The vector field equation of
motion is given by

$\nabla_{\mu}F^{\mu\nu}=\left(m_{\rm V}^{2}-8\pi\Lambda\,\Omega_{\rm
V}^{4}\tilde{T}\right)A^{\nu}\equiv\hat{z}m_{\rm V}^{2}A^{\nu}\ ,$ (5)

where

$\hat{z}=1-\frac{8\pi}{m_{\rm V}^{2}}\frac{d\ln\Omega_{\rm
V}}{d(A_{\mu}A^{\mu})}\Omega_{\rm V}^{4}\tilde{T}\ .$ (6)

$\tilde{T}$ is the trace of the stress energy tensor with respect to the
metric $\tilde{g}_{\mu\nu}=\Omega_{\rm V}^{2}g_{\mu\nu}$. In the rest of the
paper we will use

$\Omega_{\rm V}=e^{\beta A_{\mu}A^{\mu}/2}\ ,$ (7)

following the literature.

When $\hat{z}<0$ we have a naive tachyon as described above, which led to
expectations that the theory would have vectorized neutron star solutions
similar to the scalarized solutions of the DEF model [11]. However, in Ref.
[15] it was shown that action (4) also features ghost instabilities for a
neutron star that is a solution of GR when $\hat{z}$ is negative. Furthermore,
$\hat{z}$ also changes sign in many cases, which leads to other serious
problems due to divergent terms of the form $\hat{z}^{-1}$.

Despite the problems around a GR background, neutron star solutions with non-
trivial vector configurations have been numerically calculated [11]. If these
vectorized solutions are stable, one could hope that the problematic GR
backgrounds are dynamically avoided and these vectorized configurations
develop without issue. Here, we will show that this is not the case. That is,
arbitrarily small perturbations from vectorized solutions calculated in Ref.
[11] grow in an unbounded manner, meaning they are unstable. We also explore
parts of the $(\beta,m_{\rm V})$ parameter space that have not been
investigated so far, and show that there is no stable spherically symmetric
vectorized star in any sector we studied. This is in stark contrast to the
case of the DEF model where typical scalarized solutions are known to be
stable in various models [20, 21, 22, 23, 8, 24].

Despite these negative results, we shall see that the naive picture of a
tachyonic degree of freedom can apply to the axial part of the vector field,
and so, in the regions of parameter space where this holds, one should likely
be looking at axisymmetric stars for healthy end states, if any exist.

In Sec. II, we develop the perturbative framework to analyze the instabilities
of vectorized stars by generalizing the approach of Ref. [15]. In Sec. III, we
numerically solve for vectorized stars in the theory of action (4), and show
that all solutions are unstable to linearized pertubations. Not all solutions
carry all the different forms of instabilities we study, but they all carry at
least one. In Sec. IV, we will discuss our results, and argue how our methods
can be applied to other vectorization theories and other extensions of the DEF
model.

We use the geometric units $G=c=1$ and the metric signature $(-1,1,1,1)$. We
use Einstein summation conventions for greek indices $\mu,\nu,\dots$, but not
for the latin indices for spatial variables such as $i$, unless specifically
stated otherwise.

##  II Instabilities of vectorized stars

###  II.1 Spherically symmetric vectorized stars

The metric of a static and spherically symmetric vectorized star can be
represented by the ansatz

${\rm d}s^{2}=-e^{\nu(r)}{\rm d}t^{2}+\frac{1}{1-\frac{2\mu(r)}{r}}{\rm
d}r^{2}+r^{2}({\rm d}\theta^{2}+\sin^{2}\theta\,{\rm d}\phi)\ .$ (8)

The vector field also has the simple form

$\displaystyle A^{\mu}=\left(A^{0}(r),0,0,0\right)=A^{0}(r)\partial_{0}\ .$
(9)

That is, only the time component of the vector field is nonvanishing in a
static and spherically symmetric spacetime. This was “assumed” in Ref. [11]
for ease of calculation. Here, we examine the most general case, and show that
there is no need for an assumption.

The fact that spherical symmetry implies that the angular components of the
vector field, $A^{\theta,\phi}$, vanish is not a trivial one. For example,
minimally coupled massless vector fields can have nonzero angular components
under spherical symmetry such as the electromagnetic field vector $V$ of a
Reissner-Nordstrom black hole with magnetic charge. This is possible due to
the fact that the stress-energy tensor can still be spherically symmetric in
this case, even though $V^{\phi}\neq 0$ and $F^{\theta\phi}\neq 0$. However,
this is not the case for a massive field.

The metric equation of motion arising from the action (4) takes the following
form outside the star

$\displaystyle
R_{\mu\nu}=2F_{\mu\rho}F_{\nu}{}^{\rho}-\frac{1}{2}F_{\rho\sigma}F^{\rho\sigma}g_{\mu\nu}+2m_{\rm
V}^{2}A_{\mu}A_{\nu}\ .$ (10)

The line element (8) has direct implications for the Ricci and Einstein
tensors. Firstly, $R_{\mu\nu}$ has to be diagonal, hence only one component of
$A_{\mu}$ can be nonvanishing due to the $m_{\rm V}^{2}$ term. Secondly, the
surviving component cannot be $A^{\theta}$ or $A^{\phi}$, since Eq. (8) also
implies that $R_{\phi\phi}=\sin^{2}\theta R_{\theta\theta}$, which cannot be
satisfied due to the same term.

Finally, vanishing of $A^{r}$ is a direct consequence of the $\alpha=r$ case
of Eq. (5)

$\displaystyle A^{r}=\frac{\partial_{\rho}(\sqrt{|g|}F^{\rho
r})}{\sqrt{|g|}(-4\pi\Omega_{\rm V}^{4}\beta\tilde{T}+m^{2})}=0\ ,$ (11)

where we used the partial derivative formula for the divergence of an
antisymmetric tensor, and the fact that the only nonzero component of the
differentiated tensor, $F^{tr}$, is time independent for a static solution.
Overall, the only surviving component is $A^{0}$.

In the following, we will assume that the nuclear matter behaves as a perfect
fluid with stress-energy tensor

$\tilde{T}^{\mu\nu}=(\tilde{\rho}+\tilde{p})\tilde{u}^{\mu}\tilde{u}^{\nu}+\tilde{p}\tilde{g}^{\mu\nu}\
,$ (12)

where the total energy density $\tilde{\rho}$, pressure $\tilde{p}$ and the
components of the fluid 4-velocity $\tilde{u}^{\alpha}$ only depend on the
radial coordinate $r$.

###  II.2 Basics of tachyons, ghosts and gradient instabilities

Although they are related to each other, we will categorize our instabilities
into three groups: tachyons, ghosts and gradient instabilities. Their basic
behavior can be demonstrated for a classical massive scalar field theory in
$1+1D$

$\displaystyle S_{\varphi}=\int{\rm d}t{\rm d}x\
\sqrt{-g}\left[-g^{tt}(\partial_{t}\varphi)^{2}-g^{xx}(\partial_{x}\varphi)^{2}-m^{2}\varphi^{2}\right]$
(13)

whose equation of motion is

$\displaystyle
g^{tt}\partial_{t}^{2}\varphi+g^{xx}\partial_{x}^{2}\varphi=m^{2}\varphi\ $
(14)

for a metric $g_{\mu\nu}$ with constant diagonal components.

A Fourier mode $\varphi(t,x)=e^{i\left[\omega(k)t-kx\right]}$ has the
dispersion relation

$\displaystyle\omega(k)=\sqrt{\frac{m^{2}+g^{xx}k^{2}}{-g^{tt}}}\ .$ (15)

If $\omega(k)$ is real, the mode has an oscillatory behavior in time. This
happens in the “usual” case where $g^{tt}<0$, $g^{xx}>0$ and $m^{2}>0$. If
$\omega(k)$ becomes imaginary, the mode grows exponentially, which is an
instability. Let us consider the cases where one of the constants change its
sign while the other two have the usual one.

We call the $m^{2}<0$ case a “tachyon.” The modes with

$\displaystyle k^{2}<\frac{-m^{2}}{g^{xx}}$ (16)

are unstable in this case, hence it is an infrared instability. Note that high
wave number modes are stable, and the fastest growing unstable mode behaves as
$\sim e^{\sqrt{m^{2}/g^{tt}}\,t}$, hence, there is an upper limit to the
growth rate. This is the instability that underlies the spontaneous
scalarization mechanism of the DEF model [4].

We call the $g^{tt}>0$ case a “ghost.” All ghost modes exponentially grow, the
growth rate of the modes diverging at high wave numbers as $\sim
e^{\sqrt{g^{xx}/g^{tt}}k\,t}$ . Thus, unlike a tachyon, a ghost, as described
here, has infinite growth rate for generic perturbations.

We call the $g^{xx}<0$ case a “gradient instability.” The modes with

$\displaystyle k^{2}>\frac{m^{2}}{-g^{xx}}$ (17)

are unstable, and low wave number modes are stable. However, high wave number
modes can grow arbitrarily fast like the ghost as $\sim
e^{\sqrt{g^{xx}/g^{tt}}k\,t}$. Hence, the growth rate of this instability also
diverges.

In our simple model, a coexistence of two of the instabilities is equivalent
to the third instability. For example, a ghost and a gradient with $m^{2}>0$
behaves exactly as a tachyon, which is the simple mathematical result of the
fact that we can multiply the action (13) or the equation of motion (14) with
$-1$, and have the same theory. Similarly, $g^{tt}>0$, $g^{xx}<0$ and
$m^{2}<0$ can be considered as a completely stable theory. However, the scalar
field couples to other fields, and the overall sign of the action is
meaningful in most situations, e.g. in spontaneous scalarization in Eq. (1).
Hence, we will categorize each sign change as a separate instability.

Lastly, we will see that the instabilities we encounter will have $g_{\mu\nu}$
and $m^{2}$ values that are functions of spacetime coordinates, and they can
have different signs in different regions. In the case of the tachyon, having
$m^{2}<0$ only in a finite region of space effectively puts a lower bound on
the $k$ values that can be tachyonic, which might mean no mode is tachyonic
[8]. On the other hand, having different regions of spacetime with positive
and negative $g^{tt}$ values means the coefficient of the principal part of
the equation of motion (14) vanishes at the transition points between regions.
Alternatively, if we divide the whole equation by $g^{tt}$, some of the
coefficients diverge at such points. In any case, perturbative treatment of
such an equation becomes impossible, which we will see to be the case for
vectorization.

###  II.3 Ghost and gradient instabilities in vectorization

The stability analysis for any putative solution to the vectorization theory
can start at the generalized Lorenz condition [15]

$\displaystyle\nabla_{\mu}(\hat{z}A^{\mu})=0$ (18)

which is obtained from Eq. (5) by realizing that
$\nabla_{\mu}\nabla_{\nu}F^{\mu\nu}=0$ due to the antisymmetry of
$F^{\mu\nu}$. Let us fix the metric of the vectorized solution as the
background metric, and consider perturbative deviations of the vector field

$\displaystyle A^{\mu}=\bar{A}^{\mu}+\delta A^{\mu}\ .$ (19)

Here, and from now on, an overbar denotes the value of a function or operator
for the exact vectorized solution. The linearized Lorenz condition becomes

$\displaystyle\overline{\nabla}_{\mu}\left(\Xi^{\mu}{}_{\nu}\delta
A^{\nu}\right)=0\ ,$ (20)

where we define

$\displaystyle\Xi^{\mu}{}_{\nu}\equiv\bar{z}\ \delta^{\mu}{}_{\nu}+2\
\overline{\frac{d\hat{z}}{d(A^{\rho}A_{\rho})}}\ \bar{A}^{\mu}\bar{A}_{\nu}\
.$ (21)

For the exponential conformal factor $\Omega_{\rm V}=e^{\beta
A^{\mu}A_{\mu}/2}$ of Ref. [11]

$\displaystyle m_{\rm V}^{2}\ \Xi^{\mu}{}_{\nu}$
$\displaystyle=-16\pi\beta^{2}\tilde{T}_{\rm
bg}e^{2\beta\bar{A}^{\rho}\bar{A}_{\rho}}\ \bar{A}^{\mu}\bar{A}_{\nu}$
$\displaystyle+\left(m_{\rm V}^{2}-4\pi\beta\tilde{T}_{\rm
bg}e^{2\beta\bar{A}^{\rho}\bar{A}_{\rho}}\right)\ \delta^{\mu}{}_{\nu}\ ,$
(22)

where $\tilde{T}_{\rm bg}$ is the value of the trace of the stress-energy
tensor for the vectorized solution. More explicitly

$\displaystyle\Xi^{0}{}_{0}$ $\displaystyle=1-\frac{4\pi\beta\tilde{T}_{\rm
bg}e^{2\beta\bar{A}^{0}\bar{A}_{0}}\left(1+4\beta\bar{A}^{0}\bar{A}_{0}\right)}{m_{\rm
V}^{2}}$ (23) $\displaystyle\Xi^{i}{}_{i}$
$\displaystyle=1-\frac{4\pi\beta\tilde{T}_{\rm
bg}e^{2\beta\bar{A}^{0}\bar{A}_{0}}}{m_{\rm V}^{2}}\ \ \ \textrm{(no sum)}\ ,$
(24)

where there is no summation over the repeated spatial index $i$ as is the case
for most expressions in this paper. These reduce to the corresponding
expressions in Ref. [15] when $\bar{A}^{\mu}=0$, as expected.

We can already see some of the problematic aspects of the vectorized solutions
directly using the Lorenz condition. The form of the metric and the vector
field for the vectorized star imply that $\Xi^{\mu}{}_{\nu}$ is diagonal on
this background, and Eq. (20) gives the time evolution for $\delta A^{0}$

$\displaystyle\partial_{0}\left(\sqrt{-\bar{g}}\ \Xi^{0}{}_{0}\delta
A^{0}\right)=-\sum_{i=1}^{3}\partial_{i}\left(\sqrt{-\bar{g}}\
\Xi^{i}{}_{i}\delta A^{i}\right)\ ,$ (25)

where the summation over $i$ is shown explicitly.333$A^{0}$ is not a dynamical
degree of freedom since its time derivative does not appear in the action (4),
but the dynamics of the other components imply this time evolution indirectly
[25]. This is an advection equation for the vector density $\delta
B^{\mu}\equiv\sqrt{-\bar{g}}\ \Xi^{\mu}{}_{\nu}\delta A^{\nu}$. Consider a
point $p$ where $\Xi^{0}{}_{0}=0$. Note that the right hand side of Eq. (25)
is generically nonzero at $p$ since there is no symmetry principle to ensure
otherwise. This means even if the initial data is such that $\delta B(t=0)=0$
at $p$ , it will immediately evolve to nonzero values, assuming a meaningful
time evolution is possible. However, this implies that the time derivative of
$\delta A^{0}$ diverges at $p$. In summary, $\Xi^{0}{}_{0}=0$ implies
divergences in the time evolution of generic pertubations around a vectorized
star.

Eq. (20) is first order in time, hence the nature of the divergence is not
apparent. We can show that this is related to a ghost instability by following
Ref.[15] and using the Stueckelberg trick [26], which makes the following
substitution in the action (4)

$A_{\mu}\to A_{\mu}+\frac{1}{m_{\rm V}}\partial_{\mu}\psi\ .$ (26)

This separates out an explicit scalar field, whose equation of motion can be
obtained by varying the action, or equivalently by using the Stueckelberg
trick directly in the Lorenz condition (Eq. (18)),

$\nabla_{\mu}\left[\hat{z}\left(A_{\mu}+\frac{1}{m_{\rm
V}}\nabla_{\mu}\psi\right)\right]=0\ .$ (27)

For perturbative deviations from a vectorized star solution, the linearized
equation is

$\overline{\nabla}_{\mu}\left[\Xi^{\mu\nu}\left(\delta A_{\nu}+\frac{1}{m_{\rm
V}}\overline{\nabla}_{\nu}\delta\psi\right)\right]\ ,$ (28)

where

$\displaystyle\Xi^{\mu\nu}=\overline{g}^{\nu\rho}\Xi^{\mu}{}_{\rho}=\bar{z}\
\overline{g}^{\mu\nu}+2\ \overline{\frac{d\hat{z}}{d(A^{\rho}A_{\rho})}}\
\bar{A}^{\mu}\bar{A}^{\nu}\ .$ (29)

The significance of $\Xi^{\mu\nu}$ becomes clearer when we consider the
principal part, the highest derivative term, of this equation:
${\Xi}^{\mu\nu}\partial_{\mu}\partial_{\nu}\delta\psi$. This means that the
evolution of $\delta\psi$ is governed by the effective inverse metric

$\breve{g}^{\mu\nu}=\Xi^{\mu\nu}\ .$ (30)

In particular, if $\Xi^{00}$ is positive $\delta\psi$ will be a ghost, and if
$\Xi^{ii}$ is negative for any $i$, $\delta\psi$ will have a gradient
instability for large wave numbers. Using this effective metric the full
equation of motion can be expressed as

$\breve{\Box}\delta\psi=-m_{\rm V}\breve{\nabla}\cdot\delta
A+\Delta\Gamma^{\beta}_{\alpha\beta}\Xi^{\mu\alpha}\left(m_{\rm V}\delta
A_{\mu}+\breve{\nabla}_{\mu}\delta\psi\right)$ (31)

where
$\Delta\Gamma^{\mu}_{\alpha\beta}=\breve{\Gamma}^{\mu}_{\alpha\beta}-\Gamma^{\mu}_{\alpha\beta}$,
and both $\breve{\Box}$ and $\breve{\nabla}\cdot$ imply contraction with
$\breve{g}$, rather than $g$. Note that, $\breve{g}^{00}$ or equivalently
$\Xi^{00}$ changes sign exactly when $\Xi^{0}{}_{0}$ changes sign, since this
tensor and the metric are both diagonal in our cases of interest. Hence, the
divergence in the Lorenz condition Eq. (25) occurs exactly when $\delta\psi$
becomes a ghost.

Even though we considered a specific model, Ref. [11], so far, note that our
results are valid for any symmetric $\Xi$, hence applies to other
vectorization models studied in detail in Refs. [12] and [13]. It is, however,
worth noting that non-diagonal terms, which one would expect to find outside
of spherical symmetry, complicate the analysis somewhat, again see Ref. [18].

Lastly, we can also observe the ghost directly in the vector field components
when we recast the field equation (5) in the form

$\Box
A_{\alpha}+\left(\nabla_{\beta}\ln\hat{z}\right)\nabla_{\alpha}A^{\beta}=\mathcal{M}_{\alpha\beta}A^{\beta}\
,$ (32)

where we used the Lorenz constraint and the commutation rules for covariant
derivatives [15]. The mass squared tensor is defined as

$\mathcal{M}_{\alpha\beta}=\hat{z}m^{2}g_{\alpha\beta}+R_{\alpha\beta}-\nabla_{\alpha}\nabla_{\beta}\ln\hat{z}\
.$ (33)

The principle part of Eq. (32) solely consists of $\Box A_{\alpha}$ when it is
linearized around the GR solution for which the vector field vanishes. On the
other hand, for a spherically symmetric vectorized star background, the
following term also contributes to the linearized equation

$\overline{\nabla}_{\alpha}\overline{\nabla}_{\beta}\ln\hat{z}=\frac{2}{\bar{z}}\
\overline{\frac{d\hat{z}}{d(A^{\rho}A_{\rho})}}\bar{A}^{\sigma}\overline{\nabla}_{\alpha}\overline{\nabla}_{\beta}\delta
A_{\sigma}+\dots$ (34)

This means the principal part can not be expressed as a wave operator in
general, however, for $\alpha=0$ in Eq.(32) it can be written as

$\breve{g}^{\mu\nu}\partial_{\mu}\partial_{\nu}\delta A_{0}\ ,$ (35)

where we used the fact that only the $0$ component of $\bar{A}_{\rho}$ is
nonvanishing. Hence, $\breve{g}^{\mu\nu}=\Xi^{\mu\nu}$ is also the effective
metric that governs the evolution of $\delta A_{0}$, which becomes a ghost if
$\breve{g}^{00}$ changes sign, and has a gradient instability if
$\breve{g}^{ii}$ changes sign for any $i$. These are the same conditions for
the existence of instabilities in the scalar field $\psi$ arising from the
Stueckelberg trick.

To summarize, we have shown that the existence of ghost or gradient
instabilities depends on the behavior of $\Xi^{\mu\nu}$, which acts as an
effective (inverse) metric for some degrees of freedom in the theory of action
(4). Since $g_{\mu\nu}$ always has the signature $[-1,1,1,1]$, a change of
sign in $\Xi^{\mu}{}_{\nu}$ occurs exactly under the same conditions of a
change of sign of $\Xi^{\mu\nu}$ as we mentioned. Thus, we will compute
$\Xi^{\mu}{}_{\nu}$ profiles of vectorized neutron stars in Section III in
order to assess their stability.

###  II.4 Tachyons in vectorization

The original aim of vectorization as introduced in Eq. (4) was to replicate
the tachyon-based spontaneous growth of scalar fields arising from action (1)
[11]. In the standard scenario, the stability of the final solution is
realized by nonlinear effects which suppress the instability as the scalar
field grows [8].

The existence of a ghost instability for perturbations around a vectorized
star is sufficient to demonstrate that such a solution is unstable, however,
it is interesting to see if the suppression mechanism of the tachyonic
instability also fails. Namely, it is possible that the nonlinear effects that
are expected to quench the tachyon once a vectorized solution forms fail.

GR backgrounds are known to be susceptible to tachyon instabilities in
vectorization [15]. The existence of ghosts breaks the perturbative approach
due to divergent terms in the time evolution, but let us try to see the
behavior of the tachyonic modes assuming a background solution exists and we
only perturb the degrees of freedom which are tachyonic. We will repeat the
discussion of Ref. [15] for a vectorized background, starting with a vector
spherical harmonic decomposition of the perturbative deviations from the exact
solution

$\delta A_{\alpha}=\frac{1}{r}\sum^{4}_{i=1}\sum_{\ell m}c_{i}u^{\ell
m}_{(i)}(t,r)Z^{(i){\ell m}}_{\alpha}(\theta,\phi),$ (36)

where $c_{1}=c_{2}=1$, $c_{3}=c_{4}=1/\sqrt{\ell(\ell+1)}$, and

$\displaystyle Z^{(1){\ell m}}_{\alpha}$ $\displaystyle=[1,0,0,0]Y_{\ell m},$
(37) $\displaystyle Z^{(2){\ell m}}_{\alpha}$ $\displaystyle=[0,1,0,0]Y_{\ell
m},$ (38) $\displaystyle Z^{(3){\ell m}}_{\alpha}$
$\displaystyle=\frac{r}{\sqrt{\ell(\ell+1)}}[0,0,\partial_{\theta},\partial_{\phi}]Y_{\ell
m},$ (39) $\displaystyle Z^{(4){\ell m}}_{\alpha}$
$\displaystyle=\frac{r}{\sqrt{\ell(\ell+1)}}[0,0,\csc\theta\partial_{\phi},-\sin\theta\partial_{\theta}]Y_{\ell
m}\ ,$ (40)

$Y_{\ell m}$ being the scalar spherical harmonics. Details of this
decomposition can be found in Refs. [15] and [27].

As the background vector field $\bar{A}$ has no $\theta$ or $\phi$ component,
the $u_{4}$ degree of freedom, the so-called axial mode, behaves as in the
case with vanishing vector field, with a slightly different expression for the
effective mass. In particular, the perturbation Lagrangian up to quadratic
terms for any $u_{4}\equiv u_{(4)}^{\ell 0}$ becomes

$\displaystyle S_{\rm t}$
$\displaystyle=\int\mathrm{d}t\mathrm{d}r\frac{e^{\bar{\nu}/2}\left(1-2\bar{\mu}/r\right)^{-1/2}}{4\pi\ell(\ell+1)}\left\\{\vphantom{\frac{1}{2}}e^{-\bar{\nu}}(\dot{u}_{4})^{2}\right.$
$\displaystyle\quad\left.-\left(1-\frac{2\bar{\mu}}{r}\right)(u^{\prime}_{4})^{2}-\left(\frac{\ell(\ell+1)}{r^{2}}+\bar{z}m_{\rm
V}^{2}\right)u_{4}^{2}\right\\}$

after integration over the angular coordinates. This leads to the field
equation

$\displaystyle\bigg{[}$
$\displaystyle-e^{-\bar{\nu}}\partial_{t}^{2}+e^{-\bar{\nu}/2}\sqrt{1-\frac{2\bar{\mu}}{r}}\partial_{r}\left(e^{\bar{\nu}/2}\sqrt{1-\frac{2\bar{\mu}}{r}}\partial_{r}\right)$
$\displaystyle-\frac{\ell(\ell+1)}{r^{2}}-\bar{z}m_{\rm V}^{2}\bigg{]}u_{4}=0\
.$ (42)

If we perform a separation of variables $u_{4}(t,r)=e^{i\omega t}u(r)$, we
obtain

$\displaystyle\left[-\frac{d^{2}}{dr_{*}^{2}}+V_{\rm
eff}(r)\right]u=\omega^{2}\ u\ ,$ (43)

where

$\displaystyle\frac{dr_{*}}{dr}$
$\displaystyle=e^{-\bar{\nu}/2}\left(1-\frac{2\bar{\mu}}{r}\right)^{-1/2}$
(44) $\displaystyle V_{\rm eff}(r)$
$\displaystyle=e^{\bar{\nu}}\left(\frac{\ell(\ell+1)}{r^{2}}+\bar{z}m_{\rm
V}^{2}\right)\ ,$ (45)

and $V_{\rm eff}$ can be considered a function of $r_{*}$ through $r(r_{*})$.
A tachyon exists if there is a solution to this eigenfunction problem that
satisfies $\omega^{2}<0$. Note that the axial mode is only defined for
$\ell\geq 1$.

Increasing $\ell$ makes the modes strictly less tachyonic due to its positive
contribution to $V_{\rm eff}$, which means it is sufficient to check $\ell=1$
for the existence of a tachyon. Non-existence of linearized tachyons in the
axial sector is guaranteed if $V_{\rm eff}>0$ everywhere, since the
expectation value of the differential operator on the left hand side of Eq.
(43) for any trial function is positive in this case, which ensures
$\omega^{2}>0$ for any eigenfunction. We will employ this shortcut when
analyzing our vectorized stars.

The above analysis is only for a single vector spherical harmonic, and does
not have any direct implications for the even modes $u_{1,2,3}$. This means,
we can conclude the existence of tachyons from Eq. (43), but not rule them out
completely even if $u_{4}$ is not tachyonic. We will not investigate the more
complicated coupled equations for $u_{1,2,3}$ in this study, as we will soon
see that the even sector is already plagued by ghost and gradient
instabilities.

##  III Results

Figure 1: Deviations from GR measured by $\overline{\Omega}_{\rm V}(r=0)-1$ as
a function of the ADM mass of vectorized stars for various values of $\beta<0$
and $m_{\rm V}$. This is a slightly updated version of Fig. 1 of Ref. [11].
There are solutions where the vector field or the matter density does not
monotonically decrease with radius (dashed lines), and the general dependence
of deviations from GR on the neutron star mass is qualitatively different from
the case of scalarization. e.g. in Ref. [8].  Figure 2: $\Xi^{0}{}_{0}$,
$\Xi^{r}{}_{r}$, $V_{\rm eff}$, $\Phi=\sqrt{-A_{0}A^{0}}$ and the matter
energy density $\tilde{\rho}$ of a non-monotonic (left) and monotonic (right)
star for $\beta<0$. The radii where $\Xi^{0}{}_{0}$ and $\Xi^{r}{}_{r}$ cut
through $0$ can be seen as the divergences in the logarithmic scale, and the
lines are dashed when the functions attain negative values. $V_{\rm eff}$ is
monotonically decreasing with radius in both cases and is positive everywhere.
Note that the regions where $\Xi^{0}{}_{0}$ and $\Xi^{r}{}_{r}$ become
negative overlap closely, but not exactly.  Figure 3: Maximum values of
$-\Xi^{0}{}_{0}$ and $-\Xi^{r}{}_{r}$ for the neutron stars in Fig. 1 (recall
from Eq. (24) that $\Xi^{i}{}_{i}$ is the same for all $i$). Negative values
are attained for both functions in all cases, hence all solutions have ghost
and gradient instabilities.

For our purposes, we re-computed the static and spherically symmetric
vectorized neutron star solutions of Ref. [11] using the same numerical
methods, and constructed other ones in a broader part of the $(\beta,m_{\rm
V})$ parameter space. In addition, we computed $\Xi^{\mu}{}_{\nu}$ (Eqs. (23)
and (24)) and $V_{\rm eff}$ (Eq. (45)) for each solution to assess stability.
We checked the numerical convergence of our results, and also performed
independent residual analysis. All the neutron star solutions are for the
piecewise polytropic HB equation of state for nuclear matter defined in Ref.
[28].

We categorize the solutions into two groups. The first group has strictly
decreasing $A_{0}$ and matter density $\tilde{\rho}$, so there is no immediate
sign of an instability when one looks at the solutions, hereafter “monotonic
stars”. The second group have visible indications of instability, particularly
$A_{0}$, $\tilde{\rho}$ or both increase at some point within the star,
hereafter “non-monotonic stars”. However, note that non-monotonicity does not
necessarily imply instability, nor does monotonicity imply stability, hence we
examined the perturbations around all solutions.

Summarizing the previous section, we deduce the instabilities of linearized
perturbations of vectorized stars as follows:

  * •

There is a ghost instability if $\Xi^{0}{}_{0}<0$ anywhere.

  * •

There is a gradient instability if $\Xi^{i}{}_{i}<0$ anywhere.

  * •

There is a tachyonic instability if the eigenfunction problem in Eq. (43) has
a negative eigenvalue, $\omega^{2}<0$, for $V_{\rm eff}$. As a shortcut,
$u_{4}$ does not have a tachyon if $V_{\rm eff}>0$ everywhere, hence such a
case is an indication of stability in this mode. However, recall that this
does not rule out the existence of tachyons in the $u_{1,2,3}$ modes.

As we mentioned before, having all instabilities together would mean stability
if the vector field was not coupled to any other field, which is not the case
in our action (4). Moreover, issues with $\hat{z}$ cutting $0$ would remain,
and the instabilities would catch up with us at higher orders in perturbation
theory. Hence the above classification.

The general behavior of vectorized stars in terms of their deviation from GR
can be seen in Fig. 1, some sample star profiles can be seen in Fig. 2.
Applying the above instability criteria to the neutron stars in Fig. 1, we
found that all vectorized solutions have both ghost and gradient
instabilities. Therefore neither type of star is stable, which is the main
result of this study. The “strength” of the ghost and gradient instabilities
of the stars as measured by the most negative values attained by
$\Xi^{\mu}{}_{\nu}$ can be seen in Fig. 3. $\Xi^{\mu}{}_{\nu}$ is
significantly negative, so the existence of the ghost can not be avoided by
small changes to the stellar structure. The behavior of $\Xi^{\mu}{}_{\nu}$
and $V_{\rm eff}$ as functions of radius for two sample vectorized stars can
be seen in Fig. 2.

In contrast to ghosts, we find that the monotonic stars have
$V_{\mathrm{eff}}>0$ everywhere, hence the $u_{4}$ mode is not tachyonic in
this case. The $u_{4}$ tachyon is also stabilized in the majority of the non-
monotonic stars due to the same criterion. There are some non-monotonic stars
where $V_{\mathrm{eff}}<0$ in small regions, but we did not solve the
eigenfunction problem, so these stars may or may not have tachyonic
instabilities.

The naive expectation of a tachyon first growing and eventually getting
quenched in vectorization seems to be partially realized for many, possibly
all, solutions for the $u_{4}$ modes. However, note that $u_{4}=0$ for the
vectorized star, hence, it is not the $u_{4}$ modes themselves that quench the
tachyon, but possibly the indirect nonlinear effects of the other modes
through their coupling. We reiterate that regardless of the fate of the
tachyon(s), the ghost and gradient instabilities not considered in the
original work [11] are present in all solutions, which is enough to render the
vectorized stars unstable.

Figure 4: Deviations from GR measured by $1-\overline{\Omega}_{\rm V}(r=0)$ as
a function of the ADM mass of vectorized stars for various values of $\beta>0$
and $m_{\rm V}$. The behavior of the curves are relatively simpler compared to
the $\beta<0$ case in Fig. 1, which we discuss in Sec. IV. There is a second
branch of solutions appearing at high $\beta$ which consists of non-monotonic
solutions. All curves possibly continue to higher $M_{ADM}$ values than shown,
however the numerical methods we employ are ineffective in this regime. Figure
5: $\Xi^{0}{}_{0}$, $\Xi^{r}{}_{r}$, $V_{\rm eff}$, $\Phi=\sqrt{-A_{0}A^{0}}$
and the energy density $\tilde{\rho}$ of a star with (left) and without
(right) ghosts for $\beta>0$. The radii where $\Xi^{0}{}_{0}$ and
$\Xi^{r}{}_{r}$ cut through $0$ can be seen as the divergence in the
logarithmic scale, and the lines are dashed when the functions attain negative
values. Note that $\Xi^{0}{}_{0}$ and $\Xi^{r}{}_{r}$ can cut through zero
twice, being negative between these radii. $V_{\rm eff}$ is not monotonically
decreasing with radius in these examples, but it is positive everywhere.
Figure 6: Maximum values of $-\Xi^{0}{}_{0}$ and $-\Xi^{r}{}_{r}$ for the
neutron stars in Fig. 4 (recall from Eq. (24) that $\Xi^{i}{}_{i}$ is the same
for all $i$). While negative values are attained for $\Xi^{r}{}_{r}$ in all
cases, $\Xi^{0}{}_{0}>0$ everywhere for some solutions. Thus, all solutions
are unstable due to gradient instabilities, but some do not carry ghosts.

The solutions found in Ref. [11] all had $\beta<0$, however it is known, from
the scalarization case, that it is possible to destabilize massive GR stars
for $\beta>0$ [23]. Therefore we have investigated a part of the $\beta>0$
parameter space for the first time as well. The deviation from GR in this case
can be seen in Fig. 4, monotonic and non-monotonic sample star profiles can be
seen in Fig. 5.

The vectorized stars remained unstable in all cases for $\beta>0$ as well,
which can be seen in Fig 6, though the picture is slightly different. All
solutions still have gradient instabilities hence the major problem of the
$\beta<0$ solutions persists. The ghost instability remains active in the
majority of the solutions, though there do exist monotonic solutions with
$\Xi^{0}{}_{0}>0$ everywhere. As for the $u_{4}$ tachyons, there are solutions
with $V_{\rm eff}>0$ and $V_{\rm eff}<0$ among both monotonic and non-
monotonic stars, but we did not explore whether any of the latter cases indeed
have tachyonic instabilities, since $\Xi^{i}{}_{i}<0$ already ensured that all
solutions are unstable.

The shooting method we utilized to calculate the solutions does not work
reliably for high Arnowitt-Deser-Misner (ADM) mass neutron stars, especially
for higher values of $m_{\rm V}$, hence the ADM masses in our figures do not
necessarily reflect the most massive neutron stars in the theory. There are
also solutions in parts of the parameter space where $A^{0}(r)$ has nodes,
which are not shown in our figures. These are effectively “higher harmonics”
of the fundamental vectorized solution where the vector field does not have
nodes, and are always unstable. The solutions presented so far for very
different $(\beta,m_{\rm V})$ values and ADM masses are all unstable, which
strongly suggests that all static, spherically symmetric vectorized solutions
are unstable, however we cannot conclusively rule out exceptions since the
parameter space is vast. Further regions of the parameter space might be
explored in the future using alternative numerical methods [29].

One final possibility we want to mention is that, however likely or not, the
unstable part of the spherically symmetric regime might be dynamically avoided
in vectorization, and one might have to look for the stable solutions under
less stringent conditions, e.g. axisymmetry. This is possible as the non-
perturbative behaviour is triggered as $\hat{z}$ gets sufficiently small,
before it crosses $0$. Consider the concrete example of a star, far from the
regime where $\hat{z}$ crosses $0$, collapsing to form something compact
enough that, if it were to follow collapse as in GR, $\hat{z}$ would cross
$0$. Initially this should proceed as a GR collapse with some small
perturbations, however as $\hat{z}$ gets small, and before it crosses $0$,
perturbation theory will break down due to the $\hat{z}^{-1}$ terms. From that
point on we cannot say what will happen. It is merely plausible that the non-
perturbative evolution avoids the catastrophic regime we have discussed here
and in Ref. [15], ending at some stable axisymmetric star. Investigation of
such a scenario requires the simulation of collapsing stars, which have been
performed for spontaneous scalarization [30].

##  IV Conclusions

Spontaneous vectorization was motivated by the fact that scalar nature of the
field does not play any direct role in the DEF model. This naively presented
the possibility that scalarization is not an isolated mechanism for a specific
theory, but has analogs in other alternative theories of gravity. While this
may be the case for more general scalar-tensor theories [31], all indications
are that the mechanism cannot be generalized to non-scalar fields.
Specifically, the recent discovery of the existence of ghost and gradient
instabilities around GR solutions of generic spontaneous vectorization
theories showed that vectorization is radically different from scalarization
[16, 15]. Even though vectorization superficially seems to be about updating
the type of the field in spontaneous scalarization, the essence of the theory
is about changing the type of the instability carried by the field, albeit
unintentionally.

In this work, we showed that the ghost and gradient instabilities that are
present for unvectorized (GR) objects plague the static and spherically
symmetric vectorized neutron star solutions as well. This is in stark contrast
to the standard case of scalarized stars in the DEF model and other
scalarization theories inspired by it, where the tachyonic instability is
eventually quenched, leading to stability [20, 21, 22, 23, 8, 24]. Our results
were achieved by linearizing the vector field equations on a fixed vectorized
neutron star background, and showing that the existence of instabilities is
related to components of a tensor, $\Xi^{\mu}{}_{\nu}$, changing sign.

The differences between the vectorized neutron star solutions shown in Fig. 1,
and their counterparts in scalarization had already cast doubts on the
stability of the former [11]. We have confirmed these suspicions, and also
identified the nature of the instabilities. One of the potentially tachyonic
modes is suppressed in many of the vectorized stars in line with the original
expectation, but the more problematic ghost and gradient instabilities are
still present, something not foreseen when vectorization was introduced.

Another striking feature of vectorization is that the deviations from GR in
the $\beta>0$ solutions in Fig. 4 which are studied for the first time have a
much “simpler” structure compared to the $\beta<0$ case in Fig. 1. We believe
this is related to the large field behavior of the conformal factor
$\Omega_{\rm V}$ for the vectorization of spherically symmetric stars behaving
differently from that of $\Omega_{\phi}$ of scalarization. The effective mass
of both the scalar (Eq. (2)) and vector (Eq. (5)) fields depend on the
conformal factors $\Omega_{\phi}$, $\Omega_{\rm V}$ in their respective
theories, since $T=\Omega_{\phi}^{4}\tilde{T}$ in Eq. (2). In the most
familiar case of scalarization with $\beta<0$, the tachyonic nature is
suppressed as the scalar grows since $\Omega_{\phi}=e^{\beta\phi^{2}/2}$
becomes strictly smaller. The opposite is true for scalarization with
$\beta>0$, which means $\Omega_{\phi}$ makes the scalar even more tachyonic as
it grows. This does not necessarily lead to an instability since the nonlinear
changes in $\tilde{T}$ can still quench the tachyon for appropriate coupling
functions, but one can encounter unstable scalarized stars as well [23]. In
short, typically the $\beta<0$ case of scalarization is simpler than that of
$\beta>0$ [23].

The opposite picture is true for vectorization in a spherically symmetric
spacetime where $A^{0}$ is the only non-vanishing component, and
$A_{\mu}A^{\mu}$ is negative. Thus, the simpler case where $\Omega_{\rm
V}=e^{\beta A_{\mu}A^{\mu}/2}$ is suppressed with growing $|A^{0}|$ is that of
$\beta>0$. Similarly, $\Xi^{\mu}{}_{\nu}$ also becomes less unstable, i.e.
less negative, with growing $|A^{0}|$ for $\beta>0$ and $\tilde{T}_{\rm
bg}>0$. This means we might expect the $\beta>0$ solutions of vectorization to
behave more closely to the more familiar $\beta<0$ case of scalarization. This
expectation seems to be partially realized in the simplicity of Fig. 4
compared to Fig. 1, and in the fact that ghost instabilities are suppressed
($\Xi^{0}{}_{0}>0$ everywhere) for some of the vectorized stars with
$\beta>0$, unlike their $\beta<0$ counterparts (see Fig. 3 vs Fig. 6).
Nevertheless, all vectorized solutions we computed are ultimately unstable for
both positive and negative $\beta$.

Even though we have considered a single theory, that of action (4), the
structure of vectorized objects has been studied in other models as well [12,
32, 14]. Extended vector-Gauss-Bonnet theories present an especially
interesting case where spherically symmetric vectorized black holes are known
to be entropically disfavored to Schwarzschild black holes of GR, i.e. the GR
solutions have higher entropy than vectorized black holes for the same ADM
mass [32]. This result is likely due to the vector field and not due to
changing the coupling from matter to the Gauss-Bonnet term, since scalarized
black holes in scalar-Gauss-Bonnet theories are entropically favored to those
of GR [33], and the stability of at least some of these objects has been
confirmed by numerical time evolution [24]. Hence, we conjecture that
vectorized black hole solutions of Ref. [32] are also unstable, which can be
investigated by the methods we developed here.

Our approach can be adapted more broadly to any theory where the linearized
vector field equation is of the form (3). Generalizations of scalarization
beyond vectors also exist [34, 35, 36, 13, 37, 38], and are known to suffer
from similar problems [15], the only known exception being a specific form of
spontaneous spinorization [39]. Developing tools for the stability analysis of
solutions in these theories is a future topic of interest.

Finally, we should reiterate that the problem of the action (4) is not merely
the instability of vectorized neutron stars arising from it. Recall that the
neutron star solutions with vanishing vector fields, which are solutions of
GR, also suffer from ghost and gradient instabilities [16, 15]. Hence, even if
the vectorized solutions were stable, the question of how the time evolution
from a neutron star with no vector field to a vectorized one occurs, if it can
be defined in a meaningful way at all, would remain unanswered. On the other
hand, our results show that if it is possible to have an interpretation of
spontaneous vectorization of action (4) as a well-posed theory, the vectorized
stars that are computed so far are not the astrophysical objects one would
want to study further.

###### Acknowledgements.

AC acknowledges financial support from the European Commission and TÜBİTAK
under the CO-FUNDED Brain Circulation Scheme 2, Project No. 120C081. FMR is
supported by a Young Scientist (BAGEP) Award of Bilim Akademisi of Turkey and
by TÜBİTAK Grant No 117F295.

## References

  * Freire  _et al._ [2012] P. C. C. Freire, N. Wex, G. Esposito-Farese, J. P. W. Verbiest, M. Bailes, B. A. Jacoby, M. Kramer, I. H. Stairs, J. Antoniadis, and G. H. Janssen, The relativistic pulsar-white dwarf binary PSR J1738+0333 II. The most stringent test of scalar-tensor gravity, Mon. Not. Roy. Astron. Soc. 423, 3328 (2012), arXiv:1205.1450 [astro-ph.GA] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021), arXiv:2010.14529 [gr-qc] .
  * Bertotti  _et al._ [2003] B. Bertotti, L. Iess, and P. Tortora, A test of general relativity using radio links with the Cassini spacecraft, Nature 425, 374 (2003).
  * Damour and Esposito-Farèse [1993] T. Damour and G. Esposito-Farèse, Nonperturbative strong field effects in tensor-scalar theories of gravitation, Phys. Rev. Lett. 70, 2220 (1993).
  * Damour and Esposito-Farèse [1996] T. Damour and G. Esposito-Farèse, Tensor-scalar gravity and binary pulsar experiments, Phys. Rev. D 54, 1474 (1996), arXiv:gr-qc/9602056 .
  * Vainshtein [1972] A. I. Vainshtein, To the problem of nonvanishing gravitation mass, Phys. Lett. B 39, 393 (1972).
  * Babichev and Deffayet [2013] E. Babichev and C. Deffayet, An introduction to the Vainshtein mechanism, Class. Quant. Grav. 30, 184001 (2013), arXiv:1304.7240 [gr-qc] .
  * Ramazanoğlu and Pretorius [2016] F. M. Ramazanoğlu and F. Pretorius, Spontaneous Scalarization with Massive Fields, Phys. Rev. D 93, 064005 (2016), arXiv:1601.07475 [gr-qc] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , A Massive Pulsar in a Compact Relativistic Binary, Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Barack  _et al._ [2019] L. Barack _et al._ , Black holes, gravitational waves and fundamental physics: a roadmap, Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Ramazanoğlu [2017] F. M. Ramazanoğlu, Spontaneous growth of vector fields in gravity, Phys. Rev. D 96, 064009 (2017), arXiv:1706.01056 [gr-qc] .
  * Annulli  _et al._ [2019] L. Annulli, V. Cardoso, and L. Gualtieri, Electromagnetism and hidden vector fields in modified gravity theories: spontaneous and induced vectorization, Phys. Rev. D 99, 044038 (2019), arXiv:1901.02461 [gr-qc] .
  * Ramazanoğlu [2019a] F. M. Ramazanoğlu, Spontaneous tensorization from curvature coupling and beyond, Phys. Rev. D 99, 084015 (2019a), arXiv:1901.10009 [gr-qc] .
  * Oliveira and Pombo [2021] J. a. M. S. Oliveira and A. M. Pombo, Spontaneous vectorization of electrically charged black holes, Phys. Rev. D 103, 044004 (2021), arXiv:2012.07869 [gr-qc] .
  * Silva  _et al._ [2021] H. O. Silva, A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou, The ghost of vector fields in compact stars (2021), arXiv:2110.04594 [gr-qc] .
  * Garcia-Saenz  _et al._ [2021] S. Garcia-Saenz, A. Held, and J. Zhang, Destabilization of Black Holes and Stars by Generalized Proca Fields, Phys. Rev. Lett. 127, 131104 (2021), arXiv:2104.08049 [gr-qc] .
  * Esposito-Farèse  _et al._ [2010] G. Esposito-Farèse, C. Pitrou, and J.-P. Uzan, Vector theories in cosmology, Phys. Rev. D 81, 063519 (2010), arXiv:0912.0481 [gr-qc] .
  * Babichev  _et al._ [2018] E. Babichev, C. Charmousis, G. Esposito-Farèse, and A. Lehébel, Hamiltonian unboundedness vs stability with an application to Horndeski theory, Phys. Rev. D 98, 104050 (2018), arXiv:1803.11444 [gr-qc] .
  * Beltran Jimenez  _et al._ [2013] J. Beltran Jimenez, A. L. Delvas Froes, and D. F. Mota, Screening Vector Field Modifications of General Relativity, Phys. Lett. B725, 212 (2013), arXiv:1212.1923 [astro-ph.CO] .
  * Harada [1998] T. Harada, Neutron stars in scalar-tensor theories of gravity and catastrophe theory, Phys. Rev. D 57, 4802 (1998), arXiv:gr-qc/9801049 .
  * Salgado  _et al._ [1998] M. Salgado, D. Sudarsky, and U. Nucamendi, On spontaneous scalarization, Phys. Rev. D 58, 124003 (1998), arXiv:gr-qc/9806070 .
  * Novak [1998] J. Novak, Neutron star transition to strong-scalar-field state in tensor scalar gravity, Phys. Rev. D 58, 064019 (1998), arXiv:gr-qc/9806022 .
  * Mendes and Ortiz [2016] R. F. P. Mendes and N. Ortiz, Highly compact neutron stars in scalar-tensor theories of gravity: Spontaneous scalarization versus gravitational collapse, Phys. Rev. D 93, 124035 (2016), arXiv:1604.04175 [gr-qc] .
  * East and Ripley [2021] W. E. East and J. L. Ripley, Dynamics of Spontaneous Black Hole Scalarization and Mergers in Einstein-Scalar-Gauss-Bonnet Gravity, Phys. Rev. Lett. 127, 101102 (2021), arXiv:2105.08571 [gr-qc] .
  * Heisenberg [2014] L. Heisenberg, Generalization of the Proca Action, JCAP 05, 015, arXiv:1402.7026 [hep-th] .
  * Ruegg and Ruiz-Altaba [2004] H. Ruegg and M. Ruiz-Altaba, The Stueckelberg field, Int. J. Mod. Phys. A 19, 3265 (2004), arXiv:hep-th/0304245 .
  * Rosa and Dolan [2012] J. G. Rosa and S. R. Dolan, Massive vector fields on the Schwarzschild spacetime: quasi-normal modes and bound states, Phys. Rev. D 85, 044043 (2012), arXiv:1110.4494 [hep-th] .
  * Read  _et al._ [2009] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parametrized neutron-star equation of state, Phys. Rev. D 79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Rosca-Mead  _et al._ [2020] R. Rosca-Mead, C. J. Moore, U. Sperhake, M. Agathos, and D. Gerosa, Structure of neutron stars in massive scalar-tensor gravity, Symmetry 12, 1384 (2020), arXiv:2007.14429 [gr-qc] .
  * Sperhake  _et al._ [2017] U. Sperhake, C. J. Moore, R. Rosca, M. Agathos, D. Gerosa, and C. D. Ott, Long-lived inverse chirp signals from core collapse in massive scalar-tensor gravity, Phys. Rev. Lett. 119, 201103 (2017), arXiv:1708.03651 [gr-qc] .
  * Andreou  _et al._ [2019] N. Andreou, N. Franchini, G. Ventagli, and T. P. Sotiriou, Spontaneous scalarization in generalised scalar-tensor theory, Phys. Rev. D 99, 124022 (2019), [Erratum: Phys.Rev.D 101, 109903 (2020)], arXiv:1904.06365 [gr-qc] .
  * Barton  _et al._ [2021] S. Barton, B. Hartmann, B. Kleihaus, and J. Kunz, Spontaneously vectorized Einstein-Gauss-Bonnet black holes, Phys. Lett. B 817, 136336 (2021), arXiv:2103.01651 [gr-qc] .
  * Doneva and Yazadjiev [2018] D. D. Doneva and S. S. Yazadjiev, New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories, Phys. Rev. Lett. 120, 131103 (2018), arXiv:1711.01187 [gr-qc] .
  * Ramazanoğlu [2018a] F. M. Ramazanoğlu, Regularization of instabilities in gravity theories, Phys. Rev. D 97, 024008 (2018a), [Erratum: Phys.Rev.D 99, 069905 (2019)], arXiv:1710.00863 [gr-qc] .
  * Ramazanoğlu [2018b] F. M. Ramazanoğlu, Spontaneous growth of spinor fields in gravity, Phys. Rev. D 98, 044011 (2018b), [Erratum: Phys.Rev.D 100, 029903 (2019)], arXiv:1804.00594 [gr-qc] .
  * Ramazanoğlu [2018c] F. M. Ramazanoğlu, Spontaneous growth of gauge fields in gravity through the Higgs mechanism, Phys. Rev. D 98, 044013 (2018c), arXiv:1804.03158 [gr-qc] .
  * Ramazanoğlu [2019b] F. M. Ramazanoğlu, Various paths to the spontaneous growth of $p$-form fields, Turk. J. Phys. 43, 586 (2019b).
  * Ramazanoğlu and Ünlütürk [2019] F. M. Ramazanoğlu and K. I. Ünlütürk, Generalized disformal coupling leads to spontaneous tensorization, Phys. Rev. D 100, 084026 (2019), arXiv:1910.02801 [gr-qc] .
  * Minamitsuji [2020] M. Minamitsuji, Stealth spontaneous spinorization of relativistic stars, Phys. Rev. D 102, 044048 (2020), arXiv:2008.12758 [gr-qc] .