Black holes in quartic-order beyond-generalized Proca theories

The generalized Proca theories with second-order equations of motion can be healthily extended to a more general framework in which the number of propagating degrees of freedom remains unchanged. In the presence of a quartic-order nonminimal coupling to gravity arising in beyond-generalized Proca theories, the speed of gravitational waves $c_t$ on the Friedmann-Lemaitre-Robertson-Walker (FLRW) cosmological background can be equal to that of light $c$ under a certain condition. By using this condition alone, we show that the speed of gravitational waves in the vicinity of static and spherically symmetric black holes is also equivalent to $c$ for the propagation of odd-parity perturbations along both radial and angular directions. As a by-product, the black holes arising in our beyond-generalized Proca theories are plagued by neither ghost nor Laplacian instabilities against odd-parity perturbations. We show the existence of both exact and numerical black hole solutions endowed with vector hairs induced by the quartic-order coupling.


I. INTRODUCTION
The constantly accumulating observational evidence of dark energy and dark matter implies the existence of additional degrees of freedom (DOFs) beyond those appearing in standard model of physics or General Relativity (GR) [1]. One of the candidates for such extra DOFs is a spin-0 scalar field φ. If the scalar field is nonminimally coupled to gravity, Horndeski theories [2] are the most general scalar-tensor theories with second-order equations of motion [3]. It is also possible to perform a healthy extension of Horndeski theories without increasing the propagating DOFs (one scalar and two tensor polarizations) [4][5][6].
The other candidate for extra DOFs is a spin-1 vector field A µ . A massless vector field respects the U (1) gauge symmetry in Minkowski spacetime, but the gauge invariance is explicitly broken by introducing a vector-field mass or by considering derivative and nonminimal couplings. Most general U (1)-broken vector-tensor theories with secondorder equations of motion are known as generalized Proca (GP) theories [7][8][9][10], which contain five propagating DOFs (one longitudinal scalar, two transverse vectors, and two tensor polarizations). If we apply GP theories to cosmology, there exists an interesting de Sitter attractor responsible for the late-time cosmic acceleration [11]. The dark energy models in the framework of GP theories are observationally distinguished from the cosmological constant due to different cosmic expansion and growth histories [12,13]. One can extend GP theories to the domain of beyondgeneralized Proca (BGP) theories [14][15][16] in which the propagating DOFs remain five.
The recent gravitational-wave (GW) event GW170817 [17] from a neutron star merger, together with the gammaray burst GRB 170817A [18], showed that the speed of gravitational waves c t traveling over a cosmological distance (the redshift z < 0.009) is very close to that of light c with the difference less than the order of 10 −15 . If we demand that c t is strictly equivalent to c, neither quartic-order nor quintic-order nonminimal derivative couplings appearing in Horndeski and GP theories are allowed [19][20][21][22] (see also Refs. [23,24]). In scalar-tensor theories beyond Horndeski, it is possible to realize c t = c on the FLRW cosmological background even in the presence of quartic-order nonminimal derivative couplings [25,26]. This is also the case for quartic-order BGP theories [14].
After the detection of GWs from a black hole (BH) merger [27], we are now entering an era in which the physics of BHs can be probed from precise GW measurements in nonlinear regimes of gravity. In theories beyond GR, the existence of extra DOFs can leave imprints on BH solutions as new "hairs". In Horndeski theories, for example, there are several hairy BH solutions on a static and spherically symmetric background for a radial-dependent scalar φ = φ(r) [28][29][30][31][32][33][34][35][36] or a linearly time-dependent scalar φ = qt + ψ(r) [37,38] (see also Ref. [39] and references therein). In the latter configuration there exists a stealth Schwarzschild solution for the quartic coupling G 4 containing a linear term of ∂ µ φ∂ µ φ and the reduced Planck mass squared M 2 pl , in which case the GW speed differs from c on the cosmological background. The quartic-order beyond-Horndeski interaction allows for the realization of a model in which c t is equivalent to c [40].
In GP theories, the existence of a temporal vector component A 0 besides a longitudinal component A 1 gives rise to a wide variety of hairy BH solutions [41][42][43][44][45][46][47][48][49][50]. For example, there is a stealth Schwarzschild solution with A 1 = 0 for the specific quartic coupling G 4 (X) = M 2 pl /2 + X/4, where X = −A µ A µ /2. Recently, it was shown that this BH solution is unstable against odd-parity perturbations in the vicinity of the event horizon [51]. The point is that, under the absence of ghosts, one of the propagation speed squares along the angular direction is negative. This instability problem is intrinsically related to the fact that the speed of GWs around BHs is different from c for quartic couplings G 4 (X). There is also the branch with A 1 = 0, but the model given by the coupling G 4 (X) = M 2 pl /2 + β 4 M 2 pl (X/M 2 pl ) n with n ≥ 1 also leads to the radial and angular propagation speeds whose deviations from c approach nonvanishing constants at spatial infinity [51]. Unless the coupling β 4 is very small, this behavior is at odds with the observed speed of GWs. The extension to BGP theories can give rise to the exact value c t = c, so there is a possibility for overcoming the above mentioned problems.
In this paper, we focus on quartic-order BGP theories and study whether the condition imposed for obtaining the value c t = c on the FLRW cosmological background is sufficient for realizing the same speed of GWs in the vicinity of BHs. In Sec. II, we derive the equations of motion in quartic-order BGP theories on a static and spherically symmetric background. In Sec. III, we obtain the propagation speeds of GW and vector-field perturbation in the vicinity of BHs by considering odd-parity perturbations. We show that the condition for realizing the cosmological value c t = c is sufficient to obtain the same propagation speed around BHs. In Sec. IV, we search for exact and numerical BH solutions with vector hairs in BGP theories satisfying c t = c. As a result, our new hairy BH solutions are affected by neither ghost nor Laplacian instabilities against odd-parity perturbations. In the rest of sections, we choose the natural unit c = 1.

