# Coupled vector Gauss-Bonnet theories and hairy black holes

Katsuki Aoki1 and Shinji Tsujikawa2 1Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, 606-8502, Kyoto, Japan
2Department of Physics, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo
169-8555, Japan

###### Abstract

We study vector-tensor theories in which a 4-dimensional vector field
$A_{\mu}$ is coupled to a vector quantity ${\cal J}^{\mu}$, which is expressed
in terms of $A_{\mu}$ and a metric tensor $g_{\mu\nu}$. The divergence of
${\cal J}^{\mu}$ is equivalent to a Gauss-Bonnet (GB) term. We show that an
interacting Lagrangian of the form $f(X)A_{\mu}{\cal J}^{\mu}$, where $f$ is
an arbitrary function of $X=-(1/2)A_{\mu}A^{\mu}$, belongs to a scheme of
beyond generalized Proca theories. For $f(X)=\alpha={\rm constant}$, this
interacting Lagrangian reduces to a particular class of generalized Proca
theories. We apply the latter coupling to a static and spherically symmetric
vacuum configuration by incorporating the Einstein-Hilbert term, Maxwell
scalar, and vector mass term $\eta X$ ($\eta$ is a constant). Under an
expansion of the small coupling constant $\alpha$ with $\eta\neq 0$, we derive
hairy black hole solutions endowed with nonvanishing temporal and radial
vector field profiles. The asymptotic properties of solutions around the
horizon and at spatial infinity are different from those of hairy black holes
present in scalar-GB theories. We also show that black hole solutions without
the vector mass term, i.e., $\eta=0$, are prone to ghost instability of odd-
parity perturbations.

††preprint: YITP-23-38, WUCG-23-03

##  I Introduction

General Relativity (GR) has been well tested in solar-system experiments [1]
and submillimeter laboratory tests of gravity [2, 3]. If we go beyond the
solar-system scales, however, there are several unsolved problems such as the
origins of dark energy and dark matter [4, 5]. At very high energy close to
the Planck scale, we also believe that GR should be replaced by a more
fundamental theory with an ultraviolet completion. In such extreme gravity
regimes, it is expected that some higher-order curvature corrections to the
Einstein-Hilbert action come into play. These curvature corrections can
potentially modify the physics of highly compact objects like black holes
(BHs) and neutron stars. After the dawn of gravitational wave astronomy [6],
one can now probe signatures for the possible deviation from GR in strong
gravity regimes [7, 8, 9].

For the construction of healthy gravitational theories, it is desirable to
keep the field equations of motion up to second order in the metric tensor
$g_{\mu\nu}$. In this case, one can avoid so-called Ostrogradski instability
[10, 11] arising from higher-order derivative terms. Using a general class of
Lagrangians containing polynomial functions of Riemann curvature tensors,
Lanczos [12] and Lovelock [13] constructed gravitational theories with second-
order field equations of motion. In the 4-dimensional spacetime, the field
equations uniquely reduce to those in GR. In spacetime dimensions $D$ higher
than 4, Lanczos and Lovelock theories differ from GR and have richer
structures. In particular, there is a quadratic-order curvature scalar known
as a Gauss-Bonnet (GB) term modifying the spacetime dynamics in $D>4$
dimensions [14].

When $D=4$, the GB term is a topological surface term which does not
contribute to the field equations of motion. If there is a scalar field $\phi$
coupled to the GB term $\mathcal{G}$ of the form $\mu(\phi)\mathcal{G}$, where
$\mu$ is a function of $\phi$, the 4-dimensional spacetime dynamics is
modified by the scalar-GB coupling. In string theory, for example, the low
energy effective action contains a coupling between the dilaton field $\phi$
and $\mathcal{G}$ of the form $e^{-\lambda\phi}\mathcal{G}$ [15, 16, 17]. If
we apply the scalar-GB coupling $\mu(\phi)\mathcal{G}$ to a spherically
symmetric configuration, it is known that there are hairy BH and neutron star
solutions with nontrivial scalar field profiles [18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. The role of the same
scalar-GB coupling in cosmology has been also extensively studied in the
literature [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64]. We note that the extension to
more general scalar-GB couplings $f(\phi,\mathcal{G})$ containing nonlinear
functions of $\mathcal{G}$ leads to instabilities of scalar perturbations
associated with the nonlinear GB term during decelerating cosmological epochs
[65] (see also Refs. [66, 67, 68, 69, 70, 71, 72, 69]).

The linear scalar-GB coupling $\phi\mathcal{G}$ has a peculiar property in
four dimensions. The action $\int{\rm d}^{4}x\sqrt{-g}\,\phi\mathcal{G}$ is
invariant under the global shift of the scalar field $\phi\to\phi+c$ with $c$
being constant, as the integral $c\int{\rm d}^{4}x\sqrt{-g}\,\mathcal{G}$ is a
boundary term. When a theory enjoys the global symmetry, one may promote it to
a local symmetry by introducing a gauge field $A_{\mu}$ and replacing
derivatives with covariant derivatives. When the shift symmetry is localized,
the scalar field $\phi$ corresponds to a gauge mode and can be eliminated by
fixing the gauge. Then, the resulting theory is a vector-tensor theory having
three dynamical degrees of freedom (DOFs) on top of the DOFs of spacetime
metric $g_{\mu\nu}$ (cf., Refs. [73, 74, 75]). In this paper, we shall apply
this idea to the linear scalar-GB coupling and find a vector-tensor theory
analogous to the scalar-GB theory.

In practice, the above procedure requires to find a vector quantity ${\cal
J}^{\mu}$ whose divergence is equivalent to $\mathcal{G}$, i.e.,
$\nabla_{\mu}{\cal J}^{\mu}=\mathcal{G}$, by which the linear scalar-GB
coupling can be recast in the form $-\int{\rm
d}^{4}x\sqrt{-g}\,\nabla_{\mu}\phi{\cal J}^{\mu}$ via integration by parts.
The integral ${\cal J}^{\mu}$ may not be unique since it is defined only
through the differential equation $\nabla_{\mu}{\cal J}^{\mu}=\mathcal{G}$.
One form of ${\cal J}^{\mu}$, which is expressed in terms of a scalar field
and Riemann tensor, is found in Ref. [76]. On using this expression of ${\cal
J}^{\mu}={\cal J}^{\mu}[\phi,g]$ and the property
$\mathcal{G}=\nabla_{\mu}{\cal J}^{\mu}$, it is possible to prove that the
scalar-GB coupling $\mu(\phi)\mathcal{G}$ belongs to a subclass of Horndeski
theories [77] after the integration by parts [34]. We note that the
equivalence between the scalar-GB coupling and Horndeski theories was
originally shown in Ref. [78] by taking the approach of field equations of
motion.

We will find an alternative expression of ${\cal J}^{\mu}$ by using a vector
field $A_{\mu}$, where ${\cal J}^{\mu}={\cal J}^{\mu}[A,g]$ satisfies the same
relation $\nabla_{\mu}{\cal J}^{\mu}=\mathcal{G}$. As a candidate for a
Lorentz-invariant scalar characterizing the coupling between $A_{\mu}$ and the
integrated GB term in vector-tensor theories, we propose the Lagrangian
$A_{\mu}{\cal J}^{\mu}$. We will show that the interacting Lagrangian
$A_{\mu}{\cal J}^{\mu}$ is equivalent to a subclass of generalized Proca (GP)
theories with second-order field equations of motion [79, 80, 81, 82, 83], and
by construction, it reduces to a linear scalar-GB coupling in a certain limit.
We will also extend the analysis to a more general Lagrangian
$f(X)A_{\mu}{\cal J}^{\mu}$, where $f$ is an arbitrary function of
$X=-(1/2)A_{\mu}A^{\mu}$. In this case, the resulting vector-tensor theory is
shown to be equivalent to a class of beyond GP theories originally proposed in
Ref. [84] (see also Ref. [85]). Since beyond GP theories correspond to a
healthy extension of GP theories with the same dynamical DOFs, we are now able
to construct healthy theories of a vector field coupled to the integrated GB
term.

We will also apply the interacting Lagrangian $\alpha A_{\mu}{\cal J}^{\mu}$
($\alpha$ is a coupling constant) to the search for hairy BH solutions on a
static and spherically symmetric background. For this purpose, we also take
into account the Einstein-Hilbert term, Maxwell term, and vector mass term
$\eta X$, where $\eta$ is a constant. The coupling $\alpha A_{\mu}{\cal
J}^{\mu}$ is equivalent to a Lagrangian of the quintic-order coupling function
$G_{5}(X)=4\alpha\ln|X|$ in GP theories. In Refs. [86, 87], it was shown that
there are no hairy BH solutions with regular vector field profiles for the
positive power-law quintic functions $G_{5}\propto X^{n}$ with $n\geq 1$.
However, we will show that this is not the case for $G_{5}(X)=4\alpha\ln|X|$.
Under an expansion of the small coupling constant $\alpha$, we derive
solutions to the temporal and radial vector components around the BH horizon
and at spatial infinity. Numerically it is challenging to perform accurate
integrations due to the existence of a rapidly growing mode arising from the
mass term $\eta X$, but we are able to find out solutions that mimic the
asymptotic behavior in some large-distance regions. In comparison to hairy BHs
present for the linear scalar-GB coupling, the behavior of hairy BH solutions
for $\eta\neq 0$ is different both around the horizon and at large distances.
We will also show that BH solutions for $\eta=0$ suffer from ghost instability
of odd-parity perturbations.

This paper is organized as follows. In Sec. II, we first review the
correspondence between the scalar-GB coupling $\mu(\phi)\mathcal{G}$ and the
Horndeski Lagrangian. We then introduce a vector field ${\cal J}^{\mu}$ whose
divergence $\nabla_{\mu}{\cal J}^{\mu}$ is equivalent to the GB term and show
that the Lagrangian $f(X)A_{\mu}{\cal J}^{\mu}$ belongs to a subclass of
beyond GP theories. In Sec. III, we study static and spherically symmetric BH
solutions for the coupling $\alpha A_{\mu}{\cal J}^{\mu}$ and derive
perturbative solutions to $A_{\mu}$ with respect to the small coupling
$\alpha$. We also numerically confirm the existence of vector field profiles
connecting two asymptotic regimes and analytically estimate corrections to the
gravitational potentials arising from the vector-GB couplings. Sec. IV is
devoted to conclusions.

##  II Coupled vector Gauss-Bonnet theories

The GB curvature invariant is a specific combination of Lanczos [12] and
Lovelock [13] scalars. In 4-dimensional spacetime, the GB term is given by
[88, 89]

$\mathcal{G}=\frac{1}{4}\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}R^{\alpha\beta}{}_{\mu\nu}R^{\gamma\delta}{}_{\rho\sigma}\,,$
(1)

where
$\delta^{\mu_{1}\cdots\mu_{k}}_{\nu_{1}\cdots\nu_{k}}=k!\delta^{[\mu_{1}}_{\nu_{1}}\cdots\delta^{\mu_{k}]}_{\nu_{k}}$
is the generalized Kronecker delta and $R^{\alpha\beta}{}_{\mu\nu}$ is the
Riemann tensor. More explicitly, Eq. (1) can be expressed as

$\mathcal{G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}\,,$
(2)

where $R$ is the scalar curvature and $R_{\mu\nu}$ is the Ricci tensor. In 4
dimensions, the GB term is a total derivative and does not contribute to the
equations of motion while the 4-dimensional spacetime dynamics is modified in
the presence of a scalar or vector field coupled to $\mathcal{G}$ or its
associated vector.

###  II.1 Scalar field coupled to the GB term

Let us first briefly revisit the case in which there is a scalar field $\phi$
coupled to the GB term of the form $\mu(\phi)\mathcal{G}$. Because of the
antisymmetric property of
$\delta^{\mu_{1}\cdots\mu_{k}}_{\nu_{1}\cdots\nu_{k}}$, the field equations of
motion following from the coupling $\mu(\phi)\mathcal{G}$ are of second order
in the metric tensor $g_{\mu\nu}$ and the scalar field $\phi$. On using the
Riemann tensor and covariant derivatives of $\phi$, the GB term can be
expressed in the form [76, 34]

$\mathcal{G}=\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\nabla^{\delta}\left[\frac{\nabla^{\gamma}\nabla_{\rho}\phi\nabla_{\sigma}\phi}{X_{s}}\left({R^{\alpha\beta}}_{\mu\nu}-\frac{2}{3X_{s}}\nabla^{\alpha}\nabla_{\mu}\phi\nabla^{\beta}\nabla_{\nu}\phi\right)\right]\,,$
(3)

where $\nabla^{\delta}$ is a covariant derivative operator and
$X_{s}=-(1/2)\nabla_{\mu}\phi\nabla^{\mu}\phi$. Substituting this expression
of $\mathcal{G}$ into the action

