# Dynamic instability analysis for bumblebee black holes: the odd parity

Zhan-Feng Mai<EMAIL_ADDRESS>Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, China  Rui Xu
<EMAIL_ADDRESS>Department of Astronomy, Tsinghua University, Beijing 100084,
China  Kavli Institute for Astronomy and Astrophysics, Peking University,
Beijing 100871, China  Dicong Liang  Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, China  Lijing Shao
<EMAIL_ADDRESS>Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China

###### Abstract

Spherical black-hole (BH) solutions have been found in the bumblebee gravity
where a vector field nonminimally couples to the Ricci tensor. We study
dynamic (in)stability associated with the gravitational and vector
perturbations of odd parity against these bumblebee BHs. We find that
bumblebee BHs do not suffer ghost instability, but gradient instability and
tachyonic instability exist when the bumblebee charge exceeds certain values.
The existence of the instabilities also depends on the nonminimal coupling
constant $\xi$ that, there is a minimal value $\xi\sim 4\pi G$ with $G$ the
gravitational constant for the instabilities to happen. The theoretical
consideration for bumblebee BH stability turns out to place stronger
constraints on the parameter space than those from the recent observations of
supermassive BH shadows by the Event Horizon Telescope Collaboration. It is
also reminiscent of Penrose’s cosmic censorship conjecture since the charge of
bumblebee BHs cannot be too large due to the dynamic instabilities.
Specifically, for $\xi(\xi-16\pi G)>0$, we find that the charge of a bumblebee
BH cannot be larger than its mass.

##  I Introduction

A framework of Lorentz-symmetry violation using effective field theory in low-
energy approximation was proposed by Kostelecký and Samuel (1989), which is
called the Standard-Model Extension (SME). The SME serves to systematically
study theoretical effects of all possible Lorentz-symmetry violation terms as
extensions of the Standard Model of particle physics and the general
relativity (GR) Kostelecký (2004); Kostelecký and Mewes (2009); Bailey and
Kostelecký (2006); Kostelecký and Mewes (2012, 2013); Kostelecký and Li
(2019). One of the primary Lorentz-violating terms in the gravitational sector
of the SME takes the form $s^{\mu\nu}R_{\mu\nu}$, where $R_{\mu\nu}$ denotes
the Ricci tensor and $s^{\mu\nu}$ is an extra field that breaks Lorentz
symmetry when acquiring a nonzero background value Kostelecký (2004). To
investigate the connection of SME and specific vector-tensor gravitational
theories with Lorentz-symmetry violation, an action was studied by Kostelecký
(2004),

$\displaystyle
I=\int\mathrm{d}^{4}x\sqrt{-g}\bigg{[}\frac{1}{2\kappa}R+\frac{\xi}{2\kappa}B^{\mu}B^{\nu}R_{\mu\nu}-\frac{1}{4}B^{\mu\nu}B_{\mu\nu}-V\big{(}B^{\mu}B_{\mu}\pm
b^{2}\big{)}\bigg{]}\,,$ (1)

where $\kappa=8\pi G$ is related to the Newtonian constant $G$. The theory (1)
is called the bumblebee gravity. Here we use $B_{\mu}$ to denote the bumblebee
vector field. The dynamical term of the bumblebee field in the action is
generated by the tensor
$B_{\mu\nu}\equiv\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$, analogous to
Maxwell’s theory of electromagnetism. The bumblebee theory contains a
nonminimal coupling term, $\sim B^{\mu}B^{\nu}R_{\mu\nu}$, between the Ricci
tensor and the bumblebee field, resembling the SME Lorentz-violating term
$s^{\mu\nu}R_{\mu\nu}$. This term is controlled by a coordinate-independent
parameter $\xi$. When $\xi=0$ and $V$ globally vanishes, the bumblebee theory
reduces to the Einstein-Maxwell theory, identifying $B_{\mu}$ with the four-
dimensional electromagnetic potential. Compared with the Einstein-Maxwell
theory, the vector field in the bumblebee theory has a self-interaction
potential, $V(B^{\mu}B_{\mu}\pm b^{2})$. For a stable vacuum of spacetime, we
require that the potential $V$ is minimized when $B^{\mu}=b^{\mu}$ and
$b^{\mu}b_{\mu}=\mp b^{2}$, implying that the bumblebee vector field has a
nonzero background for a preferred frame and it violates Lorentz symmetry in
the stable vacuum, analogous to the Higgs mechanism for the Higgs scalar
field. In general, a nonvanishing minimum of the potential $V$ is equivalent
to the cosmological constant. For an unknown expression of $V$, it is an often
practice to probe Lorentz-symmetry violation in the asymptotically flat
spacetime under the consideration of a vanishing minimum of $V$. We thus
consider

$\displaystyle V(B^{\mu}B_{\mu}\pm b^{2})\Big{|}_{B^{\mu}=b^{\mu}}=0\,,$ (2)
(3) $\displaystyle V^{\prime}(B^{\mu}B_{\mu}\pm
b^{2})\Big{|}_{B^{\mu}=b^{\mu}}=0\,,$

where $V^{\prime}(x)\equiv{\mathrm{d}V}/{\mathrm{d}x}$.

For the potential with a vanishing minimum, the bumblebee theory is consistent
with the SME in the weak field limit Bailey and Kostelecký (2006); Kostelecký
and Mewes (2009, 2012). Furthermore, to probe Lorentz-symmetry violation in
the strong-field region where relativistic effects of gravity become
important, the properties and associated applications to compact objects, such
as black holes (BHs) and neutron stars, are thus interesting topics in the
bumblebee theory. Recently, several spherical and slowly rotating BH solutions
with bumblebee vector hair and their properties in the bumblebee theory have
been studied Bertolami and Paramos (2005); Casana  _et al._ (2018); Güllü and
Övgün (2022); Izmailov and Nandi (2022); Liu  _et al._ (2023); Xu  _et al._
(2023a); Xu (2023). In particular, Xu  _et al._ (2023a) constructed a class of
spherical BHs with nonvanishing temporal component for the bumblebee field.
These spherical bumblebee BHs only have two degrees of freedom, the mass $M$
and the vector charge $Q$. Specifically, there are two special analytical
solutions in the bumblebee theory: (i) the Reissner-Nordström (RN) BH when
$\xi=0$, and (ii) a stealth Schwarzschild BH carrying a nontrivial bumblebee
hair when $\xi=2\kappa$ Xu  _et al._ (2023a). The current observations of the
shadows of the supermassive BHs, M87∗ and Sgr ${\rm A}^{*}$ from the Event
Horizon Telescope (EHT) already placed some interesting bounds on $Q/M$ for
various values of $\xi$ Xu  _et al._ (2023b, a). These constraints show that
there is still a large parameter space for the bumblebee BHs yet to be
excluded observationally. Liang  _et al._ (2023) illustrated that it is
possible to probe the vector charge-to-mass ratio up to $Q/M\sim{\cal
O}(10^{-3})$ with future observations of millihertz gravitational waves from
extreme-mass-ratio inspirals, so that we might distinguish bumblebee BHs from
the Schwarzschild BH. A recent study by Hu  _et al._ (2023) shows that timing
a radio pulsar orbiting around Sgr ${\rm A}^{*}$ in a close orbit can as well
probe a vector charge-to-mass ratio as small as $Q/M\sim{\cal O}(10^{-3})$.

Regardless of the constraints on the vector charge from observations, as one
of the important intrinsic properties of BHs, BH stability answers whether
such a BH can exist in our Universe from the theoretical perspective. It might
give extra constraints on the parameter space of the bumblebee BHs. In a
previous study, the extended thermodynamics and associated local thermodynamic
stability of bumblebee BHs have been studied by us Mai  _et al._ (2023).
Unlike the RN BH which has one phase transition point separating locally
thermodynamic stability and instability, the bumblebee BHs can have zero to
two phase transition points, depending on the value of the coupling constant
$\xi$ Mai  _et al._ (2023).

In this work, we study the dynamic stability of bumblebee BHs. As shown by the
action (1), the bumblebee theory can be considered as a class of vector-tensor
theories. For the special case of a vanishing $\xi$, the bumblebee theory
recovers the Einstein-Maxwell theory and the RN BH solution has been proven to
be stable against generic perturbations Moncrief (1974); Regge and Wheeler
(1957). For general bumblebee BHs, we are going to focus on the gradient
instability and the ghost instability. These two types of instabilities have
been widely investigated for hairy BHs in various modified gravity theories,
such as the $f(R)$ gravity De Felice  _et al._ (2011); Motohashi and Suyama
(2011), scalar-tensor theories Kobayashi  _et al._ (2012, 2014); Ogawa  _et
al._ (2016); Ganguly  _et al._ (2018), vector-tensor theories Heisenberg  _et
al._ (2017); Kase  _et al._ (2018a, b), and scalar-vector-tensor theories
Gannouji and Baez (2022); Kase and Tsujikawa (2023). It seems that compact
objects suffer the gradient and ghost instabilities due to nonminimal
couplings between the gravitational field and extra fields. Specifically, Kase
_et al._ (2018a) considered vectorized BHs in the vector-tensor theory
including a Horndeski-like nonminimal coupling term
$G_{\mu\nu}A^{\mu}A^{\nu}$, finding that there are parameter-dependent ghost
and gradient instabilities associated with the odd-parity gravitational and
vector perturbations. In addition, BHs in the Einstein-æther theory, which is
one of the Lorentz-violating gravity theories, suffer the ghost and gradient
instabilities associated with the odd-parity perturbations as well Tsujikawa
_et al._ (2021). The gradient and ghost instabilities of bumblebee BHs have
not been studied yet. Our work complements to the existing knowledge of BH
instabilities in vector-tensor theories and provides concrete examples with
detailed numerical results on the conditions for the gradient and ghost
instabilities of bumblebee BHs.