II. QUARTIC-ORDER BEYOND-GENERALIZED PROCA THEORIES
We consider quartic-order BGP theories [14] with the vector field A µ and the field strength where ∇ µ is the covariant derivative operator. The corresponding action is given by where g is the determinant of four-dimensional metric tensor g µν , R is the Ricci scalar, and G 4 is a function of is a new term appearing beyond the domain of second-order GP theories, which is given by where f 4 is a function of X, and E α1α2γ3γ4 is the Levi-Civita tensor satisfying the normalization E α1α2γ3γ4 E α1α2γ3γ4 = −4!. We note that, by taking the scalar limit A µ → ∇ µ φ, the action (2.1) reduces to that of quartic-order shiftsymmetric Horndeski theories and its Gleyzes-Langlois-Piazza-Vernizzi (GLPV) extension [4]. We study BH solutions on a static and spherically symmetric background described by the line element where t, r and (θ, ϕ) represent the time, radial, and angular coordinates, respectively, and f, h are functions of r. The vector-field profile compatible with the background (2.3) is [52] A µ = (A 0 (r), A 1 (r), 0, 0) , where A 0 and A 1 are functions of r. The quantity X is expressed in the form We compute the action (2.1) on the background (2.3) and vary it with respect to f, h, A 0 , A 1 . The resulting equations of motion are given by ) where a prime represents the derivative with respect to r. The coefficients c 1 , · · · , c 7 and d 1 , · · · , d 10 are given in Appendix A.
On the FLRW cosmological background, the propagation speed c t of tensor perturbations was computed in Ref. [14]. For the theories given by the action (2.1), we have The condition for realizing the value c 2 t = 1 translates to where X = 0. In Sec. III, we show that, under the condition (2.11), the propagation speed squared of gravitational waves in the odd-parity sector around the static and spherically symmetric background (2.3) is also equivalent to 1. In Sec. IV, we search for hairy BH solutions by imposing the condition (2.11).