${\cal S}_{\rm sGB}=\int{\rm d}^{4}x\sqrt{-g}\,\mu(\phi)\mathcal{G}\,,$ (4)

and integrating (4) by parts, it follows that

${\cal S}_{\rm sGB}=-\int{\rm
d}^{4}x\sqrt{-g}\,\mu_{,\phi}(\phi)\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{1}{X_{s}}\nabla^{\gamma}\nabla_{\rho}\phi\nabla_{\sigma}\phi\nabla^{\delta}\phi\left({R^{\alpha\beta}}_{\mu\nu}-\frac{2}{3X_{s}}\nabla^{\alpha}\nabla_{\mu}\phi\nabla^{\beta}\nabla_{\nu}\phi\right)\,,$
(5)

where $g$ is the determinant of $g_{\mu\nu}$, and we use the notations
$F_{,\phi}\equiv\partial F/\partial\phi$ and $F_{,X_{s}}\equiv\partial
F/\partial X_{s}$ for any arbitrary function $F$. We will expand the
generalized Kronecker delta, integrate the action (5) by parts, and exploit
the relation
$[\nabla_{\mu},\nabla_{\nu}]\nabla^{\alpha}\phi={R^{\alpha}}_{\lambda\mu\nu}\nabla^{\lambda}\phi$
to eliminate contractions of the Riemann tensors. Up to boundary terms, the
action (5) is equivalent to [34]

$\displaystyle{\cal S}_{\rm sGB}$ $\displaystyle=$ $\displaystyle\int{\rm
d}^{4}x\sqrt{-g}\,\biggl{[}G_{2s}-G_{3s}\square\phi+G_{4s}R+G_{4s,X_{s}}\left\\{(\square\phi)^{2}-(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)\right\\}$
(6)
$\displaystyle+G_{5s}G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi-\frac{1}{6}G_{5s,X_{s}}\left\\{(\square\phi)^{3}-3(\square\phi)\,(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)+2(\nabla^{\mu}\nabla_{\alpha}\phi)(\nabla^{\alpha}\nabla_{\beta}\phi)(\nabla^{\beta}\nabla_{\mu}\phi)\right\\}\biggr{]}\,,$

where $G_{\mu\nu}=R_{\mu\nu}-(1/2)g_{\mu\nu}R$ is the Einstein tensor, and

$\displaystyle
G_{2s}=-8\mu_{,\phi\phi\phi\phi}(\phi)X_{s}^{2}(3-\ln|X_{s}|)\,,\qquad
G_{3s}=4\mu_{,\phi\phi\phi}(\phi)X_{s}(7-3\ln|X_{s}|)\,,$ $\displaystyle
G_{4s}=4\mu_{,\phi\phi}(\phi)X_{s}(2-\ln|X_{s}|)\,,\qquad
G_{5s}=-4\mu_{,\phi}(\phi)\ln|X_{s}|\,.$ (7)

The action (6) belongs to a subclass of scalar Horndeski theories [77] with
second-order Euler equations of motion. Originally, the equivalence of scalar-
GB theories with Horndeski theories given by the coupling functions (7) was
shown in Ref. [78] by using the field equations of motion. In Ref. [34], the
same equivalence was proven at the level of the action (as explained above).

For the linear scalar-GB coupling $\mu(\phi)=-\alpha\phi$, where $\alpha$ is a
constant, we have $G_{5s}=4\alpha\ln|X_{s}|$ and $G_{2s}=G_{3s}=G_{4s}=0$.
This falls into a subclass of shift-symmetric Horndeski theories where the
field equations of motion are invariant under the shift $\phi\to\phi+c$. In
the original form of the GB coupling $\phi\mathcal{G}$, the action is quasi-
invariant under the shift $\phi\to\phi+c$, i.e., invariant up to a total
derivative, while the Lagrangian in the Horndeski form is manifestly invariant
under the shift. For $\mu(\phi)$ containing nonlinear functions of $\phi$, we
generally have the $\phi$ dependence in $G_{2s},G_{3s},G_{4s},G_{5s}$. As we
mentioned in Introduction, there are BH and neutron star solutions endowed
with scalar hairs for such scalar-GB couplings.

###  II.2 Vector field coupled to the integrated GB term

If we want to construct theories in which a vector field $A_{\mu}$ is coupled
to the GB term in some way, we need to construct a Lorentz-invariant scalar
appearing in the Lagrangian. For instance, one may consider the coupling
$A^{\mu}\nabla_{\mu}\mathcal{G}$. However, the equations of motion associated
with this coupling contain derivatives higher than second order and hence such
theories are generally prone to Ostrogradski instability. Another possible
coupling would be $\mu_{\rm v}(X)\mathcal{G}$ where

$X\equiv-\frac{1}{2}A_{\mu}A^{\mu}\,.$ (8)

Again, this interaction may summon a ghostly DOF in the longitudinal sector of
$A_{\mu}$ which can be understood by taking the decoupling of the longitudinal
and transverse DOFs. The longitudinal mode becomes manifest by introducing the
Stückelberg field according to the replacement $A_{\mu}\to g_{\rm
v}A_{\mu}+\nabla_{\mu}\phi$. Then, the decoupling limit $g_{\rm v}\to 0$ gives
$X\to X_{s}$. The interacting Lagrangian $\mu_{\rm v}(X)\mathcal{G}$ reduces
to a coupling between the GB term and the derivative of $\phi$, not the scalar
field itself, which should yield equations of motion with derivatives higher
than second order.

As we already explained, the linear coupling $-\phi\mathcal{G}$ has a global
shift symmetry and the vector-tensor theory can be obtained by localizing this
global symmetry. This requires finding a vector field $\mathcal{J}^{\mu}$
whose divergence agrees with the GB term,
$\nabla_{\mu}\mathcal{J}^{\mu}=\mathcal{G}$. After integration by parts, the
coupling $-\phi\mathcal{G}$ becomes
$\nabla_{\mu}\phi\mathcal{J}^{\mu}[\phi,g]$. As shown in Eq. (5), the action
contains the derivatives of $\phi$ but not the field itself, for the linear
coupling $\mu_{,\phi}=-\alpha$. The global shift symmetry can be localized by
the replacement $\nabla_{\mu}\phi\to\nabla_{\mu}\phi+g_{\rm v}A_{\mu}$ with
the help of the vector field $A_{\mu}$. The symmetry transformation is now
$\phi\to\phi+\chi(x),A_{\mu}\to A_{\mu}-\nabla_{\mu}\chi(x)/g_{\rm v}$. The
scalar field $\phi$ can be eliminated by setting the unitary gauge $\phi=0$.
All in all, what we need is the replacement $\nabla_{\mu}\phi\to A_{\mu}$,
where the gauge coupling $g_{\rm v}$ is absorbed into the definition of
$A_{\mu}$. In this way, we can obtain a vector-tensor theory coupled to the
integrated GB term.

However, there is an ambiguity in the above procedure. The second derivative
$\nabla_{\mu}\nabla_{\nu}\phi$ is symmetric in its indices while the replaced
quantity $\nabla_{\mu}A_{\nu}$ is not symmetric. We resolve this ambiguity by
imposing the condition $\nabla_{\mu}\mathcal{J}^{\mu}=\mathcal{G}$ even after
the replacement $\nabla_{\mu}\phi\to A_{\mu}$. The expression of
$\mathcal{J}^{\mu}$ consistent with such requirement is given by

$\displaystyle\mathcal{J}^{\mu}\equiv\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\left[\frac{A^{\alpha}\nabla_{\nu}A^{\beta}}{X}\left(R^{\gamma\delta}{}_{\rho\sigma}-\frac{2}{3X}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}\right)\right]\,.$
(9)

In the following, we will show that this vector field satisfies the relation
$\nabla_{\mu}\mathcal{J}^{\mu}=\mathcal{G}$. In doing so, we use the relation
$\nabla_{[\mu}R^{\alpha\beta}{}_{\nu\rho]}=0$ and the antisymmetric property
of the generalized Kronecker delta. We also exploit the following equality

$\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}\nabla_{\mu_{1}}\left(\frac{A^{\nu_{1}}}{X^{2}}\right)\nabla_{\mu_{2}}A^{\nu_{2}}\nabla_{\mu_{3}}A^{\nu_{3}}\nabla_{\mu_{4}}A^{\nu_{4}}=-\frac{1}{2X^{3}}\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}\nu_{5}}A_{\mu_{1}}A^{\nu_{1}}\nabla_{\mu_{2}}A^{\nu_{2}}\nabla_{\mu_{3}}A^{\nu_{3}}\nabla_{\mu_{4}}A^{\nu_{4}}\nabla_{\mu_{5}}A^{\nu_{5}}=0\,,$
(10)

together with the expansion of
$\delta^{\mu_{1}\cdots\mu_{d}}_{\nu_{1}\cdots\nu_{d}}$ in $d$ dimensions:

$\delta^{\mu_{1}\cdots\mu_{d}}_{\nu_{1}\cdots\nu_{d}}=\sum_{k=1}^{d}(-1)^{d+k}\delta^{\mu_{d}}_{\nu_{k}}\delta^{\mu_{1}\cdots\mu_{k}\cdots\mu_{d-1}}_{\nu_{1}\cdots\bar{\nu}_{k}\cdots\nu_{d}}\,,$
(11)

where $\bar{\nu}_{k}$ means that this index is omitted. Notice that the
5-dimensional generalized Kronecker delta
$\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}\nu_{5}}$
vanishes in 4-dimensional spacetime, whose property was used in the second
equality of Eq. (10). Then, it follows that

$\nabla_{\mu}\mathcal{J}^{\mu}={\cal F}_{1}+{\cal F}_{2}+{\cal F}_{3}\,,$ (12)

where

$\displaystyle{\cal F}_{1}$ $\displaystyle=$
$\displaystyle\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{A^{\alpha}}{X}\nabla_{\mu}\nabla_{\nu}A^{\beta}R^{\gamma\delta}{}_{\rho\sigma}\,,$
(13) $\displaystyle{\cal F}_{2}$ $\displaystyle=$
$\displaystyle\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\nabla_{\mu}\left(\frac{A^{\alpha}}{X}\right)\nabla_{\nu}A^{\beta}R^{\gamma\delta}{}_{\rho\sigma}\,,$
(14) $\displaystyle{\cal F}_{3}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{A^{\alpha}}{X^{2}}\nabla_{\mu}\left(\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}\right)\,.$
(15)

Since
$\nabla_{\mu}\nabla_{\nu}A^{\beta}-\nabla_{\nu}\nabla_{\mu}A^{\beta}={R^{\beta}}_{\lambda\mu\nu}A^{\lambda}$,
Eq. (13) reduces to

${\cal
F}_{1}=\frac{1}{2}\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{A^{\alpha}A^{\lambda}}{X}{R^{\beta}}_{\lambda\mu\nu}R^{\gamma\delta}{}_{\rho\sigma}=\frac{1}{4}\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}{R^{\alpha\beta}}_{\mu\nu}R^{\gamma\delta}{}_{\rho\sigma}=\mathcal{G}\,,$
(16)

where, in the second equality, we used the relation

$\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}\left(\frac{2A^{\nu_{1}}A^{\lambda}}{X}R^{\nu_{2}}{}_{\lambda\mu_{1}\mu_{2}}R^{\nu_{3}\nu_{4}}{}_{\mu_{3}\mu_{4}}-R^{\nu_{1}\nu_{2}}{}_{\mu_{1}\mu_{2}}R^{\nu_{3}\nu_{4}}{}_{\mu_{3}\mu_{4}}\right)=\frac{1}{2X}\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}\nu_{5}}A_{\mu_{1}}A^{\nu_{1}}R^{\nu_{2}\nu_{3}}{}_{\mu_{2}\mu_{3}}R^{\nu_{4}\nu_{5}}{}_{\mu_{4}\mu_{5}}=0\,.$
(17)

For the computation of ${\cal F}_{2}$, we exploit the following property

$\displaystyle\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}\left[\nabla_{\mu_{1}}\left(\frac{A^{\nu_{1}}}{X}\right)\nabla_{\mu_{2}}A^{\nu_{2}}R^{\nu_{3}\nu_{4}}{}_{\mu_{3}\mu_{4}}-\frac{A^{\nu_{1}}A^{\alpha}}{X^{2}}R^{\nu_{2}}{}_{\alpha\mu_{1}\mu_{2}}\nabla_{\mu_{3}}A^{\nu_{3}}\nabla_{\mu_{4}}A^{\nu_{4}}\right]$
$\displaystyle=-\frac{1}{2X^{2}}\delta^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}\nu_{5}}A_{\mu_{1}}A^{\nu_{1}}R^{\nu_{2}\nu_{3}}{}_{\mu_{2}\mu_{3}}\nabla_{\mu_{4}}A^{\nu_{4}}\nabla_{\mu_{5}}A^{\nu_{5}}=0\,.$
(18)

