# New Non-Abelian Reissner-Nordström Black Hole Solutions in the Generalized
SU(2) Proca Theory And Some Astrophysical Implications

Gabriel Gómez<EMAIL_ADDRESS>Departamento de Física, Universidad de
Santiago de Chile,
Avenida Víctor Jara 3493, Estación Central, 9170124, Santiago, Chile  José F.
Rodríguez<EMAIL_ADDRESS>Escuela de Física, Universidad
Industrial de Santander, Ciudad Universitaria, Bucaramanga 680002, Colombia
ICRANet, Piazza della Repubblica 10, 65122, Pescara PE, Italy

###### Abstract

The Generalized SU(2) Proca theory is a vector-tensor theory of gravity whose
action is invariant under global transformations of the SU(2) group and
includes second-order derivative self-interactions of the vector field beyond
the massive Yang-Mills theory. We find, in particular, that the presence of
two Lagrangian pieces consisting of four gauge fields minimally coupled to
gravity gives rise to an exact Reissner-Nordström black hole solution endowed
with two different non-Abelian effective charges that depend on the specific
combination, $\chi=2\chi_{1}+\chi_{2}$, of the respective coupling constants.
After studying the spacetime structure of the black hole, which allows us to
characterize the parameter space that preserves the weak cosmic censorship
conjecture, some astrophysical implications of the black hole solutions are
investigated. First, joint analysis of observations of the EHT’s first images
of Sagittarius A⋆ of our Galaxy and the Keck telescope set the first serious
constraint on the free parameters of the theory beyond the theoretical bounds
found. Second, we investigate the accretion properties of spherical steady
flows around this class of non-Abelian Reissner-Nordström black hole.
Specifically, we examine the general conditions under which transonic flow is
allowed. An analytical solution for critical accretion is found in terms of
the coupling constant. In addition, we explore the effect of changing $\chi$
on the radial velocity and mass density numerically and show how the extremal
Reissner-Nordström and the standard Schwarzschild solutions as limit cases are
achieved. Finally, working in the fully relativistic regime, an analytical
expression for the critical mass accretion rate of a polytropic fluid onto a
black hole is derived. As a main result, we find that the critical accretion
rate efficiency can be noticeably improved compared to the Schwarzschild case
for a specific region of the parameter space where the non-Abelian charge
becomes imaginary.

##  I Introduction

Black holes (BHs) are among the most fascinating objects in the Universe,
resulting from the gravitational collapse of massive objects, as predicted by
general relativity (GR), see e.g. Poisson (2009). Apart from the fundamental
conceptions and interesting properties they harbor, BHs are ideal
laboratories, due to their intense gravitational fields, to study high-energy
astrophysical processes occurring in their vicinity Shapiro and Teukolsky
(1983). Furthermore, observations of BHs in the strong field regime provide a
unique opportunity to study the properties of spacetime and understand the
nature of gravity in extreme environments. This is the primary program of
theories beyond GR, which aim to predict deviations from Einstein’s theory.
Motivated by these concerns, BHs have been the central target of current
astrophysical experiments, including the Event Horizon Telescope (EHT) and the
Very Large Telescope global networks Akiyama  _et al._ (2019, 2022), GRAVITY
collaboration Abuter  _et al._ (2020), and the LIGO-Virgo collaboration Abbott
_et al._ (2017a, b) among others. While the predictions of GR are well
consistent with all available observational data within the current
uncertainties Will (2014), some theories beyond Einstein’s theory can also
explain the observed phenomena Abbott  _et al._ (2021); Cardoso and Pani
(2019); Psaltis  _et al._ (2020); Vagnozzi  _et al._ (2022)). Therefore,
current and future measurements of BHs at the event-horizon scale pose a
challenge for theories beyond GR. Such observations, however, do provide
strong evidence about the existence of BHs.

BHs can possess electric and magnetic charges as described by the Reissner-
Nordström (RN) solution in the Einstein-Maxwell theory Reissner (1916);
Nordström (1918). However, it is widely believed that astrophysical BHs are
electrically neutral due to charge neutralization by astrophysical plasma,
among other suitable physical mechanisms. Alternatively, it is possible for
BHs to carry, for instance, U(1) charge instead of electromagnetic charge due
to mechanisms in the early Universe within the dark (hidden) sector with no
coupling to Standard-Model particles (see, e.g., De Rujula  _et al._ (1990);
Cardoso  _et al._ (2016)). See also Zajaček and Tursunov (2019) for a
discussion from an astrophysical point of view. Regardless of the underlying
physical process behind the charge mechanism, this remains an open issue that
has recently received significant attention after measurements made by the EHT
of the supermassive BH M87⋆ shadow size, and the detection of gravitational
waves from compact object binaries Kocherlakota  _et al._ (2021); Bozzola and
Paschalidis (2021); Liu  _et al._ (2020); Christiansen  _et al._ (2021); Wang
_et al._ (2021); Benavides-Gallego and Han (2022). These discoveries demand a
careful examination of BH charges beyond academic considerations. See also
Zakharov (2014); Zajaček  _et al._ (2018) for observational limits on the
charge of the Galactic Center BH.

Moreover, several observations have confirmed that BHs must rotate to account
for various astrophysical phenomena, such as X-rays streaming off material
near BHs due to the formation of an accretion disk (see e.g., Ref. Done  _et
al._ (2007)). On the other hand, the RN solution provides a useful first
approximation for studying realistic and complex phenomena in the presence of
electric charge, such as accretion and shadow phenomena. The charge of the RN
black hole plays a significant role in the event horizon structure, similar to
the role of spin. However, it is important to note that RN and Kerr black
holes are fundamentally different scenarios.

The electric and magnetic charges of the RN black hole can have similar
effects as the spin parameter of a rotating Kerr black hole, which has been
observed in the magnetar J1745-2900 orbiting around the supermassive BH
Sagittarius A⋆ Juraeva  _et al._ (2021). In addition, recent observations on
the motion of S-stars have constrained the spin of the central object in the
Milky Way to be relatively small, with $a/M\lesssim 0.1$ Fragione and Loeb
(2020). This suggests that even though BHs are expected to have significant
spin, _there exist objects with small angular momentum_ that can be described
by a stationary and spherically symmetric spacetime Misner  _et al._ (2017).

After the discovery of (purely magnetic) static spherically symmetric non-
Abelian BH solutions in the Einstein SU(2) Yang-Mills (EYM) model Bizon
(1990); Volkov and Galtsov (1989); Kuenzle and Masood-ul Alam (1990), it was
soon demonstrated that they are perturbatively unstable Lee  _et al._ (1992).
To address this problem, higher-order curvature terms of the gauge field Radu
and Tchrakian (2012); Mazharimousavi and Halilsoy (2009), as well as non-
trivial combinations with other theories, have been introduced into the
gravitational sector (see, e.g., Ref. Volkov and Gal’tsov (1999) for previous
proposals). Furthermore, some BH solutions with non-Abelian hairs have been
found in theories beyond the canonical Yang-Mills theory but still within the
framework of GR Herdeiro  _et al._ (2017); Radu and Tchrakian (2012);
Mazharimousavi and Halilsoy (2009). In addition, higher curvature terms of the
metric tensor, such as $f(R)$ gravity coupled to the Yang-Mills field, also
admit BH solutions with single or double horizons Mazharimousavi and Halilsoy
(2011).

The EYM case is interesting as there exists a RN solution Harnad  _et al._
(1980); Bartnik and McKinnon (1988). However, unlike the Einstein-Maxwell
case, the EYM RN solution is unstable, indicating that although both models
share the same spacetime configuration, they are perturbatively different. The
same is true in the Einstein-Yang-Mills-Higgs case, where a RN solution exists
but can be stabilized by the Higgs mechanism Breitenlohner  _et al._ (1995).

The quest to construct classical theories beyond Einstein’s theory is an
ongoing and active area of research, fueled by the need to address
longstanding issues related to singularities Penrose (1965); Hawking and
Penrose (1970) and renormalization Deser  _et al._ (1975) in GR. While
Einstein’s theory of gravity provides an effective description of the
gravitational interaction, it is only valid up to a certain cutoff scale
before it loses its regime of validity Donoghue (1994); Burgess (2004). If
this breakdown occurs, for instance, in the strong gravity regime, then
modified gravity theories may play a significant role in describing the
behavior of compact objects such as BHs and neutron stars. More importantly,
these theories have the potential to provide a more complete and accurate
description of the gravitational interaction at these scales, which is crucial
for a deeper understanding of the nature of gravity and the behavior of
astrophysical objects.

One promising approach to modify Einstein’s theory is to introduce new
gravitational degrees of freedom (see e.g., Heisenberg (2019)). The simplest
example is the Horndeski theory Horndeski (1974), which introduces a scalar
field and yields field equations that are, at most, of second order. This is
crucial for avoiding the Ostrogradski ghost Ostrogradsky (1850). The
Generalized Proca theory Heisenberg (2014); Allys  _et al._ (2016a, b);
Gallego Cadavid and Rodriguez (2019), is a vector-tensor version of the
Horndeski theory, where the internal gauge symmetry of the vector field has
been abandoned to allow for the existence of extra terms Allys  _et al._
(2016c); Gómez and Rodríguez (2019); Gallego Cadavid  _et al._ (2020, 2022).
Adding a global SU(2) internal symmetry to the Generalized Proca leads to what
we called the Generalized SU(2) Proca (GSU2P) theory Allys  _et al._ (2016c);
Gómez and Rodríguez (2019); Gallego Cadavid  _et al._ (2020, 2022), whose
action is invariant under diffeomorphisms and globally invariant under the
SU(2) group transformations. However, the equivalence between these two
vector-tensor theories is not straightforward due to the non-Abelian nature of
the GSU2P theory, leading to the presence or absence of new terms, which can
potentially result in new phenomenology Garnica  _et al._ (2022). Although the
GSU2P theory was first formulated in Allys  _et al._ (2016c) by imposing a
primary constraint-enforcing relation to eliminate the non-physical degree of
freedom from the vector field, a secondary constraint-enforcing relation was
required to close the constraint algebra Gallego Cadavid  _et al._ (2020).
Therefore, all Lagrangian building blocks in the GSU2P theory were constructed
to ensure the propagation of the correct number of physical degrees of
freedom, thereby avoiding the Ostrogradski instability.

The aim of this paper is to explore the astrophysical implications of the
GSU2P theory, which provides a useful framework for investigating the theory
at that relevant scale. The field equations in the GSU2P theory are highly
non-trivial, which makes it challenging to obtain analytical solutions for the
entire theory. Therefore, our focus in this study is to examine the
phenomenological aspects of individual Lagrangian pieces, building upon our
previous work on particle-like solutions Martinez  _et al._ (2023).
Specifically, we delve into the analysis of two Lagrangian terms that involve
quartic order non-derivative self-interactions of the gauge field. We choose
these Lagrangian terms because they enable the derivation of exact analytical
solutions for BH, unlike the other Lagrangian terms of the theory. Moreover,
these solutions are distinguished by possessing an effective global charge, as
was unveiled in our recent paper Martinez  _et al._ (2023). Interestingly, the
solutions we obtained correspond to a Reissner-Nordström solution with a non-
Abelian magnetic charge. This magnetic charge adds an intriguing aspect to the
traditional Reissner-Nordström solution, highlighting the presence of non-
Abelian gauge fields and their impact on the properties of BHs. We plan to
report on numerical BH solutions for the derivative self-interaction terms in
a separate work.