This paper is organized as follows. In Sec. II, we review the ghost and
gradient instabilities briefly. In Sec. III.1, we review the basic properties
of the vectorized bumblebee BHs given by Xu  _et al._ (2023a). In Sec. III.2,
we first introduce the gravitational and the bumblebee vector perturbations of
odd parity. From the perturbed action up to the second order, we find that
bumblebee BHs have no ghost instability, and then we give the generic
conditions of the gradient instability. In Sec. III.3, we analyze the gradient
instability in detail with the numerical BH solutions, finding that in the
absence of the gradient instability, bumblebee BHs cannot carry very large
values of the bumblebee charge. In Sec. IV, we summarize our investigation and
give some discussions for further studies. In this paper, we adopt $(-,+,+,+)$
as the metric convention and the Planck natural units, namely
$G=c=\hbar=4\pi\epsilon_{0}=k_{B}=1$.

##  II Gradient and ghost instabilities in a toy field theory

We aim to investigate the dynamic stability of bumblebee BHs. To do this, we
first generally consider a field perturbation under the BH background and
study whether it will grow with time. Let us introduce various stabilities
using a toy model Demirboğa  _et al._ (2022); Doneva  _et al._ (2022). As a
simple example, we first consider a single scalar field
$\psi(t,r,\theta,\varphi)$ in four-dimensional spherical coordinate. In
general, the action of a free scalar field in Minkowski spacetime is

$\displaystyle S$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int\mathrm{d}^{4}x\left(\partial_{\mu}\psi\partial^{\mu}\psi\right)$
(4) (5) $\displaystyle=$ $\displaystyle\frac{1}{2}\int
r^{2}\mathrm{d}t\mathrm{d}r\mathrm{d}\Omega\left[\dot{\psi}^{2}-(\bm{\nabla}\psi)^{2}\right]\,,$
(6)

where $\dot{\psi}$ denotes the time derivative and $\bm{\nabla}$ denotes the
gradient operator in three spatial dimensions. To investigate various
instabilities, we introduce two coefficients, $c_{1}$ and $c_{2}$, and
construct the following action

$S=\frac{1}{2}\int
r^{2}\mathrm{d}t\mathrm{d}r\mathrm{d}\Omega\left[c_{1}\dot{\psi}^{2}-c_{2}(\bm{\nabla}\psi)^{2}\right]\,.$
(7)

Due to the spherical symmetry, one can set

$\psi(t,r,\theta,\varphi)=\sum_{\ell,m}\frac{\phi_{\ell m}(t,r)}{r}{\rm
Y}_{\ell m}(\theta,\varphi)\,,$ (8)

where $\ell$ and $m$ are the total and azimuthal angular momentum numbers
respectively, satisfying $\ell=0,1,2,\cdots$ and $m=0,\pm 1,\pm
2,\cdots,\pm\ell$. After integrating over the angular coordinates, the action
(7) yields

$\displaystyle S=\frac{1}{2}\sum_{\ell,m}\int\mathrm{d}t\mathrm{d}r\bigg{[}$
$\displaystyle c_{1}|\dot{\phi}_{\ell m}|^{2}-c_{2}|\phi^{\prime}_{\ell
m}|^{2}-\frac{c_{2}\ell(\ell+1)}{r^{2}}|\phi_{\ell m}|^{2}\bigg{]}\,.$ (9)

Here we have performed the integration by parts with respect to $r$.

In the following, we drop the $\ell,m$ indices of the scalar mode for
simplicity. Considering a Fourier mode $\phi(t,r)=\phi_{k}\mathrm{e}^{i(\omega
t-kr)}$, the dispersion relation from Eq. (9) is approximately

$\omega^{2}c_{1}-k^{2}c_{2}-\frac{c_{2}\ell(\ell+1)}{r^{2}}=0.$ (10)

Here we have taken the high-energy limit so that the derivatives of
$\phi_{k},\,\omega$ and $k$ are neglected. Depending on the relative size of
$k$ and $\ell/r$, there are two simple cases useful to consider, namely

$\displaystyle k\gg\ell/r:~{}~{}\quad\quad$
$\displaystyle\omega^{2}=\frac{c_{2}}{c_{1}}k^{2}\,,$ (11) $\displaystyle
k\ll\ell/r:~{}~{}\quad\quad$
$\displaystyle\omega^{2}=\frac{c_{2}}{c_{1}}\frac{\ell^{2}}{r^{2}}\,.$ (12)

We can further define the propagating speed along the radial direction as

$c_{r}\equiv\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{\omega}{k}=\sqrt{\frac{c_{2}}{c_{1}}}\,,$
(13)

and the propagating speed along the angular direction as

$\quad
c_{\Omega}\equiv\frac{r\mathrm{d}\theta}{\mathrm{d}t}=\frac{r\omega}{\ell}=\sqrt{\frac{c_{2}}{c_{1}}}\,.$
(14)

In this simple model (7), $c_{r}$ and $c_{\Omega}$ are the same.

Obviously, $\phi$ is stable when $c_{1}$ and $c_{2}$ have the same sign,
implying that all modes of $\phi$ are propagating modes. On the contrary,
$\phi$ is unstable when $c_{1}$ and $c_{2}$ have different signs, giving an
imaginary propagating speed, implying that $\phi$ grows exponentially with
time. Furthermore, the instability is usually classified into two types. When
$c_{1}<0$ and $c_{2}>0$, $\phi$ is said to have the _ghost instability_. When
$c_{1}>0$ and $c_{2}<0$, $\phi$ is said to have the _gradient instability_. In
Table 1, we list case by case the conditions for ghost and gradient
instabilities for different $c_{1}$ and $c_{2}$.

Table 1: The conditions of ghost instability and gradient instability for the action (7). Instability | Conditions
---|---
Propagating mode |  $c_{1}>0$ |  $c_{2}>0$ |  $c_{r}^{2}\,,c_{\Omega}^{2}>0$
Ghost instability |  $c_{1}<0$ |  $c_{2}>0$ |  $c_{r}^{2}\,,c_{\Omega}^{2}<0$
Gradient instability |  $c_{1}>0$ |  $c_{2}<0$ |  $c_{r}^{2}\,,c_{\Omega}^{2}<0$
Propagating mode |  $c_{1}<0$ |  $c_{2}<0$ |  $c_{r}^{2}\,,c_{\Omega}^{2}>0$

Before we generalize the conditions to the case of multiple scalar fields, we
need to point out that the instability analysis becomes less straightforward
when $c_{1}$ and $c_{2}$ depend on $r$. One can only conclude that when
$c_{1}$ and $c_{2}$ have the same sign for any $r$, there are no
instabilities. If $c_{1}$ or $c_{2}$ becomes negative in some intervals of
$r$, the existence of ghost instability or gradient instability depends on the
specific expressions of $c_{1}$ and $c_{2}$. In conclusion, for $r$-dependent
$c_{1}$ and $c_{2}$, we only have a necessary condition for ghost or gradient
instability: $c_{1}<0$ or $c_{2}<0$ for some intervals of $r$, or
equivalently, a sufficient condition for $\phi$ being stable: $c_{1}c_{2}>0$
for any $r$.

Now we discuss the ghost and gradient instabilities generalized to a toy model
with multiple scalar fields. Consider $n$ coupled scalar fields
($\psi_{1},\psi_{2},\cdots,\psi_{n}$) in spherical coordinates, which can be
written in the vector form

$\vv{\psi}=(\psi_{1},\psi_{2},\cdots\psi_{n})^{\intercal}\,.$ (15)

Expanding with the spherical harmonics,

$\psi_{i}=\sum_{\ell,m}\frac{\phi_{i}(t,r)}{r}{\rm Y}_{\ell
m}(\theta,\varphi)$ (16)

and after integrating over the angular coordinates, the action of the model
reduces to

$S=\sum_{\ell,m}\int\mathrm{d}t\mathrm{d}r\Big{(}{\dot{\vv{\phi}}}^{\intercal}\bm{K}\dot{\vv{\phi}}+{\vv{\phi}^{\prime}}^{\intercal}\bm{G}\vv{\phi}^{\prime}+\vv{\phi}^{\intercal}\bm{M}\vv{\phi}\Big{)}\,,$
(17)

where $\bm{K},\bm{G},\bm{M}$ are $n\times n$ matrices that may depend on $r$.
In addition, we assume that $\bm{M}$ depends on $\ell$. Similarly, we consider
a Fourier mode $\vv\phi(t,r)=\vv{\phi}_{k}\mathrm{e}^{i(\omega t-kr)}$, and
take the high-energy limit to get the dispersion relation