Then, we have

${\cal
F}_{2}=\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{A^{\alpha}A^{\lambda}}{X^{2}}{R^{\beta}}_{\lambda\mu\nu}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}\,.$
(19)

Expanding the covariant derivative $\nabla_{\mu}$ in ${\cal F}_{3}$ and using
the commutation relation of $\nabla_{\mu}\nabla_{\nu}A^{\beta}$, we obtain

${\cal
F}_{3}=-\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{A^{\alpha}A^{\lambda}}{X^{2}}{R^{\beta}}_{\lambda\mu\nu}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}=-{\cal
F}_{2}\,.$ (20)

On using Eqs. (16) and (20), Eq. (12) yields

$\nabla_{\mu}\mathcal{J}^{\mu}=\mathcal{G}\,.$ (21)

Thus, the divergence of $\mathcal{J}^{\mu}$ is equivalent to the GB term.

We consider a scalar quantity $A_{\mu}\mathcal{J}^{\mu}$ as a candidate for
the ghost-free Lagrangian of the vector field $A_{\mu}$ coupled to the
integrated GB term. As a generalization, we also multiply an arbitrary
function $f(X)$ of $X$ with the scalar product $A_{\mu}\mathcal{J}^{\mu}$. The
interaction part of the action in such theories is given by

$\displaystyle{\cal S}_{\rm vGB}=\int{\rm
d}^{4}x\sqrt{-g}\,f(X)A_{\mu}\mathcal{J}^{\mu}\,,$ (22)

which is composed of two parts:

$\displaystyle{\cal S}_{{\rm vGB}1}$ $\displaystyle=$ $\displaystyle\int{\rm
d}^{4}x\sqrt{-g}\,\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{f(X)}{X}A_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}{R^{\gamma\delta}}_{\rho\sigma}\,,$
(23) $\displaystyle{\cal S}_{{\rm vGB}2}$ $\displaystyle=$
$\displaystyle-\int{\rm
d}^{4}x\sqrt{-g}\,\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\frac{2f(X)}{3X^{2}}A_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}\,.$
(24)

On using the properties $\nabla_{\mu}X=-A_{\nu}\nabla_{\mu}A^{\nu}$ and
$\nabla_{\mu}\ln X=-A_{\nu}\nabla_{\mu}A^{\nu}/X$, we find that Eqs. (23) and
(24) reduce, respectively, to

$\displaystyle{\cal S}_{{\rm vGB}1}$ $\displaystyle=$ $\displaystyle\int{\rm
d}^{4}x\sqrt{-g}\,\biggl{[}G_{5}(X)G_{\mu\nu}\nabla^{\mu}A^{\nu}$ (25)
$\displaystyle\qquad\qquad\qquad-\frac{4f(X)}{X}\left(A^{\rho}\nabla_{\sigma}A^{\sigma}-A^{\sigma}\nabla_{\sigma}A^{\rho}\right)\nabla_{\nu}\nabla_{\rho}A^{\nu}-\frac{4f(X)}{X}\left(A^{\sigma}\nabla_{\sigma}A^{\nu}-A^{\nu}\nabla_{\sigma}A^{\sigma}\right)\nabla_{\nu}\nabla_{\rho}A^{\rho}\biggr{]}\,,$
$\displaystyle{\cal S}_{{\rm vGB}2}$ $\displaystyle=$ $\displaystyle\int{\rm
d}^{4}x\sqrt{-g}\,\biggl{[}-\frac{2f(X)}{3X}\delta^{\mu\nu\rho}_{\alpha\beta\gamma}\nabla_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}-3f_{5}(X)\delta^{\mu\nu\rho}_{\alpha\beta\gamma}A^{\alpha}A_{\lambda}\nabla_{\mu}A^{\lambda}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}$
(26)
$\displaystyle\qquad\qquad\qquad+\frac{4f(X)}{X}\left(A^{\rho}\nabla_{\sigma}A^{\sigma}-A^{\sigma}\nabla_{\sigma}A^{\rho}\right)\nabla_{\nu}\nabla_{\rho}A^{\nu}+\frac{4f(X)}{X}\left(A^{\sigma}\nabla_{\sigma}A^{\nu}-A^{\nu}\nabla_{\sigma}A^{\sigma}\right)\nabla_{\nu}\nabla_{\rho}A^{\rho}\biggr{]}\,,$

where

$G_{5}(X)\equiv 4\int{\rm
d}\tilde{X}\frac{f(\tilde{X})}{\tilde{X}}+8f(X)\,,\qquad
f_{5}(X)\equiv-\frac{2f_{,X}}{3X}\,.$ (27)

Taking the sum of ${\cal S}_{{\rm vGB}1}$ and ${\cal S}_{{\rm vGB}2}$, terms
on the second lines of Eqs. (25) and (26) cancel each other. On using the
property

$f_{5}\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}A_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}=\frac{4}{3}f_{,X}\delta^{\mu\nu\rho}_{\alpha\beta\gamma}\nabla_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}-3f_{5}\delta^{\mu\nu\rho}_{\alpha\beta\gamma}A^{\alpha}A_{\lambda}\nabla_{\mu}A^{\lambda}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\,,$
(28)

we obtain the reduced action

${\cal S}_{\rm vGB}=\int{\rm
d}^{4}x\sqrt{-g}\left[G_{5}(X)G_{\mu\nu}\nabla^{\mu}A^{\nu}-\frac{1}{6}G_{5,X}(X)\delta^{\mu\nu\rho}_{\alpha\beta\gamma}\nabla_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}+f_{5}(X)\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}A_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\nabla_{\sigma}A^{\delta}\right]\,.$
(29)

This belongs to a subclass of beyond GP theories proposed in Ref. [84] and
thus it is free from the Ostrogradski ghost.

In theories where the function $f(X)$ is constant, i.e.,

$f(X)=\alpha={\rm constant}\,,$ (30)

we have

$G_{5}(X)=4\alpha\ln|X|\,,\qquad f_{5}(X)=0\,.$ (31)

Note that we dropped a constant term $8\alpha$ in $G_{5}(X)$, as it does not
contribute to the field equations of motion. Since the last term on the right
hand side of Eq. (29) vanishes, the action ${\cal S}_{\rm vGB}$ belongs to a
subclass of GP theories [79, 80, 82]. We recall that, from Eq. (7), the linear
scalar-GB coupling $\mu(\phi)\mathcal{G}$ with $\mu(\phi)=-\alpha\phi$
corresponds to $G_{5s}(X_{s})=4\alpha\ln|X_{s}|$ with
$G_{2s}(X_{s})=G_{3s}(X_{s})=G_{4s}(X_{s})=0$. The coupling (31) is the
vector-tensor analogue to the linear scalar-GB coupling in Horndeski theories.

In the following, we will focus on GP theories given by the functions (31).
Varying the action (29) with respect to $A_{\mu}$, it follows that

$\displaystyle\frac{\delta{\cal S}_{\rm vGB}}{\delta
A_{\mu}}=\alpha\left(\mathcal{J}^{\mu}+\mathcal{J}_{F}^{\mu}\right)\,,$ (32)

where $\mathcal{J}^{\mu}$ is defined by Eq. (9), and

$\displaystyle\mathcal{J}_{F}^{\mu}\equiv\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}A_{\nu}F^{\alpha\beta}\left(\frac{1}{X^{2}}\nabla^{\gamma}A_{\rho}\nabla^{\delta}A_{\sigma}-\frac{1}{2X}R^{\gamma\delta}{}_{\rho\sigma}\right)\,,$
(33)

with $F^{\alpha\beta}\equiv\nabla^{\alpha}A^{\beta}-\nabla^{\beta}A^{\alpha}$.
The sum of $\mathcal{J}^{\mu}$ and $\mathcal{J}_{F}^{\mu}$ can be expressed in
a compact form

$\mathcal{J}^{\mu}+\mathcal{J}_{F}^{\mu}=\delta^{\mu\nu\rho\sigma}_{\alpha\beta\gamma\delta}\left[\frac{A^{\alpha}\nabla^{\beta}A_{\nu}}{X}\left(R^{\gamma\delta}{}_{\rho\sigma}-\frac{2}{3X}\nabla^{\gamma}A_{\rho}\nabla^{\delta}A_{\sigma}\right)\right]\,.$
(34)

Let us consider the case in which the Maxwell term
$F\equiv-(1/4)F^{\alpha\beta}F_{\alpha\beta}$ and the mass term $X$ are
present in addition to the action (29) with $f(X)=\alpha$. The action in such
a subclass of GP theories is given by

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{1}{g_{\rm v}^{2}}F+\eta
X+G_{5}(X)G_{\mu\nu}\nabla^{\mu}A^{\nu}-\frac{1}{6}G_{5,X}(X)\delta^{\mu\nu\rho}_{\alpha\beta\gamma}\nabla_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\right]\,,$
(35)

where $g_{\rm v}$ and $\eta$ are constants and $G_{5}(X)=4\alpha\ln|X|$.
Varying this action with respect to $A_{\mu}$, it follows that

$\frac{1}{g_{\rm v}^{2}}\nabla_{\nu}F^{\mu\nu}+\eta
A^{\mu}=\alpha(\mathcal{J}^{\mu}+\mathcal{J}_{F}^{\mu})\,.$ (36)

While the GB term does not explicitly show up in Eq. (36), it appears by
taking the divergence of Eq. (36) as

$\eta\nabla_{\mu}A^{\mu}=\alpha\mathcal{G}+\alpha\nabla_{\mu}\mathcal{J}_{F}^{\mu}\,,$
(37)

where we used the relation $\nabla_{\mu}\nabla_{\nu}F^{\mu\nu}=0$ and Eq.
(21).

Taking the decoupling limit $g_{\rm v}\to 0$ with the replacement $A^{\mu}\to
g_{\rm v}A^{\mu}+\nabla^{\mu}\phi$, we have $\mathcal{J}_{F}^{\mu}\to 0$ and
hence Eqs. (36) and (37) reduce to the Maxwell equation,
$\nabla_{\nu}F^{\mu\nu}=0$, and the equation of motion for the scalar field,
$\eta\nabla_{\mu}\nabla^{\mu}\phi=\alpha\mathcal{G}$, respectively. This
latter scalar field equation also follows by varying the Lagrangian $L=\eta
X_{s}-\alpha\phi\mathcal{G}$ with respect to $\phi$, so the linearly coupled
scalar-GB theory is recovered by taking the above decoupling limit. In shift-
symmetric Horndeski theories, using the expression of $\mathcal{G}$ in Eq. (3)
shows that the scalar field equation can be expressed in the form
$\nabla_{\mu}j^{\mu}=0$, where $j^{\mu}$ is a conserved current. When this
equation is solved for $j^{\mu}$, there is an integration constant
corresponding to boundary/initial conditions of the system.

In vector-tensor theories the equation of motion for the vector field
$A^{\mu}$ is given by Eq. (36), which does not contain an integration
constant. Although Eq. (37) corresponds to the differential version of Eq.
(36), one cannot choose an arbitrary integration constant when integrating Eq.
(37). This property is different from that in shift-symmetric Horndeski
theories discussed above. If we apply GP theories to the isotropic and
homogeneous cosmological background, the temporal vector component $A_{0}$ is
always related to the Hubble expansion rate $H$ [90, 91, 92]. This is known as
a tracker solution, in which case we do not have a freedom of changing initial
conditions of the vector field. In scalar-tensor theories, on the other hand,
it is possible to choose initial conditions away from the tracker because of
the existence of the integration constant said above [93, 94]. As a result,
one can distinguish between shift-symmetric Horndeski theories and GP theories
from the background cosmological dynamics.

If we apply linear scalar-GB theory to the static and spherically symmetric
background in vacuum, there are hairy BH solutions satisfying the boundary
condition $X_{s}=0$ on the horizon [24, 25]. This boundary condition fixes the
integration constant mentioned above [36]. In vector-GB theory, the similar
boundary condition, like $X=0$, cannot be necessarily imposed because of the
absence of an arbitrary constant in Eq. (36). This implies that the BH
solution in vector-GB theory should be different from that in linear scalar-GB
theory. In Sec. III, we will investigate the property of hairy BH solutions in
vector-GB theory.

##  III Black hole solutions in vector-GB theory

We study BH solutions by incorporating the Einstein-Hilbert term in the action
(35). The corresponding action belongs to a subclass of GP theories given by