III. ODD-PARITY PERTURBATIONS
We study the stability of BHs against odd-parity perturbations on top of the spacetime metric (2.3) and the vectorfield profile (2.4). We decompose the metric g µν and the vector field A µ into the background and perturbed parts as g µν =ḡ µν + h µν and A µ =Ā µ + δA µ , where a bar represents the background values. The components of metric perturbations h µν in the odd-parity sector are expressed in the forms [51,[53][54][55][56][57]: where a, b represent θ or ϕ, and Q lm , W lm , U lm are functions of t and r. The tensor E ab is defined by where γ is the determinant of two-dimensional metric γ ab on the sphere and ǫ ab is the Levi-Civita symbol with ǫ θϕ = 1, and Y lm is the spherical harmonics. We choose the Regge-Wheller gauge [58,59], in which the perturbation U lm vanishes. The vector perturbation δA lm for odd-parity modes is given by where δA lm is a function of t and r. We expand the action (2.1) up to quadratic order in odd-parity perturbations and then perform the integrals with respect to θ and ϕ. Integrating the action by parts with respect to t and r, and using the background equations of motion (2.6)-(2.9), we obtain the second-order action of odd-parity perturbations in the form where L = l(l + 1), and where a dot represents the derivative with respect to t, and Since the coefficient C 13 is not needed in the following discussion, we do not write its explicit expression here. The coefficient C 6 vanishes in quartic-order BGP theories, but this is not the case in the presence of other interactions [51]. We can derive conditions for the absence of ghosts and Laplacian instabilities by following the procedure given in Ref. [51]. There are two dynamically propagating modes: for l ≥ 2. For the monopole mode (l = 0), the Lagrangian (3.7) vanishes identically. For the dipole mode (l = 1), the perturbation χ becomes non-dynamical and the vector-field perturbation δA 1m is the only propagating DOF. As shown in Ref. [51], the mode δA 1m possesses the propagation speed same as that for δA lm (l ≥ 2) in GP theories. Hence, the perturbation δA lm corresponds to the intrinsic vector mode, and consequently the other mode χ is associated with the tensor perturbation arising from the gravity sector. Introducing χ as a Lagrange multiplier in the action and eliminating W lm and Q lm from S odd by using their perturbation equations of motion, the second-order Lagrangian is expressed in the form where X t = (χ, δA lm ), and K, R, G, M are 2 × 2 matrices. In general, there are other contributions X ′t S X anḋ X t T X to the Lagrangian L odd [51]. The diagonal components of matrices S and T can be absorbed into M after integration by parts. Moreover, the off-diagonal components of S and T vanish by using the coefficients given in Eq. (3.8). Hence the second-order Lagrangian in quartic-order BGP theories can be expressed in the form (3.10) without the contributions X ′t S X and˙ X t T X . The nonvanishing components of the kinetic matrix K are K 11 = q 1 and The sufficient conditions for the absence of ghosts correspond to q 1 > 0 and q 2 > 0. Let us first consider the radial propagation of odd-parity modes by assuming the solution of the form X t ∝ e i(ωt−kr) . In the limit of large ω and k, the dispersion relation reduces to det(ω 2 K − ωkR + k 2 G) = 0. The radial propagation speed c r in proper time is given by c r = ω/( √ f h k) [51]. On using the fact that the nonvanishing components of R and G are given by 3 )/C 1 , we obtain the two propagation speeds from the dispersion relation: The Laplacian instability along the radial direction can be avoided for c 2 r1 ≥ 0 and c 2 r2 ≥ 0. For the modes L ≫ 1, we substitute the solution X t ∝ e i(ωt−lθ) into Eq. (3.10) to derive propagation speeds along the angular direction. Then, the dispersion relation corresponds to det(ω 2 K + M ) = 0. The leading-order diagonal components of the matrix M are M 11 = −LC 1 and M 22 = L(L − 2)D 1 , where The propagation speed squared along the angular direction in proper time is given by c 2 Ω = ω 2 r 2 /(f l 2 ). Taking the limit L → ∞ in the dispersion relation, we obtain the two values: (3.14) We require the two conditions c 2 Ω1 ≥ 0 and c 2 Ω2 ≥ 0 to avoid the Laplacian instability along the angular direction. Since the matrices K, R, and G are diagonal and the matrix M also becomes diagonal in the limit L ≫ 1, the tensor mode χ and the intrinsic vector mode δA lm are orthogonal and decoupled in the high-frequency limit.
We recall that, under the condition (2.11), the cosmological value of c 2 t is equivalent to 1. We compute the quantities In Ref. [60] it was argued that, if the Lagrangian contains cross terms of both the time and spatial derivatives (˙ X t R X ′ in our theory), the positivity of kinetic matrix K is not necessarily required for the Hamiltonian bounded from below. In other words, provided that the cross terms associated with the matrix R do not vanish, the two conditions q 1 > 0 and q 2 > 0 are sufficient but not necessary for the absence of ghosts. In our BGP theory the matrix components of R vanish identically by using Eq. (3.15), so the sufficient conditions for the absence of ghosts translate to q 1 > 0 and q 2 > 0. These quantities yield (3.16) Provided that G 4 > 0, the conditions q 1 > 0 and q 2 > 0 are trivially satisfied outside the horizon. The squares of the radial propagation speeds in Eq. (3.12) are given by On using the fact that D 1 is equivalent to C 12 = −1/(2r 4 ), the propagation speed squares in the angular direction are We have thus shown that, under the condition (2.11), the propagation speeds for odd-parity perturbations on the static and spherically symmetric background are all equivalent to 1. The propagation speeds c r1 and c Ω1 can be identified with those arising from tensor perturbations. Then, under the condition (2.11), the speed of gravitational waves propagating around BHs is the same as the cosmological value c t = 1. The other speeds c r2 and c Ω2 correspond to those arising from vector-field perturbations. For quartic-order BGP theories, the propagation speed squared of vector perturbations on the FLRW cosmological background is given by [14] Under the condition (2.11), it follows that c 2 v = 1. This is consistent with the fact that both c 2 r2 and c 2 Ω2 are equivalent to 1 on the background (2.3). The coincidence of the propagation speed of the vector perturbation with that of the tensor perturbation and their coincidence with the speed of light arises from the specific choice of our theory (2.1) with the condition (2.11). For instance, if the action (2.1) contains nonlinear kinetic terms of the vector field G 2 (X, F, Y ) with F = −F µν F µν /4 and Y = A µ A ν F µρ F ν ρ , the propagation speed of vector perturbations generally differs from the speed of light, while that of tensor perturbations remains the same.
It is also natural to expect that, if the propagation speed of a mode on the cosmological background coincides with the speed of light, that on the BH background should also coincide with the speed of light, since the propagation speed of perturbations is locally fixed on scales much shorter than background curvature radii. Thus, if the propagation speed of the vector mode on the cosmological background c v is equivalent to the speed of light, those on the static and spherically symmetric background, c 2 r2 and c 2 Ω2 , also coincide with the speed of light. For the dipole perturbation (l = 1), only the vector perturbation δA lm propagates with the radial and angular speed squares c 2 r2 and c 2 Ω2 , respectively. They are equivalent to 1 under the condition (2.11). We note that the configuration of a linearly time-dependent scalar φ = qt + ψ(r) in quartic-order shift-symmetric Horndeski theories and its GLPV extension [4] can be recovered by taking the limits δA lm → 0, A 0 → q, and A 1 → ψ ′ , where q is a constant and ψ is a function of r. The fact that the condition (2.11) is sufficient to guarantee the values c 2 r1 = c 2 Ω1 = 1 in BGP theories means that the same result also holds in quartic-order shift-symmetric GLPV theories. Thus, we proved that the claim of Ref. [40] is correct for odd-parity perturbations without putting any restriction on the models.
As we mentioned in Introduction, the charged stealth Schwarzschild solution arising from the specific quartic coupling G 4 (X) = M 2 pl /2+X/4 in GP theories is unstable against odd-parity perturbations in the vicinity of the event horizon [51]. We note that, by the "charged stealth Schwarzschild" solution [41], we distinguish it from the "stealth Schwarzschild" solution obtained in Ref. [37] and its straightforward extension to the GP theory with G 4 = M 2 pl /2+βX and F µν = 0 [43], where β is an arbitrary dimensionless coupling constant. One may wonder if this instability can be alleviated according to the discussion of no-ghost criterion claimed in Ref. [60]. As we will show in Appendix B, this is not the case since the origin of this instability is not the appearance of ghosts but the propagation speed squared being negative. Thus, the conclusion of Ref. [51] was rather obtained from the same criterion as the hyperbolicity condition employed in Ref. [60]. This charged stealth Schwarzschild solution has a nonzero electric field and hence there is no counterpart solution in scalar-tensor theories obtained by the replacement of A µ with ∂ µ φ. Thus, our argument here is peculiar for vector-tensor theories.
Finally, one may concern that the instability of the stealth Schwarzschild solution stemming from the model G 4 (X) = M 2 pl /2 + X/4 in GP theories [51] would contradict with the stability of our model in BGP theories described by the action (2.1), as these two theories may be related to each other via a disformal transformation. As we show in Appendix C, however, the disformal transformation cannot exactly map the former into the latter. After the transformation, there are new interactions of the forms (C3) [15]. Hence the quadratic GP theory after the disformal transformation is not physically equivalent to our BGP theory given by the action (2.1).