$\det\big{(}\omega^{2}\bm{K}+k^{2}\bm{G}+\bm{M}\big{)}=0\,.$ (18)

In general, it is difficult to deduce necessary and sufficient conditions from
Eq. (18) for the ghost and gradient instabilities. The two simplified cases in
Eqs. (11) and (12) are useful to consider. The case where the $\ell$-terms are
dropped corresponds to

$\det\big{(}\omega^{2}\bm{K}+k^{2}\bm{G}\big{)}=0\,,$ (19)

from which the radial propagating speed can be solved,

$c_{r}^{2}=\lim_{k\rightarrow\infty}\frac{\omega^{2}}{k^{2}}\,.$ (20)

The case where $k$ is dropped corresponds to

$\det\big{(}\omega^{2}\bm{K}+\bm{M}\big{)}=0\,,$ (21)

from which the angular propagating speed can be solved,

$c_{\Omega}^{2}=\lim_{\ell\rightarrow\infty}\frac{r^{2}\omega^{2}}{\ell^{2}}.$
(22)

A set of necessary conditions for stable propagating modes then can be
obtained by requiring

$c_{r}^{2}>0\,,\quad\quad c_{\Omega}^{2}>0\,.$ (23)

Because the matrices $\bm{K},\bm{G},\bm{M}$ generally depend on $r$, the
speeds calculated using Eqs. (20) and (22) are functions of $r$. Equation (23)
requires them to be positive everywhere in the domain of interest.

If either $c_{r}^{2}$ or $c_{\Omega}^{2}$ is negative, then instability
occurs. Analogous to the single-field model, if $c_{r}^{2}<0$ or
$c_{\Omega}^{2}<0$ while $\bm{K}$ is non-positive definite, then the
instability is the ghost instability. If $c_{r}^{2}<0$ while $\bm{G}$ is non-
negative definite, then the instability is gradient instability. Besides the
ghost instability and the gradient instability, there might be another type of
instability when $c_{\Omega}^{2}<0$ and $\bm{M}$ is non-negative definite.
This is called the _tachyonic instability_. In the single-field model, it
happens to be the same coefficient $c_{2}$ appearing in both $c_{r}$ and
$c_{\Omega}$, so the gradient instability and the tachyonic instability occur
simultaneously. In the multi-field model, the elements of the matrix $\bm{M}$
are in general not related to the matrix $\bm{G}$, so the gradient instability
and the tachyonic instability can occur independently.

##  III Gradient and ghost instabilities of bumblebee BHs

###  III.1 BHs with bumblebee charge

We first briefly review the bumblebee BH solutions that extend the RN BH
solution, and detailed derivations can be found in Refs. Xu  _et al._ (2023a);
Mai  _et al._ (2023). The key assumptions to obtain these bumblebee BH
solutions include the following.

  * •

The potential term $V(\cdot)$ in Eq. (1) can be ignored because either the
background bumblebee field $b_{\mu}$ associated with the BH solutions
satisfies Eq. (2) or the potential has a characteristic length at the
cosmological scale so that it plays little role at the lengthscales of BHs.

  * •

The background bumblebee field has only the temporal component nonvanishing,
namely $b_{\mu}=(b_{t},0,0,0)$. In principle, all spherical, static BH
solutions with $b_{\mu}=(b_{t},b_{r},0,0)$ in the bumblebee gravity were
obtained by Xu  _et al._ (2023a). We here focus on the branch of bumblebee BHs
with temporal component only.

The field equations for the metric of bumblebee gravity and the associated
bumblebee background field $b_{\mu}$ are then Xu  _et al._ (2023a)

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\kappa
T^{b}_{\mu\nu}\,,\quad\quad\nabla_{\mu}b^{\mu\nu}=\frac{\xi}{\kappa}R^{\mu\nu}b_{\mu}\,,$
(24)

where $b_{\mu\nu}=\partial_{\mu}b_{\nu}-\partial_{\nu}b_{\mu}$, and

$\displaystyle T^{b}_{\mu\nu}=\frac{\xi}{2\kappa}$
$\displaystyle\bigg{(}g_{\mu\nu}b^{\alpha}b^{\beta}R_{\alpha\beta}-2b_{\mu}b_{\lambda}R_{\nu}^{\phantom{\nu}\lambda}-2b_{\nu}b_{\lambda}R_{\mu}^{\phantom{\mu}\lambda}-\Box_{g}(b_{\mu}b_{\nu})-g_{{\mu\nu}}\nabla_{\alpha}\nabla_{\beta}(b^{\alpha}b^{\beta})+\nabla_{\kappa}\nabla_{\mu}\left(b^{\kappa}b_{\nu}\right)+\nabla_{\kappa}\nabla_{\nu}(b_{\mu}b^{\kappa})\bigg{)}$
(25) (27)
$\displaystyle+b_{\mu\lambda}b_{\nu}{}^{\lambda}-\frac{1}{4}g_{\mu\nu}b^{\alpha\beta}b_{\alpha\beta}\,.$

Here $\Box_{g}\equiv\nabla_{\alpha}\nabla^{\alpha}$ is defined as the
d’Alembertian in the curved spacetime and $\nabla_{\alpha}$ is the covariant
derivative operator.

With the spherical ansatz for the metric

$\displaystyle\mathrm{d}s^{2}=-h(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{f(r)}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\varphi^{2}\right),$
(28)

one gets three equations from Eq. (24) to solve for the metric functions
$h(r)$, $f(r)$, and the bumblebee field component, $b_{t}(r)$. Two analytical
solutions are worth mentioning.

  1. (I)

When $\xi=0$, we have

$\displaystyle f=h=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}},\quad
b_{t}=\mu_{\infty}-\sqrt{\frac{2}{\kappa}}\frac{Q}{r}\,.$ (29)

It is the RN solution as one expects for $\xi=0$ ($\mu_{\infty}$ is a constant
usually taken to be zero).

  2. (II)

When $\xi=2\kappa$, we have

$\displaystyle
f=h=1-\frac{2M}{r}\,,~{}b_{t}=\frac{Q}{\sqrt{2\kappa}M}\left(1-\frac{2M}{r}\right)\,.$
(30)

It is an interesting case showing that the Schwarzschild metric can be
accompanied by a simple but nonzero bumblebee field.

For other arbitrary values of $\xi$, we could not find analytical solutions.
Numerical BH solutions are calculated instead. These numerical BHs have the
following asymptotic behavior,

$\displaystyle
h(r)\Big{|}_{r\to\infty}=1-\frac{2M}{r}+\frac{\tilde{h}_{2}}{r^{2}}+\cdots\,,$
(31) (33) $\displaystyle
f(r)\Big{|}_{r\to\infty}=1-\frac{2M}{r}+\frac{\tilde{f}_{2}}{r^{2}}+\cdots\,,$
$\displaystyle
b_{t}(r)\Big{|}_{r\to\infty}=\mu_{\infty}-\sqrt{\frac{2}{\kappa}}\frac{Q}{r}+\frac{\tilde{b}_{2}}{r^{2}}+\cdots\,,$

where the ADM mass $M$ and the bumblebee charge $Q$ are the only two free
parameters for the solutions, and other expansion coefficients,
$\mu_{\infty},\,\tilde{h}_{2},\,\tilde{f}_{2},\,\tilde{b}_{2},\,\cdots$, are
recursively related to $M$ and $Q$. At the event horizon $r=r_{h}$, the
functions $h,\,f,\,b_{t}$ behave as

$\displaystyle f(r)\Big{|}_{r\to
r_{h}}=f_{1}(r-r_{h})+f_{2}(r-r_{h})^{2}+\cdots\,,$ (35) (37) $\displaystyle
h(r)\Big{|}_{r\to r_{h}}=h_{1}(r-r_{h})+h_{2}(r-r_{h})^{2}+\cdots\,,$ (39)
$\displaystyle b_{t}(r)\Big{|}_{r\to
r_{h}}=b_{t1}(r-r_{h})+b_{t2}(r-r_{h})^{2}+\cdots,$

where the coefficients,
$f_{1},\,h_{1},\,b_{t1},\,f_{2},\,h_{2},\,b_{t2},\,\cdots$, and the horizon
radius $r_{h}$ can be related to $M$ and $Q$ numerically.

In Fig. 1, we plot $r_{h}$ with respect to the bumblebee charge $Q$ for
different values of $\xi$. We note that when $\xi<2\kappa$, the BHs have
maximal vaules for the bumblebee charge; a typical example is $Q_{\rm max}=M$
in the RN case when $\xi=0$. When $\xi=2\kappa$, the numerical result shows
that $Q$ becomes unrestricted, agreeing with the analytical solution in Eq.
(30). For $\xi>2\kappa$, our numerical method can produce very large values of
$Q$ with extremely small $r_{h}$, suggesting that $Q$ is also unrestricted
when $\xi>2\kappa$. But the errors of our numerical solutions increase rapidly
at large $Q$ so that the solutions quickly become untrustworthy before we can
extract further information.

To figure out whether $Q$ has maximal values when $\xi>2\kappa$, we have made
another insufficient yet plausible attempt. We find approximate solutions
around $\xi=2\kappa$,