$\displaystyle{\cal S}$ $\displaystyle=\int{\rm
d}^{4}x\sqrt{-g}\left[\frac{M_{\rm Pl}^{2}}{2}R+F+\eta X+\alpha A_{\mu}{\cal
J}^{\mu}\right]$ $\displaystyle=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm
Pl}^{2}}{2}R+F+\eta
X+(4\alpha\ln|X|)G_{\mu\nu}\nabla^{\mu}A^{\nu}-\frac{2\alpha}{3X}\delta^{\mu\nu\rho}_{\alpha\beta\gamma}\nabla_{\mu}A^{\alpha}\nabla_{\nu}A^{\beta}\nabla_{\rho}A^{\gamma}\right]\,,$
(38)

where $M_{\rm Pl}$ is the reduced Planck mass and we set $g_{\rm v}=1$. Let us
consider the static and spherically symmetric background given by the line
element

${\rm d}s^{2}=-f(r){\rm d}t^{2}+h^{-1}(r){\rm d}r^{2}+r^{2}\left({\rm
d}\theta^{2}+\sin^{2}\theta\,{\rm d}\varphi^{2}\right)\,,$ (39)

where $t$, $r$ and $(\theta,\varphi)$ represent the time, radial, and angular
coordinates, respectively, and $f$ and $h$ are functions of $r$. We will focus
on the case of positive mass squared $\eta>0$ as in standard massive Proca
theory. For the vector field, we consider the following configuration

$A_{\mu}=\left[A_{0}(r),A_{1}(r),0,0\right]\,,$ (40)

where $A_{0}$ and $A_{1}$ are functions of $r$. On the background (39), we
have

$X=\frac{A_{0}^{2}}{2f}-\frac{hA_{1}^{2}}{2}\,,\qquad F=\frac{hA_{0}^{\prime
2}}{2f}\,,$ (41)

where a prime represents the derivative with respect to $r$. Note that
$F=-F_{\mu\nu}F^{\mu\nu}/4\neq 0$ implies a nonvanishing temporal vector
component; that is, $A_{\mu}$ cannot be expressed by a scalar gradient. Hence,
the presence of a nontrivial temporal component $A_{0}(r)$ is essential to
differentiate solutions in vector-GB theory from those in scalar-GB theory.

Varying the action (38) with respect to $f$ and $h$, respectively, we obtain

$\displaystyle M_{\rm
Pl}^{2}f\left(rh^{\prime}+h-1\right)+\frac{r^{2}hA_{0}^{\prime
2}}{2}+\frac{\eta r^{2}}{2}\left(A_{0}^{2}+fhA_{1}^{2}\right)+\frac{4\alpha
f}{(A_{0}^{2}-fhA_{1}^{2})^{2}}[2A_{0}^{3}A_{0}^{\prime}A_{1}h(h-1)+2A_{0}A_{0}^{\prime}A_{1}^{3}fh^{2}(h+1)$
$\displaystyle+A_{0}^{4}\\{2A_{1}^{\prime}h(h-1)+A_{1}h^{\prime}(3h-1)\\}+A_{1}^{4}f^{2}h^{2}\\{2A_{1}^{\prime}h(h-1)+A_{1}h^{\prime}(3h-1)\\}$
$\displaystyle+2A_{0}^{2}A_{1}^{2}fh\\{A_{1}h^{\prime}(1-3h)-2A_{1}^{\prime}h(h-1)\\}]=0\,,$
(42) $\displaystyle M_{\rm
Pl}^{2}\left[rhf^{\prime}+f(h-1)\right]+\frac{r^{2}hA_{0}^{\prime
2}}{2}-\frac{\eta r^{2}}{2}\left(A_{0}^{2}+fhA_{1}^{2}\right)+\frac{4\alpha
hA_{1}}{(A_{0}^{2}-fhA_{1}^{2})^{2}}[2A_{0}^{\prime}A_{0}^{3}f(1-3h)+2A_{0}^{2}A_{1}^{2}f^{\prime}fh(1-3h)$
$\displaystyle+2A_{0}^{\prime}A_{0}A_{1}^{2}f^{2}h(h-1)+A_{0}^{4}f^{\prime}(3h-1)+A_{1}^{4}f^{\prime}f^{2}h^{2}(3h-1)]=0\,.$
(43)

Variations of the action (38) with respect to $A_{0}$ and $A_{1}$ lead,
respectively, to

$\displaystyle
A_{0}^{\prime\prime}+\left(\frac{2}{r}-\frac{f^{\prime}}{2f}+\frac{h^{\prime}}{2h}\right)A_{0}^{\prime}-\frac{\eta}{h}A_{0}-\frac{4\alpha
A_{0}}{r^{2}h(A_{0}^{2}-fhA_{1}^{2})^{2}}[A_{0}^{2}\\{A_{1}f^{\prime}h(h-1)+2A_{1}^{\prime}fh(h-1)+A_{1}fh^{\prime}(3h-1)\\}$
$\displaystyle+A_{1}^{2}fh\\{A_{1}f^{\prime}h(h+1)+2A_{1}^{\prime}fh(h+1)+A_{1}fh^{\prime}(1-h)\\}]=0\,,$
(44) $\displaystyle\eta
A_{1}-\frac{4\alpha}{r^{2}f(A_{0}^{2}-fhA_{1}^{2})^{2}}\left[(h-1)(f^{\prime}f^{2}h^{2}A_{1}^{4}-2f^{\prime}fhA_{0}^{2}A_{1}^{2}+f^{\prime}A_{0}^{4}-2fA_{0}^{\prime}A_{0}^{3})-2f^{2}h(h+1)A_{1}^{2}A_{0}^{\prime}A_{0}\right]=0\,.$
(45)

In the absence of the vector-GB coupling ($\alpha=0$) with $\eta\neq 0$, we
have $A_{1}(r)=0$ from Eq. (45). For the asymptotically flat boundary
conditions where $f$ and $h$ approach $1$ at spatial infinity, the large-
distance solution to Eq. (44) is given by
$A_{0}=C_{1}e^{-\sqrt{\eta}r}/r+C_{2}e^{\sqrt{\eta}r}/r$. To avoid the
divergence of $A_{0}$ at spatial infinity, we have to choose $C_{2}=0$ and
hence $A_{0}=C_{1}e^{-\sqrt{\eta}r}/r$ at large distances. From Eqs. (42) and
(43), we obtain the relation $(f/h)^{\prime}=\eta A_{0}^{2}r/(M_{\rm
Pl}^{2}h^{2})$. Since $f/h$ should be a finite constant on the BH horizon, we
need to require that $A_{0}=0$. Indeed, the solution consistent with all the
background equations and boundary conditions is $A_{0}(r)=0$ at any radius. In
this case, we end up with the Schwarzschild solution with the metric
components $f=h=1-r_{h}/r$, where $r_{h}$ is the horizon radius.

For $\alpha\neq 0$, Eq. (45) shows that it is possible to realize the solution
with $A_{1}\neq 0$. Moreover, the nonvanishing radial component $A_{1}$
affects the temporal component $A_{0}$ through the $\alpha$-dependent terms in
Eq. (44). For simplicity, we shall seek solutions with $A_{\mu}\neq 0$ for a
small coupling constant $\alpha$ and leave general analysis for future work.
When $\alpha=0$ we only have the trivial solution $A_{\mu}=0$, so the
solutions for a small $\alpha$ may be scaled as $A_{\mu}=\mathcal{O}(\alpha)$.
Let us express the leading-order solutions to $A_{0}$ and $A_{1}$ in the forms

$A_{0}(r)=\alpha\tilde{A}_{0}(r)\,,\qquad A_{1}(r)=\alpha\tilde{A}_{1}(r)\,,$
(46)

where $\tilde{A}_{0}$ and $\tilde{A}_{1}$ are functions of $r$. Then, from
Eqs. (44) and (45), we obtain

$\displaystyle\tilde{A}_{0}^{\prime\prime}+\left(\frac{2}{r}-\frac{f^{\prime}}{2f}+\frac{h^{\prime}}{2h}\right)\tilde{A}_{0}^{\prime}-\frac{\eta}{h}\tilde{A}_{0}-\frac{4\tilde{A}_{0}}{r^{2}h(\tilde{A}_{0}^{2}-fh\tilde{A}_{1}^{2})^{2}}[\tilde{A}_{0}^{2}\\{\tilde{A}_{1}f^{\prime}h(h-1)+2\tilde{A}_{1}^{\prime}fh(h-1)+\tilde{A}_{1}fh^{\prime}(3h-1)\\}$
$\displaystyle+\tilde{A}_{1}^{2}fh\\{\tilde{A}_{1}f^{\prime}h(h+1)+2\tilde{A}_{1}^{\prime}fh(h+1)+\tilde{A}_{1}fh^{\prime}(1-h)\\}]=0\,,$
(47)
$\displaystyle\eta\tilde{A}_{1}-\frac{4}{r^{2}f(\tilde{A}_{0}^{2}-fh\tilde{A}_{1}^{2})^{2}}\left[(h-1)(f^{\prime}f^{2}h^{2}\tilde{A}_{1}^{4}-2f^{\prime}fh\tilde{A}_{0}^{2}\tilde{A}_{1}^{2}+f^{\prime}\tilde{A}_{0}^{4}-2f\tilde{A}_{0}^{\prime}\tilde{A}_{0}^{3})-2f^{2}h(h+1)\tilde{A}_{1}^{2}\tilde{A}_{0}^{\prime}\tilde{A}_{0}\right]=0\,.$
(48)

As we observe in Eqs. (42) and (43), the vector field contributions to metric
components $f$ and $h$ arise at second order in $\alpha$. Then, up to first
order in $\alpha$, we can exploit the Schwarzschild metric components:

$f=h=1-\frac{r_{h}}{r}\,.$ (49)

We substitute Eq. (49) and its derivatives into Eqs. (47) and (48).

Note that we only need to impose two boundary conditions to solve Eqs. (47)
and (48). Although Eq. (47) contains a second-order derivative of
$\tilde{A}_{0}$, one can express $\tilde{A}_{0}^{\prime\prime}$ with respect
to first-order derivatives of $\tilde{A}_{0}$ and $\tilde{A}_{1}$ by
differentiating Eq. (48). Then, one obtains a set of first-order differential
equations of $\tilde{A}_{0}$ and $\tilde{A}_{1}$.

###  III.1 Boundary conditions

Around the horizon, we expand the temporal vector component in the form

$\tilde{A}_{0}=\sum_{i=0}a_{i}\left(\frac{r-r_{h}}{r_{h}}\right)^{i}=a_{0}+a_{1}\frac{r-r_{h}}{r_{h}}+a_{2}\frac{(r-r_{h})^{2}}{r_{h}^{2}}\cdots\,,$
(50)

where $a_{i}$’s are constants. If $\tilde{A}_{0}$ decreases around $r=r_{h}$,
$a_{1}$ is negative. We are interested in regular vector field solutions where
both $X$ and $F$ are finite on the horizon. The leading-order contribution to
$F$ at $r=r_{h}$ is $a_{1}^{2}/(2r_{h}^{2})$. To keep $X$ finite on the
horizon, we require that the leading-order radial vector component diverges as
$\tilde{A}_{1}=\tilde{A}_{0}/\sqrt{fh}=a_{0}r/(r-r_{h})$ at $r=r_{h}$. In this
case, we can expand $\tilde{A}_{1}$ in the form

$\tilde{A}_{1}=a_{0}\frac{r}{r-r_{h}}+\sum_{i=0}b_{i}\left(\frac{r-r_{h}}{r_{h}}\right)^{i}=a_{0}\frac{r}{r-r_{h}}+b_{0}+b_{1}\frac{r-r_{h}}{r_{h}}\cdots\,,$
(51)

where $b_{i}$’s are constants. Then, we have
$X=a_{0}(a_{1}-b_{0})\alpha^{2}+{\cal O}(r-r_{h})$ in the vicinity of
$r=r_{h}$. Substituting Eqs. (50)-(51) into Eqs. (47)-(48), we find that
$b_{0}$ and $b_{1}$ are related to $a_{0}$, $a_{1}$, and $a_{2}$ according to

$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle\frac{2r_{m}^{2}a_{1}+r_{h}^{3}a_{0}a_{1}\pm
2\sqrt{r_{m}^{2}a_{1}[(a_{1}-4a_{0})r_{m}^{2}-r_{h}^{3}a_{0}^{2}]}}{4r_{m}^{2}+r_{h}^{3}a_{0}}\,,$
(52) $\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle[2(a_{1}-b_{0})^{3}a_{1}r_{m}^{2}+2a_{0}^{2}(6a_{1}^{2}-2a_{1}b_{0}+4a_{2}b_{0})r_{m}^{2}+4a_{0}(a_{1}-b_{0})\\{a_{1}^{2}+a_{1}(a_{2}-7b_{0})-2(a_{2}-2b_{0})b_{0}\\}r_{m}^{2}$
(53) $\displaystyle-
a_{0}(a_{1}-b_{0})^{3}(a_{0}+b_{0})r_{h}^{3}]/[4a_{0}a_{1}(2a_{0}-a_{1}+b_{0})r_{m}^{2}]\,,$