IV. HAIRY BH SOLUTIONS
In this section, we derive hairy BH solutions in quartic-order BGP theories. The background equations of motion (2.6)-(2.9) can be expressed in the form  Hence we cannot solve Eq. (4.1) for x to derive closed-form differential equations. This property generally holds in quartic-order BGP theories on the static and spherically symmetric background without imposing the condition (2.11). We note that the determinant also vanishes for the dynamics of anisotropic cosmology in quartic-order BGP theories [61]. Then, the property of vanishing determinant arises for vector-tensor theories with the equations of motion higher than second order under the metric ansatz with maximally-symmetric two-dimensional space. It is an open question whether such behavior generally occurs in the spacetime with the two-dimensional maximally-symmetric space for other gravitational theories beyond second order (e.g., GLPV theories), which we would like to address in a future publication.
The fact that the background equations of motion are not closed means that we need additional conditions to close the system. From Eq. (2.9), there are in general two branches: (a) A 1 = 0, or (b) A 1 = 0.
For the branch (a), Eq. (2.9) is redundant, so the differential equations (4.1) reduce to the system of the 3 × 3 matrix Z with x = t (f ′ , h ′ , A ′′ 0 ). In this case, the determinant of Z is given by which does not generally vanish. Then, we can solve Eq. (4.1) for the variables f, h, A 0 . In Sec. IV A, we will obtain numerical BH solutions for the branch A 1 = 0 by considering quartic-order power-law couplings. For the branch (b), we need to impose at least one condition to close the system (4.1). In Refs. [47], the authors found exact BH solutions in GP theories by imposing the two conditions f = h and X = constant. In Sec. IV B, we will find exact BH solutions in quartic-order BGP theories by imposing the same conditions.
A. Numerical solutions for the branch A1 = 0 In this subsection, we will focus on the branch and numerically obtain hairy BH solutions for power-law couplings where n ≥ 1 is an integer and β 4 is a constant. We also impose the condition (2.11), under which the function f 4 is given by Around the event horizon characterized by the distance r h , we iteratively derive the solutions to Eqs. (2.6)-(2.8) by using the expansions: where f i , h i , a 0 are constants. The coupling β 4 works as corrections to the metric components of the Reissner-Nordström (RN) solution: where µ is a constant in the range 0 < µ < 1. Substituting Eq. (4.7) into Eqs. (2.6)-(2.8) for the branch A 1 = 0, the leading-order coefficients are given by where we have assumed f 1 = h 1 . The result (4.8) holds irrespective of the values of n, but the next-to-leading order coefficients depend on the power n. For n = 1, the nontrivial β 4 dependence appears at the order of O((r − r h ) 2 ), as (4.9) For n = 2, the coupling β 4 appears at the order of O((r − r h ) 3 ), as (4.10) For n > 2, the nontrivial β 4 dependence around the horizon appears at the order of O((r − r h ) n+1 ). Thus, the regularity of f, h, A 0 at the horizon is ensured for general n (≥ 1).
At large distances (r ≫ r h ), the iterative solutions for general n are given by (4.12) (4.13) The coupling β 4 works as corrections to the RN solution with A 0 = P + Q/r. In Fig. 1, we plot the numerically integrated solutions of f, h, A 0 , f − h for n = 2, β 4 = 0.49, and µ = 0.3. We employ the iterative solutions (4.7) up to third order as boundary conditions in the vicinity of the horizon and solve Eqs. (2.6)-(2.8) for the branch A 1 = 0. The two asymptotic solutions in the regimes r ≃ r h and r ≫ r h smoothly connect to each other without any discontinuity. As estimated above, the temporal vector component A 0 is close to 0 around the horizon and then it increases toward the asymptotic value P as r → ∞.
We also numerically confirmed that the curvature invariants such as R, R µν R µν , and R µναβ R µναβ (where R µν is the Ricci tensor and R µναβ is the Riemann tensor) are regular at/outside the horizon and hence there is no curvature singularity. Using the iterative solutions (4.7) for n = 2 and picking up the dominant contributions around the BH event horizon, these quantities reduce to R → [20µ 2 β 4 /{(1 − µ)r 3 h }](r − r h ), R µν R µν → 4µ 2 /r 4 h , and R µναβ R µναβ → 4(5µ 2 − 6µ + 3)/r 4 h , while they converge to 0 at spatial infinity. We note that for general n (≥ 1), R ∼ (r − r h ) n−1 , while R µν R µν and R µναβ R µναβ approach constant as r → r h .
In our numerical simulation, we have shifted the value of f to 1 at the distance r = 10 7 r h by using the freedom of time rescaling. In Fig. 1, we observe that the difference between f and h induced by the coupling β 4 is most significant in the vicinity of the horizon (f − h ≃ 0.1 around r ≃ 3r h ). This difference may be potentially probed in future high-precision GW measurements in nonlinear regimes of gravity.
We have thus shown the existence of hairy BH solutions regular throughout the horizon exterior for n = 2. Numerically, we have also confirmed that the two asymptotic solutions (4.7) and (4.13) are smoothly joined each other for general powers of n (≥ 1). Since there are two independent parameters r h and µ for the near-horizon solutions, the charge P generally depends on M and Q. Hence the Proca hair P is of the secondary type.