In this paper, we report an analytical exact Reissner Nordström solution
within the GSU2P theory. The solutions are characterized by two non-Abelian
effective charges that depend on the corresponding coupling constants of the
Lagrangian pieces. While this solution shares the same RN spacetime structure
as the standard solution, there is a crucial difference: the charge in our
solution is not of electromagnetic origin and can even be imaginary, resulting
in negative energy density. This leads to interesting astrophysical
implications that we uncover in this work. It is worthwhile mentioning that
our findings are rooted in modified theories of gravity, although they were
particularly derived from Lagrangian pieces minimally coupled to gravity.

Accretion processes of ideal and polytropic fluids onto black holes have been
an area of intense study in astrophysics (see e.g Bondi (1952); Michel (1972);
Richards  _et al._ (2021); Aguayo-Ortiz  _et al._ (2021)), serving as a probe
of concept in both the context of GR and more general frameworks. In
particular, accretion flows in an arbitrary spacetime have been extensively
studied as valuable astrophysical probes for detecting any deviations from GR
and testing alternative theories of gravity.Bauer  _et al._ (2022); Salahshoor
and Nozari (2018); Feng  _et al._ (2022); Zuluaga and Sánchez (2021); Ditta
and Abbas (2020); Uniyal  _et al._ (2022); Chakhchi  _et al._ (2022); John
(2019); Pérez  _et al._ (2017); Liu  _et al._ (2022); Heydari-Fard  _et al._
(2021); Stashko  _et al._ (2021); Shaikh and Joshi (2019); Van Aelst  _et al._
(2021).

The paper is structured as follows. In Section II, we introduce the model and
derive in detail an exact non-Abelian RN BH solution in terms of the coupling
constants of the theory. Some properties of the BH solutions, such as the
event horizon, photon sphere, and shadow, are studied. In particular,
observational data of the EHT’s first images of Sagittarius A⋆ are used to
infer the first constraints on the effective coupling constant. In Section
III, we present a general description of the hydrodynamics equations of the
accretion flow, and in Section IV, we calculate the critical accretion rate
for isothermal and polytropic fluids using analytical and numerical
computations. In Section V, we discuss the main findings and possible
extensions of the work, along with further observational constraints on the
theory that can be used in the future. Throughout the manuscript, Latin
indices are internal SU(2) group indices and run from 1 to 3, while Greek
indices stand for spacetime indices and run from 0 to 3. We use geometrized
units with $c=G=1$.

##  II Reissner Nordstrom black hole with non-Abelian charge

The action of the model, which corresponds to some Lagrangian pieces of the
GSU2P theory Gallego Cadavid  _et al._ (2020), includes quartic order self-
interactions of the vector field111As the inclusion of a mass term
$\mu^{2}B_{a\alpha}B^{a\alpha}$ spoils the existence of the solution, it has
been taken away from the model. This result is similar to the classical
massive vector field, where the mass needs to vanish to guarantee regularity
of the solution and to allow, therefore, a vector hair to exist Bekenstein
(1972).,

$S=\frac{1}{16\pi}\int\sqrt{-g}\,d^{4}x[R-F_{a\mu\nu}F^{a\mu\nu}\\\
+\chi_{1}B_{a\mu}B^{a\mu}B_{b\nu}B^{b\nu}+\chi_{2}B_{a\mu}B^{a}{}_{\nu}B_{b}{}^{\mu}B^{b\nu}],$
(1)

where $R$ is the Ricci scalar, $B_{a\mu}$ represents the vector fields,
$F_{a\mu\nu}=\partial_{\mu}B_{a\nu}-\partial_{\nu}B_{a\mu}+\tilde{g}\epsilon_{abc}B^{b}{}_{\mu}B^{c}{}_{\nu}$
is the field strength, $\tilde{g}$ is the gauge coupling constant and
$\epsilon_{abc}$ is the structure constant tensor of the SU(2) group. In
geometrized units $\tilde{g}$ has units of inverse length, and the free
parameters $\chi_{1}$ and $\chi_{2}$ have units of inverse square length.

The line element in a stationary and spherical symmetric spacetime has the
following form,

$\displaystyle ds^{2}$
$\displaystyle=g_{tt}(r)dt^{2}+g_{rr}(r)dr^{2}+r^{2}d\Omega^{2}$
$\displaystyle=-e^{-2\delta}Ndt^{2}+N^{-1}dr^{2}+r^{2}d\Omega^{2},$ (2)

where $N=1-2m/r$, $\delta$ and $m$ are functions of the coordinate $r$, and
$d\Omega$ is the solid angle element. Regarding the vector fields we chose the
Wu-Yang monopole, given by,

$\mathbf{B}=(w/v+1)\,\mathbf{t}_{\phi}d\theta+(v-w)\sin\theta\,\mathbf{t}_{\theta}d\phi$
(3)

where,

$\displaystyle\mathbf{t}_{\theta}$
$\displaystyle=\cos\theta\cos\phi\,\mathbf{t}_{1}+\cos\theta\sin\phi\,\mathbf{t}_{2}-\sin\theta\,\mathbf{t}_{3},$
(4) $\displaystyle\mathbf{t}_{\phi}$
$\displaystyle=-\sin\phi\,\mathbf{t}_{1}+\cos\phi\,\mathbf{t}_{2},$ (5)

in which $\mathbf{t}_{i}=-i\sigma_{i}/2$ correspond to the vector basis of the
SU(2) algebra with $\sigma_{i}$ being the Pauli matrices, $w$ is constant and
$v$ is an integer denoting the azimuthal winding number. We use the coupling
constant $\tilde{g}$ to define the normalized variables, $\hat{r}=r\tilde{g}$,
$\hat{m}=m\tilde{g}$, $\hat{\chi}_{1}=\chi_{1}/\tilde{g}^{2}$ and
$\hat{\chi}_{2}=\chi_{2}/\tilde{g}^{2}$. The form of the equations in the
normalized variables can be obtained effectively by setting $\tilde{g}=1$.
Hereafter, all the equations are normalized, but we drop the hat to ease the
notation.

The field equation obtained after varying the action with respect to
$B_{a\mu}$ is,

$(v+w)\left[\left(v^{2}\chi_{1}+v^{2}\chi_{2}+\chi_{1}\right)(v+w)^{2}+w(v-w)\right]=0.$
(6)

The solutions of this last equation are,

$w_{\rm schw}=-v,$ (7) $w_{\rm
I,II}=\frac{v+2v\chi_{1}+2v^{3}\left(\chi_{1}+\chi_{2}\right)\pm\sqrt{v^{2}\left[8v^{2}\chi_{2}+8\left(v^{2}+1\right)\chi_{1}+1\right]}}{2\left(1-v^{2}\chi_{1}-v^{2}\chi_{2}-\chi_{1}\right)}.$
(8)

The first solution (7) is the trivial solution with vanishing vector field
which corresponds to the Schwarzschild spacetime222To fully understand this
solution, we need to consider Eq. (9). If we set $w_{\rm sch}=-v$ we find
$m^{\prime}=0$, implying that the mass function takes the form of the
Schwarzschild solution, where $m=M$ and $\delta=0$. Here $M$ is the total
gravitational mass. This, in turn, is consistent with an asymptotically flat
solution.. Instead, the latter solutions (8), with two branches I and II,
allow the existence of a non-trivial vector field we shall focus on, and
constitutes, therefore, an important outcome of this work, as will be
described in detail below.

The field equations obtained by varying the action (1) with respect to metric
are given by,

$m^{\prime}-\frac{\left(v^{3}-vw^{2}\right)^{2}}{2r^{2}v^{4}}\\\
+\frac{(v+w)^{4}\left[v^{4}\chi_{2}+\left(v^{2}+1\right)^{2}\chi_{1}+\chi_{2}\right]}{{4r^{2}v^{4}}}=0$
(9) $\delta^{\prime}=0$ (10)
$m^{\prime\prime}+\frac{v^{2}(v-w)(v+3w)(v+w)^{2}}{r^{3}v^{4}}\\\
+\frac{\left[\left(v^{4}+5\right)\chi_{2}+\left(v^{4}+6v^{2}+5\right)\chi_{1}\right](v+w)^{4}}{2r^{3}v^{4}}=0.$
(11)

We look for asymptotically flat solutions, i.e. the components of the metric
have the following behavior $g_{\mu\nu}\to\eta_{\mu\nu}+\mathcal{O}(r^{-1})$
when $r\to\infty$, where $\eta_{\mu\nu}$ are the components of the Minkowski
metric (see e.g. Wald (1984)). This implies that the functions tend
asymptotically to $\delta\to 0$, $m\to M\equiv{\rm finite}$ Wald (1984). The
Arnowitt-Desser-Misner mass, which coincides with the Komar mass in this
spherically symmetric and stationary case, is given by the asymptotic value of
$m(r)$, thus $M$ corresponds the total gravitational mass Ashtekar and Magnon-
Ashtekar (1979). Under these conditions Eq. (10) is easily solved as
$\delta=0$. On the other hand, the solution of mass function $m$ for the cases
given by (8) has the following Reissner-Nordström solution,

$m=M-\frac{Q_{\rm NA}^{2}}{2r},$ (12)

where is a constant representing the effective charge and depends on the free
parameters of the theory.

Eqs. (9) and (11) must be consistent, which implies an additional constraint
between $\chi_{1},\chi_{2}$ and $v$ given by,

$\frac{3\left(v^{4}-1\right)\left(\chi_{1}+\chi_{2}\right)\left\\{4v^{3}\chi_{2}+3\sqrt{v^{2}\left[8v^{2}\chi_{2}+8\left(v^{2}+1\right)\chi_{1}+1\right]}+4\left(v^{3}+v\right)\chi_{1}+5v\right\\}^{2}}{\left[v^{2}\chi_{2}+\left(v^{2}+1\right)\chi_{1}-1\right]{}^{4}}=0.$
(13)

One solution fixes the winding number as,

$v=\pm 1,$ (14)

with the parameters $\chi_{1}$ and $\chi_{2}$ being independent. The other
solution gives a relation between the free parameters of the action model,

$\chi_{2}=-\chi_{1},$ (15)

with the winding number now unconstrained333It seems that there is another
solution given by the vanishing of the expression between braces in the
numerator of Eq. (13). Nevertheless, this does not constitute a solution
because it makes also the denominator to vanish, inducing a divergence..

Despite the fact that the solution set by (14) gives two possible solutions,
it represents only one since changing the sign of the winding number
interchanges the solutions $w_{\rm I}$ and $w_{\rm II}$ as can be verified in
Eq. (8). In the solution given by (14) the value of $w$ depends on the
combination $\chi=2\chi_{1}+\chi_{2}$,