where we set

$\eta\equiv\frac{1}{r_{m}^{2}}\,.$ (54)

We would like to derive asymptotic solutions of $\tilde{A}_{0}$ and
$\tilde{A}_{1}$ approaching 0 at spatial infinity. Note that the GB
corrections to (47) and (48) are “zeroth” order in $A_{\mu}$, that is, the
power of $A_{\mu}$ in the denominator and that in the numerator are the same.
Therefore, the limit $A_{\mu}\to 0$ needs to be carefully analyzed. Since the
equation of motion for $\tilde{A}_{0}$ contains a mass term
$\tilde{A}_{0}/r_{m}^{2}$, we search for solutions where $\tilde{A}_{0}$
decreases as $e^{-r/r_{m}}/r$ or faster while $\tilde{A}_{1}$ decreases slower
than $\tilde{A}_{0}$. Ignoring the $\tilde{A}_{0}$-dependent contributions to
Eq. (48), it follows that

$\frac{\tilde{A}_{1}}{r_{m}^{2}}-\frac{4(h-1)f^{\prime}}{r^{2}f}=0\,.$ (55)

Then, the radial vector component should have the asymptotic behavior

$\tilde{A}_{1}=-\frac{4r_{h}^{2}r_{m}^{2}}{r^{5}}\,,$ (56)

whose amplitude decreases in proportion to $r^{-5}$. For the last terms of Eq.
(47) containing the squared bracket, we neglect the $\tilde{A}_{0}$-dependent
contributions and substitute the solution (56) into Eq. (47). Then, at large
distances, we have

$\tilde{A}_{0}^{\prime\prime}+\frac{2}{r}\tilde{A}_{0}^{\prime}-\frac{\tilde{A}_{0}}{r_{m}^{2}}+\frac{8r}{r_{m}^{2}r_{h}}\tilde{A}_{0}-\frac{20r^{2}}{r_{h}^{2}r_{m}^{2}}\tilde{A}_{0}=0\,.$
(57)

Ignoring the third and fourth terms relative to the fifth one in the regime
$r\gg r_{h}$, we obtain the following asymptotic solution

$\tilde{A}_{0}=\frac{C_{1}}{\sqrt{r}}I_{1/4}\left(\frac{\sqrt{5}r^{2}}{r_{h}r_{m}}\right)+\frac{C_{2}}{\sqrt{r}}K_{1/4}\left(\frac{\sqrt{5}r^{2}}{r_{h}r_{m}}\right)\,,$
(58)

where $I_{1/4}$ and $K_{1/4}$ are the Bessel functions of first and second
kinds, respectively, and $C_{1}$, $C_{2}$ are integration constants. The
boundary condition avoiding the divergence of $\tilde{A}_{0}$ corresponds to
$C_{1}=0$ and hence

$\tilde{A}_{0}=\frac{C_{2}}{\sqrt{r}}K_{1/4}\left(\frac{\sqrt{5}r^{2}}{r_{h}r_{m}}\right)\,,$
(59)

which decreases as $\tilde{A}_{0}\propto
r^{-3/2}e^{-\sqrt{5}r^{2}/(r_{h}r_{m})}$. Note that this solution decreases
even faster than $e^{-r/r_{m}}/r$, but the discussion for deriving Eqs. (56)
and (59) does not lose its validity. Therefore, we have obtained a consistent
asymptotic solution (56) and (59) which contains one parameter undetermined by
the asymptotic boundary condition $A_{\mu}\to 0$.

###  III.2 Numerical solutions

In this section, we will numerically study the existence of hairy BH solutions
in theories given by the action (38). For this purpose, we introduce a new
variable $\tilde{B}_{1}$ defined by

$\displaystyle\tilde{B}_{1}\equiv\tilde{A}_{1}-\frac{\tilde{A}_{0}}{f}\,,$
(60)

with $f=1-r_{h}/r$. Around the horizon, using the expanded solutions (50) and
(51) leads to

$\tilde{B}_{1}=b_{0}-a_{1}+(b_{1}-a_{1}-a_{2})\frac{r-r_{h}}{r_{h}}+{\cal
O}\left(\frac{(r-r_{h})^{2}}{r_{h}^{2}}\right)\,.$ (61)

Unlike $\tilde{A}_{1}$, the new variable has a finite value
$\tilde{B}_{1}(r_{h})=b_{0}-a_{1}$ on the horizon. We also have

$\displaystyle
X=-\frac{\alpha^{2}}{2}\tilde{B}_{1}(2\tilde{A}_{0}+f\tilde{B}_{1})\,,$ (62)

which manifests the regularity of $X$ for finite values of $\tilde{A}_{0}$ and
$\tilde{B}_{1}$.

Taking the $r$ derivative of Eq. (48) and using Eq. (47) to eliminate
$\tilde{A}_{0}^{\prime\prime}$, we obtain the first-order differential
equation containing $\tilde{A}_{1}^{\prime}$. We also note that Eq. (48) can
be regarded as the first-order coupled differential equation for
$\tilde{A}_{0}$. Then, the background equations can be expressed in the forms

$\displaystyle\tilde{A}_{0}^{\prime}=F_{0}[\tilde{A}_{0},\tilde{B}_{1}]\,,$
(63)
$\displaystyle\tilde{B}_{1}^{\prime}=F_{1}[\tilde{A}_{0},\tilde{B}_{1}]\,,$
(64)

where the functions $F_{0}$ and $F_{1}$ depend on $\tilde{A}_{0}$ and
$\tilde{B}_{1}$. The right hand sides of Eqs. (63) and (64) are regular on the
horizon.

At large distances ($r\gg r_{h}$), $\tilde{A}_{0}$ decreases much faster than
$\tilde{A}_{1}$ and hence $X\simeq-\alpha^{2}\tilde{B}_{1}^{2}/2<0$, where
$\tilde{B}_{1}\simeq\tilde{A}_{1}=-4r_{h}^{2}r_{m}^{2}/r^{5}$. Since the
background Eqs. (47) and (48) contain terms proportional to $X$ in their
denominators, the regularity of these equations means that $X$ should not
change its sign throughout the horizon exterior, i.e., $X<0$ for $r>r_{h}$. In
the vicinity of $r=r_{h}$, we have

$X=\alpha^{2}\frac{2a_{0}[a_{1}r_{m}^{2}\mp\sqrt{a_{1}^{2}r_{m}^{4}-a_{0}a_{1}r_{m}^{2}(4r_{m}^{2}+a_{0}r_{h}^{3})}]}{4r_{m}^{2}+a_{0}r_{h}^{3}}+{\cal
O}\left(\frac{r-r_{h}}{r_{h}}\right)\,,$ (65)

where the minus and plus signs correspond to the plus and minus signs of
$b_{0}$ in Eq. (52), respectively. If we choose the minus branch of Eq. (65),
we have $X<0$ for $a_{0}>0$ and $a_{1}<0$. For the plus branch of Eq. (65), it
is possible to realize $X<0$ for $a_{0}<0$ and $a_{1}>0$, so long as the
condition $4r_{m}^{2}+a_{0}r_{h}^{3}>0$ is satisfied.

Around the horizon, we have $\tilde{A}_{1}>0$ and $\tilde{A}_{1}^{\prime}<0$
for $a_{0}>0$ from Eq. (51), whereas, for $a_{0}<0$, $\tilde{A}_{1}<0$ and
$\tilde{A}_{1}^{\prime}>0$. At spatial infinity, Eq. (56) gives
$\tilde{A}_{1}<0$ and $\tilde{A}_{1}^{\prime}>0$. Then, the branch smoothly
matching the solutions in two asymptotic regimes ($r\simeq r_{h}$ and $r\gg
r_{h}$) without the sign changes of $\tilde{A}_{1}$ and
$\tilde{A}_{1}^{\prime}$ corresponds to the minus sign of $b_{0}$. We will
search for numerical solutions in this latter regime, i.e., $a_{0}<0$,
$a_{1}>0$, and $4r_{m}^{2}+a_{0}r_{h}^{3}>0$. In this case, we have $b_{0}<0$
from Eq. (52). Around $r=r_{h}$, the variable $\tilde{B}_{1}$ behaves as Eq.
(61), whose leading-order term $b_{0}-a_{1}$ is negative. The first-order
coupled differential Eqs. (63) and (64) can be integrated outward for two
given boundary conditions of $\tilde{A}_{0}(r_{h})=a_{0}$ and
$\tilde{B}_{1}(r_{h})=b_{0}-a_{1}$.

Since $b_{0}$ is expressed by using $a_{0}$ and $a_{1}$ as Eq. (52), choosing
two boundary conditions on the horizon amounts to fixing the two constants
$a_{0}$ and $a_{1}$ in the expansion of $\tilde{A}_{0}$ given by Eq. (50),
with $\tilde{A}_{1}=a_{0}r/(r-r_{h})+b_{0}$. While we need two boundary
conditions on the horizon, there is only one undetermined integration constant
$C_{2}$ at spatial infinity. Two boundary conditions at $r=r_{h}$ may not be
uniquely fixed even if we impose the regularity on the horizon and $A_{\mu}\to
0$ at spatial infinity. This suggests the existence of a family of solutions,
which we will address in the following.

Figure 1:  Numerically derived solutions to $-\tilde{A}_{0}$ (red solid) and
$-\tilde{B}_{1}$ (blue solid) versus $r/r_{h}$ for $r_{m}/r_{h}=10$ with the
boundary conditions (i) $\tilde{A}_{0}(r_{h})=-0.1/r_{h}$,
$\tilde{B}_{1}(r_{h})=-0.17026504/r_{h}$ (left), and (ii)
$\tilde{A}_{0}(r_{h})=-0.2/r_{h}$, $\tilde{B}_{1}(r_{h})=-0.34443001/r_{h}$
(right). Both $-\tilde{A}_{0}$ and $-\tilde{B}_{1}$ are normalized by
$1/r_{h}$. The dashed red and blue curves represent the large-distance
analytic solutions to $-\tilde{A}_{0}$ and $-\tilde{B}_{1}$ determined by Eqs.
(59) and (56), respectively.

Numerically, it is difficult (if possible) to find solutions satisfying
$A_{\mu}\to 0$ at spatial infinity because the asymptotic solution generically
has a nonvanishing growing mode. Instead, we find solutions that smoothly
connect to their large-distance analytic solutions at a finite distance by
using the shooting method. For this purpose, we fix $r_{m}=10r_{h}$ in our
numerical simulations. In this case, unless the boundary conditions on the
horizon are carefully chosen, the growing-mode solution of $\tilde{A}_{0}$,
which corresponds to the first term on the right hand side of Eq. (58), should
manifest itself for the distance $r\gtrsim(10/\sqrt{5})^{1/2}=1.4r_{h}$. Under
appropriate boundary conditions, on the other hand, the solutions can match
the large-distance analytic solutions at least at certain distances, which we
shall regard as our numerical solutions.

In the left panel of Fig. 1, we plot the numerically derived values of
$-\tilde{A}_{0}$ and $-\tilde{B}_{1}$ versus $r/r_{h}$ as solid lines for
$r_{m}/r_{h}=10$ with the boundary conditions
$\tilde{A}_{0}(r_{h})=-0.1/r_{h}$ and
$\tilde{B}_{1}(r_{h})=-0.17026504/r_{h}$. The analytic solutions of
$-\tilde{A}_{0}$ and $-\tilde{B}_{1}$ at large distances are also shown as
dashed lines, which are obtained from Eqs. (59) and (56) with the choice
$C_{2}=-1.6$. We observe that both $-\tilde{A}_{0}$ and $-\tilde{B}_{1}$
approach these large-distance solutions around the distance $r\gtrsim
4.5r_{h}$. We also find that $\tilde{A}_{0}$ starts to deviate from the
decaying mode (59) for $r\gtrsim 6.5r_{h}$, which is followed by the departure
of $\tilde{B}_{1}$ from the large-distance solution
$-4r_{h}^{2}r_{m}^{2}/r^{5}$. This behavior is attributed to the presence of a
growing mode $(C_{1}/\sqrt{r})I_{1/4}(\sqrt{5}r^{2}/(r_{h}r_{m}))$ in
$\tilde{A}_{0}$. Since such a rapidly growing mode is very sensitive to the
accumulation of tiny numerical errors, it is challenging to find the exact
boundary conditions at $r=r_{h}$ realizing $C_{1}=0$ at spatial infinity.