B. Exact BH solutions
The exact BH solution found for the specific coupling G 4 (X) = M 2 pl /2 + X/4 in GP theories [41] satisfies the two relations where X c is a constant. In the following, we will search for exact BH solutions in quartic-order BGP theories by imposing the two conditions (4.14). The second condition gives the relation A 2 1 = (A 2 0 − 2f X c )/(f h) between A 1 and A 0 1 . From Eq. (2.9), it follows that so there are three branches satisfying (i) For this branch, the derivative f ′ is given by Substituting this relation into Eq. (2.8), we obtain whose integrated solution is The second condition is satisfied for either (A) G 4 (X c ) = X c /4 = P 2 /8, or (B) Q = 0. In the case (A), Eq. (4.16) reduces to f ′ = [P 2 (1 − f )r 2 − 2Q 2 ]/(P 2 r 3 ), which is integrated to give where M is an integration constant. If we identify the constant P as 2M pl , the metric components in Eq. (4.21) reduce to the RN solution with G 4 (X c ) = M 2 pl /2. The difference from the RN solution in GR is that there is a nonvanishing longitudinal mode A 1 . We require that 2P (M P + Q)r > Q 2 for the existence of the exact solution (4.21). At spatial infinity, the longitudinal mode decreases as A 1 ∝ 1/ √ r. The solution (4.21) exists for the couplings where b n are arbitrary constants. The case (B) corresponds to the special case of (A), i.e., Q = 0 in Eq. (4.21). Namely, this is the stealth Schwarzschild solution f = h = 1 − 2M/r with A 0 = P and A 1 = P 2M/r/f . This solution exists for arbitrary regular functions G 4 (X).