$w_{\rm I,II}=\frac{1+2\chi\pm\sqrt{1+8\chi}}{2-2\chi}\,.$ (16)

In the other solution corresponding to (15), the value of $w/v$ has the same
functional form of (16) after making $\chi\mapsto\chi_{1}$. Therefore, the
physical behavior of both cases can be analyzed by means of (16). Hence, all
subsequent analysis will be carried out in terms of the new effective coupling
constant $\chi$. Notice that once the value of $\chi$ is determined, the value
of $w$ is fully specified. This is consistent with the Wu-Yang monopole
Ansatz, which assumes that $w$ is a constant, as stated in Eq. (3).

Consequently, the value of the effective charge $Q^{2}_{\rm NA}$ is given by,

$Q^{2}_{\rm
NA,I,II}=\frac{1-4\chi(5+2\chi)\mp(1+8\chi)^{3/2}}{2(1-\chi)^{3}}\,,$ (17)

where the $-$ sign corresponds to the branch I, and the + sign to the branch
II. The dependence of both quantities on the coupling constant is displayed in
Figure 1.

There exist solutions for the interval $-1/8<\chi<1\cup\chi>1$. When $\chi=0$
the branch I solution corresponds to the Schwarzschild solution, and the
branch II solution corresponds to the EYM charged solution with $Q^{2}=1$.
When $\chi\to 1$ the branch I has a divergence, in contrast, in the same
limit, the branch II is finite and tends to $w_{\rm II}\to-1/3$, $Q^{2}_{\rm
NA,II}\to 16/27$.

On the hand, if we assume that the total mass of the black hole in normalized
units is $M=1$, we can find regions where $Q^{2}_{\rm NA}>1$, corresponding to
naked singularities. These regions are given by
$-1/8<\chi<(11-5\sqrt{5})/2\cup 1<\chi<(11+5\sqrt{5})/2$, for the branch I,
and $-1/8<\chi<0$, for the branch II (see the gray regions in Figure 1).
Notice that we conjecture the existence of a naked singularity based solely on
the vanishing of the event horizon. However, this criterion alone is
insufficient to ensure the formation of such a spacetime. Therefore, a formal
stability analysis must be carried out to determine whether a naked
singularity can indeed form in this theory444There are in fact physical
reasons to believe that the formation of a naked singularity, at least in a
Kerr BH spacetime, is highly unlikely in any realistic astrophysical collapse
scenario Visser (2007)..

Finally, it is worthwhile to mention that the energy density associated with
the vector fields is $\rho=Q^{2}_{\rm NA}/(8\pi r^{4})$. Thus, when the non-
Abelian charge is imaginary the energy density is negative, and this happens
only in the branch I in the interval $0<\chi<1$.

Figure 1: Values of the vector field and effective charge as a function of
$\chi$ for the branch I (left panel) and for the branch II (right panel). In
both cases when $\chi\to\infty$ the charge tends to zero, thus the solution
becomes the Schwarzschild spacetime. If the mass of the black hole in
normalized units is $M=1$, values of $Q^{2}_{\rm NA}>1$ represent a naked
singularity. These last cases are shown as gray regions. In the branch II the
charge is always real. In the branch I for $0<\chi<1$ the charge is imaginary,
which implies that the energy density is negative. Notice also that the shift
$\chi+1/8$ in the abscissa has been done for convenience.

Figure 2: Left panel: position of the event horizon $x_{H}$ as a function of
the coupling constant for both charges as specified by the legend. Right
panel: innermost stable circular orbit $x_{\rm ISCO}$ as a function of the
coupling parameter $\chi$ as well for both non-Abelian charges. In both plots
red points indicate the match with the extremal RN solution. The two-fold
degeneracy of $q_{\rm NA,I}$ and the convergence to the Schwarzschild
solutions when $\chi\to\infty$ for both charges is evident here. Note that the
parameter $\chi$ explicitly appears in the expression for $q_{NA}$ as given by
Eq. (17). Therefore, the event horizon and ISCO radii, which are determined by
Eqs. (18) and (22), respectively, are dependent on $\chi$ as well.

###  II.1 Event horizon

In the stationary and spherically symmetric case, the vanishing of the metric
function $g_{tt}$ defines unequivocally the horizon (see e.g. Vishveshwara
(1968)). When multiple solutions exist the greatest positive solution is
identified with the event horizon of the black hole $r_{+}=r_{H}$. In
particular, the metric function associated with the Reissner-Nordström
solution has two distinct roots

$r\pm=M\pm\sqrt{M^{2}-Q_{\rm NA}^{2}}.$ (18)

The internal solution $r_{-}$ is an apparent horizon and the external solution
corresponds to the event horizon. Hence, any observer outside the black hole
(in asymptotically flat spacetimes), or on the event horizon itself, cannot
see any singularity because they are protected by an event horizon. Otherwise
it is said to posses a naked singularity at $r=0$. We do not mention here all
the minor details and conditions about the precise formulation of what is
called the weak cosmic censorship conjecture555We refer reader to Wald (1984)
for a general and robust formulation of the cosmic censor conjecture..

The structure of the RN solution has been studied extensively whereby we do
not pretend to make here a detailed examination on this subject. Nevertheless,
an intriguing query arises when one asks about the implications of the
coupling constant $\chi$ on the charge and, therefore, on the event horizon.
In particular, we are interested in finding which values of the coupling
constant account for the convergence to both the Schwarzschild and the
extremal RN black holes, as limit cases of our solution.

When performing numerical analysis and presenting the corresponding general
discussion, we shall work with dimensionless variables by normalizing all
physical quantities by the black hole mass $M$, unless otherwise said.
Accordingly, we introduce, as usual, the charge to mass ratio $q_{\rm
NA}=Q_{\rm NA}/M$ and the dimensionless radial coordinate $x=r/M$.

The event horizon is depicted in left panel of Figure 2 as a function of the
coupling constant for both charges666When plotting, the shift $\chi+1/8$ in
the abscissa is done for convenience.. For $q_{\rm NA,II}$ (purple solid line)
the event horizon is a well-behaved function of the coupling constant. Indeed,
it covers continually the full range $\chi\in(0,\infty)$ where the finite
extreme value corresponds to the extremal RN case, in which case both
solutions of Eq. (18) meet at $x_{H}=1$ (red point on the purple curve); while
large values lead to the uncharged solution where $x_{H}=2$ and the apparent
horizon (dashed curve) coincides with the singularity. This is hence a quite
normal behavior that reproduces plainly the standard RN solution.

On the contrary, the event horizon for the case $q_{\rm NA,I}$ (blue solid
curves) exhibits a peculiar structure: there is a two-fold degeneracy with
respect to the constant coupling. It means that the limit cases, that is, the
Schwarzschild solution and the extremal RN black hole, can be described in two
distinct regions of the parameter space. The same degeneracy is also presented
for the apparent horizon (blue dashed curves). This is clearly appreciated in
the regions defined by the ranges $\chi\in(-0.0901,0)$ and
$\chi\in(11.0902,\infty)$ of the left panel of Figure 2. There is also a
special region of the parameter space $\chi\in(0,1)$, where the square of the
effective charge becomes negative and gravity repulsive. Out of the mentioned
ranges results in $q_{\rm NA}^{2}>1$, a case that describes a naked
singularity.

###  II.2 Inner most stable circular orbit

Next, we turn our attention to the motion of massive test particles in the
equatorial plane, as usual, to find the position of the so-called innermost
stable circular orbits (ISCO). Following the standard procedure (see e.g.
Pugliese  _et al._ (2011)), one finds that the geodesic equation can be
written as

$\left(\frac{dr}{ds}\right)^{2}=\left[E^{2}-N\left(1+\frac{L^{2}}{r^{2}}\right)\right],$
(19)

where $E$ and $L$ are identified, respectively, as the energy and angular
momentum which correspond to conserved quantities as in the classical
Keplerian motion. From the above equation, one can also identify the effective
potential

$V_{\rm eff}=N\left(1+\frac{L^{2}}{r^{2}}\right).$ (20)

We remind the reader that $N=1-2m/r$, the mass function $m(r)$ has the form
given by Eq. (12), and the (squared) angular momentum is defined as

$L^{2}=\frac{N^{\prime}(r)r^{3}}{2N(r)-rN^{\prime}(r)}.$ (21)

In marginally stable circular orbits, the local extremum of the effective
potential, i.e. $V_{\rm eff}^{\prime\prime}=0$, determines the position of the
ISCO. It yields

$r_{\rm ISCO}=2M+\frac{4M^{4}-3M^{2}Q_{\rm NA}^{2}+\left(8M^{6}-9M^{4}Q_{\rm
NA}^{2}+2M^{2}Q_{\rm NA}^{4}+\sqrt{5M^{8}Q_{\rm NA}^{4}-9M^{6}Q_{\rm
NA}^{6}+4M^{4}Q_{\rm NA}^{8}}\right)^{2/3}}{M\left(8M^{6}-9M^{4}Q_{\rm
NA}^{2}+2M^{2}Q_{\rm NA}^{4}+\sqrt{5M^{8}Q_{\rm NA}^{4}-9M^{6}Q_{\rm
NA}^{6}+4M^{4}Q_{\rm NA}^{8}}\right)^{1/3}}.$ (22)

This expression of course equals the corresponding standard RN case Pugliese
_et al._ (2011). At this point we are able to compute the location of the ISCO
in terms of the coupling constant $\chi$ with the aid of Eq. (17). The result
is shown in right panel of Figure 2. Clearly the effect of the coupling
constant is replicated on the ISCO structure in the same fashion as it is on
the event horizon. Hence the same parameter space region establishes the
corresponding Schwarzschild solution $x_{\rm ISCO}=6$, and the extremal RN
case $x_{\rm ISCO}=4$ (red point on all curves) for both positive square of
the non-Abelian charges as can be read from the plot. Thus, all the discussed
properties for the event horizon are preserved for the ISCO location. It is
interesting to notice, however, that for the charge $q_{\rm NA,I}$, there
appears another region where gravity is repulsive in the interval
$\chi\in(1,11.092)$, in comparison with the ones seen for the event horizon.
For this charge it is admissible then to have stable circular orbits for
almost all the parameter space with the exception of its extreme values even
though there does not exist an event horizon. This particular behavior is
reminiscent of having stable accretion flows onto a RN black hole, even
thought there is not a physical event horizon. Hence, in a naked singularity
situation it is (mathematically) possible to have stable circular orbits. We
will not discuss this point in detail, however, since it is not of physical
interest for the present work.

###  II.3 Photon sphere and shadow

A black hole has a central dark area called the shadow. This shadow is not
delimited by the event horizon, but by the photon sphere, which is made of
circular photon orbits. The radius of the photon sphere $r_{\rm ph}$ is given
by solving,

$r_{\rm ph}g_{tt}^{\prime}(r_{\rm ph})-2g_{tt}(r_{\rm ph})=0.$ (23)

However, the observed shadow radius, $r_{\rm sh}$ is given by the lensed image
of this surface Psaltis (2008),

$r_{\rm sh}=\frac{r_{\rm ph}}{\sqrt{g_{tt}(r_{\rm ph})}}.$ (24)