However, as we have mentioned, the fact that there are regions of the distance
in which $\tilde{A}_{0}$ and $\tilde{B}_{1}$ can be well approximated by their
large-distance analytic expressions means that solutions with the proper
asymptotic behavior should exist for the boundary conditions close to
$\tilde{A}_{0}(r_{h})=-0.1/r_{h}$ and
$\tilde{B}_{1}(r_{h})=-0.17026504/r_{h}$. If we fix
$\tilde{A}_{0}(r_{h})=-0.1/r_{h}$ and vary $\tilde{B}_{1}(r_{h})$, we have not
found other ranges of $\tilde{B}_{1}(r_{h})$ in which $-\tilde{A}_{0}$ and
$-\tilde{B}_{1}$ temporally approach their large-distance solutions like the
left panel of Fig. 1. If we consider other boundary conditions of
$\tilde{A}_{0}(r_{h})$ around $-0.1/r_{h}$, there are solutions which can be
well approximated by the large-distance solutions for some ranges of $r$. The
right panel of Fig. 1 corresponds to such an example, in which case
$\tilde{A}_{0}(r_{h})=-0.2/r_{h}$ and
$\tilde{B}_{1}(r_{h})=-0.34443001/r_{h}$. We also found similar cases for
$\tilde{A}_{0}(r_{h})$ larger than $-0.1/r_{h}$, say
$\tilde{A}_{0}(r_{h})=-0.09/r_{h}$, by choosing the values of
$\tilde{B}_{1}(r_{h})$ properly. These facts show that there are appropriate
solutions of $\tilde{A}_{0}$ and $\tilde{B}_{1}$ connecting two asymptotic
regimes ($r\simeq r_{h}$ and $r\gg r_{h}$) for some ranges of
$\tilde{A}_{0}(r_{h})$ around $-0.1/r_{h}$. If $\tilde{A}_{0}(r_{h})$ is far
away from the order $-0.1/r_{h}$, it is typically difficult to find the
parameter spaces of $\tilde{A}_{0}(r_{h})$ and $\tilde{B}_{1}(r_{h})$ in which
both $\tilde{A}_{0}$ and $\tilde{B}_{1}$ approach their asymptotic solutions.

We thus showed that, for $r_{m}=10r_{h}$, there are some ranges of
$\tilde{A}_{0}(r_{h})$ and $\tilde{B}_{1}(r_{h})$ in which the solutions in
two asymptotic regimes can be smoothly connected. This means that the
solutions are not uniquely fixed even for a fixed vector mass term. From Eq.
(51), the radial vector component $\tilde{A}_{1}$ diverges at $r=r_{h}$. In
massless scalar-tensor theories with the linear scalar-GB coupling
$-\alpha\phi{\cal G}$, the field derivative $\phi^{\prime}$ for linearly
stable hairy BH solutions is finite on the horizon and hence
$X_{s}=-(1/2)h\phi^{\prime 2}$ vanishes there [36]. In the latter case, the BH
solution with $X_{s}(r_{h})=0$ is uniquely fixed by performing the expansions
of scalar field and metrics with respect to the small coupling $\alpha$ [24,
25, 36]. It is also possible to consider the boundary condition where
$X_{s}(r_{h})$ is a nonvanishing constant, but in such cases the hairy BHs in
scalar-tensor theories are subject to instabilities of even-parity
perturbations in the vicinity of the horizon [36, 35] (see also Refs. [95, 96,
97, 98] for the general formulation of BH perturbations in Horndeski
theories).

In vector-GB theories discussed above, Eq. (65) shows that $X$ is a
nonvanishing constant at $r=r_{h}$. Unlike scalar-tensor theories, however, we
have to caution that $X$ contains the temporal vector component $A_{0}$
besides the radial component $A_{1}$. Hence the instability argument performed
in Refs. [36, 35] for scalar-tensor theories cannot be applied to vector-
tensor theories. To address this issue, we need to derive linear stability
conditions of even-parity perturbations in the subclass of GP theories. In the
most general class of GP theories the stability of BHs against odd-parity
perturbations was addressed in Ref. [99], but the analysis in the even-parity
sector was not done yet. We also note that we only considered the case
$r_{m}=10r_{h}$ in our numerical simulations, but there should be appropriate
solutions to $A_{\mu}$ for other values of $r_{m}$, too. It is beyond the
scope of this paper to scrutinize all the parameter spaces of boundary
conditions as well as to study linear stabilities of BHs.

###  III.3 Corrections to gravitational potentials

Let us estimate vector field corrections to the metric components $f$ and $h$
both around $r=r_{h}$ and at spatial infinity. In the vicinity of $r=r_{h}$,
we substitute the expanded solutions (50) and (51) into the gravitational Eqs.
(42) and (43). We exploit the leading-order solutions (49) for computing
corrections to $f$ and $h$ of order $\alpha^{2}$. On using the relations (52)
and (53), the differential equations for $h$ and $f$, up to the order of
$\alpha^{2}$, are given by

$\displaystyle-\frac{M_{\rm
Pl}^{2}}{r^{2}}\left(rh^{\prime}+h-1\right)-\mu_{h}\alpha^{2}=0\,,$ (66)
$\displaystyle\frac{M_{\rm
Pl}^{2}}{r^{2}}\left(h-1+\frac{rhf^{\prime}}{f}\right)+\mu_{h}\alpha^{2}=0\,,$
(67)

where $\mu_{h}$ is a constant defined by

$\mu_{h}\equiv\frac{a_{1}^{2}}{2r_{h}^{2}}-\frac{a_{0}a_{1}}{r_{m}^{2}}-\frac{4b_{0}[3a_{0}a_{1}-(2a_{0}+a_{1})b_{0}+b_{0}^{2}]}{(a_{1}-b_{0})^{2}r_{h}^{3}}\,.$
(68)

We have ignored the corrections of order ${\cal O}(r-r_{h})$ for the
derivation of $\mu_{h}$.

Now, we search for solutions of the forms

$\displaystyle f$ $\displaystyle=$
$\displaystyle\left(1-\frac{r_{h}}{r}\right)\left[1+\alpha^{2}F(r)\right]\,,$
(69) $\displaystyle h$ $\displaystyle=$
$\displaystyle\left(1-\frac{r_{h}}{r}\right)\left[1+\alpha^{2}H(r)\right]\,,$
(70)

where $F$ and $H$ are functions of $r$. Substituting Eq. (70) into Eq. (66),
we obtain the integrated solution

$H(r)=\frac{1}{r-r_{h}}\left(-\frac{\mu_{h}r^{3}}{3M_{\rm Pl}^{2}}+{\cal
C}_{1}\right)\,.$ (71)

The integration constant ${\cal C}_{1}$ should be chosen to avoid the
divergence of $H(r)$ at $r=r_{h}$, such that ${\cal
C}_{1}=\mu_{h}r_{h}^{3}/(3M_{\rm Pl}^{2})$. Then, the solution to $h$ up to
the order of $\alpha^{2}$ is given by

$h=\left(1-\frac{r_{h}}{r}\right)\left[1-\alpha^{2}\frac{\mu_{h}(r^{2}+r_{h}r+r_{h}^{2})}{3M_{\rm
Pl}^{2}}\right]\qquad{\rm for}\quad r-r_{h}\ll r_{h}\,.$ (72)

In the limit that $r\to r_{h}$, the $\alpha^{2}$-correction in the square
bracket of Eq. (72) approaches a constant value
$-\alpha^{2}\mu_{h}r_{h}^{2}/M_{\rm Pl}^{2}$. Integrating Eq. (67) after the
substitution of Eq. (72), we obtain

$F(r)=-\frac{\mu_{h}r(r+r_{h})}{3M_{\rm Pl}^{2}}+{\cal C}_{2}\,.$ (73)

Setting ${\cal C}_{2}=0$ by a suitable time reparametrization, it follows that

$f=\left(1-\frac{r_{h}}{r}\right)\left[1-\alpha^{2}\frac{\mu_{h}r(r+r_{h})}{3M_{\rm
Pl}^{2}}\right]\qquad{\rm for}\quad r-r_{h}\ll r_{h}\,,$ (74)

whose $\alpha^{2}$-order correction is finite around $r=r_{h}$. The fact that
the finite corrections to $f$ and $h$ arise at the order of $\alpha^{2}$ is
analogous to the case of linearly coupled scalar-GB theory. We recall however
that $\tilde{A}_{1}$ is divergent on the horizon in vector-GB theory, while
this is not the case for the field derivative $\phi^{\prime}$ in scalar-GB
theory.

At large distances, $\tilde{A}_{0}$ decreases much faster than
$\tilde{A}_{1}$. Then, we can ignore contributions of the
$\tilde{A}_{0}$-dependent terms to the equations of motion of $f$ and $h$. On
using the leading-order solution (49) for the computation of
$\alpha^{2}$-order corrections arising from the vector field, we obtain the
following differential equations

$\displaystyle-\frac{M_{\rm
Pl}^{2}}{r^{2}}\left(rh^{\prime}+h-1\right)+\frac{192r_{h}^{3}r_{m}^{2}}{r^{9}}\alpha^{2}=0\,,$
(75) $\displaystyle\frac{M_{\rm
Pl}^{2}}{r^{2}}\left(h-1+\frac{rhf^{\prime}}{f}\right)-\frac{32r_{h}^{3}r_{m}^{2}}{r^{9}}\alpha^{2}=0\,.$
(76)

We substitute Eqs. (69)-(70) and their $r$ derivatives into Eqs. (75) and (76)
to obtain the differential equations for $H(r)$ and $F(r)$. The integrated
solution to Eq. (75), which is derived by setting the integration constant 0
(whose contribution can be absorbed into $r_{h}$), is given by

$h=\left(1-\frac{r_{h}}{r}\right)\left(1-\alpha^{2}\frac{32r_{h}^{3}r_{m}^{2}}{M_{\rm
Pl}^{2}r^{7}}\right)\qquad{\rm for}\quad r\gg r_{h}\,.$ (77)

Thus, the $\alpha^{2}$-correction term rapidly decreases at large distances.
On using Eq. (77) for Eq. (76) with Eq. (69), we obtain the integrated
solution to $F(r)$. Setting the integration constant 0, the resulting solution
to $f$ is

$f=\left(1-\frac{r_{h}}{r}\right)\left(1-\alpha^{2}\frac{64r_{h}^{3}r_{m}^{2}}{7M_{\rm
Pl}^{2}r^{7}}\right)\qquad{\rm for}\quad r\gg r_{h}\,.$ (78)

The $\alpha^{2}$-order corrections to $f$ and $h$ decrease much faster in
comparison to linearly coupled scalar-GB theory where $F(r)$ and $H(r)$ are
proportional to $1/r$ at large distances [24, 25, 36].

More precisely, when one chooses a specific integration constant to integrate
the scalar field equation in scalar-GB theory, the integrated equations of
motion are the same as Eqs. (42)-(45) with $A_{0}=0$ and replacing $A_{1}$
with $\phi^{\prime}$. In scalar-GB theory, the integration constant $Q$
arising from the scalar field equation is uniquely fixed to a nonzero value
[24, 25] to realize a finite value of $\phi^{\prime}(r)$ on the horizon [36,
35]. In vector-tensor theory, the existence of $A_{0}$ besides $A_{1}$ allows
us to satisfy the field equations around the horizon even with a divergent
value of $A_{1}$. In this case, two boundary conditions of $A_{0}$ and $A_{1}$
at $r=r_{h}$ are not uniquely fixed in general. Since $A_{0}$ and the field
strength $F=-F^{\mu\nu}F_{\mu\nu}/4=hA_{0}^{\prime}{}^{2}/(2f)$ exponentially
decrease in the regime $r\gg r_{h}$, the large-distance solution is dominated
by the longitudinal mode $A_{1}$ proportional to $r^{-5}$. In scalar-tensor
theory, the field derivative has the large-distance behavior
$\phi^{\prime}(r)\propto r^{-2}$ and hence its radial dependence is different
from that of $A_{1}$. In the vicinity of the horizon, the BH solution in
vector-GB theory is particularly different from that in scalar-tensor theory
due to the interplay of both temporal and longitudinal vector components.

###  III.4 BH solutions with $\eta=0$

Finally, we consider the massless vector field case, i.e.,

$\eta=0\,.$ (79)

In this case, we can solve Eq. (45) for the radial vector component $A_{1}$ as

$A_{1}^{2}=\frac{A_{0}[(A_{0}f^{\prime}+A_{0}^{\prime}f)h+fA_{0}^{\prime}-A_{0}f^{\prime}\pm\sqrt{A_{0}^{\prime}f\\{A_{0}^{\prime}f(h+1)^{2}+4A_{0}f^{\prime}h(h-1)\\}}]}{f^{\prime}fh(h-1)}\,.$
(80)

Around the BH horizon characterized by the radius $r_{h}$, we expand $f$, $h$,
and $A_{0}$ in the forms