Branch (ii)
We proceed to the second branch characterized by G 4,X (X c ) = 0. In this case, Eq. (2.8) reduces to (4.17), so the solution to A 0 is given by Eq. (4.18). From Eq. (2.6), we obtain under which Eq. (2.7) is also satisfied. This gives the following integrated solution with A 0 = P + Q/r. Provided that X c = P 2 /2, the longitudinal mode A 1 approaches a constant for r → ∞. This behavior is different from the branch (i) in which A 1 decreases toward 0 due to the condition X c = P 2 /2. The exact solution (4.24) can be realized for the couplings If we choose G 4 (X c ) = M 2 pl /2, the metric components f and h in Eq. (4.24) are the same as those of the RN solution.

Branch (iii)
Let us finally discuss exact solutions for the branch (iii) satisfying A 1 = 0. In this case, the two conditions (4.14) give A 0 = √ 2f X c , where we have chosen the branch A 0 > 0. Multiplying Eqs. (2.6) and (2.7) by G 4 (X c ) and G 4 (X c ) − 2X c G 4,X (X c ), respectively, and taking their sums, it follows that Since we are considering the case X c = 0, we obtain which is integrated to give where C 1 and M are constants. From Eq. (2.6), we obtain which also follows from Eq. (2.7). This relation is satisfied for Then, the resulting solution is which corresponds to the extremal RN solution. The above exact solution can be realized by the couplings The solution (4.32) is the special case of Eq. (4.24) with the correspondence under which A 1 vanishes.