For the Reissner-Nordström solution the shadow radius is,

$\frac{r_{\rm sh}}{M}=x_{\rm sh}=\frac{\sqrt{2}\left(\sqrt{9-8q_{\rm
NA}^{2}}+3\right)}{\sqrt{\frac{4q_{\rm NA}^{2}+\sqrt{9-8q_{\rm
NA}^{2}}-3}{q_{\rm NA}^{2}}}}\,.$ (25)

This last equation can be inverted such that, for a given observed shadow
radius, we can constrain the charge $q_{\rm NA}$. From this value we can
constrain in turn the free parameters of the theory $\tilde{g}$ and $\chi$.
The observations of Sagittarius A∗ made by the EHT collaboration give the
following constraints on the size of the shadow, which depend on the mass-to-
distance ratio Event Horizon Telescope Collaboration  _et al._ (2022),

$4.5\lesssim x_{\rm sh}\lesssim 5.5,$ (26)

with Keck and

$4.3\lesssim x_{\rm sh}\lesssim 5.3,$ (27)

with VLTI. With this information, we can place constraints on the parameters
$\tilde{g}$ and $\chi$, shown in Figure 3. Since (17) is normalized but
depends physically on the $\tilde{g}$ and $\chi$, only in this part of the
analysis do we put back the units. The first conclusion is that an extremal RN
black hole is not consistent with the current observations.

If the non-Abelian charge is real, the maximum value of the shadow is $r_{\rm
sh}=3\sqrt{3}M$, obtained when $Q_{\rm NA}=0$, i.e., when the solution is the
Schwarzschild spacetime. This last case corresponds to the limit
$\chi\to\infty$ in both branches I and II. The shadow of a Schwarzschild black
hole is inside the intervals (26) and (27). Since the charge is always real
for _branch II_ , we can give only lower limits on $\chi M^{2}$ for this
branch (see the blue and orange curves in the lower plot of Figure 3). For
example, if $\tilde{g}M=1.01$, the minimum values of $\chi M^{2}$ are
$0.608985$ and $0.226688$, for Keck and VLTI, respectively.

In contrast, the _branch I_ exhibits imaginary values of the charge, allowing
a greater shadow than the one of a Schwarzschild black hole. Also, in this
branch the allowed region of values for $\tilde{g}M$ and $\chi M^{2}$ are
greater than in branch II. For instance, when $\tilde{g}M=0.9$ the minimum
values of $\chi M^{2}$ are $-0.0326136$ and $-0.0414562$, for mass-to-distance
ratio with Keck and VLTI, respectively. This implies a less stringent
constraint in the parameters than in the branch discussed above. More exactly,
if $\tilde{g}M=1.01$, the constraints are $(-0.0545704\lesssim\chi
M^{2}\lesssim 0.0223263)\,\cup\,14.3271\lesssim\chi M^{2}$, for Keck and
$(-0.0708403\lesssim\chi M^{2}\lesssim 0.00774184)\,\cup\,12.6237\lesssim\chi
M^{2}$, for VLTI, see the upper plots in 3. More cases for branch I can be
inferred from the Figure 3.

Figure 3: Constraints on the parameters $\tilde{g}$, and $\chi$ obtained from
the observation of the shadow of the object located at the Galactic Center of
the Milky Way, Sagittarius A⋆. The upper plots correspond to the branch I
which has been split into two ranges of $\chi$ for better representation, and
the lower plot corresponds to the branch II. The different colored curves show
the possible values of $g_{c}$ and $\chi$ for the lower and upper limits on
$x_{\rm sh}$ given by (26) and (27) as indicated by the legend. For
comparison, we have added the curves for $x_{\rm sh}=4,5.19,3\sqrt{3}$. The
regions with slanted orange lines represent the constraints obtained from the
mass-to-distance ratio with Keck, and the blue regions represent the
constraints with VLTI. It can be seen that a naked singularity is not
consistent with the observational data.

The main purpose now is to figure out how the structure of the non-Abelian RN
black hole can impact the transonic properties of accretion flows. After
formulating the basic hydrodynamics equations, this subject will be
investigated by employing both numerical and analytical treatments in the
subsequent sections.

##  III Spherical steady accretion flows in a spherically symmetric spacetime

Bondi accretion process in a spherically symmetric spacetime is briefly
described in this section following the general prescription presented in Ref.
Bauer  _et al._ (2022). Accordingly, we consider a spacetime with line element
given by (2), with $m$ and $Q^{2}_{\rm NA}$ given, respectively, by (12) and
(17), describing a black hole of mass $M$ and non-Abelian charge $Q_{\rm
NA}^{2}$ in Schwarzschild coordinates. Although the non-Abelian RN solution
Eq. (12) is, in principle, globally indistinguishable for the standard RN BH
solution, it was found, formally speaking, in the framework of modified
gravity. So, we keep as much as possible the generality in the
description777This general treatment serves also as a starting point to other
(non-analytical) BH solutions found in the GSU2P theory Martinez  _et al._
(2023).. On the other hand, we consider a steady fluid with total density
$\rho$, mass density $\rho_{0}$ and internal energy density $\epsilon$, such
that $\rho=\rho_{0}+\epsilon$. For isentropic fluids the pressure can be
defined as $P=k\rho^{\gamma}$ where $k$ is a constant and $\gamma$ is the
adiabatic index. For perfect fluids, the stress energy momentum tensor is
given by

$T^{\mu\nu}=(\rho+P)u^{\mu}u^{\nu}+Pg^{\mu\nu},$ (28)

where $u^{\mu}=(u^{t},u^{r},0,0)$ is the four velocity of the fluid
characterized by infall radial flow. The normalization condition allows ones
to obtain the relation between the components
$u^{t}=\sqrt{\frac{g_{rr}(u^{2}+g_{tt})}{g_{tt}}}$, where we have defined for
abbreviation $u\equiv u^{r}$. From the baryon conservation and energy momentum
conservation

$\displaystyle\nabla_{\mu}(\rho_{0}u^{\mu})$ $\displaystyle=0,$ (29)
$\displaystyle\nabla_{\mu}T^{\mu\nu}$ $\displaystyle=0,$ (30)

one obtains two master equations, respectively

$\displaystyle\frac{\rho^{\prime}_{0}}{\rho_{0}}+\frac{u^{\prime}}{u}+$
$\displaystyle\Sigma=0,$ (31) $\displaystyle
uu^{\prime}+\frac{g_{tt}\prime}{2g_{rr}g_{tt}}(1+g_{rr}u^{2})+\frac{g_{rr}^{\prime}}{2g_{rr}}u^{2}$
$\displaystyle+\frac{c_{s}^{2}}{g_{rr}}(1+g_{rr}u^{2})\frac{\rho^{\prime}_{0}}{\rho_{0}}=0,$
(32)

where prime denotes radial derivative and the quantity
$\Sigma\equiv\frac{(\sqrt{-g})^{\prime}}{\sqrt{-g}}$ has been introduced for
brevity as in Bauer  _et al._ (2022). These equations reduce to the
Schwarzschild case when $-g_{tt}=\frac{1}{g_{rr}}=1-\frac{2M}{r}$. In finding
Eq. (32), we have used the definition of the sound speed of a medium at
constant entropy $c_{s}^{2}\equiv\frac{dP}{d\rho}$, and the useful relation
$P^{\prime}=\frac{(\rho+P)}{\rho_{0}}c_{s}^{2}\;\rho_{0}^{\prime}$ derived
from the first law of thermodynamics in the form
$\frac{d\rho}{d\rho_{0}}=\frac{\rho+P}{\rho_{0}}$. Integration of Eqs. (31)
and (32) gives, respectively, the mass accretion rate

$\dot{M}=4\pi r^{2}u\rho_{0},$ (33)

and, after some algebraic manipulations, the relativistic version of the
Bernoulli equation (see Ref. Bauer  _et al._ (2022). for more details)

$g_{tt}(1+g_{rr}u^{2})\left(\frac{\rho+P}{\rho_{0}}\right)=C,$ (34)

where $C$ is an integration constant. Once $C$ is defined by the boundary
conditions at infinity for instance, and the metric functions are specified,
the inward radial velocity can be computed for a given equation of state
$P=P(\rho_{0})$. Our primary concern is to solve this equation to determine,
in turn, the accretion rate given by Eq. (33). Before computing this, we find
critical values at which accretion flow is regular and causality is
guaranteed.

##  IV Accretion in a non-Abelian Reissner-Nordström black hole

###  IV.1 Critical accretion

We study in this part general conditions under which transonic flow can take
place in the vicinity of black holes. Next, we shall describe both the
spacetime geometry and the fluid nature. Let us first write Eqs. (31) and (32)
in the more convenient form

$\displaystyle\frac{u^{\prime}}{u}$
$\displaystyle=\frac{g_{tt}(2c_{s}^{2}(1+u^{2}g_{rr})\Sigma-u^{2}g_{rr}^{\prime})-(1+u^{2}g_{rr})g_{tt}^{\prime}}{2g_{rr}g_{tt}(u^{2}-c_{s}^{2}(g^{rr}u^{2}))},$
(35) $\displaystyle\frac{\rho_{0}^{\prime}}{\rho_{0}}$
$\displaystyle=-\frac{u^{2}g_{tt}(-2g_{rr}\Sigma+g_{rr}^{\prime})+(1+u^{2}g_{rr})g_{tt}^{\prime}}{2g_{rr}g_{tt}(u^{2}-c_{s}^{2}(g^{rr}u^{2}))}.$
(36)

Imposing regular condition in both equations, implies that both numerators
must vanish simultaneously at some critical point $r_{c}$, resulting in

$\displaystyle u_{c}^{2}$
$\displaystyle=-\frac{g_{tt}^{\prime}}{g_{tt}g_{rr}^{\prime}+g_{rr}(-2g_{tt}\Sigma+g_{tt}^{\prime})},$
(37) $\displaystyle c_{s,c}^{2}$
$\displaystyle=\frac{g_{rr}g_{tt}^{\prime}}{2g_{rr}g_{tt}\Sigma-
g_{tt}g_{rr}^{\prime}}.$ (38)

Causality constraint $c_{s}^{2}<1$ in the flow sets a special point ($r_{c}$)
by which the flow must pass. This physical requirement leads to the relation

$u_{c}^{2}=-\frac{c_{s,c}^{2}\;g_{rr}}{-1+c_{s,c}^{2}},$ (39)

between the radial velocity and the sound speed. At large radius, the flow is
in the subsonic regime $u^{2}<c_{s}^{2}$. So far these results are general in
the sense that can be applied to any spherically symmetric black hole
solution. At this point we must specify necessarily the metric functions to
obtain the exact forms for the critical velocity and sound speed Eqs. (37) and
(38), respectively. Thus, considering the non-Abelian Reissner Nordström BH
solution Eq. (12), leads to the critical values

$u_{c}^{2}=-\frac{Q_{\rm NA}^{2}-Mr}{2r^{2}},\;c_{s,c}^{2}=\frac{-Q_{\rm
NA}^{2}+Mr}{Q_{NA}^{2}+r(-3M+2r)}.$ (40)