$\displaystyle
f=h=1-\frac{2M}{r}-\frac{\delta\xi}{2\kappa}\frac{Q^{2}}{r^{2}}+O\big{(}\delta\xi^{2}\big{)},$
$\displaystyle
b_{t}=\frac{Q}{\sqrt{2\kappa}M}\left(1-\frac{2M}{r}-\frac{\delta\xi}{2\kappa}\frac{Q^{2}}{r^{2}}+O\big{(}\delta\xi^{2}\big{)}\right),$
(40)

where $\delta\xi=\xi-2\kappa$. Then the radius of horizon has an approximation

$\displaystyle
r_{h}=2M+\frac{\delta\xi}{4\kappa}\frac{Q^{2}}{M}+O\big{(}\delta\xi^{2}\big{)}.$
(41)

For the approximation to work at a small enough $|\delta\xi|$, we have an
estimation

$\displaystyle\frac{|\delta\xi|}{4\kappa}\frac{Q^{2}}{M}\lesssim 2M,$ (42)

which tells us that $Q/M$ is bounded at least at a small enough $|\delta\xi|$
via

$\displaystyle\frac{Q}{M}\lesssim\sqrt{\frac{8\kappa}{|\delta\xi|}}.$ (43)

We cannot think of a reason for $Q$ being unbounded when $\delta\xi$ is large
if it already has a bound around $\xi=2\kappa$. So we tentatively conclude
that $Q$ may also have maximal values when $\xi>2\kappa$. The corresponding
boundary is currently represented by the largest values of $Q$ that we can
find using our numerical code while limiting the relative error in the radius
of horizon to $5\%$.

Figure 1: BH horizon $r_{h}$ versus BH charge $|Q|$ for different $\xi$.

###  III.2 Field perturbation: odd parity

With the bumblebee BH solutions in hand, we are ready to investigate the
dynamical stability of bumblebee BHs. We consider a gravitational perturbation
$\epsilon h_{\mu\nu}$ and a vector perturbation $\epsilon\delta b_{\mu}$ in
the bumblebee BH background with $\epsilon\ll 1$. Then, including the
bumblebee BH background $g_{\mu\nu}$ and bumblebee vector $b_{\mu}$, the
perturbed metric and bumblebee vector field are given by

$g_{\mu\nu}\to g_{\mu\nu}+\epsilon h_{\mu\nu}\,,\quad\quad b_{\mu}\to
b_{\mu}+\epsilon\delta b_{\mu}\,.$ (44)

In general, on the spherical ansatz background, the perturbations can be
decomposed into odd parity modes and even parity modes based on the following
rules. For the gravitational perturbation $h_{\mu\nu}$, the components
$h_{tt}$, $h_{tr}$ and $h_{rr}$ transform as scalars under rotation on the
two-dimensional sphere of $\theta$ and $\varphi$. For example, an arbitrary
scalar function, $\Psi(t,r,\theta,\varphi)$, can be written as summation of
spherical harmonic function, ${\rm Y}_{\ell m}(\theta,\varphi)$, as

$\Psi(t,r,\theta,\varphi)=\sum_{\ell,m}\Psi_{\ell m}(t,r){\rm Y}_{\ell
m}(\theta,\varphi)\,.$ (45)

Since the spherical harmonic function transforms as

${\rm Y}_{\ell m}(\pi-\theta,\varphi+\pi)=(-1)^{\ell}{\rm Y}_{\ell
m}(\theta,\varphi)$

under reflection on the two-dimensional sphere, one finds

$\Psi(t,r,\pi-\theta,\varphi+\pi)=(-1)^{\ell}\Psi(t,r,\theta,\varphi)\,.$

It is called the even parity mode. On the other hand, with
$a,b=\theta,\varphi$, $h_{ta}$ and $h_{ra}$ transform as vector, while
$h_{ab}$ transforms as tensor under two-dimensional reflection. For the vector
or tensor components, the mentioned even parity modes transform with the
$(-1)^{\ell}$ factor, while the others transform with the $(-1)^{\ell+1}$
factor and are called odd parity modes Regge and Wheeler (1957). Given an
arbitrary vector or symmetric tensor on the two-dimensional sphere, $V_{a}$
and ${\cal T}_{ab}$, they can be decomposed into,

$\displaystyle V_{a}=$ $\displaystyle
D_{a}\Psi_{1}+\sqrt{\gamma}\varepsilon^{b}{}_{a}D_{b}\Psi_{2}\,,$ (46)
$\displaystyle{\cal T}_{ab}=$ $\displaystyle
D_{a}D_{b}\Psi_{3}+\gamma_{ab}\Psi_{4}+\frac{\sqrt{\gamma}}{2}\big{(}\varepsilon^{c}{}_{a}D_{c}D_{b}+\varepsilon^{c}{}_{b}D_{c}D_{a}\big{)}\Psi_{5}\,,$

where the anti-symmetric tensor $\varepsilon_{ab}$ is defined on the two-
dimensional sphere with $\varepsilon_{\theta\varphi}=1$, and $\Psi_{i}$ are
scalar functions. We denote $\gamma_{ab}$ as the induced metric on this
sphere, with $\gamma$ its determinant and $D_{a}$ its corresponding covariant
derivative.

In this paper, we focus on the bumblebee BHs’ dynamical stability associated
with the odd parity perturbations. Following the above illustrations, the
gravitational perturbation $h_{\mu\nu}$ of odd parity can be written as

$\displaystyle h_{tt}$ $\displaystyle=h_{rr}=h_{rt}=0\,,$ $\displaystyle
h_{ta}$ $\displaystyle=\sum_{\ell,m}h_{0\ell
m}(t,r)\sqrt{\gamma}\varepsilon_{a}{}^{b}D_{b}{\rm Y}_{\ell
m}(\theta,\varphi)\,,$ $\displaystyle h_{ra}$
$\displaystyle=\sum_{\ell,m}h_{1\ell
m}(t,r)\sqrt{\gamma}\varepsilon_{a}{}^{b}D_{b}{\rm Y}_{\ell
m}(\theta,\varphi)\,,$ $\displaystyle h_{ab}$
$\displaystyle=\frac{1}{2}\sum_{\ell,m}h_{2\ell
m}(t,r)\sqrt{\gamma}\big{(}\varepsilon^{c}{}_{a}D_{c}D_{b}+\varepsilon^{c}{}_{b}D_{c}D_{a}\big{)}{\rm
Y}_{\ell m}(\theta,\varphi)\,.$

The odd parity perturbation of the bumblebee vector can be written as

$\displaystyle\delta b_{t}$ $\displaystyle=\delta b_{r}=0\,,$ (47)
$\displaystyle\delta b_{a}$ $\displaystyle=\sum_{\ell,m}\beta_{\ell
m}(t,r)\sqrt{\gamma}\varepsilon_{a}{}^{b}\partial_{b}{\rm Y}_{\ell m}.$ (48)

It is worthwhile to point out that not all perturbations are physical because
of the gauge degrees of freedom. Note that the diffeomorphism still holds in
the bumblebee theory, and one can consider an infinitesimal diffeomorphism
transformation, $x^{\mu}\to x^{\mu}+\lambda^{\mu}$, yielding

$h^{\prime}_{\mu\nu}\to
h_{\mu\nu}+\nabla_{\nu}\lambda_{\mu}+\nabla_{\mu}\lambda_{\nu}\,.$ (49)

It indicates that one can always find four scalar functions to simplify the
perturbations. Based on the decomposition rules for odd parity modes, we
consider

$\displaystyle\lambda_{t}$ $\displaystyle=\lambda_{r}=0\,,$ (50)
$\displaystyle\lambda_{a}$ $\displaystyle=\sum_{\ell,m}\Lambda_{\ell
m}(t,r)\sqrt{\gamma}\varepsilon_{a}{}^{b}\partial_{b}{\rm Y}_{\ell m}\,,$ (51)

such that the gravitational perturbations, $h_{0\ell m}$, $h_{1\ell m}$, and
$h_{2\ell m}$, respectively transform as

$\displaystyle h_{0\ell m}\to h_{0\ell m}+\frac{\mathrm{d}\Lambda_{\ell
m}}{\mathrm{d}t}\,,$ (52) (54) $\displaystyle h_{1\ell m}\to h_{1\ell
m}+\frac{\mathrm{d}\Lambda_{\ell m}}{\mathrm{d}r}-\frac{2}{r}\Lambda_{\ell
m}\,,$ $\displaystyle h_{2\ell m}\to h_{2\ell m}+2\Lambda_{\ell m}\,.$

Specifically, the monopole with $\ell=0$ does not exist in odd-parity modes.
As for the dipole mode with $\ell=1$, $h_{ab}$ vanishes identically. We
therefore discuss this case separately. For $\ell\geq 2$, we shall choose the
Regge-Wheeler gauge such that $h_{2\ell m}=0$ and two physical components,
$h_{0\ell m}$ and $h_{1\ell m}$, are left Regge and Wheeler (1957). Following
the above arguments, the gravitational perturbation, together with the
bumblebee vector perturbation, can be rewritten in the following matrix form,