$f=\sum_{i=1}f_{i}\left(\frac{r-r_{h}}{r_{h}}\right)^{i}\,,\qquad
h=\sum_{i=1}h_{i}\left(\frac{r-r_{h}}{r_{h}}\right)^{i}\,,\qquad
A_{0}=\sum_{i=0}a_{i}\left(\frac{r-r_{h}}{r_{h}}\right)^{i}\,.$ (81)

Substituting Eq. (81) with Eq. (80) into the background Eqs. (42)-(44), the
coefficients $f_{i}$, $h_{i}$, and $a_{i}$ can be obtained at each order. On
using the relation (80), we have

$X=\frac{a_{0}[a_{1}\pm\sqrt{a_{1}(a_{1}-4a_{0}h_{1})}]}{2f_{1}}+{\cal
O}\left(\frac{r-r_{h}}{r_{h}}\right)\,,$ (82)

which is finite on the horizon. At spatial infinity, one can show that there
are also solutions respecting the asymptotic flatness.

Even if there were BH solutions connecting the solutions around $r=r_{h}$ and
$r\gg r_{h}$, the absence of a mass term $\eta X$ may induce some
instabilities of perturbations on the static and spherically symmetric
background. To address this issue, we study the BH stability against odd-
parity perturbations by using linear stability conditions derived in Ref.
[99]. For the action (38) with $\eta=0$, we compute the quantity $q_{2}$
defined in Eq. (3.24) of Ref. [99]. Exploiting Eq. (80) togehter with the
expansions (81), it follows that

$q_{2}=-\frac{16h_{1}a_{0}^{2}\alpha^{2}}{r_{h}(M_{\rm
Pl}^{2}r_{h}f_{1}+8a_{0}\sqrt{f_{1}h_{1}}\alpha)[\sqrt{a_{1}(a_{1}-4a_{0}h_{1})}\pm
a_{1}]^{2}(r-r_{h})^{2}}+{\cal O}\left(\frac{r_{h}}{r-r_{h}}\right)\,.$ (83)

For small $\alpha$, the leading-order term of $q_{2}$ is given by

$q_{2}=-\frac{16a_{0}^{2}\alpha^{2}}{r_{h}^{2}M_{\rm
Pl}^{2}[\sqrt{a_{1}(a_{1}-4a_{0}h_{1})}\pm a_{1}]^{2}(r-r_{h})^{2}}\,,$ (84)

where we used the fact that $f_{1}$ is equivalent to $h_{1}$ in the limit that
$\alpha\to 0$. The absence of ghosts for vector field perturbations requires
that $q_{2}>0$, but we have $q_{2}<0$ from Eq. (84). Unless $\alpha$ is
strictly 0, there is ghost instability for small values of $\alpha$. This
problem arises only for the theories with $\eta=0$.111The existence of a
problem in the theory with $\eta=0$ can be also seen in the absence of a
kinetic term of the longitudinal mode in the decoupling limit $g_{\rm v}\to
0$. The longitudinal mode, which is a part of even-parity perturbations in the
context of BH perturbations, would be pathological at least in the asymptotic
region.

##  IV Conclusions

We constructed a class of vector-tensor theories in which a vector field
$A_{\mu}$ is coupled to the vector ${\cal J}^{\mu}[A,g]$ whose divergence
corresponds to the GB term ${\cal G}$. We showed that the interacting
Lagrangian $\alpha A_{\mu}{\cal J}^{\mu}$ is equivalent to a subclass of GP
theories with the quintic coupling function $G_{5}=4\alpha\ln|X|$, where
$X=-(1/2)A_{\mu}A^{\mu}$. This is analogous to the fact that a linear scalar-
GB coupling of the form $-\alpha\phi\mathcal{G}$ falls into a subclass of
Horndeski theories with the coupling function $G_{5s}=4\alpha\ln|X_{s}|$,
where $X_{s}=-(1/2)\nabla_{\mu}\phi\nabla^{\mu}\phi$. We also extended the
analysis to a more general Lagrangian $f(X)A_{\mu}{\cal J}^{\mu}$ and found
that this belongs to a subclass of beyond GP theories given by the action
(29). Since beyond GP theories are the healthy extension of GP theories
without modifying the dynamical DOFs, our vector-GB theories are free from
Ostrogradski-type instabilities.

Even though the Lagrangian $\alpha A_{\mu}{\cal J}^{\mu}$ has correspondence
with the scalar-GB coupling $-\alpha\phi\mathcal{G}$, the fact that the vector
field $A_{\mu}$ has a longitudinal component besides transverse components
generally gives rise to spacetime dynamics different from those in scalar-
tensor theories. We applied the vector-GB coupling $\alpha A_{\mu}{\cal
J}^{\mu}$ to the search for static and spherically symmetric BHs with vector
hairs by incorporating the Einstein-Hilbert term, Maxwell scalar, and vector
mass term $\eta X$ with $\eta>0$. Under an expansion of the small coupling
$\alpha$, we derived solutions to the temporal and radial vector components
both around the horizon ($r=r_{h}$) and at spatial infinity ($r\gg r_{h}$).
Unless the boundary conditions on the horizon are chosen very accurately, it
is difficult to numerically integrate the vector field differential equations
up to sufficiently large distances due to the existence of a rapidly growing
mode. Nevertheless, we confirmed the existence of regular BH solutions
approaching the large-distance solution for some ranges of $r$.

We also computed corrections to the Schwarzschild metric arising from the
vector field coupled to the GB term. At large distances, these corrections
rapidly decrease in comparison to linear scalar-GB theory. On the horizon the
radial vector component diverges for the coordinate (39), but this is just a
coordinate singularity. In fact, scalar quantities such as $A_{\mu}A^{\mu}$
and $F_{\mu\nu}F^{\mu\nu}$ and the backreaction to the spacetime metric remain
finite around $r=r_{h}$. In linear scalar-GB theory, the field kinetic term
$X_{s}$ needs to vanish on the horizon to avoid linear instability of even-
parity perturbations. In this case, the perturbatively derived BH solution is
uniquely determined for a given small coupling $\alpha$. In our vector-GB
theory, we generally have a nonvanishing value of $X$ on the horizon. For
small couplings $\alpha$, we numerically found that there are some ranges of
boundary conditions of the vector field in which analytic solutions in two
asymptotic regimes are joined each other. This suggests that, unlike linear
scalar-GB theory, hairy BH solutions in vector-GB theory for given $\alpha$
are not unique.

The fact that $A^{\mu}$ cannot be expressed in terms of a scalar gradient due
to the nonvanishing field strength $F=hA_{0}^{\prime}{}^{2}/(2f)$
differentiates our BH solution from that in scalar-GB theories. At large
distances, $A_{0}$ decays rapidly relative to the longitudinal mode $A_{1}$.
The latter asymptotic solution is given by $A_{1}\propto r^{-5}$, whose radial
dependence is different from a large-distance solution of the field derivative
($\phi^{\prime}\propto r^{-2}$) in scalar-GB theory. In the vicinity of the
horizon, the contribution of $A_{0}$ to the vector-field equation of motion is
as important as that of $A_{1}$. Thus, the boundary conditions on the horizon
are different from those in scalar-tensor theories. This means that taking the
scalar limit $A^{\mu}\to\nabla^{\mu}\phi$ for our BH solution does not recover
the hairy BH in scalar-GB theory.

It will be of interest to search for the parameter space of boundary
conditions in more detail for broader ranges of the vector field mass. We
would also like to formulate BH perturbations in GP theories especially for
the even-parity sector to study the linear stability of hairy BHs. As a
byproduct, it will be possible to compute the quasi-normal modes of BHs which
can be probed by the gravitational wave observations. These issues are left
for future works.

## Acknowledgements

The work of K.A. was supported in part by Grants-in-Aid from the Scientific
Research Fund of the Japan Society for the Promotion of Science, No. 20K14468.
S.T. was supported by the Grant-in-Aid for Scientific Research Fund of the
JSPS Nos. 19K03854 and 22K03642.