From the transonic condition that $u_{c}^{2}=c_{s,c}^{2}$ at the critical
radius, and considering Eq. (39), one obtains unequivocally the critical
radius

$r_{c}=\frac{M+3c_{s,c}^{2}M\pm\sqrt{(M+3c_{s,c}^{2}M)^{2}-8c_{s,c}^{2}(1+c_{s,c}^{2})Q_{\rm
NA}^{2}}}{4c_{s,c}^{2}},$ (41)

where the non-trivial contribution of the non-Abelian charge to the
Schwarzschild solution, $r_{c,\rm Sch}=M(1+3c_{s,c}^{2})/2c_{s,c}^{2}$ and
$u_{c,\rm Sch}^{2}=M/2r_{c}$, is clearly manifested. From Eq. (41) we can see
that there exist, mathematically speaking, two distinct critical points but
the positive branch corresponds only to the physical solution which resides
outside the event horizon, and therefore, it is the one to which we shall pay
our attention. Of course, the negative branch is in a region observationally
inaccessible since it is delimited by the event horizon. Existence of the
critical radius demands

$\frac{Q_{\rm
NA}^{2}}{M^{2}}<\frac{1+6c_{s,c}^{2}+9c_{s,c}^{4}}{8c_{s,c}^{2}+8c_{s,c}^{4}},$
(42)

which in turn puts constraints on the coupling parameter for a given non-
Abelian charge. Nevertheless, the resulting expressions are very lengthy (and
not illuminating) to be reported here. We shall illustrate below this aspect
numerically.

In order to understand the striking structure of the critical radius we must
necessarily specify the non-Abelian charge and the nature of the fluid that
characterizes the sound speed. For simplicity in the former analysis, we
assume an isothermal fluid888Another possibility is to take a polytropic fluid
but it introduces an extra parameter that overclouds the present analysis
about the structure of the critical radius. Nevertheless, it will be
successfully addressed in the next part for the sake of completeness when
looking for the critical mass accretion rate. so that the sound speed equals
its equation of state. This simple choice will give us, however, a profound
insight about the behavior of the critical radius in terms of the model
parameters. So the present case suffices to prove the rich structure of the
critical radius due to the non-Abelian charge. It should be noticed, on the
other hand, that for the (physical) critical radius there are two solutions
due to the existence of the non-Abelian charges $Q_{\rm NA,I,II}$.

Let us first explore the effect of changing the sound speed on the critical
radius for some specific values of the coupling constant that cover mostly the
(physical) BH solutions of interest. The behavior of the critical radius
depends, however, on the non-Abelian charge chosen. This is shown in Fig. 4,
where the left and right panels correspond, respectively, to the non-Abelian
charges $q_{\rm NA,I}$ and $q_{\rm NA,II}$. In particular, for $q_{\rm NA,I}$
and $\chi=0$, the critical radius matches exactly the Schwarzschild solution,
whereas for $q_{\rm NA,II}$ the convergence is possible provided that
$\chi\to\infty$. This is the reason why we have displayed the exact
Schwarzschild solution (green curve) in right panel as an asymptotic limit of
large coupling constant values. Accordingly, for the larger case shown,
$\chi=20$, the solution is barely distinguishable from the uncharged case but
physically different since the former is endowed with a charge $q_{\rm
NA,I}=0.113438$. Hence, despite taking the same values of $\chi$ for both non-
Abelian charges, the physical implications do not hold for the critical radius
as it also happens for the event horizon and ISCO.

To better illustrate the role of $\chi$ in the critical radius, we take the
same values of the coupling constant, as described below, for both charges.
For instance, taking $q_{\rm NA,I}$ (see left panel) and a given coupling
parameter lead to the following physical situations: the extremal case
$\chi=-0.0901$, the Schwarzschild solution $\chi=0$, the charged solution
$\chi=20$ ($q_{\rm NA,I}=0.4112$) and a naked singularity $\chi=11.0222$, that
affect in a way different the behavior of the critical radius. Now, taking the
same values for $\chi$ as before but for $q_{\rm NA,II}$, the previous
physical meaning is lost. Naked singularities can take place, for instance,
for $\chi=-0.0901$ which in the $q_{\rm NA,I}$ case occurs for $\chi=11.0222$
. Notice that for the former case accretion is possible only for (subsonic)
sound speeds $c_{s,c}^{2}<0.3478$ while for the latter it occurs out of the
range $0.6930<c_{s,c}^{2}<1.5677$. This explains the discontinuity of those
curves. Hence accretion may be possible even though the event horizon
vanishes. An interesting discussion about how to distinguish a BH from a naked
singularity spacetime by using the image of thin accretion disks is addressed
in Shaikh and Joshi (2019).

As a general trend, for subsonic sound speeds the critical radius can be
located far away from the event horizon, and for supersonic sound speeds it
can be accommodated, on the contrary, between the event horizon and the
apparent horizon. For $c_{s}^{2}=1$ all critical points coincide with their
respective position of the event horizon as can be easily verified. Another
interesting feature is that for the Schwarzschild case and subsonic speeds
$c_{s,c}^{2}<1$ the critical radius is outside the event horizon while for
$c_{s,c}^{2}>1$ the critical radius is located behind it. See purple and green
curves in the left and right panels, respectively. This general discussion is
in agreement with former studies about the accretion of perfect fluids in the
RN metric Babichev  _et al._ (2011).

In contrast, in Fig. 5, the coupling parameter is fixed while the sound speed
varies. Notice that for the case $q_{\rm NA,I}$ the critical radius exhibits a
discontinuity as it was also perceived for the event horizon structure. It
should be noticed that as $c_{s,c}^{2}$ decreases (see for instance
$c_{s,c}^{2}$=1/4) the critical radius can be in a region where naked
singularity takes place. As to the case $q_{\rm NA,II}$, the critical radius
is also a monotonically increasing function of $\chi$, but it is well-behaved
until the condition Eq. (42) is broken for small values of $\chi$ and given
squared sound speeds. This discussion provides, in complement to the one made
around Fig. 4, a full picture of the general conditions under which the
critical radius exists in terms of the coupling constant through Eq. (42).

Now, we are in a more grounded position to investigate the effect of changing
the coupling parameter on the accretion properties, concretely on the radial
velocity and on the mass density around this class of black hole. This will be
carried out first numerically for an isothermal fluid, and later with the aid
of some analytical treatments for polytropic fluids, to allow a more robust
and complete exploration of the involved parameters. These two cases will be
then treated separately in the next sections.

Figure 4: Left panel: critical radius $x_{c}\equiv r_{c}/M$ (41) for the non-
abelian charge $q_{\rm NA,I}$ for a given coupling parameters that correspond
to the extremal case $\chi=-0.0901$, Schwarzschild solution $\chi=0$, charged
solution $\chi=20$ ($q_{\rm NA,I}=0.4112$) and a naked singularity
$\chi=11.0222$ ($q_{\rm NA,I}=1.01$). Right panel: critical radius for the
non-abelian charge $q_{\rm NA,II}$ for the same values of parameters as left
panel. The exact Schwarzschild solution (green curve) has been included here
for comparison. Notice that, on the contrary, the solutions $\chi=-0.0901$,
$\chi=0$, $\chi=20$ and $\chi=11.0222$ correspond, respectively, to a naked
singularity ($q_{\rm NA,II}=1.11352$), a extremal case and charge BH cases
$q_{\rm NA,II}=0.113438$ and $q_{\rm NA,II}=0.17345$.

Figure 5: Left panel: critical radius for the non-abelian charge $q_{\rm
NA,I}$ for certain critical sound speeds, as described in the legend, as a
function of the coupling parameter. Right panel: critical radius for the non-
abelian charge $q_{\rm NA,II}$ for the same parameters as left panel.

###  IV.2 Isothermal test fluid

Before describing accretion flows for a more general fluid, let us first
consider a simplistic but useful isothermal test fluid. It will provide us
some physical insights on how the coupling constant influences the behavior of
the infall radial velocity and mass density. This inquiry is complementary to
the discussion on the critical point realized previously.

We focus for comparison reasons in a range of the coupling constant that
resembles the RN BH and the Schwarzschild solution and leave the (allowed)
range that provides $q_{\rm NA,I}^{2}<0$ out of this analysis. It implies that
the corresponding parameter space of $\chi$ for a given non-Abelian charge,
provides the same physics whereby we focus on $q_{\rm NA,I}$. To illustrate
this point we show only a limit case for $q_{\rm NA,II}$ to see the
convergence. The full range of $\chi$ will be, however, considered in the
calculation of the critical accretion rate for a polytropic fluid. The above
are also advantageous for numerical facilities since we have many variables
involved.

Accordingly, the equation of state is of the form $P=\kappa\rho$, with
$\kappa$ being a constant, from which the simple relation for the sound speed
$c_{s,c}^{2}=\kappa$ is derived.

Let us start our analysis by considering a stiff fluid $\kappa=1$. As we
already discussed, for this case ($c_{s,c}^{2}=1$) critical points coincide
with the event horizons no matter the value of the coupling constant. The
latter spans the allowed region of the parameter space $\chi\in(-0.09,0)$ for
the charge $q_{\rm NA,I}$, as can be seen in the bar legend of Figure. 6. We
are not considering the other possible range of values
$\chi\in(11.0902,\infty)$ because the same physical properties are replicated
in the already shown range. For the other cases describing an ultra-
relativistic fluid $\kappa=1/2$, a radiation fluid $\kappa=1/3$ and a sub-
relativistic fluid $\kappa=1/4$; all critical points move out of the BH as
$\kappa$ decreases according to Eq. (41). As a general trend, all transonic
solutions are delimited from below to the extremal RN solution $\chi=-0.0901$
and from above to the Schwarzschild case. This latter case is attained
whenever the coupling constant increases until it reaches the maximum value
showed here ($\chi=0$). We have also included the solution $\chi=11.09$
(dashed magenta curve) for the charge $q_{\rm NA,II}$ which matches the RN
solution of the $q_{\rm NA,I}$ case. What is of physical interest here is that
all infall radial velocities pass through the corresponding critical points
marked as points on the curves and computed from Eq. (41), guaranteeing thus
the transonic flow of all solutions, otherwise a stellar wind is generated.
This selects a unique solution with constant inward mass flux that connects
the subsonic with the supersonic regimes, as expected in Bondi-type accretion
Bondi (1952); Michel (1972). Far beyond the BH influence, that is, in the non-
relativistic regime, the radial velocities of the particles are too low such
that the inflow rate is decreased compared with particles with high-speed
velocities.

For the mass density distribution due to the BH gravitational potential, the
RN solution now delimits all solutions from above, as seen in Figure. 7. As in
the radial velocity case, $\chi$ covers the same range of values for each
value of the constant $\kappa$ considered. We can observe that as $\kappa$
reduces, the mass density is more spread out along the radial coordinate.