$\displaystyle h_{\mu\nu}$
$\displaystyle=\sum_{\ell,m}\begin{pmatrix}0&&0&&-{h_{0}(t,r)}\sin^{-1}\theta\,\partial_{\varphi}&&\sin\theta\,h_{0}(t,r)\partial_{\theta}\\\
0&&0&&-{h_{1}(t,r)}\sin^{-1}\theta\,\partial_{\varphi}&&\sin\theta\,h_{1}(t,r)\partial_{\theta}\\\
*&&*&&0&&0\\\ *&&*&&0&&0\\\ \end{pmatrix}{\rm Y}_{\ell m}(\theta,\varphi)\,,$
(56) $\displaystyle\delta b_{\mu}(t,r)$
$\displaystyle=\sum_{\ell,m}\beta(t,r)\begin{pmatrix}0\\\ 0\\\
-\sin^{-1}\theta\,\partial_{\varphi}\\\
\sin\theta\,\partial_{\theta}\end{pmatrix}{\rm Y}_{\ell m}(\theta,\varphi)\,,$
(57)

where the “*” represents the symmetric part of the matrix. For convenience, in
the following, we drop the $\ell,m$ indices. With Eqs. (56) and (57), we
expand action (1) up to ${\cal O}(\epsilon^{2})$. Then, after performing
integration with respect to $\theta$ and $\varphi$, as well as integrating by
parts with respect to $t$ and $r$, the resulting perturbed action of odd
parity modes reads

$I=I_{0}+\epsilon^{2}\sum_{\ell,m}\frac{\ell(\ell+1)}{4\kappa}I_{\rm odd}\,,$
(58)

where $I_{0}$ denotes the on-shell action satisfying the equations of motion
(24). Moreover, we have

$\displaystyle I_{\rm odd}=\int\mathrm{d}t\mathrm{d}r\Bigg{[}$ $\displaystyle
C_{1}\left(\dot{h}_{1}-h_{0}^{\prime}+\frac{2}{r}h_{0}\right)^{2}+2(C_{2}\beta^{\prime}+C_{3}\beta)\left(\dot{h}_{1}-h_{0}^{\prime}+\frac{2}{r}h_{0}\right)+C_{4}\dot{\beta}^{2}+C_{5}\beta^{\prime
2}$ (59) (61) $\displaystyle+\left(C_{6}h_{1}^{2}+C_{7}h_{0}^{2}+C_{8}\beta
h_{0}\right)+C_{9}\beta^{2}\Bigg{]}\,,$

where the coefficients $C_{i}$ $(i=1,2,\cdots,9)$ are listed as follows,

$\displaystyle C_{1}$ $\displaystyle=\frac{\sqrt{f}(h-\xi
b_{t}^{2})}{h^{3/2}}\,,$ $\displaystyle C_{2}$
$\displaystyle=\frac{\sqrt{f}\xi b_{t}}{\sqrt{h}}\,,$ $\displaystyle C_{3}$
$\displaystyle=\sqrt{\frac{f}{h}}\left(\frac{(\xi-2\kappa)b_{t}^{\prime}}{r}-\frac{2\xi
b_{t}}{r}\right)\,,$ $\displaystyle C_{4}$
$\displaystyle=\frac{2\kappa}{\sqrt{fh}}\,,$ $\displaystyle C_{5}$
$\displaystyle=-2\kappa\sqrt{fh}\,,$ $\displaystyle C_{6}$
$\displaystyle=-\frac{(\ell-1)(\ell+2)\sqrt{fh}}{r^{2}}\,,$ $\displaystyle
C_{7}$ $\displaystyle=\frac{(\ell-1)(\ell+2)(h-\xi
b_{t}^{2})}{\sqrt{fh}hr^{2}}\,,$ $\displaystyle C_{8}$
$\displaystyle=-\frac{(\ell-1)(\ell+2)2\xi b_{t}}{r^{2}\sqrt{fh}}\,,$
$\displaystyle C_{9}$
$\displaystyle=-\frac{2\kappa\ell(\ell+1)}{r^{2}}\sqrt{\frac{h}{f}}-\frac{\xi\big{(}h(-2+2f+rf^{\prime})+rfh^{\prime}\big{)}}{r^{2}\sqrt{fh}}\,.$

####  III.2.1 The $\ell\geq 2$ modes

In this subsection, we consider the odd parity modes with $\ell\geq 2$. Note
that the perturbed action (59) does not contain $\dot{h}_{0}$ term, thus there
are only two dynamical fields, $h_{1}$ and $\beta$. However, since Eq. (59)
involves $h_{0}^{\prime}$, $h_{0}$ cannot be solved directly. By introducing a
Lagrangian multiplier $\chi(t,r)$, Eq. (59) yields

$\displaystyle I^{\prime}_{\rm odd}=\int\mathrm{d}t\mathrm{d}r$
$\displaystyle\Bigg{\\{}C_{1}\bigg{[}-\chi^{2}+2\chi\left(\dot{h}_{1}-h_{0}^{\prime}+\frac{2}{r}h_{0}+\frac{C_{2}\beta^{\prime}+C_{3}\beta}{C_{1}}\right)\bigg{]}-\frac{(C_{2}\beta^{\prime}+C_{3}\beta)^{2}}{C_{1}}+C_{4}\dot{\beta}^{2}+C_{5}\beta^{\prime
2}$ (62) (64) $\displaystyle+\left(C_{6}h_{1}^{2}+C_{7}h_{0}^{2}+C_{8}\beta
h_{0}\right)+C_{9}\beta^{2}\Bigg{\\}}\,.$

It is not difficult to find that the former action (59) can be recovered by
varying $\chi(t,r)$ in the new action (62). Performing variation on
$I^{\prime}_{\rm odd}$ with respect to $h_{0}$ and $h_{1}$ gives

$\displaystyle C_{1}\dot{\chi}-C_{6}h_{1}=0\,,$ (65) (67) $\displaystyle
2C_{1}(\chi^{\prime}-\chi)-\frac{4}{r}\chi-(2C_{7}h_{0}+C_{8}\beta)=0\,.$

Then $I^{\prime}_{\rm odd}$ can be rewritten in the quadratic form

$\displaystyle I^{\prime}_{\rm
odd}=\int\mathrm{d}t\mathrm{d}r\left[-\frac{C_{1}^{2}}{C_{6}}\dot{\chi}^{2}-\frac{C_{1}^{2}}{C_{7}}\chi^{\prime
2}+U_{\chi}\chi^{2}+U_{\chi\beta}\chi\beta+C_{4}\dot{\beta}^{2}+\left(C_{5}-\frac{C_{2}^{2}}{C_{1}}\right)\beta^{\prime
2}+U_{\beta}\beta^{2}\right]\,,$ (68)

or in the matrix form

$I^{\prime}_{\rm
odd}=\int\mathrm{d}t\mathrm{d}r\left(\dot{\vv{\chi}}^{\intercal}\bm{K}\dot{\vv{\chi}}+\vv{\chi}^{\prime\intercal}\bm{G}\vv{\chi}^{\prime}+\vv{\chi}^{\intercal}\bm{M}\vv{\chi}\right)\,.$
(69)

Here $\vv{\chi}$ = $(\chi,\beta)^{\intercal}$ and $\bm{K},\bm{G},\bm{M}$ are
$2\times 2$ matrices

$\displaystyle\bm{K}$
$\displaystyle=\begin{pmatrix}-\frac{C_{1}^{2}}{C_{6}}&&0\\\
0&&C_{4}\end{pmatrix}\,,$ (70) $\displaystyle\bm{G}$
$\displaystyle=\begin{pmatrix}-\frac{C_{1}^{2}}{C_{7}}&&0\\\
0&&C_{5}-\frac{C_{2}^{2}}{C_{1}}\end{pmatrix}\,,$ (71) $\displaystyle\bm{M}$
$\displaystyle=\begin{pmatrix}U_{\chi}&&U_{\chi\beta}\\\
0&&U_{\beta}\end{pmatrix}\,,$ (72)

where

$\displaystyle
U_{\chi}=-C_{1}-\frac{6C_{1}^{2}-r^{2}C_{1}C_{1}^{\prime\prime}}{r^{2}C_{7}}-\frac{rC_{1}C_{1}^{\prime}C_{7}^{\prime}+2C_{1}^{2}C_{7}^{\prime}}{rC_{7}^{2}}\,,$
(73) (75) $\displaystyle
U_{\beta}=C_{9}+\left(\frac{C_{2}C_{3}}{C_{1}}\right)^{\prime}-\frac{C_{3}^{2}}{C_{1}}-\frac{C_{8}^{2}}{4C_{7}}\,,$
(77) $\displaystyle
U_{\chi\beta}=2C_{3}-\frac{2C_{1}C_{8}}{rC_{7}}+C_{1}\left(\frac{C_{8}}{C_{7}}\right)^{\prime}\,.$

We next shall analyze the ghost and gradient (in)stability of $\vv{\chi}$
following the scheme introduced in Sec. II. One can find that the coefficients
of kinetic terms of $\dot{\chi}^{2}$ and $\dot{\beta}^{2}$,

$\displaystyle C_{4}$ $\displaystyle=\frac{2\kappa}{\sqrt{fh}}>0\,,$ (78)
$\displaystyle\frac{C_{1}^{2}}{C_{6}}$
$\displaystyle=-\sqrt{\frac{f}{h}}\frac{r^{2}(h-\xi
b_{t})^{2}}{(\ell+2)(\ell-1)h^{3}}<0\,,$ (79)

