Coupled vector Gauss-Bonnet theories and hairy black holes

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.

parts. The integral J µ may not be unique since it is defined only through the differential equation ∇ µ J µ = G. One form of J µ , which is expressed in terms of a scalar field and Riemann tensor, is found in Ref. [76]. On using this expression of J µ = J µ [ϕ, g] and the property G = ∇ µ J µ , it is possible to prove that the scalar-GB coupling µ(ϕ)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 J µ by using a vector field A µ , where J µ = J µ [A, g] satisfies the same relation ∇ µ J µ = G. As a candidate for a Lorentz-invariant scalar characterizing the coupling between A µ and the integrated GB term in vector-tensor theories, we propose the Lagrangian A µ J µ . We will show that the interacting Lagrangian A µ J µ 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 µ J µ , where f is an arbitrary function of X = −(1/2)A µ A µ . 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 αA µ J µ (α 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 ηX, where η is a constant. The coupling αA µ J µ is equivalent to a Lagrangian of the quintic-order coupling function G 5 (X) = 4α 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 ∝ X n with n ≥ 1. However, we will show that this is not the case for G 5 (X) = 4α ln |X|. Under an expansion of the small coupling constant α, 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 η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 η ̸ = 0 is different both around the horizon and at large distances. We will also show that BH solutions for η = 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 µ(ϕ)G and the Horndeski Lagrangian. We then introduce a vector field J µ whose divergence ∇ µ J µ is equivalent to the GB term and show that the Lagrangian f (X)A µ J µ belongs to a subclass of beyond GP theories. In Sec. III, we study static and spherically symmetric BH solutions for the coupling αA µ J µ and derive perturbative solutions to A µ with respect to the small coupling α. 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] ν k is the generalized Kronecker delta and R αβ µν is the Riemann tensor. More explicitly, Eq. (2.1) can be expressed as where R is the scalar curvature and R µν 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 G or its associated vector.
A. Scalar field coupled to the GB term Let us first briefly revisit the case in which there is a scalar field ϕ coupled to the GB term of the form µ(ϕ)G. Because of the antisymmetric property of δ µ1···µ k ν1···ν k , the field equations of motion following from the coupling µ(ϕ)G are of second order in the metric tensor g µν and the scalar field ϕ. On using the Riemann tensor and covariant derivatives of ϕ, the GB term can be expressed in the form [34,76] where ∇ δ is a covariant derivative operator and X s = −(1/2)∇ µ ϕ∇ µ ϕ. Substituting this expression of G into the action and integrating (2.4) by parts, it follows that where g is the determinant of g µν , and we use the notations F ,ϕ ≡ ∂F/∂ϕ and F ,Xs ≡ ∂F/∂X s for any arbitrary function F . We will expand the generalized Kronecker delta, integrate the action (2.5) by parts, and exploit the relation [∇ µ , ∇ ν ]∇ α ϕ = R α λµν ∇ λ ϕ to eliminate contractions of the Riemann tensors. Up to boundary terms, the action (2.5) is equivalent to [34] where G µν = R µν − (1/2)g µν R is the Einstein tensor, and The action (2.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 (2.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 µ(ϕ) = −αϕ, where α is a constant, we have G 5s = 4α 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 ϕ → ϕ + c. In the original form of the GB coupling ϕG, the action is quasi-invariant under the shift ϕ → ϕ + c, i.e., invariant up to a total derivative, while the Lagrangian in the Horndeski form is manifestly invariant under the shift. For µ(ϕ) containing nonlinear functions of ϕ, we generally have the ϕ 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.

B. Vector field coupled to the integrated GB term
If we want to construct theories in which a vector field A µ 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 µ ∇ µ 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 µ v (X)G where Again, this interaction may summon a ghostly DOF in the longitudinal sector of A µ 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 µ → g v A µ + ∇ µ ϕ. Then, the decoupling limit g v → 0 gives X → X s .
The interacting Lagrangian µ v (X)G reduces to a coupling between the GB term and the derivative of ϕ, not the scalar field itself, which should yield equations of motion with derivatives higher than second order. As we already explained, the linear coupling −ϕ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 J µ whose divergence agrees with the GB term, ∇ µ J µ = G. After integration by parts, the coupling −ϕG becomes ∇ µ ϕJ µ [ϕ, g]. As shown in Eq. (2.5), the action contains the derivatives of ϕ but not the field itself, for the linear coupling µ ,ϕ = −α. The global shift symmetry can be localized by the replacement ∇ µ ϕ → ∇ µ ϕ+g v A µ with the help of the vector field A µ . The symmetry transformation is now ϕ → ϕ + χ(x), A µ → A µ − ∇ µ χ(x)/g v . The scalar field ϕ can be eliminated by setting the unitary gauge ϕ = 0. All in all, what we need is the replacement ∇ µ ϕ → A µ , where the gauge coupling g v is absorbed into the definition of A µ . 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 ∇ µ ∇ ν ϕ is symmetric in its indices while the replaced quantity ∇ µ A ν is not symmetric. We resolve this ambiguity by imposing the condition ∇ µ J µ = G even after the replacement ∇ µ ϕ → A µ . The expression of J µ consistent with such requirement is given by In the following, we will show that this vector field satisfies the relation ∇ µ J µ = G. In doing so, we use the relation ∇ [µ R αβ νρ] = 0 and the antisymmetric property of the generalized Kronecker delta. We also exploit the following equality whereν k means that this index is omitted. Notice that the 5-dimensional generalized Kronecker delta δ µ1µ2µ3µ4µ5 ν1ν2ν3ν4ν5 vanishes in 4-dimensional spacetime, whose property was used in the second equality of Eq. (2.10). Then, it follows that where 14) where, in the second equality, we used the relation (2.17) For the computation of F 2 , we exploit the following property Then, we have Expanding the covariant derivative ∇ µ in F 3 and using the commutation relation of ∇ µ ∇ ν A β , we obtain Thus, the divergence of J µ is equivalent to the GB term. We consider a scalar quantity A µ J µ as a candidate for the ghost-free Lagrangian of the vector field A µ coupled to the integrated GB term. As a generalization, we also multiply an arbitrary function f (X) of X with the scalar product A µ J µ . The interaction part of the action in such theories is given by which is composed of two parts: Taking the sum of S vGB1 and S vGB2 , terms on the second lines of Eqs. (2.25) and (2.26) cancel each other. On using the property we obtain the reduced action (2.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., we have Note that we dropped a constant term 8α 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. (2.29) vanishes, the action S vGB belongs to a subclass of GP theories [79,80,82]. We recall that, from Eq. (2.7), the linear scalar-GB coupling µ(ϕ)G with µ(ϕ) = −αϕ corresponds to G 5s (X s ) = 4α ln |X s | with G 2s (X s ) = G 3s (X s ) = G 4s (X s ) = 0. The coupling (2.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 (2.31). Varying the action (2.29) with respect to A µ , it follows that where J µ is defined by Eq. (2.9), and The sum of J µ and J µ F can be expressed in a compact form Let us consider the case in which the Maxwell term F ≡ −(1/4)F αβ F αβ and the mass term X are present in addition to the action (2.29) with f (X) = α. The action in such a subclass of GP theories is given by where g v and η are constants and G 5 (X) = 4α ln |X|. Varying this action with respect to A µ , it follows that where we used the relation ∇ µ ∇ ν F µν = 0 and Eq. (2.21).
Taking the decoupling limit g v → 0 with the replacement A µ → g v A µ +∇ µ ϕ, we have J µ F → 0 and hence Eqs. (2.36) and (2.37) reduce to the Maxwell equation, ∇ ν F µν = 0, and the equation of motion for the scalar field, η∇ µ ∇ µ ϕ = αG, respectively. This latter scalar field equation also follows by varying the Lagrangian L = ηX s − αϕG with respect to ϕ, so the linearly coupled scalar-GB theory is recovered by taking the above decoupling limit. In shift-symmetric Horndeski theories, using the expression of G in Eq. (2.3) shows that the scalar field equation can be expressed in the form ∇ µ j µ = 0, where j µ is a conserved current. When this equation is solved for j µ , 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 µ is given by Eq. (2.36), which does not contain an integration constant. Although Eq. (2.37) corresponds to the differential version of Eq. (2.36), one cannot choose an arbitrary integration constant when integrating Eq. (2.37). This property is different from that in shiftsymmetric 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. (2.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 (2.35). The corresponding action belongs to a subclass of GP theories given by where M Pl is the reduced Planck mass and we set g v = 1. Let us consider the static and spherically symmetric background given by the line element where t, r and (θ, φ) 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 η > 0 as in standard massive Proca theory. For the vector field, we consider the following configuration where A 0 and A 1 are functions of r. On the background (3.2), we have where a prime represents the derivative with respect to r. Note that F = −F µν F µν /4 ̸ = 0 implies a nonvanishing temporal vector component; that is, A µ 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 (3.1) with respect to f and h, respectively, we obtain Variations of the action (3.1) with respect to A 0 and A 1 lead, respectively, to In the absence of the vector-GB coupling (α = 0) with η ̸ = 0, we have A 1 (r) = 0 from Eq. (3.8). For the asymptotically flat boundary conditions where f and h approach 1 at spatial infinity, the large-distance solution to Eq. (3.7) is given by A 0 = C 1 e − √ ηr /r + C 2 e √ η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 − √ ηr /r at large distances. From Eqs. (3.5) and (3.6), we obtain the relation (f /h) ′ = ηA 2 0 r/(M 2 Pl 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 α ̸ = 0, Eq. (3.8) shows that it is possible to realize the solution with A 1 ̸ = 0. Moreover, the nonvanishing radial component A 1 affects the temporal component A 0 through the α-dependent terms in Eq. (3.7). For simplicity, we shall seek solutions with A µ ̸ = 0 for a small coupling constant α and leave general analysis for future work. When α = 0 we only have the trivial solution A µ = 0, so the solutions for a small α may be scaled as A µ = O(α). Let us express the leading-order solutions to A 0 and A 1 in the forms A 0 (r) = αÃ 0 (r) , A 1 (r) = αÃ 1 (r) , (3.9) whereÃ 0 andÃ 1 are functions of r. Then, from Eqs. (3.7) and (3.8), we obtaiñ As we observe in Eqs. (3.5) and (3.6), the vector field contributions to metric components f and h arise at second order in α. Then, up to first order in α, we can exploit the Schwarzschild metric components: We substitute Eq. (3.12) and its derivatives into Eqs. (3.10) and (3.11).

A. Boundary conditions
Around the horizon, we expand the temporal vector component in the form where a i 's are constants. IfÃ 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 2 1 /(2r 2 h ). To keep X finite on the horizon, we require that the leading-order radial vector component diverges asÃ 1 =Ã 0 / √ f h = a 0 r/(r − r h ) at r = r h . In this case, we can expandÃ 1 in the form where b i 's are constants. Then, we have X = a 0 (a 1 − b 0 )α 2 + O(r − r h ) in the vicinity of r = r h . Substituting Eqs. (3.13)-(3.14) into Eqs. (3.10)-(3.11), we find that b 0 and b 1 are related to a 0 , a 1 , and a 2 according to where we set (3.17) We would like to derive asymptotic solutions ofÃ 0 andÃ 1 approaching 0 at spatial infinity. Note that the GB corrections to (3.10) and (3.11) are "zeroth" order in A µ , that is, the power of A µ in the denominator and that in the numerator are the same. Therefore, the limit A µ → 0 needs to be carefully analyzed. Since the equation of motion for A 0 contains a mass termÃ 0 /r 2 m , we search for solutions whereÃ 0 decreases as e −r/rm /r or faster whileÃ 1 decreases slower thanÃ 0 . Ignoring theÃ 0 -dependent contributions to Eq. (3.11), it follows that Then, the radial vector component should have the asymptotic behavior whose amplitude decreases in proportion to r −5 . For the last terms of Eq. (3.10) containing the squared bracket, we neglect theÃ 0 -dependent contributions and substitute the solution (3.19) into Eq. (3.10). Then, at large distances, we haveÃ (3.20) Ignoring the third and fourth terms relative to the fifth one in the regime r ≫ r h , we obtain the following asymptotic solutionÃ 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Ã 0 corresponds to C 1 = 0 and hencẽ which decreases asÃ 0 ∝ r −3/2 e − √ 5r 2 /(r h rm) . Note that this solution decreases even faster than e −r/rm /r, but the discussion for deriving Eqs. (3.19) and (3.22) does not lose its validity. Therefore, we have obtained a consistent asymptotic solution (3.19) and (3.22) which contains one parameter undetermined by the asymptotic boundary condition A µ → 0.

B. Numerical solutions
In this section, we will numerically study the existence of hairy BH solutions in theories given by the action (3.1). For this purpose, we introduce a new variableB 1 defined bỹ with f = 1 − r h /r. Around the horizon, using the expanded solutions (3.13) and (3.14) leads tõ UnlikeÃ 1 , the new variable has a finite valueB 1 (r h ) = b 0 − a 1 on the horizon. We also have which manifests the regularity of X for finite values ofÃ 0 andB 1 .
Taking the r derivative of Eq. (3.11) and using Eq. (3.10) to eliminateÃ ′′ 0 , we obtain the first-order differential equation containingÃ ′ 1 . We also note that Eq. (3.11) can be regarded as the first-order coupled differential equation forÃ 0 . Then, the background equations can be expressed in the forms (3.27) where the functions F 0 and F 1 depend onÃ 0 andB 1 . The right hand sides of Eqs. (3.26) and (3.27) are regular on the horizon. At large distances (r ≫ r h ),Ã 0 decreases much faster thanÃ 1 and hence X ≃ −α 2B2 1 /2 < 0, whereB 1 ≃Ã 1 = −4r 2 h r 2 m /r 5 . Since the background Eqs. (3.10) and (3.11) 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 where the minus and plus signs correspond to the plus and minus signs of b 0 in Eq. (3.15), respectively. If we choose the minus branch of Eq. (3.28), we have X < 0 for a 0 > 0 and a 1 < 0. For the plus branch of Eq. (3.28), it is possible to realize X < 0 for a 0 < 0 and a 1 > 0, so long as the condition 4r 2 m + a 0 r 3 h > 0 is satisfied. Around the horizon, we haveÃ 1 > 0 andÃ ′ 1 < 0 for a 0 > 0 from Eq. (3.14), whereas, for a 0 < 0,Ã 1 < 0 and A ′ 1 > 0. At spatial infinity, Eq. (3.19) givesÃ 1 < 0 andÃ ′ 1 > 0. Then, the branch smoothly matching the solutions in two asymptotic regimes (r ≃ r h and r ≫ r h ) without the sign changes ofÃ 1 andÃ ′ 1 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 2 m + a 0 r 3 h > 0. In this case, we have b 0 < 0 from Eq. (3.15). Around r = r h , the variableB 1 behaves as Eq. (3.24), whose leading-order term b 0 − a 1 is negative. The first-order coupled differential Eqs. (3.26) and (3.27) can be integrated outward for two given boundary conditions ofÃ 0 (r h ) = a 0 andB 1 (r h ) = b 0 − a 1 .
Since b 0 is expressed by using a 0 and a 1 as Eq. (3.15), choosing two boundary conditions on the horizon amounts to fixing the two constants a 0 and a 1 in the expansion ofÃ 0 given by Eq. (3.13), withÃ 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 µ → 0 at spatial infinity. This suggests the existence of a family of solutions, which we will address in the following. Numerically, it is difficult (if possible) to find solutions satisfying A µ → 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Ã 0 , which corresponds to the first term on the right hand side of Eq. (3.21), should manifest itself for the distance r ≳ (10/ √ 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 −Ã 0 and −B 1 versus r/r h as solid lines for r m /r h = 10 with the boundary conditionsÃ 0 (r h ) = −0.1/r h andB 1 (r h ) = −0.17026504/r h . The analytic solutions of −Ã 0 and −B 1 at large distances are also shown as dashed lines, which are obtained from Eqs. (3.22) and (3.19) with the choice C 2 = −1.6. We observe that both −Ã 0 and −B 1 approach these large-distance solutions around the distance r ≳ 4.5r h . We also find thatÃ 0 starts to deviate from the decaying mode (3.22) for r ≳ 6.5r h , which is followed by the departure ofB 1 from the large-distance solution −4r 2 h r 2 m /r 5 . This behavior is attributed to the presence of a growing mode (C 1 / √ r)I 1/4 ( √ 5r 2 /(r h r m )) inÃ 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Ã 0 andB 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Ã 0 (r h ) = −0.1/r h andB 1 (r h ) = −0.17026504/r h . If we fixÃ 0 (r h ) = −0.1/r h and varyB 1 (r h ), we have not found other ranges ofB 1 (r h ) in which −Ã 0 and −B 1 temporally approach their large-distance solutions like the left panel of Fig. 1. If we consider other boundary conditions of 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Ã 0 (r h ) = −0.2/r h and B 1 (r h ) = −0.34443001/r h . We also found similar cases forÃ 0 (r h ) larger than −0.1/r h , sayÃ 0 (r h ) = −0.09/r h , by choosing the values ofB 1 (r h ) properly. These facts show that there are appropriate solutions ofÃ 0 andB 1 connecting two asymptotic regimes (r ≃ r h and r ≫ r h ) for some ranges ofÃ 0 (r h ) around −0.1/r h . IfÃ 0 (r h ) is far away from the order −0.1/r h , it is typically difficult to find the parameter spaces ofÃ 0 (r h ) andB 1 (r h ) in which bothÃ 0 and B 1 approach their asymptotic solutions.
We thus showed that, for r m = 10r h , there are some ranges ofÃ 0 (r h ) andB 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. (3.14), the radial vector componentÃ 1 diverges at r = r h . In massless scalar-tensor theories with the linear scalar-GB coupling −αϕG, the field derivative ϕ ′ for linearly stable hairy BH solutions is finite on the horizon and hence X s = −(1/2)hϕ ′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 α [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 [35,36] (see also Refs. [95][96][97][98] for the general formulation of BH perturbations in Horndeski theories).
In vector-GB theories discussed above, Eq. (3.28) 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. [35,36] 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 µ 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.

C. 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 (3.13) and (3.14) into the gravitational Eqs. (3.5) and (3.6). We exploit the leading-order solutions (3.12) for computing corrections to f and h of order α 2 . On using the relations (3.15) and (3.16), the differential equations for h and f , up to the order of α 2 , are given by where µ h is a constant defined by (3.31) We have ignored the corrections of order O(r − r h ) for the derivation of µ h . Now, we search for solutions of the forms (3.33) where F and H are functions of r. Substituting Eq. (3.33) into Eq. (3.29), we obtain the integrated solution The integration constant C 1 should be chosen to avoid the divergence of H(r) at r = r h , such that C 1 = µ h r 3 h /(3M 2 Pl ). Then, the solution to h up to the order of α 2 is given by (3. 35) In the limit that r → r h , the α 2 -correction in the square bracket of Eq.
Setting C 2 = 0 by a suitable time reparametrization, it follows that whose α 2 -order correction is finite around r = r h . The fact that the finite corrections to f and h arise at the order of α 2 is analogous to the case of linearly coupled scalar-GB theory. We recall however thatÃ 1 is divergent on the horizon in vector-GB theory, while this is not the case for the field derivative ϕ ′ in scalar-GB theory. At large distances,Ã 0 decreases much faster thanÃ 1 . Then, we can ignore contributions of theÃ 0 -dependent terms to the equations of motion of f and h. On using the leading-order solution (3.12) for the computation of α 2 -order corrections arising from the vector field, we obtain the following differential equations We substitute Eqs. (3.32)-(3.33) and their r derivatives into Eqs. (3.38) and (3.39) to obtain the differential equations for H(r) and F (r). The integrated solution to Eq. (3.38), which is derived by setting the integration constant 0 (whose contribution can be absorbed into r h ), is given by Thus, the α 2 -correction term rapidly decreases at large distances. On using Eq.
The α 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. (3.5)-(3.8) with A 0 = 0 and replacing A 1 with ϕ ′ . 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 ϕ ′ (r) on the horizon [35,36]. 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 µν F µν /4 = hA ′ 0 2 /(2f ) exponentially decrease in the regime r ≫ 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 ϕ ′ (r) ∝ 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.
D. BH solutions with η = 0 Finally, we consider the massless vector field case, i.e., η = 0 . (3.42) In this case, we can solve Eq. (3.8) for the radial vector component A 1 as .
Around the BH horizon characterized by the radius r h , we expand f , h, and A 0 in the forms , the coefficients f i , h i , and a i can be obtained at each order. On using the relation (3.43), we have 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 ≫ r h , the absence of a mass term η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 (3.1) with η = 0, we compute the quantity q 2 defined in Eq. (3.24) of Ref. [99]. Exploiting Eq. (3.43) togehter with the expansions (3.44), it follows that For small α, the leading-order term of q 2 is given by where we used the fact that f 1 is equivalent to h 1 in the limit that α → 0. The absence of ghosts for vector field perturbations requires that q 2 > 0, but we have q 2 < 0 from Eq. (3.47). Unless α is strictly 0, there is ghost instability for small values of α. This problem arises only for the theories with η = 0. 1

IV. CONCLUSIONS
We constructed a class of vector-tensor theories in which a vector field A µ is coupled to the vector J µ [A, g] whose divergence corresponds to the GB term G. We showed that the interacting Lagrangian αA µ J µ is equivalent to a subclass of GP theories with the quintic coupling function G 5 = 4α ln |X|, where X = −(1/2)A µ A µ . This is analogous to the fact that a linear scalar-GB coupling of the form −αϕG falls into a subclass of Horndeski theories with the coupling function G 5s = 4α ln |X s |, where X s = −(1/2)∇ µ ϕ∇ µ ϕ. We also extended the analysis to a more general Lagrangian f (X)A µ J µ and found that this belongs to a subclass of beyond GP theories given by the action (2.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 αA µ J µ has correspondence with the scalar-GB coupling −αϕG, the fact that the vector field A µ 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 αA µ J µ to the search for static and spherically symmetric BHs with vector hairs by incorporating the Einstein-Hilbert term, Maxwell scalar, and vector mass term ηX with η > 0. Under an expansion of the small coupling α, we derived solutions to the temporal and radial vector components both around the horizon (r = r h ) and at spatial infinity (r ≫ 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 (3.2), but this is just a coordinate singularity. In fact, scalar quantities such as A µ A µ and F µν F µν 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 α. In our vector-GB theory, we generally have a nonvanishing value of X on the horizon. For small couplings α, 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 α are not unique.
The fact that A µ cannot be expressed in terms of a scalar gradient due to the nonvanishing field strength F = hA ′ 0 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 ∝ r −5 , whose radial dependence is different from a large-distance solution of the field derivative (ϕ ′ ∝ 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 µ → ∇ µ ϕ 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.