So far we have obtained the expected behavior for the radial velocity and mass
density within the steady-state, spherical accretion scenario for the
discussed range of $\chi$-values. As in previous sections, this has been very
useful to understand the role of the coupling constant on the transonic flows
and how our solutions approach to the extremal RN and Schwarzschild solutions.
Making a more robust description of the fluid, in next section we will perform
analytical calculations of the critical accretion rate in the fully
relativistic regime along with numerical computations for the entire range of
values of $\chi$ and both non-Abelian charges.

Figure 6: Infall radial velocity for some specific values of $\kappa$
describing a stiff fluid ($\kappa=1$), ultra-relativistic fluid
($\kappa=1/2$), radiation fluid ($\kappa=1/3$) and sub-relativistic fluid
($\kappa=1/4$) while the coupling parameter $\chi$ spans the allowed range as
depicted by the bar legend. The dotted magenta curve that delimits all the
possible transonic solutions from below corresponds to the extremal case
$\chi=11.09$ for $q_{\rm NA,I}$ (or $\chi=-0.09$ for $q_{\rm NA,II}$), while
$\chi\rightarrow 0$ resembles the Schwarzschild solution from above. It is
worth noting that the event horizon is always located in the interval
$x\in[1,2]$, with the lower and upper limits corresponding to the extremal
Reissner-Nordström and Schwarzschild solutions, respectively. In the specific
case of stiff fluid $\kappa=1$ (top left panel), the position of the critical
points coincide with the location of the event horizons. However, for smaller
and larger values of $\kappa$ and $\chi$, all critical points shift toward
larger values of the corresponding event horizons. Note also that the
divergences observed are not physical, but they are related, instead, to the
coordinate singularity problem associated with the Schwarzschild coordinates.
This singularity can be removed by moving to non-singular coordinates at the
event horizon as the Eddington-Finkelstein coordinates.

Figure 7: Mass density distribution for some specific values of $\kappa$
describing a stiff fluid ($\kappa=1$), ultra-relativistic fluid
($\kappa=1/2$), radiation fluid ($\kappa=1/3$) and sub-relativistic fluid
($\kappa=1/4$) while the coupling parameter $\chi$ spans the allowed range as
depicted by the bar legend. The dotted magenta curve that delimits all the
possible transonic solutions from below corresponds to the extremal case
$\chi=11.09$ for $q_{\rm NA,I}$ (or $\chi=-0.09$ for $q_{\rm NA,II}$), while
$\chi\rightarrow 0$ (or $\chi\rightarrow\infty$) resembles the Schwarzschild
solution from above.

###  IV.3 Polytropic fluid

Once the properties of the steady flows are known at the sonic point, this is,
radial velocity, sound speed and critical radius, we can proceed to express
such quantities in terms of the boundary conditions, with the help of the
Bernoulli equation, as commonly done. The purpose of doing so is to calculate
the accretion rate explicitly. It is necessary also to adopt an equation of
state for the gas. So we study a non-relativistic baryonic gas with a
polytropic equation

$P=Kn^{\gamma},$ (43)

where $\gamma$ is the adiabatic index and $K$ is a constant. With this and
from the first of law of thermodynamics one can get Shapiro and Teukolsky
(1983)

$\rho=mn+\frac{K}{\gamma-1}n^{\gamma}.$ (44)

####  IV.3.1 Relativistic regime

Surprisingly, the fully relativistic accretion rate for the Reissner Nordström
solution has not been thoroughly treated in the literature, where most of
existing works have focused on the non-relativistic limit. It leaves
undoubtedly an incomplete comprehension of the full picture. We first derive
some analytical expressions and show some numerical examples to illustrate
better the role of the coupling constant on the mass accretion rate. We follow
closely Ref. Richards  _et al._ (2021), where Bondi accretion of steady
spherical gas flow onto a Schwarzschild black hole has been studied. We extend
this work to the charged case.

With the aid of the polytropic equation, it is possible to relate the sound
speed with the mass density

$c_{s}^{2}=\frac{\gamma k\rho_{0}^{\gamma-1}}{1+\gamma
k\rho_{0}^{\gamma-1}/(\gamma-1)}.$ (45)

This expression can be evaluated at the critical point and in the asymptotic
region to provide the useful relation

$\rho_{0,s}=\rho_{0,\infty}\left(\frac{c_{s,c}^{2}}{c_{s,\infty}^{2}}\right)^{\frac{1}{\gamma-1}}\left(\frac{\gamma-1-c_{s,\infty}^{2}}{\gamma-1-c_{s,c}^{2}}\right)^{\frac{1}{\gamma-1}}.$
(46)

This closed expression requires the knowledge of the critical sound speed that
can be extracted from the relativistic Bernoulli equation

$(1+3c_{s,c}^{2})\left(1-\frac{c_{s,c}^{2}}{\gamma-1}\right)^{2}=\left(1-\frac{c_{s,\infty}^{2}}{\gamma-1}\right)^{2}.$
(47)

So, once the sound speed at infinity is specified, the sound speed and the
mass density at the critical point are uniquely determined. The Bernoulli
equation is actually a cubic equation for $c_{s,c}^{2}$ with one real solution
for the range $1<\gamma<5/3$. There are some procedures to solve analytically
this equations as, for instance, a standard root-finding schema which we
implement to.

Having expressed all quantities at the critical radius in terms of the
boundary conditions, the critical accretion rate

$\dot{M}=4\pi\rho_{0,s}u_{s}r_{s}^{2},$ (48)

can be computed easily

$\dot{M}=4\pi\left(\frac{M}{c_{s,\infty}^{2}}\right)^{2}c_{s,\infty}^{2}\;\rho_{0,\infty}\;\lambda_{\rm
RN}^{\rm NA},$ (49)

with the accretion rate eigenvalue

$\lambda_{\rm RN}^{\rm
NA}\equiv\left(\frac{c_{s,c}^{2}}{c_{s,\infty}^{2}}\right)^{\frac{5-3\gamma}{\gamma-1}}\left(\frac{\gamma-1-c_{s,\infty}^{2}}{\gamma-1-c_{s,c}^{2}}\right)^{\frac{1}{\gamma-1}}\frac{(1+3c_{s,c}^{2})^{3/2}}{4}\beta,$
(50)

and the $\beta$ factor, containing information of the non-Abelian charge, is

$\beta=\frac{1}{4}\left[1+\sqrt{1-\frac{8c_{s,c}^{2}(1+c_{s,c}^{2})q_{\rm
NA}^{2}}{(1+3c_{s,c}^{2})^{2}}}\right]^{2}.$ (51)

This quantity clearly accounts for the deviation from the Schwarzschild case.
In what follows, we quantify such a deviation by computing the ratio of both
accretion rates

$\frac{\dot{M}_{\rm RN}^{\rm NA}}{\dot{M}_{\rm Sch}}=\frac{\lambda_{\rm
RN}^{\rm NA}}{\lambda_{\rm Sch}}=\beta.$ (52)

As it is known, the electric charge of the RN black hole reduces the accretion
rate compared to the Schwarzschild black hole. Our case may be, however,
different if the imaginary charge of the black hole is allowed, i.e. when
$q_{NA}^{2}=q_{NA,I}^{2}<0$ and $\chi\in(0,1)$. Under this choice,
$\dot{M}_{\rm RN}^{\rm NA}>\dot{M}_{\rm Sch}$ as can be verified in all left
panels of Figure 8 where the accretion rate has been plotted as a function of
the coupling constant $\chi$ for different adiabatic indices as indicated.
Notice that all curves meet in the corresponding position of the event horizon
$\chi\to 0,\infty$. Out of the mentioned range, the expected behavior of the
Reissner Nordström black hole is displayed just as the case $q_{\rm NA}=q_{\rm
NA,II}$ (right panels). Even though the non-Abelian charge case is effectively
distinguishable from its electric counterpart in the small range
$\chi\in(0,1)$ for $q_{\rm NA,I}$, these numerical computations allow ones, in
addition, to understand better the multiplicity of the non-Abelian Reissner
Nordström black solution and its implications in the accretion rate, in
particular for the range of values of $\chi$ that leads to a significant
enhancement of the accretion flow.

As a last remark about the boundary conditions, the larger the boundary sound
speed, the shorter the accretion rate variations are among the polytropic
fluid considered. The stiff case $\gamma=1$ is more sensible to the increase
of the boundary sound speed: notice how the dashed black curves in right
panels meet the other curves for low $\chi$ or, which is equivalent, in the
non-Abelian extremal case $q_{\rm NA}\to 1$.

Figure 8: Ratio of mass accretion rate in the non-Abelian RN BH to the
accretion rate in the Schwarzschild BH for different adiabatic indices as
indicated in the legend. Left and right panels correspond to the cases $q_{\rm
NA}=q_{\rm NA,I}$ and $q_{\rm NA}=q_{\rm NA,II}$, respectively. In tops panels
the boundary condition approaching to the non-relativistic regime
$c_{s,\infty}^{2}=0.001$ has been taken, whereas in middle and bottom panels a
relativistic boundary sound speed $c_{s,\infty}^{2}=0.1$ and
$c_{s,\infty}^{2}=0.5$ have been chosen, respectively, in contrasting. In all
cases the condition $u_{c}>c_{s,\infty}$ is guaranteed, ensuring thus the
transonic flow.

####  IV.3.2 Asymptotic limit

We do not describe here in detail all derivations concerning the asymptotic
limit of the mass accretion rate since it can be found, for instance, in
reference Shapiro and Teukolsky (1983). The purpose of this part is to check
the consistency with the well-known Newtonian limit and the corresponding
dependence on the non-Abelian charge. From the Bernoulli equation, one can
derive the useful relation

$c_{s,c}^{2}\approx\frac{2c_{s,\infty}^{2}}{(5-3\gamma)},$ (53)

under the non-relativistic condition $c_{s,c}\ll 1$ which holds for reasonable
large radius, $r\gg r_{c}$, far away from the BH gravitational influence. The
same condition leads to the simple relation

$c_{s,c}^{2}\approx K\gamma\;\rho_{0,c}^{\gamma-1},$ (54)

between the sound speed and the mass density at the critical point. This
implies that the mass density can be expressed in terms of the sound speed at
the infinity in view of Eq. (53) to yield

$\rho_{0,c}\approx\rho_{0,\infty}\left(\frac{c_{s,c}^{2}}{c_{s,\infty}^{2}}\right)^{\frac{1}{\gamma-1}}\approx\left(\frac{2}{5-3\gamma}\right)^{\frac{1}{\gamma-1}}.$
(55)

As in the uncharged case, physical solutions require $\gamma<5/3$ in this. At
this point, all of the above are standard and these quantities do not receive
contributions from the effective charge $Q_{\rm NA}$ at lowest order in
$c_{s,c}$. This is not the case however for the critical radius Eq. (41)
where, at leading order, $Q_{\rm NA}$ already appears explicitly as second
order power

$r_{c}\approx\frac{M}{2c_{s,c}^{2}}+\frac{3M^{2}-2Q_{\rm
NA}^{2}}{2M}+\frac{2(M^{2}Q_{\rm NA}^{2}-Q_{\rm
NA}^{4})}{M^{3}}c_{s,c}^{2}+\mathcal{O}(c_{s,c}^{4}).$ (56)