yielding

$\bm{K}_{11}>0\,,\quad\quad\det\bm{K}>0\,.$ (80)

It, therefore, shows that for the $\ell\geq 2$ case of the odd parity modes,
the bumblebee BHs have _no ghost instability_.

In the following, we shall investigate the gradient instability of bumblebee
BHs following the scheme introduced in Sec. II. We assume that $\vv{\chi}$ has
a wave solution along the radial direction,

$\vv{\chi}=\vv{\chi}_{k}\mathrm{e}^{{i}(\omega t-kr)}\,,$ (81)

where $\vv{\chi}_{k}$ is a constant vector. We first consider the high-energy
limit along the radial direction. In other words, we take limits
$\omega\to\infty$ and $k\to\infty$, but keeping $\omega/k$ a finite value. In
order for $\vv{\chi}_{k}$ to have a nontrivial solution, we have

$\det\big{(}\omega^{2}\bm{K}+k^{2}\bm{G}\big{)}=0\,.$ (82)

We then define the locally propagating speed of $\vv{\chi}$ along the radial
direction,

$c_{r}\equiv\frac{\mathrm{d}r_{*}}{\mathrm{d}\tau}=\frac{\omega}{k\sqrt{fh}}\,,$

where the proper time $\tau\equiv\int h\mathrm{d}t$ and the tortoise
coordinate satisfies ${\mathrm{d}r_{*}}/{\mathrm{d}r}=f$. Then, Eq. (82) gives
two solutions for the radial speed,

$c_{r_{1}}^{2}=\frac{C_{2}^{2}-C_{1}C_{5}}{C_{1}C_{4}fh}\,,\quad\quad
c_{r_{2}}^{2}=-\frac{C_{6}}{C_{7}fh}\,.$ (83)

Recall that the no-ghost condition implies $C_{4}>0$ and $C_{6}<0$. For the
absence of gradient instability along the radial direction, namely

$c_{r}^{2}>0\,,$ (84)

we have the following constraints of gradient stability for bumblebee BHs,

$\displaystyle\frac{C_{2}^{2}-C_{1}C_{5}}{C_{1}}=$
$\displaystyle\frac{\xi^{2}b_{t}^{2}+2\kappa(h-\xi b_{t}^{2})}{h-\xi
b_{t}^{2}}$ $\displaystyle>0\,,$ (85) (86) $\displaystyle C_{7}=$
$\displaystyle\frac{(\ell-1)(\ell+2)(h-\xi b_{t}^{2})}{r^{2}\sqrt{f}h^{3/2}}$
$\displaystyle>0\,.$ (87)

It is easy to verify that when $h(r)>\xi b_{t}(r)^{2}$, namely that
$c_{r_{1}}^{2}>0$ is satisfied, $c_{r_{2}}^{2}>0$ is also satisfied.
Therefore, for bumblebee BHs, the condition for the absence of gradient
instability associated with the radial direction is

$h(r)-\xi b_{t}(r)^{2}>0\,.$ (88)

In addition to the radial direction, instability may also arise along the
angular direction, depending on the propagating speed of $\vv{\chi}$. Now we
consider the high-energy limit but along the angular direction. In other
words, we take $\omega\to\infty$ and $\ell\to\infty$, but keep $\omega/\ell$
being finite. This limit indicates that the perturbation field, $\vv{\chi}$,
propagates with large energy and angular momentum but the speed of propagation
is finite. For nontrivial $\vv{\chi}_{\ell}$ solution, we have

$\det\big{(}\omega^{2}\bm{K}+\bm{M}\big{)}=0\,.$ (89)

In the large-$\ell$ limit, the speed of propagation along the angular
direction is defined as

$c_{\Omega}\equiv
r\frac{\mathrm{d}\theta}{\mathrm{d}\tau}=r\frac{\mathrm{d}\theta}{\sqrt{h}\mathrm{d}t}=\frac{r\omega}{\sqrt{h}\ell}\,.$

With this definition and in the large-$\ell$ limit, $c_{\Omega}^{2}$ can be
solved from Eq. (89) as

$\displaystyle c_{\Omega_{1}}^{2}$ $\displaystyle=\frac{h}{(h-\xi
b_{t}^{2})}\,,$ (90) $\displaystyle c_{\Omega_{2}}^{2}$
$\displaystyle=\frac{\xi^{2}b_{t}^{2}+2\kappa(h-\xi b_{t}^{2})}{2\kappa(h-\xi
b_{t}^{2})}\,.$ (91)

We can see that the no tachyonic instability condition associated with the
angular speed, namely,

$c_{\Omega}^{2}>0\,,$ (92)

shares the same constraints given by the no gradient instability condition
associated with the radial speed in Eq. (88).

####  III.2.2 The $\ell=1$ modes

We have discussed the instabilities of bumblebee BHs for $\ell\geq 2$ odd-
parity modes. Now we separately discuss the dipole mode, namely, $\ell=1$. As
we discussed previously, when $\ell=1$, the gauge degree of freedom is not
fixed in Eq. (58) since $h_{ab}$ identically vanishes.

Recalling the gauge transformation rules given in Eq. (52), one can simplify
the gravitational perturbation by choosing a gauge with $h_{1}=0$, implying
that

$\Lambda(t,r)=-r^{2}\int\frac{h_{1}(t,r^{\prime})}{r^{\prime
2}}\mathrm{d}r^{\prime}+r^{2}{\cal F}(t)\,,$ (93)

where ${\cal F}(t)$ is an arbitrary function with respect to $t$. Performing
variation on the action (59) with respect to $\beta$ and $h_{0}$, and setting
$h_{1}=0$, we have

$\dot{\Xi}=0\,,\quad\quad\Xi^{\prime}=0\,,$ (94)

where

$\Xi\equiv
2C_{1}\left(h_{0}^{\prime}-\frac{2}{r}h_{0}\right)-2(C_{2}\beta^{\prime}+C_{3}\beta)\,.$
(95)

It yields

$\Xi={\cal C}_{1}\,,$ (96)

where ${\cal C}_{1}$ is an integration constant. Then $h_{0}$ can be solved as

$h_{0}=r^{2}{\cal F}_{2}(t)+r^{2}\int\mathrm{d}\tilde{r}\frac{{\cal
C}_{1}\tilde{r}+2C_{4}\beta+2C_{3}\tilde{r}\beta^{\prime}}{2C_{1}\tilde{r}^{3}}\,,$
(97)

where ${\cal F}_{2}(t)$ is a gauge mode which can be eliminated by setting
${\cal F}(t)=\int\mathrm{d}t{\cal F}_{2}(t)$. In addition, if shutting down
the vector perturbation with $\beta=0$, namely that $h_{0}$ does not depend on
time, then ${\cal C}_{1}$ is related to the angular momentum of a slowly
rotating BH Kobayashi  _et al._ (2012); Ogawa  _et al._ (2016).

After integrating by parts, action (59) then reduces to

$\displaystyle I^{\prime}_{\rm odd}=$
$\displaystyle\int\mathrm{d}t\mathrm{d}r\left\\{C_{4}\dot{\beta}^{2}-\left(C_{5}-\frac{C_{2}^{2}}{C_{1}}\right)\beta^{\prime
2}+\left[C_{9}+\left(\frac{C_{2}C_{3}}{C_{1}}\right)^{\prime}-\frac{C_{3}^{2}}{C_{1}}\right]\beta^{2}\right\\}\,,$
(98)

where we set ${\cal C}_{1}=0$ for convenience. The reduced action (98) shows
that only one dynamical field, $\beta(t,r)$, propagates. It is obvious that
for $\ell=1$, a bumblebee BH has _n_ o ghost instability because

$C_{4}=\frac{2\kappa}{\sqrt{fh}}>0\,.$ (99)

Since $\ell$ has been fixed, one can only read $c_{r}^{2}$ in the high-energy
limit along the radial direction and it is easy to find that the _no gradient
instability_ condition is

$C_{5}-\frac{C_{2}^{2}}{C_{1}}=\frac{\xi^{2}b_{t}^{2}+2\kappa(h-\xi
b_{t}^{2})}{h-\xi b_{t}^{2}}>0\,,$ (100)

which is consistent with condition (85). To summarize our results, for the
odd-parity modes, the condition of _no gradient instability_ for bumblebee BHs
is

$h(r)-\xi b_{t}(r)^{2}>0\,.$ (101)

In addition to the ghost instability, gradient instability, and tachyonic
instability, instabilities associated with the quasinormal modes (QNMs) of BHs
may arise when solving the equations of motion of $\vv{\chi}$ under
appropriate boundary conditions. However, it is out of the scope of this work.
We leave it to our future work to analyze the QNMs of bumblebee BHs.

###  III.3 Numerical results

In the previous section, we find that a bumblebee BH has no ghost instability
but may have gradient instability and tachyonic instability, and obtain the
condition (101) for avoiding them. In the following, we numerically
investigate the condition and give constraints on the vector charge of
bumblebee BHs.

Based on Eq. (101), a direct conclusion is that for $\xi\leq 0$,