## References

  * Will [2014] C. M. Will,  _Living Rev. Rel._ 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Hoyle  _et al._ [2001] C. D. Hoyle, U. Schmidt, B. R. Heckel, E. G. Adelberger, J. H. Gundlach, D. J. Kapner, and H. E. Swanson,  _Phys. Rev. Lett._ 86, 1418 (2001), arXiv:hep-ph/0011014 .
  * Adelberger  _et al._ [2003] E. G. Adelberger, B. R. Heckel, and A. E. Nelson,  _Ann. Rev. Nucl. Part. Sci._ 53, 77 (2003), arXiv:hep-ph/0307284 .
  * Copeland  _et al._ [2006] E. J. Copeland, M. Sami, and S. Tsujikawa,  _Int. J. Mod. Phys. D_ 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk,  _Phys. Rept._ 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. Lett._ 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Berti  _et al._ [2015] E. Berti _et al._ ,  _Class. Quant. Grav._ 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Barack  _et al._ [2019] L. Barack _et al._ ,  _Class. Quant. Grav._ 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Berti  _et al._ [2018] E. Berti, K. Yagi, and N. Yunes,  _Gen. Rel. Grav._ 50, 46 (2018), arXiv:1801.03208 [gr-qc] .
  * Ostrogradsky [1850] M. Ostrogradsky,  _Mem. Acad. St. Petersbourg_ 6, 385 (1850).
  * Woodard [2015] R. P. Woodard,  _Scholarpedia_ 10, 32243 (2015), arXiv:1506.02210 [hep-th] .
  * Lanczos [1938] C. Lanczos,  _Annals Math._ 39, 842 (1938).
  * Lovelock [1971] D. Lovelock,  _J. Math. Phys._ 12, 498 (1971).
  * Stelle [1978] K. S. Stelle,  _Gen. Rel. Grav._ 9, 353 (1978).
  * Zwiebach [1985] B. Zwiebach,  _Phys. Lett. B_ 156, 315 (1985).
  * Gross and Sloan [1987] D. J. Gross and J. H. Sloan,  _Nucl. Phys. B_ 291, 41 (1987).
  * Metsaev and Tseytlin [1987] R. R. Metsaev and A. A. Tseytlin,  _Nucl. Phys. B_ 293, 385 (1987).
  * Kanti  _et al._ [1996] P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley,  _Phys. Rev. D_ 54, 5049 (1996), arXiv:hep-th/9511071 .
  * Torii  _et al._ [1997] T. Torii, H. Yajima, and K.-i. Maeda,  _Phys. Rev. D_ 55, 739 (1997), arXiv:gr-qc/9606034 .
  * Chen  _et al._ [2007] C.-M. Chen, D. V. Gal’tsov, and D. G. Orlov,  _Phys. Rev. D_ 75, 084030 (2007), arXiv:hep-th/0701004 .
  * Guo  _et al._ [2008] Z.-K. Guo, N. Ohta, and T. Torii,  _Prog. Theor. Phys._ 120, 581 (2008), arXiv:0806.2481 [gr-qc] .
  * Pani and Cardoso [2009] P. Pani and V. Cardoso,  _Phys. Rev. D_ 79, 084031 (2009), arXiv:0902.1569 [gr-qc] .
  * Kleihaus  _et al._ [2011] B. Kleihaus, J. Kunz, and E. Radu,  _Phys. Rev. Lett._ 106, 151104 (2011), arXiv:1101.2868 [gr-qc] .
  * Sotiriou and Zhou [2014a] T. P. Sotiriou and S.-Y. Zhou,  _Phys. Rev. Lett._ 112, 251102 (2014a), arXiv:1312.3622 [gr-qc] .
  * Sotiriou and Zhou [2014b] T. P. Sotiriou and S.-Y. Zhou,  _Phys. Rev. D_ 90, 124063 (2014b), arXiv:1408.1698 [gr-qc] .
  * Ayzenberg and Yunes [2014] D. Ayzenberg and N. Yunes,  _Phys. Rev. D_ 90, 044066 (2014), [Erratum: Phys.Rev.D 91, 069905 (2015)], arXiv:1405.2133 [gr-qc] .
  * Maselli  _et al._ [2015] A. Maselli, P. Pani, L. Gualtieri, and V. Ferrari,  _Phys. Rev. D_ 92, 083014 (2015), arXiv:1507.00680 [gr-qc] .
  * Kleihaus  _et al._ [2016] B. Kleihaus, J. Kunz, S. Mojica, and E. Radu,  _Phys. Rev. D_ 93, 044047 (2016), arXiv:1511.05513 [gr-qc] .
  * Doneva and Yazadjiev [2018] D. D. Doneva and S. S. Yazadjiev,  _Phys. Rev. Lett._ 120, 131103 (2018), arXiv:1711.01187 [gr-qc] .
  * Silva  _et al._ [2018] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti,  _Phys. Rev. Lett._ 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Antoniou  _et al._ [2018] G. Antoniou, A. Bakopoulos, and P. Kanti,  _Phys. Rev. Lett._ 120, 131102 (2018), arXiv:1711.03390 [hep-th] .
  * Minamitsuji and Ikeda [2019] M. Minamitsuji and T. Ikeda,  _Phys. Rev. D_ 99, 044017 (2019), arXiv:1812.03551 [gr-qc] .
  * Silva  _et al._ [2019] H. O. Silva, C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti,  _Phys. Rev. D_ 99, 064011 (2019), arXiv:1812.05590 [gr-qc] .
  * Langlois  _et al._ [2022] D. Langlois, K. Noui, and H. Roussille,  _JCAP_ 09, 019 (2022), arXiv:2204.04107 [gr-qc] .
  * Minamitsuji  _et al._ [2022a] M. Minamitsuji, K. Takahashi, and S. Tsujikawa,  _Phys. Rev. D_ 106, 044003 (2022a), arXiv:2204.13837 [gr-qc] .
  * Minamitsuji  _et al._ [2022b] M. Minamitsuji, K. Takahashi, and S. Tsujikawa,  _Phys. Rev. D_ 105, 104001 (2022b), arXiv:2201.09687 [gr-qc] .
  * Minamitsuji and Tsujikawa [2022] M. Minamitsuji and S. Tsujikawa,  _Phys. Rev. D_ 106, 064008 (2022), arXiv:2207.04461 [gr-qc] .
  * Antoniadis  _et al._ [1994] I. Antoniadis, J. Rizos, and K. Tamvakis,  _Nucl. Phys. B_ 415, 497 (1994), arXiv:hep-th/9305025 .
  * Gasperini  _et al._ [1997] M. Gasperini, M. Maggiore, and G. Veneziano,  _Nucl. Phys. B_ 494, 315 (1997), arXiv:hep-th/9611039 .
  * Kawai  _et al._ [1998] S. Kawai, M.-a. Sakagami, and J. Soda,  _Phys. Lett. B_ 437, 284 (1998), arXiv:gr-qc/9802033 .
  * Cartier  _et al._ [2000] C. Cartier, E. J. Copeland, and R. Madden,  _JHEP_ 01, 035 (2000), arXiv:hep-th/9910169 .
  * Cartier  _et al._ [2001] C. Cartier, J.-c. Hwang, and E. J. Copeland,  _Phys. Rev. D_ 64, 103504 (2001), arXiv:astro-ph/0106197 .
  * Tsujikawa  _et al._ [2002] S. Tsujikawa, R. Brandenberger, and F. Finelli,  _Phys. Rev. D_ 66, 083513 (2002), arXiv:hep-th/0207228 .
  * Toporensky and Tsujikawa [2002] A. Toporensky and S. Tsujikawa,  _Phys. Rev. D_ 65, 123509 (2002), arXiv:gr-qc/0202067 .
  * Amendola  _et al._ [2006] L. Amendola, C. Charmousis, and S. C. Davis,  _JCAP_ 12, 020 (2006), arXiv:hep-th/0506137 .
  * Nojiri  _et al._ [2005] S. Nojiri, S. D. Odintsov, and M. Sasaki,  _Phys. Rev. D_ 71, 123509 (2005), arXiv:hep-th/0504052 .
  * Calcagni  _et al._ [2005] G. Calcagni, S. Tsujikawa, and M. Sami,  _Class. Quant. Grav._ 22, 3977 (2005), arXiv:hep-th/0505193 .
  * Calcagni  _et al._ [2006] G. Calcagni, B. de Carlos, and A. De Felice,  _Nucl. Phys. B_ 752, 404 (2006), arXiv:hep-th/0604201 .
  * Koivisto and Mota [2007a] T. Koivisto and D. F. Mota,  _Phys. Lett. B_ 644, 104 (2007a), arXiv:astro-ph/0606078 .
  * Koivisto and Mota [2007b] T. Koivisto and D. F. Mota,  _Phys. Rev. D_ 75, 023518 (2007b), arXiv:hep-th/0609155 .
  * Tsujikawa and Sami [2007] S. Tsujikawa and M. Sami,  _JCAP_ 01, 006 (2007), arXiv:hep-th/0608178 .
  * Guo  _et al._ [2007] Z.-K. Guo, N. Ohta, and S. Tsujikawa,  _Phys. Rev. D_ 75, 023520 (2007), arXiv:hep-th/0610336 .
  * Amendola  _et al._ [2007] L. Amendola, C. Charmousis, and S. C. Davis,  _JCAP_ 10, 004 (2007), arXiv:0704.0175 [astro-ph] .
  * Satoh and Soda [2008] M. Satoh and J. Soda,  _JCAP_ 09, 019 (2008), arXiv:0806.4594 [astro-ph] .
  * Guo and Schwarz [2009] Z.-K. Guo and D. J. Schwarz,  _Phys. Rev. D_ 80, 063523 (2009), arXiv:0907.0427 [hep-th] .
  * Kanti  _et al._ [2015] P. Kanti, R. Gannouji, and N. Dadhich,  _Phys. Rev. D_ 92, 041302 (2015), arXiv:1503.01579 [hep-th] .
  * Hikmawan  _et al._ [2016] G. Hikmawan, J. Soda, A. Suroso, and F. P. Zen,  _Phys. Rev. D_ 93, 068301 (2016), arXiv:1512.00222 [hep-th] .
  * Kawai and Kim [2019] S. Kawai and J. Kim,  _Phys. Lett. B_ 789, 145 (2019), arXiv:1702.07689 [hep-th] .
  * Yi  _et al._ [2018] Z. Yi, Y. Gong, and M. Sabir,  _Phys. Rev. D_ 98, 083521 (2018), arXiv:1804.09116 [gr-qc] .
  * Barton  _et al._ [2021] S. Barton, B. Hartmann, B. Kleihaus, and J. Kunz,  _Phys. Lett. B_ 817, 136336 (2021), arXiv:2103.01651 [gr-qc] .
  * Kawai and Kim [2021a] S. Kawai and J. Kim,  _Phys. Rev. D_ 104, 083545 (2021a), arXiv:2108.01340 [astro-ph.CO] .
  * Zhang [2022] F. Zhang,  _Phys. Rev. D_ 105, 063539 (2022), arXiv:2112.10516 [gr-qc] .
  * Kawai and Kim [2021b] S. Kawai and J. Kim,  _Phys. Rev. D_ 104, 043525 (2021b), arXiv:2105.04386 [hep-ph] .
  * Kawaguchi and Tsujikawa [2023] R. Kawaguchi and S. Tsujikawa,  _Phys. Rev. D_ 107, 063508 (2023), arXiv:2211.13364 [astro-ph.CO] .
  * Tsujikawa [2023] S. Tsujikawa,  _Phys. Lett. B_ 838, 137751 (2023), arXiv:2212.10022 [gr-qc] .
  * Carroll  _et al._ [2005] S. M. Carroll, A. De Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner,  _Phys. Rev. D_ 71, 063513 (2005), arXiv:astro-ph/0410031 .
  * Chiba [2005] T. Chiba,  _JCAP_ 03, 008 (2005), arXiv:gr-qc/0502070 .
  * De Felice  _et al._ [2006] A. De Felice, M. Hindmarsh, and M. Trodden,  _JCAP_ 08, 005 (2006), arXiv:astro-ph/0604154 .
  * De Felice  _et al._ [2010a] A. De Felice, D. F. Mota, and S. Tsujikawa,  _Phys. Rev. D_ 81, 023532 (2010a), arXiv:0911.1811 [gr-qc] .
  * De Felice  _et al._ [2010b] A. De Felice, J.-M. Gerard, and T. Suyama,  _Phys. Rev. D_ 82, 063526 (2010b), arXiv:1005.1958 [astro-ph.CO] .
  * De Felice and Suyama [2009a] A. De Felice and T. Suyama,  _Phys. Rev. D_ 80, 083523 (2009a), arXiv:0907.5378 [astro-ph.CO] .
  * De Felice and Suyama [2009b] A. De Felice and T. Suyama,  _JCAP_ 06, 034 (2009b), arXiv:0904.2092 [astro-ph.CO] .
  * Cheng  _et al._ [2006] H.-C. Cheng, M. A. Luty, S. Mukohyama, and J. Thaler,  _JHEP_ 05, 076 (2006), arXiv:hep-th/0603010 .
  * Mukohyama [2007] S. Mukohyama,  _JHEP_ 05, 048 (2007), arXiv:hep-th/0610254 .
  * Aoki  _et al._ [2022] K. Aoki, M. A. Gorji, S. Mukohyama, and K. Takahashi,  _JCAP_ 01, 059 (2022), arXiv:2111.08119 [hep-th] .
  * Colleaux [2019] A. Colleaux,  _Regular black hole and cosmological spacetimes in Non-Polynomial Gravity theories_ , Ph.D. thesis, Trento U. (2019).
  * Horndeski [1974] G. W. Horndeski,  _Int. J. Theor. Phys._ 10, 363 (1974).
  * Kobayashi  _et al._ [2011] T. Kobayashi, M. Yamaguchi, and J. Yokoyama,  _Prog. Theor. Phys._ 126, 511 (2011), arXiv:1105.5723 [hep-th] .
  * Heisenberg [2014] L. Heisenberg,  _JCAP_ 05, 015 (2014), arXiv:1402.7026 [hep-th] .
  * Tasinato [2014] G. Tasinato,  _JHEP_ 04, 067 (2014), arXiv:1402.6450 [hep-th] .
  * Allys  _et al._ [2016] E. Allys, P. Peter, and Y. Rodriguez,  _JCAP_ 02, 004 (2016), arXiv:1511.03101 [hep-th] .
  * Beltran Jimenez and Heisenberg [2016] J. Beltran Jimenez and L. Heisenberg,  _Phys. Lett. B_ 757, 405 (2016), arXiv:1602.03410 [hep-th] .
  * Gallego Cadavid and Rodriguez [2019] A. Gallego Cadavid and Y. Rodriguez,  _Phys. Lett. B_ 798, 134958 (2019), arXiv:1905.10664 [hep-th] .
  * Heisenberg  _et al._ [2016] L. Heisenberg, R. Kase, and S. Tsujikawa,  _Phys. Lett. B_ 760, 617 (2016), arXiv:1605.05565 [hep-th] .
  * Kimura  _et al._ [2017] R. Kimura, A. Naruko, and D. Yoshida,  _JCAP_ 01, 002 (2017), arXiv:1608.07066 [gr-qc] .
  * Heisenberg  _et al._ [2017a] L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa,  _Phys. Rev. D_ 96, 084049 (2017a), arXiv:1705.09662 [gr-qc] .
  * Heisenberg  _et al._ [2017b] L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa,  _JCAP_ 08, 024 (2017b), arXiv:1706.05115 [gr-qc] .
  * Padmanabhan and Kothawala [2013] T. Padmanabhan and D. Kothawala,  _Phys. Rept._ 531, 115 (2013), arXiv:1302.2151 [gr-qc] .
  * Fernandes  _et al._ [2022] P. G. S. Fernandes, P. Carrilho, T. Clifton, and D. J. Mulryne,  _Class. Quant. Grav._ 39, 063001 (2022), arXiv:2202.13908 [gr-qc] .
  * De Felice  _et al._ [2016a] A. De Felice, L. Heisenberg, R. Kase, S. Mukohyama, S. Tsujikawa, and Y.-l. Zhang,  _JCAP_ 06, 048 (2016a), arXiv:1603.05806 [gr-qc] .
  * De Felice  _et al._ [2016b] A. De Felice, L. Heisenberg, R. Kase, S. Mukohyama, S. Tsujikawa, and Y.-l. Zhang,  _Phys. Rev. D_ 94, 044024 (2016b), arXiv:1605.05066 [gr-qc] .
  * De Felice  _et al._ [2017] A. De Felice, L. Heisenberg, and S. Tsujikawa,  _Phys. Rev. D_ 95, 123540 (2017), arXiv:1703.09573 [astro-ph.CO] .
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa,  _Phys. Rev. Lett._ 105, 111301 (2010), arXiv:1007.2700 [astro-ph.CO] .
  * De Felice and Tsujikawa [2012] A. De Felice and S. Tsujikawa,  _JCAP_ 02, 007 (2012), arXiv:1110.3878 [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] .
  * Kase and Tsujikawa [2022] R. Kase and S. Tsujikawa,  _Phys. Rev. D_ 105, 024059 (2022), arXiv:2110.12728 [gr-qc] .
  * Kase and Tsujikawa [2023] R. Kase and S. Tsujikawa,  _Phys. Rev. D_ 107, 104045 (2023), arXiv:2301.10362 [gr-qc] .
  * Kase  _et al._ [2018] R. Kase, M. Minamitsuji, S. Tsujikawa, and Y.-L. Zhang,  _JCAP_ 02, 048 (2018), arXiv:1801.01787 [gr-qc] .