Keeping only higher order contributions in $c_{s,c}$ and using Eq. (53), the
critical radius can be approximated to

$r_{c}\approx\frac{(5-3\gamma)}{4c_{s,\infty}^{2}}M\eta(Q_{\rm NA}),$ (57)

where the dimensionless correction factor due to the effective charge has been
defined as

$\eta(Q_{\rm NA})=1+\frac{c_{s,\infty}^{2}}{(5-3\gamma)}\left(6-4\frac{Q_{\rm
NA}^{2}}{M^{2}}\right).$ (58)

With all this, the critical accretion mass rate can be written in the familiar
form

$\dot{M}\approx
4\pi\rho_{0,\infty}M^{2}c_{s,\infty}^{-3}\left(\frac{1}{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}\left(\frac{5-3\gamma}{4}\right)^{\frac{3\gamma-5}{2(\gamma-1)}}\eta(Q_{\rm
NA})^{2},$ (59)

which agrees with Freitas Pacheco (2011). Finally, the uncharged Newtonian
limit is straightforwardly achieved in the limit $Q_{\rm NA}\to 0$ and,
appropriately, $c_{c,\infty}^{2}\ll 1$.

##  V Discussion and Conclusions

Within the framework of modified theories of gravity, in particular in the
context of vector-tensor theories that follow the same spirit of Horndeski’s
theory, a kind of RN black hole solution with two different non-Abelian
effective charges has been found in terms of the coupling constants of the
involved Lagrangian pieces of the generalized SU(2) Proca theory. Such new
solutions correspond to genuine non-Abelian RN BH solutions in the sense they
were derived from a theory where the vector fields belong to the Lie algebra
of the SU(2) group. These objects then carry a dark charge because the
magnetic charge is not of electromagnetic origin. There do exist other
solutions coming from other more involved Lagrangians of the theory Gallego
Cadavid  _et al._ (2020), but they do not possess effective charges in the
asymptotic limit, which means that they converge, formally speaking, to the
Schwarzschild solution. This will be reported in a separate paper.

Even though these solutions retain the main global properties of the standard
RN BH derived from the Einstein-Maxwell theory, the solutions found exhibit an
appealing structure due to the non-trivial dependence on the coupling
constants. Interestingly, the solution with $q_{\rm NA,I}$ is characterized by
having a negative square of the effective charge for a certain range of values
of the coupling constant $\chi$, similar to the tidal charge of brane BHs
Dadhich  _et al._ (2000); Zakharov (2022); Neves (2020) and BH solutions in
Horndeski theory Babichev  _et al._ (2017). Preventing the naked singularity
from forming in both non-Abelian RN BHs leads to discard a small region of the
parameter space. This happens particularly for the solutions with $q_{\rm
NA,I}$ since the solution with $q_{\rm NA,II}$ is always real and well-behaved
in the sense that no divergences are present in the event horizon and ISCO
structures. Studying the RN solution within the framework of the GSU2P theory
allowed us to infer constraints on the free parameters of the model. This
finding is of particular physical interest as it can be compared with future
astrophysical constraints to check the consistency of the theory at the
relevant scales.

Some phenomenological implications of the BH solutions were also investigated
in the astrophysical setting in order to constrain the parameter space in a
joint way with the aforementioned theoretical considerations. They are
summarized as follows.

  * •

Observations of the EHT’s first images of Sagittarius A⋆ of our Galaxy, along
with Keck telescope results, set the first serious constraint on the free
parameters of the theory ($\tilde{g},\chi$), leaving almost all the available
parameter space of the non-Abelian RN BHs basically unconstrained. For a given
$\tilde{g}$, a lower limit on $\chi$ is determined as can be inferred from the
parameter space Fig. 3. As in the electric RN BH case, these observational
constraints also rule out regions of the parameter space of the BH solutions
associated with naked singularities and with the extremal BH. On the contrary,
the corresponding region of the parameter space for which $q_{\rm
NA,I}^{2}<0$, i.e., imaginary non-Abelian charge, is allowed.

  * •

As a first step toward a more realistic and elaborate description of accretion
processes, a fully relativistic treatment of spherical accretion of isothermal
and polytropic fluids onto this class of BHs has been performed to quantify
the effect of the non-Abelian charge and, therefore, of the coupling constant,
on the critical accretion rate. Interestingly, we have found some
dissimilarities in the accretion process with respect to the standard electric
RN case that can serve as a potential observational signature to test the
theory. Concretely, the critical accretion rate efficiency can be noticeably
improved compared to the Schwarzschild case (and also to the electric RN
case), that is, $\dot{M}_{\rm RN}^{\rm NA}>\dot{M}_{\rm Sch}$, provided that
$\chi\in(0,1)$ and $q_{\rm NA}=q_{\rm NA,I}^{2}<0$ which is, as discussed
above, admitted form the observational side. In this regard, we have examined
carefully, with the aid of numerical computations for different adiabatic
indices of a polytropic fluid, the role of the coupling parameter on the
transonic properties of steady flows.

  * •

As a way of probing the consistency of the non-Abelian BH solutions, the
Schwarzschild solution and the extremal RN BH solution are recovered in our
solutions, as limit cases of the theory, for certain values of $\chi$. This is
a probe of concept of how the BH solutions found behave in extreme regimes of
the parameter space.

An immediate theoretical extension of this work is to implement the Newman-
Janis algorithm to find rotating non-Abelian charge BH solutions. Effectively,
a Kerr-(non-Abelian) Newman BH solution is naturally expected. Although the
applicability of this algorithm must be taken with great care Hansen and Yunes
(2013), the absence of direct couplings of the gauge fields to curvature terms
guarantee the viability of this future work. Then, we plan to study the main
properties of the image of the resulting BHs, such as the shadows and photon
rings, surrounded by an optically and geometrically thin accretion disk and
the subsequent comparison with current observations. In this regard, it is
imperative to use observational constraints from the shadow of the
supermassive BH galaxy M87⋆, as was recently done for the tidal charge Neves
(2020). Gravitational and electromagnetic waveforms for charged black hole
binaries can be used to estimate the charges of BHs in current and future
gravitational wave experiments as has been discussed recently Christiansen
_et al._ (2021); Benavides-Gallego and Han (2022); Bozzola and Paschalidis
(2021); Liu  _et al._ (2020); Wang  _et al._ (2021). This is another
interesting way to assess the effect of the coupling constants in the strong-
field regime in the vicinity of BHs. Hence, gravitational wave observations
have also the potential of putting constraints on the coupling constants of
the theory.

## Acknowledgments

G. G. acknowledges financial support from Agencia Nacional de Investigación y
Desarrollo (ANID) through the FONDECYT postdoctoral Grant No. 3210417. J.F.R.
is thankful for financial support from the Patrimonio Autónomo - Fondo
Nacional de Financiamiento para la Ciencia, la Tecnología y la Innovación
Francisco José de Caldas (MINCIENCIAS - COLOMBIA) under the grant No.
110685269447 RC-80740–465–2020, project 69553.

## References

  * Poisson (2009) E. Poisson,  _A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics_ (Cambridge University Press, 2009).
  * Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky,  _Black holes, white dwarfs, and neutron stars: The physics of compact objects_ (1983).
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ (2022) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022).
  * Abuter  _et al._ (2020) R. Abuter _et al._ (GRAVITY), Astron. Astrophys. 636, L5 (2020), arXiv:2004.07187 [astro-ph.GA] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017a), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017b), arXiv:1710.05832 [gr-qc] .
  * Will (2014) C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021), arXiv:2010.14529 [gr-qc] .
  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Psaltis  _et al._ (2020) D. Psaltis _et al._ (Event Horizon Telescope), Phys. Rev. Lett. 125, 141104 (2020), arXiv:2010.01055 [gr-qc] .
  * Vagnozzi  _et al._ (2022) S. Vagnozzi _et al._ , (2022), arXiv:2205.07787 [gr-qc] .
  * Reissner (1916) H. Reissner, Annalen der Physik 355, 106 (1916).
  * Nordström (1918) G. Nordström, Koninklijke Nederlandse Akademie van Wetenschappen Proceedings Series B Physical Sciences 20, 1238 (1918).
  * De Rujula  _et al._ (1990) A. De Rujula, S. L. Glashow, and U. Sarid, Nucl. Phys. B 333, 173 (1990).
  * Cardoso  _et al._ (2016) V. Cardoso, C. F. B. Macedo, P. Pani, and V. Ferrari, JCAP 05, 054 (2016), [Erratum: JCAP 04, E01 (2020)], arXiv:1604.07845 [hep-ph] .
  * Zajaček and Tursunov (2019) M. Zajaček and A. Tursunov, (2019), arXiv:1904.04654 [astro-ph.GA] .
  * Kocherlakota  _et al._ (2021) P. Kocherlakota _et al._ (Event Horizon Telescope), Phys. Rev. D 103, 104047 (2021), arXiv:2105.09343 [gr-qc] .
  * Bozzola and Paschalidis (2021) G. Bozzola and V. Paschalidis, Phys. Rev. Lett. 126, 041103 (2021), arXiv:2006.15764 [gr-qc] .
  * Liu  _et al._ (2020) L. Liu, O. Christiansen, Z.-K. Guo, R.-G. Cai, and S. P. Kim, Phys. Rev. D 102, 103520 (2020), arXiv:2008.02326 [gr-qc] .
  * Christiansen  _et al._ (2021) O. Christiansen, J. Beltrán Jiménez, and D. F. Mota, Class. Quant. Grav. 38, 075017 (2021), arXiv:2003.11452 [gr-qc] .
  * Wang  _et al._ (2021) H.-T. Wang, P.-C. Li, J.-L. Jiang, G.-W. Yuan, Y.-M. Hu, and Y.-Z. Fan, Eur. Phys. J. C 81, 769 (2021).
  * Benavides-Gallego and Han (2022) C. A. Benavides-Gallego and W.-B. Han, (2022), arXiv:2209.00874 [gr-qc] .
  * Zakharov (2014) A. F. Zakharov, Phys. Rev. D 90, 062007 (2014), arXiv:1407.7457 [gr-qc] .
  * Zajaček  _et al._ (2018) M. Zajaček, A. Tursunov, A. Eckart, and S. Britzen, Mon. Not. Roy. Astron. Soc. 480, 4408 (2018), arXiv:1808.07327 [astro-ph.GA] .
  * Done  _et al._ (2007) C. Done, M. Gierlinski, and A. Kubota, Astron. Astrophys. Rev. 15, 1 (2007), arXiv:0708.0148 [astro-ph] .
  * Juraeva  _et al._ (2021) N. Juraeva, J. Rayimbaev, A. Abdujabbarov, B. Ahmedov, and S. Palvanov, Eur. Phys. J. C 81, 70 (2021).
  * Fragione and Loeb (2020) G. Fragione and A. Loeb, ApJ 901, L32 (2020).
  * Misner  _et al._ (2017) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (2017).
  * Bizon (1990) P. Bizon, Phys. Rev. Lett. 64, 2844 (1990).
  * Volkov and Galtsov (1989) M. S. Volkov and D. V. Galtsov, JETP Lett. 50, 346 (1989).
  * Kuenzle and Masood-ul Alam (1990) H. P. Kuenzle and A. K. M. Masood-ul Alam, J. Math. Phys. 31, 928 (1990).
  * Lee  _et al._ (1992) K.-M. Lee, V. P. Nair, and E. J. Weinberg, Phys. Rev. Lett. 68, 1100 (1992), arXiv:hep-th/9111045 .
  * Radu and Tchrakian (2012) E. Radu and D. H. Tchrakian, Phys. Rev. D 85, 084022 (2012), arXiv:1111.0418 [gr-qc] .
  * Mazharimousavi and Halilsoy (2009) S. H. Mazharimousavi and M. Halilsoy, Phys. Lett. B 681, 190 (2009), arXiv:0908.0308 [gr-qc] .
  * Volkov and Gal’tsov (1999) M. S. Volkov and D. V. Gal’tsov, Phys. Rept. 319, 1 (1999), arXiv:hep-th/9810070 .
  * Herdeiro  _et al._ (2017) C. Herdeiro, V. Paturyan, E. Radu, and D. H. Tchrakian, Phys. Lett. B 772, 63 (2017), arXiv:1705.07979 [gr-qc] .
  * Mazharimousavi and Halilsoy (2011) S. H. Mazharimousavi and M. Halilsoy, Phys. Rev. D 84, 064032 (2011), arXiv:1105.3659 [gr-qc] .
  * Harnad  _et al._ (1980) J. Harnad, J. Tafel, and S. Shnider, Journal of Mathematical Physics 21, 2236 (1980).
  * Bartnik and McKinnon (1988) R. Bartnik and J. McKinnon, Phys. Rev. Lett. 61, 141 (1988).
  * Breitenlohner  _et al._ (1995) P. Breitenlohner, P. Forgács, and D. Maison, Nuclear Physics B 442, 126 (1995).
  * Penrose (1965) R. Penrose, Phys. Rev. Lett. 14, 57 (1965).
  * Hawking and Penrose (1970) S. W. Hawking and R. Penrose, Proc. Roy. Soc. Lond. A 314, 529 (1970).
  * Deser  _et al._ (1975) S. Deser, P. Van Nieuwenhuizen, and D. Boulware, in  _7th International Conference on Gravitation and Relativity_ (1975) pp. 1–18.
  * Donoghue (1994) J. F. Donoghue, Phys. Rev. D 50, 3874 (1994), arXiv:gr-qc/9405057 .
  * Burgess (2004) C. P. Burgess, Living Rev. Rel. 7, 5 (2004), arXiv:gr-qc/0311082 .
  * Heisenberg (2019) L. Heisenberg, Phys. Rept. 796, 1 (2019), arXiv:1807.01725 [gr-qc] .
  * Horndeski (1974) G. W. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * Ostrogradsky (1850) M. Ostrogradsky, Mem. Acad. St. Petersbourg 6, 385 (1850).
  * Heisenberg (2014) L. Heisenberg, JCAP 05, 015 (2014), arXiv:1402.7026 [hep-th] .
  * Allys  _et al._ (2016a) E. Allys, P. Peter, and Y. Rodriguez, JCAP 02, 004 (2016a), arXiv:1511.03101 [hep-th] .
  * Allys  _et al._ (2016b) E. Allys, J. P. Beltran Almeida, P. Peter, and Y. Rodríguez, JCAP 09, 026 (2016b), arXiv:1605.08355 [hep-th] .
  * Gallego Cadavid and Rodriguez (2019) A. Gallego Cadavid and Y. Rodriguez, Phys. Lett. B 798, 134958 (2019), arXiv:1905.10664 [hep-th] .
  * Allys  _et al._ (2016c) E. Allys, P. Peter, and Y. Rodriguez, Phys. Rev. D 94, 084041 (2016c), arXiv:1609.05870 [hep-th] .
  * Gómez and Rodríguez (2019) L. G. Gómez and Y. Rodríguez, Phys. Rev. D 100, 084048 (2019), arXiv:1907.07961 [gr-qc] .
  * Gallego Cadavid  _et al._ (2020) A. Gallego Cadavid, Y. Rodriguez, and L. G. Gómez, Phys. Rev. D 102, 104066 (2020), arXiv:2009.03241 [hep-th] .
  * Gallego Cadavid  _et al._ (2022) A. Gallego Cadavid, C. M. Nieto, and Y. Rodriguez, Phys. Rev. D 105, 104051 (2022), arXiv:2204.04328 [hep-th] .
  * Garnica  _et al._ (2022) J. C. Garnica, L. G. Gomez, A. A. Navarro, and Y. Rodriguez, Annalen Phys. 534, 2100453 (2022), arXiv:2109.10154 [gr-qc] .
  * Martinez  _et al._ (2023) J. N. Martinez, J. F. Rodriguez, Y. Rodriguez, and G. Gomez, JCAP 04, 032 (2023), arXiv:2212.13832 [gr-qc] .
  * Bondi (1952) H. Bondi, Mon. Not. Roy. Astron. Soc. 112, 195 (1952).
  * Michel (1972) F. C. Michel, Ap&SS 15, 153 (1972).
  * Richards  _et al._ (2021) C. B. Richards, T. W. Baumgarte, and S. L. Shapiro, Mon. Not. Roy. Astron. Soc. 502, 3003 (2021), [Erratum: Mon.Not.Roy.Astron.Soc. 506, 3935 (2021)], arXiv:2101.08797 [astro-ph.HE] .
  * Aguayo-Ortiz  _et al._ (2021) A. Aguayo-Ortiz, E. Tejeda, O. Sarbach, and D. López-Cámara,  (2021), 10.1093/mnras/stab1127, arXiv:2102.12529 [astro-ph.HE] .
  * Bauer  _et al._ (2022) A. M. Bauer, A. Cárdenas-Avendaño, C. F. Gammie, and N. Yunes, Astrophys. J. 925, 119 (2022), arXiv:2111.02178 [gr-qc] .
  * Salahshoor and Nozari (2018) K. Salahshoor and K. Nozari, Eur. Phys. J. C 78, 486 (2018), arXiv:1806.08949 [gr-qc] .
  * Feng  _et al._ (2022) H. Feng, M. Li, G.-R. Liang, and R.-J. Yang, JCAP 04, 027 (2022), arXiv:2203.02924 [gr-qc] .
  * Zuluaga and Sánchez (2021) F. H. Zuluaga and L. A. Sánchez, Eur. Phys. J. C 81, 840 (2021), arXiv:2106.03140 [gr-qc] .
  * Ditta and Abbas (2020) A. Ditta and G. Abbas, Gen. Rel. Grav. 52, 77 (2020).
  * Uniyal  _et al._ (2022) A. Uniyal, R. C. Pantig, and A. Övgün, (2022), arXiv:2205.11072 [gr-qc] .
  * Chakhchi  _et al._ (2022) L. Chakhchi, H. El Moumni, and K. Masmar, Phys. Rev. D 105, 064031 (2022).
  * John (2019) A. J. John, Mon. Not. Roy. Astron. Soc. 490, 3824 (2019), arXiv:1603.09425 [gr-qc] .
  * Pérez  _et al._ (2017) D. Pérez, F. G. Lopez Armengol, and G. E. Romero, Phys. Rev. D 95, 104047 (2017), arXiv:1705.02713 [astro-ph.HE] .
  * Liu  _et al._ (2022) C. Liu, S. Yang, Q. Wu, and T. Zhu, JCAP 02, 034 (2022), arXiv:2107.04811 [gr-qc] .
  * Heydari-Fard  _et al._ (2021) M. Heydari-Fard, M. Heydari-Fard, and H. R. Sepangi, Eur. Phys. J. C 81, 473 (2021), arXiv:2105.09192 [gr-qc] .
  * Stashko  _et al._ (2021) O. S. Stashko, V. I. Zhdanov, and A. N. Alexandrov, Phys. Rev. D 104, 104055 (2021), arXiv:2107.05111 [gr-qc] .
  * Shaikh and Joshi (2019) R. Shaikh and P. S. Joshi, JCAP 10, 064 (2019), arXiv:1909.10322 [gr-qc] .
  * Van Aelst  _et al._ (2021) K. Van Aelst, E. Gourgoulhon, and F. H. Vincent, Phys. Rev. D 104, 124034 (2021), arXiv:2103.01827 [gr-qc] .
  * Bekenstein (1972) J. D. Bekenstein, Phys. Rev. D 5, 1239 (1972).
  * Wald (1984) R. M. Wald,  _General Relativity_ (1984).
  * Ashtekar and Magnon-Ashtekar (1979) A. Ashtekar and A. Magnon-Ashtekar, Journal of Mathematical Physics 20, 793 (1979).
  * Visser (2007) M. Visser, in  _Kerr Fest: Black Holes in Astrophysics, General Relativity and Quantum Gravity_ (2007) arXiv:0706.0622 [gr-qc] .
  * Vishveshwara (1968) C. V. Vishveshwara, Journal of Mathematical Physics 9, 1319 (1968).
  * Wald (1984) R. M. Wald,  _General relativity_ (Chicago Univ. Press, Chicago, IL, 1984).
  * Pugliese  _et al._ (2011) D. Pugliese, H. Quevedo, and R. Ruffini, Phys. Rev. D 83, 024021 (2011), arXiv:1012.5411 [astro-ph.HE] .
  * Psaltis (2008) D. Psaltis, Phys. Rev. D 77, 064006 (2008), arXiv:0704.2426 [astro-ph] .
  * Event Horizon Telescope Collaboration  _et al._ (2022) Event Horizon Telescope Collaboration, K. Akiyama, A. Alberdi, W. Alef, J. C. Algaba, R. Anantua, K. Asada, R. Azulay, U. Bach, A.-K. Baczko, D. Ball, M. Baloković, and Barrett, ApJ 930, L17 (2022).
  * Babichev  _et al._ (2011) E. Babichev, S. Chernov, V. Dokuchaev, and Y. Eroshenko, J. Exp. Theor. Phys. 112, 784 (2011), arXiv:0806.0916 [gr-qc] .
  * Freitas Pacheco (2011) J. A. d. Freitas Pacheco, (2011), arXiv:1109.6798 [gr-qc] .
  * Dadhich  _et al._ (2000) N. Dadhich, R. Maartens, P. Papadopoulos, and V. Rezania, Phys. Lett. B 487, 1 (2000), arXiv:hep-th/0003061 .
  * Zakharov (2022) A. F. Zakharov, Universe 8, 141 (2022), arXiv:2108.01533 [gr-qc] .
  * Neves (2020) J. C. S. Neves, Eur. Phys. J. C 80, 717 (2020), arXiv:2005.00483 [gr-qc] .
  * Babichev  _et al._ (2017) E. Babichev, C. Charmousis, and A. Lehébel, JCAP 04, 027 (2017), arXiv:1702.01938 [gr-qc] .
  * Hansen and Yunes (2013) D. Hansen and N. Yunes, Phys. Rev. D 88, 104020 (2013), arXiv:1308.6631 [gr-qc] .