$h(r)-\xi b_{t}(r)^{2}>0$ (102)

is always satisfied, indicating that these bumblebee BHs (including the RN BH)
have no gradient instability nor tachyonic instability.

For $\xi>0$, the condition becomes complicated so we shall discuss it case by
case. We begin with $\xi=2\kappa$ since then the bumblebee BH has an
analytical solution, the Schwarzschild metric with a nontrivial vector field
given in Eq. (30). Equation (101) follows that

$h-\xi
b_{t}^{2}=\frac{(r-2M)\big{[}2MQ^{2}+r(M^{2}-Q^{2})\big{]}}{r^{2}M^{2}}>0\,,$
(103)

implying that when $M>Q$, $h-\xi b_{t}^{2}$ is always positive, while $Q>M$,
$h-\xi b_{t}^{2}$ is positive only when $r$ is small. However, there must
exist a critical radius,

$r_{c}=\frac{2MQ^{2}}{Q^{2}-M^{2}}\,,$ (104)

where $h-\xi b_{t}^{2}$ changes the sign, implying that instabilities occur.
If $Q=M$, $r_{c}$ is located at the infinity. We thus easily draw a simple
conclusion that when $\xi=2\kappa$, the condition for no gradient instability
nor tachyonic instability for a stealth Schwarzschild BH is

$M\geq Q\,.$ (105)

It is a conclusion that, while concise, is rich in physical significance. It
implies that even though stealth Schwarzschild BHs admit solutions with $Q>M$,
the conditions of no instabilities indicate that such solutions cannot exist
stably. Similar phenomena where Schwarzschild BHs with nontrivial vector field
suffer the gradient instability also appear in the Einstein-aether theory and
other vector-tensor theories Tsujikawa  _et al._ (2021); Kase  _et al._
(2018a).

For a general $\xi$, there is yet no analytical but numerical BH solutions.
Following the hints given by the case of $\xi=2\kappa$, we thus choose a
gradient stability indicator to show whether these numerical bumblebee BHs
satisfy the condition (101). Note that

$h(r)\big{|}_{r\to\infty}=1\,,\quad\quad
b_{t}(r)\big{|}_{r\to\infty}=\mu_{\infty}\,,$ (106)

we thus choose the gradient stability indicator

$\left.h(r)-\xi b_{t}^{2}\right|_{r\to\infty}=1-\xi\mu_{\infty}^{2}$ (107)

to show that when $1-\xi\mu_{\infty}^{2}>0$, bumblebee BHs have no gradient
instability. Recall the example of $\xi=2\kappa$, we have

$1-\xi\mu_{\infty}^{2}=\frac{M^{2}-Q^{2}}{M^{2}}$ (108)

showing that the result of gradient stability for stealth Schwarzschild BHs is
consistent with Eq. (105). In addition to the case of $\xi=2\kappa$, we find
that $1-\xi\mu_{\infty}^{2}$ changes sign when $\xi\geq\xi_{c}$, where

$0.5\kappa<\xi_{c}<0.51\kappa,$ (109)

indicating that bumblebee BHs begin to suffer gradient instability when
$\xi\geq\xi_{c}$ (see Table 2).

Table 2: Values of $|Q|/M$ and $1-\xi\mu_{\infty}^{2}$ of bumblebee BHs for $\xi\simeq 0.51\kappa$. It shows that $1-\xi\mu_{\infty}^{2}$ changes sign when $\xi=0.51\kappa$, indicating the critical point when the gradient instability happens.  $\xi/\kappa$ | $0.48$ | $0.49$ | $0.50$ | $0.51$ | $0.52$ | $0.53$
---|---|---|---|---|---|---
$|Q|/[M]$ | $1.00016$ | $1.00227$ | $1.00193$ | $1.00187$ | $1.00216$ | $1.00287$
$1-\xi\mu_{\infty}^{2}$ | $0.040$ | $0.024$ | $0.0039$ | $-0.016$ | $-0.036$ | $-0.055$

Here we practically assume that when $\xi<\xi_{c}$, $h-\xi b_{t}^{2}$ is a
monotonically increasing function with respect to $r$. Since when $r\to
r_{h}$, one has $h-\xi b_{t}^{2}\to 0$, such an assumption implies that $h-\xi
b_{t}^{2}$ will always be positive for $r>r_{h}$ when $\xi<\xi_{c}$. To verify
it numerically, we show $h-\xi b_{t}^{2}$ against $r$ for various $\xi$ around
$\xi_{c}$ in Fig. 2.

Figure 2: $h-\xi b_{t}^{2}$ against $r$ for various $\xi$ around $\xi_{c}$. We
verify that when $\xi<\xi_{c}$, $h-\xi b_{t}^{2}$ is a monotonically
increasing function with respect to $r$. Figure 3: The relation between the
gradient (in)stability indicator, $1-\xi\mu_{\infty}^{2}$, and the charge of
bumblebee BHs.

Recall that bumblebee BHs can be described by two parameters, $Q$ and $M$. We
plot the gradient (in)stability indicator $1-\xi\mu_{\infty}^{2}$ with respect
to $|Q|$ in Fig. 3. In the figure, we use the same $\xi$ as in Fig. 1. When
$\xi=-2\kappa$ and $0$, the gradient stability indicators are always positive
for all $|Q|$, which is consistent with our previous statement that a
bumblebee BH has no gradient instability when $\xi\leq 0$. Moreover, in Fig.
3, we denote the critical points where $1-\xi b_{t}^{2}$ changes sign with red
dots and the maximum of $|Q|$ with black dots. For $\xi<\xi_{c}$, there does
not exist any critical point (red dot), showing that a bumblebee BH does not
suffer from the gradient instability. We found that the critical points (red
dots) appear when $\xi>\xi_{c}$. When $\xi$ increases but is less than
$1.4\kappa$, there exist two bumblebee BH solutions for the same $|Q|$. When
$\xi$ is close to but larger than $\xi_{c}$, only one of the two BH solutions
suffers gradient instability (see the green solid line for $0.6\kappa$ in Fig.
3). The critical points first move to the right, and then to the left after
the critical points coincide with the maximum points (the black dots). When
$0.75\kappa<\xi<1.4\kappa$, both of the two branches of bumblebee BHs suffer
the gradient instability for larger $|Q|$ (see the orange solid line for
$\xi=\kappa$ in Fig. 3 as an example). Furthermore, when $\xi>1.4\kappa$,
there is only one bumblebee BH solution for one $|Q|$, and the bumblebee BH
suffers the gradient instability when $1-\xi\mu^{2}_{\infty}$ is negative (see
the case of $\xi=2\kappa$ in Fig. 3).

Figure 4: Constraints on the $\xi$-$|Q|$ plane from EHT observations and the
gradient instability. The red region shows the parameter space where only one
of two branches of BH solutions is unstable, while the blue region shows that
both two branches of BHs are gradient unstable. The red solid line gives the
boundary of the gradient instability. The orange and green dot-dashed lines
show the bounds given by the observations of the shadows of supermassive BHs
in M87 and Milky Way respectively Xu  _et al._ (2023b). The black lines give
the boundary of the existence domain of BH solutions. The region enclosed by
the black and blue lines indicates that there exist two branches of BHs for
the same $|Q|$ and $\xi/\kappa$.

As mentioned in the Introduction, Xu  _et al._ (2023b, a) have put
observational bounds on $|Q|$ associated with bumblebee BHs for various $\xi$
using the supermassive BH shadows of M87∗ and Sgr ${\rm A}^{*}$ from the EHT.
Based on our results of the gradient stability for bumblebee BHs, in Fig. 4 we
give an extended plot of Fig. 3 for more choices of $\xi$. The black line,
formed by connecting the maximum $|Q|$ of bumblebee BHs, serves as the
boundary of the BH solutions. The blue line serves as the boundary of the
existence of two branches of solutions for the same $|Q|$ and $\xi$. On the
other hand, the green and orange dot-dashed lines are given by the
observations of the shadows for the supermassive BHs in M87 and the Milky Way
respectively. The red line, formed by connecting the critical points for
different values of $\xi$ (some of which are denoted as red dots in Fig. 4)
serves as the boundary of the gradient stable and unstable regions. Therefore,
the colored region represents bumblebee BHs suffering the gradient
instability.

In addition, due to the bound given from observations on shadows for the
supermassive BHs, we find that when $\xi<2\kappa$, even though the bumblebee
BHs do not suffer the gradient instability, $|Q|$ of bumblebee BHs cannot be
sufficiently large since such BH solutions do not exist. Specifically, when
$\xi$ is negative, $|Q|<M$ is always satisfied. Nevertheless, when
$\xi>2\kappa$, it admits bumblebee BHs carrying a considerable vector charge
$|Q|$. However, the bumblebee BHs with $|Q|>M$ can rarely exist due to
gradient instability. Therefore, we draw a conclusion that the bumblebee BHs
could not carry considerable vector charge $|Q|$, and further when $\xi\geq
2\kappa$ or $\xi\leq 0$—equivalently $\xi(\xi-2\kappa)>0$—the vector charge of
a bumblebee BH cannot be larger than its mass,

$|Q|<M\,.$ (110)

##  IV Conclusions