V. CONCLUSIONS
The recent event GW170817 showed that the GW speed c t traveling over the cosmological distance is very close to 1. This fact put strong constraints on models of cosmic acceleration in the framework of modified gravity theories. In GP theories with second-order equations of motion, the quartic-and quintic-order interactions are not allowed, unless their coupling constants are very small. In the healthy extension of GP theories (dubbed BGP theories), the additional quartic-order interaction (2.2) gives rise to a model in which the cosmological value of c t is equivalent to 1 under the condition (2.11).
The remaining question is whether the condition (2.11) is sufficient to ensure that the speed of GWs around massive bodies like BHs is equal to 1 as well. To address this point, we considered metric and vector-field perturbations in the odd-parity sector on the static and spherically symmetric background in quartic-order BGP theories. We explicitly showed that, under the condition (2.11), the propagation speeds c r1 and c Ω1 along the radial and angular directions in the gravity sector are both equivalent to 1. Under the same condition, we also found that the speeds of vector-field perturbations in the radial and angular directions reduce to 1. The no-ghost conditions are trivially satisfied for G 4 > 0. Our result about the GW speed around BHs is also valid in quartic-order shift-symmetric Horndeski and GLPV theories with the time-dependent scalar field φ = qt + ψ(r), where r is the radial coordinate, by taking the limits δA lm → 0, A 0 → q, and A 1 → ψ ′ . Hence we proved the claim of Ref. [40] for odd-parity perturbations without restricting models.
We also searched for hairy BH solutions in quartic-order BGP theories by imposing the condition (2.11). In general, the additional interaction beyond the domain of GP theories leads to a vanishing determinant for the equations of motion on the static and spherically symmetric background. This property does not hold under additional conditions, say, by choosing a branch with the vanishing longitudinal component (A 1 = 0) or by imposing the condition f = h.
For the branch A 1 = 0, we analytically derived iterative solutions around the horizon and at spatial infinity for the quartic-order power-law model (4.5) with the BGP interaction (4.6). Numerically, we also confirmed that the solutions in two asymptotic regimes connect to each other without any discontinuity outside the horizon. The coupling β 4 works as corrections to the RN metric. As we see in Fig. 1, the difference between two metric components f and h induced by β 4 is most significant in the vicinity of the horizon.
Imposing the two conditions f = h and X = X c = constant, we also obtained three branches of exact solutions in quartic-order BGP theories satisfying the condition (2.11). The branch (i) corresponds to the RN-type solution (4.21) present for the model (4.22), in which case the longitudinal mode has the dependence A 1 ∝ 1/ √ r at spatial infinity. The branch (ii) arises for the model (4.25) with the RN-type metric given in Eq. (4.24), but A 1 approaches a constant for r → ∞. The branch (iii), which exists for the model (4.33), corresponds to A 1 = 0 with the extremal RN metric given in Eq. (4.32).
In GP theories with the quartic power-law coupling (4.5), the branch A 1 = 0 is unstable against odd-parity perturbations [51]. Moreover, the branch with A 1 = 0 gives rise to the speed of GWs approaching a constant different from 1 at spatial infinity, so this behavior can be odd with the observational bound of c t . In contrast, all the numerical and exact BH solutions derived in this paper satisfy c t = 1 even in the vicinity of BHs, so they are not prone to the instability problem against odd-parity perturbations. Thus, the extension from GP theories to BGP theories allows the possibility for realizing hairy BH solutions in which the behavior of tensor perturbations is similar to that in GR.
In this paper, we focused on perturbations in the odd-parity sector, but it is necessary to study the behavior of even-parity perturbations in order to ensure the stability of BHs in the model with c t = 1. In particular, the existence of scalar perturbations in the even-parity sector may give rise to additional constraints on the model parameters. The numerical solutions with A 1 = 0 and exact solutions with Q = 0 presented in Secs. IV A and IV B do not exist as the counterparts of shift-symmetric Horndeski theories, so it is of interest to investigate the stabilities of them against even-parity perturbations. It is also interesting to place observational constraints on dark energy models in quartic-order BGP theories satisfying the condition (2.11). These issues are left for future works. at r = 2M and picking up the leading-order contribution, we obtain c 2 Ω1 = − M (4M 2 pl + P 2 ) (2M P + Q) 2 (r − 2M ) + O((r − 2M ) 2 ) . (B3) Provided that M > 0, we have c 2 Ω1 < 0 outside the horizon. Thus, the instability of BHs (B2) arises due to the negative sound speed squared.