In this paper, we have studied the dynamic (in)stability of bumblebee BHs in
the bumblebee vector-tensor theory. Treating the BH spacetime as the
background, we have considered perturbations of the gravitational field and
the bumblebee field with odd parity and have investigated the associated
(in)stability. Our conclusions can be summarized as follows.

  1. (I)

Bumblebee BHs do not suffer the ghost instability.

  2. (II)

Bumblebee BHs suffer the gradient/tachyonic instability when $\xi>\xi_{c}$,
where the critical coupling constant $\xi_{c}$ is between $0.5\kappa$ and
$0.51\kappa$ according to our numerical results.

  3. (III)

The conditions for avoiding the gradient/tachyonic instability for bumblebee
BHs indicate a no-go theorem on the bumblebee charge. Specifically, when
$\xi(\xi-2\kappa)\geq 0$, the vector charge of a bumblebee BH can not be
larger than its mass.

As shown in Fig. 4, the theoretical condition for no gradient/tachyonic
instability gives a stronger constraint on the charge of bumblebee BHs than
the constraints given by the observed BH images from EHT when $\xi\gtrsim
1.3\kappa$. In general, Fig. 4 also shows that if $\xi(\xi-2\kappa)\geq 0$,
bumblebee BHs with $|Q|>M$ can hardly exist in our Universe due to the
gradient/tachyonic instability. This seems like a realization of the weak
cosmic censorship conjecture in the bumblebee theory. In GR, for a RN BH, the
weak cosmic censorship conjecture states that the charge cannot be larger than
the mass, to avoid the naked singularity Hawking and Penrose (1970); Penrose
(1974). Our finding suggests that the gradient/tachyonic instability might be
the physical mechanism behind the weak cosmic censorship, if we extend the
Einstein-Maxwell theory to the bumblebee vector-tensor theory, potentially
providing a physical mechanism for the weak cosmic censorship in a broader
context of gravity theories. As the gradient/tachyonic instability also exists
in other modified gravity theories, we expect there to be also such worthwhile
constraints on the corresponding additional hairs of BHs in these theories. In
addition to the ghost instability, gradient instability, and tachyonic
instability of odd-parity perturbations considered in this work, instabilities
related to the QNMs of bumblebee BHs, together with the even-parity
perturbations and their stability properties, are also worthwhile topics for
the stability analysis of bumblebee BHs. We leave them for future studies.

###### Acknowledgements.

We thank Zexin Hu, Yingli Zhang, Jinbo Yang and Run-Qiu Yang for useful
discussions. This work was supported by the National Natural Science
Foundation of China (Grants No. 12247128, No. 11991053, No. 11975027, No.
11721303), the China Postdoctoral Science Foundation (No. 2021TQ0018, No.
2023M741999), the National SKA Program of China (No. 2020SKA0120300), the Max
Planck Partner Group Program funded by the Max Planck Society, and the High-
Performance Computing Platform of Peking University.

## References

  * Kostelecký and Samuel (1989) V. A. Kostelecký and S. Samuel, Phys. Rev. Lett. 63, 224 (1989).
  * Kostelecký (2004) V. A. Kostelecký, Phys. Rev. D 69, 105009 (2004).
  * Kostelecký and Mewes (2009) V. A. Kostelecký and M. Mewes, Phys. Rev. D 80, 015020 (2009).
  * Bailey and Kostelecký (2006) Q. G. Bailey and V. A. Kostelecký, Phys. Rev. D 74, 045001 (2006), arXiv:gr-qc/0603030 .
  * Kostelecký and Mewes (2012) V. A. Kostelecký and M. Mewes, Phys. Rev. D 85, 096005 (2012).
  * Kostelecký and Mewes (2013) V. A. Kostelecký and M. Mewes, Phys. Rev. D 88, 096006 (2013).
  * Kostelecký and Li (2019) V. A. Kostelecký and Z. Li, Phys. Rev. D 99, 056016 (2019).
  * Bertolami and Paramos (2005) O. Bertolami and J. Paramos, Phys. Rev. D 72, 044001 (2005), arXiv:hep-th/0504215 .
  * Casana  _et al._ (2018) R. Casana, A. Cavalcante, F. P. Poulis, and E. B. Santos, Phys. Rev. D 97, 104001 (2018), arXiv:1711.02273 [gr-qc] .
  * Güllü and Övgün (2022) I. Güllü and A. Övgün, Annals Phys. 436, 168721 (2022), arXiv:2012.02611 [gr-qc] .
  * Izmailov and Nandi (2022) R. N. Izmailov and K. K. Nandi, Class. Quant. Grav. 39, 215006 (2022).
  * Liu  _et al._ (2023) W. Liu, X. Fang, J. Jing, and J. Wang, Eur. Phys. J. C 83, 83 (2023), arXiv:2211.03156 [gr-qc] .
  * Xu  _et al._ (2023a) R. Xu, D. Liang, and L. Shao, Phys. Rev. D 107, 024011 (2023a), arXiv:2209.02209 [gr-qc] .
  * Xu (2023) R. Xu, in  _9th Meeting on CPT and Lorentz Symmetry_ (2023) arXiv:2301.12666 [gr-qc] .
  * Xu  _et al._ (2023b) R. Xu, D. Liang, and L. Shao, Astrophys. J. 945, 148 (2023b), arXiv:2302.05671 [gr-qc] .
  * Liang  _et al._ (2023) D. Liang, R. Xu, Z.-F. Mai, and L. Shao, Phys. Rev. D 107, 044053 (2023), arXiv:2212.09346 [gr-qc] .
  * Hu  _et al._ (2023) Z. Hu, L. Shao, R. Xu, D. Liang, and Z.-F. Mai, (2023), arXiv:2312.02486 [astro-ph.HE] .
  * Mai  _et al._ (2023) Z.-F. Mai, R. Xu, D. Liang, and L. Shao, Phys. Rev. D 108, 024004 (2023), arXiv:2304.08030 [gr-qc] .
  * Moncrief (1974) V. Moncrief, Phys. Rev. D 10, 1057 (1974).
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
  * De Felice  _et al._ (2011) A. De Felice, T. Suyama, and T. Tanaka, Phys. Rev. D 83, 104035 (2011), arXiv:1102.1521 [gr-qc] .
  * Motohashi and Suyama (2011) H. Motohashi and T. Suyama, Phys. Rev. D 84, 084041 (2011), arXiv:1107.3705 [gr-qc] .
  * Kobayashi  _et al._ (2012) T. Kobayashi, H. Motohashi, and T. Suyama, Phys. Rev. D 85, 084025 (2012), [Erratum: Phys.Rev.D 96, 109903 (2017)], arXiv:1202.4893 [gr-qc] .
  * Kobayashi  _et al._ (2014) T. Kobayashi, H. Motohashi, and T. Suyama, Phys. Rev. D 89, 084042 (2014), arXiv:1402.6740 [gr-qc] .
  * Ogawa  _et al._ (2016) H. Ogawa, T. Kobayashi, and T. Suyama, Phys. Rev. D 93, 064078 (2016), arXiv:1510.07400 [gr-qc] .
  * Ganguly  _et al._ (2018) A. Ganguly, R. Gannouji, M. Gonzalez-Espinoza, and C. Pizarro-Moya, Class. Quant. Grav. 35, 145008 (2018), arXiv:1710.07669 [gr-qc] .
  * Heisenberg  _et al._ (2017) L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa, JCAP 08, 024 (2017), arXiv:1706.05115 [gr-qc] .
  * Kase  _et al._ (2018a) R. Kase, M. Minamitsuji, S. Tsujikawa, and Y.-L. Zhang, JCAP 02, 048 (2018a), arXiv:1801.01787 [gr-qc] .
  * Kase  _et al._ (2018b) R. Kase, M. Minamitsuji, and S. Tsujikawa, Phys. Lett. B 782, 541 (2018b), arXiv:1803.06335 [gr-qc] .
  * Gannouji and Baez (2022) R. Gannouji and Y. R. Baez, JHEP 02, 020 (2022), arXiv:2112.00109 [gr-qc] .
  * Kase and Tsujikawa (2023) R. Kase and S. Tsujikawa, Phys. Rev. D 107, 104045 (2023), arXiv:2301.10362 [gr-qc] .
  * Tsujikawa  _et al._ (2021) S. Tsujikawa, C. Zhang, X. Zhao, and A. Wang, Phys. Rev. D 104, 064024 (2021), arXiv:2107.08061 [gr-qc] .
  * Demirboğa  _et al._ (2022) E. S. Demirboğa, A. Coates, and F. M. Ramazanoğlu, Phys. Rev. D 105, 024057 (2022), arXiv:2112.04269 [gr-qc] .
  * Doneva  _et al._ (2022) D. D. Doneva, F. M. Ramazanoğlu, H. O. Silva, T. P. Sotiriou, and S. S. Yazadjiev, (2022), arXiv:2211.01766 [gr-qc] .
  * Hawking and Penrose (1970) S. W. Hawking and R. Penrose, Proc. Roy. Soc. Lond. A 314, 529 (1970).
  * Penrose (1974) R. Penrose, in  _Gravitational Radiation and Gravitational Collapse_ , Vol. 64, edited by C. Dewitt-Morette (1974) p. 82.