Spontaneously scalarized black holes in dynamical Chern-Simons gravity: dynamics and equilibrium solutions
SSpontaneously scalarized black holes in dynamical Chern-Simons gravity: dynamicsand equilibrium solutions
Daniela D. Doneva ∗ and Stoytcho S. Yazadjiev
1, 2, 3, † Theoretical Astrophysics, Eberhard Karls University of T ¨ubingen, T ¨ubingen 72076, Germany Department of Theoretical Physics, Faculty of Physics, Sofia University, Sofia 1164, Bulgaria Institute of Mathematics and Informatics, Bulgarian Academy of Sciences,Acad. G. Bonchev St. 8, Sofia 1113, Bulgaria
In the present paper, we construct spontaneously scalarized rotating black hole solutionsin dynamical Chern-Simons (dCS) gravity by following the scalar field evolution in the de-coupling limit. For the range of parameters where the Kerr black hole becomes unstablewithin dCS gravity the scalar field grows exponentially until it reaches an equilibrium con-figuration that is independent of the initial perturbation. Interestingly, the Z symmetryof the scalar field is broken and a strong maximum around only one of the rotational axescan be observed. The black hole scalar charge is calculated for two coupling functions sug-gesting that the main observations would remain qualitatively correct even if one considerscoupling functions/coupling parameters producing large deviations from the Kerr solutionbeyond the decoupling limit approximation. I. INTRODUCTION
An interesting class of modified gravity theories that reduce exactly to General relativity forspherically symmetric solutions and can deviate in the presence of a parity-odd source such asrotation is the dynamical Chern-Simons (dCS) gravity [1]. This modified theory possesses anadditional dynamical (pseudo) scalar field coupled non-minimally to the Pontryagin topologicalinvariant. Such a coupling arises naturally in loop quantum gravity [2–4] and in effective fieldtheories [5]. For a review on the subject, we refer the reader to [6]. DCS gravity can be viewedas a special case of the extended scalar-tensor theories where the usual Einstein-Hilbert actionis supplemented with all possible algebraic curvature invariants of second order with a dynam-ical scalar field nonminimally coupled to these invariants [7]. What makes the dCS gravity sointeresting is the fact that the theory is parity-violating – the deviations from general relativityoccur only for systems that violate parity through the presence of a preferred axis. Such systemsare for example the isolated rotating black holes. The spinning black holes within dCS gravitywere studied perturbatively in a series of papers [8–14]. The only non-perturbative study of spin-ning black holes in dCS gravity is that in [15] for a scalar field linearly coupled to the Pontryagininvariant, while dCS black holes without Z isometry were considered in [16]. The binary blackhole merger of such black holes was examined in [17].In the last years the phenomenon of spontaneous scalarization of black holes attracted a lot ofinterest [18–22]. Most of the studies of this phenomenon were performed within the Gauss-Bonnet(GB) gravity and very little was done for the investigation of the spontaneous scalarization in thecase of the dCS gravity. An exception is [23] (see also [24]) where the authors have studied the ∗ [email protected] † [email protected]fia.bg a r X i v : . [ g r- q c ] F e b tachyonic instability that triggers the spontaneous scalarization of the Kerr black hole within thequadratic dCS gravity. It was shown that the Kerr black hole becomes unstable under linear scalarperturbations in certain region of the parameter space. This is an indication that the Kerr blackhole in the dCS gravity would get scalarized giving rise to a new non-Kerr black hole solution. Inaddition, as a toy model for rotation, the scalarization of Schwarzschild-NUT spacetime in dCSgravity was examined in [25]. However, no spontaneously scalarized black hole solutions havebeen explicitly found/constructed in the dCS gravity.The purpose of the present paper is to study the very dynamics of the spontaneous scalariza-tion in dCS gravity, i.e. the process of forming scalarized black holes from Kerr black holes withindCS gravity taking into account the non-linearities in the scalar field coupling similar to [26].Since the scalarization dynamics of the rotating black holes in its full generality and nonlinearityis an extremely difficult task, in the present paper we consider the scalarization dynamics in the“decoupling limit” – we numerically evolve the nonlinear scalar field equation on the fixed geom-etry background of a Kerr black hole, i.e. we neglect the back reaction of the scalar field dynamicson the spacetime geometry. This approximation has proven to be very accurate in order cases ofblack hole scalarization, such as in GB gravity, if we keep the scalar charge small enough [26–29].As end states of the dynamics, we also obtain the stationary scalarized black hole solutionsin the dCS gravity in the “decoupling limit” and study their properties. The solution of the fullstationary field equations without approximations is very difficult since we have to deal with ro-tating solutions and field equations containing third-order derivative of the metric functions [15].That is why, despite the adopted approximation, the present studies can give us very valuableinsight into the properties of scalarized black holes and allow us to compare with gravitationaltheories without parity violation, such as GB gravity.The paper is organized as follows. In Section II we present the basic background behind dCSgravity and derive the relevant scalar field evolution equation. Section III is devoted to the ob-tained numerical results starting with a description of the scalar field time evolution followed byan exploration of the properties of the equilibrium black holes. The paper ends with Conclusions. II. DYNAMICAL CHERN-SIMONS GRAVITY AND SCALAR FIELD EVOLUTION EQUATION
The action for the dCS gravity is given by S = S = π (cid:90) d x (cid:112) − g (cid:104) R − ∇ µ ϕ ∇ µ ϕ + λ f ( ϕ ) (cid:63) RR (cid:105) , (1)where R denotes the Ricci scalar with respect to the spacetime metric g µν , ∇ µ is the covariantderivative with respect to the spacetime metric g µν and f ( ϕ ) is the coupling function for the scalarfield ϕ . The coupling constant λ has dimension of length and (cid:63) RR denotes the Pontryagin invari-ant defined by (cid:63) RR = (cid:63) R µναβ R µναβ , where R µναβ is the Riemann tensor and (cid:63) R µναβ = (cid:101) µνδγ R δγαβ is its dual with (cid:101) µνδγ being the 4-dimensional Levi-Civita tensor.The field equations derived from the action are R µν − Rg µν + λ (cid:104) ∇ α f ( ϕ ) (cid:101) αβγ ( µ ∇ γ R βν ) + ∇ α ∇ β f ( ϕ ) (cid:63) R β ( µν ) α (cid:105) = ∇ µ ϕ ∇ ν ϕ − g µν ∇ α ϕ ∇ α ϕ , (2) ∇ α ∇ α ϕ = − λ d f ( ϕ ) d ϕ (cid:63) RR . (3)In the present paper we are interested in asymptotically flat spacetimes and we shall considerthe case with ϕ ∞ =
0. Without loss of generality we can impose on the coupling function f ( ϕ ) the following conditions f ( ) = d fd ϕ ( ) = (cid:101) with (cid:101) = ±
1. In order for the spontaneousscalarization to occur we have to impose one more condition on the coupling function, namely d fd ϕ ( ) =
0. When this condition is fulfilled the Kerr solution (with mass M and angular momen-tum per unit mass a ) ds = − ∆ − a sin θ Σ dt − a sin θ r + a − ∆Σ dtd φ ++ ( r + a ) − ∆ a sin θ Σ sin θ d φ + Σ∆ dr + Σ d θ , (4)where ∆ = r − Mr + a and Σ = r + a cos θ , is also a solutions to the dCS field equationswith a trivial scalar field ϕ =
0. However, for certain range of the parameters M , a and λ theKerr solution suffers from a tachyonic instability – it becomes unstable within the dCS gravity asshown in [23].In the present paper using an approximate approach we show that the exponential growth ofthe scalar field will last until the scalar field becomes large enough so that the nonlinear terms inthe coupling function suppress the instability giving rise to a new stationary scalarized solutionwith a nontrivial scalar hair. As we have already commented the fully nonlinear dynamical prob-lem is extremely difficult and that is why we shall base our study on an approximate approachwhich is free from heavy technical complications but preserves the leading role of the nonlinear-ity associated with the coupling function. In our approach we keep the spacetime geometry fixedand the whole dynamics is governed by the nonlinear equation for the scalar field. This dynam-ical approach is a very good approximation for example in the vicinity of the bifurcation pointwhere the back reaction of the scalar field on the geometry is small, or away from the bifurcationbut for coupling functions leading to relatively weak scalar field [26]. It was successfully appliedin the case of binary black hole merger in GB gravity [28].Following this approach we consider the nonlinear wave equation (3) on the Kerr backgroundgeometry. In explicit form the equation (3) takes the form − (cid:2) ( r + a ) − ∆ a sin θ (cid:3) ∂ t ϕ + ( r + a ) ∂ x ϕ + r ∆ ∂ x ϕ − Mar ∂ t ∂ φ ∗ ϕ + a ( r + a ) ∂ x ∂ φ ∗ ϕ + ∆ (cid:20) θ ∂ θ ( sin θ∂ θ ϕ ) + θ ∂ φ ∗ ϕ (cid:21) (5) = − λ aM ∆Σ r cos θ ( r − a cos θ )( r − a cos θ ) d f ( ϕ ) d ϕ .where we have introduced the tortoise coordinate x and the new azimuthal coordinate φ ∗ definedby dx = r + a ∆ dr , d φ ∗ = d φ + a ∆ dr . (6)We have also taken into account the explicit form of the Pontryagin invariant for the metric (4),namely (cid:63) RR = aM Σ r cos θ ( r − a cos θ )( r − a cos θ ) . (7)We will conclude this section with comments on the conditions for scalarization (see also [23,24]). For this purpose it is useful to consider the linearized version of eq. (5), namely ∇ α ∇ α δϕ = − (cid:101)λ (cid:63) RR δϕ , (8)where δϕ is the scalar field perturbation and (cid:101) = d fd ϕ ( ) . The right-hand side of the equation givesrise to an effective scalar field mass that can be written in the form µ = − (cid:101)λ (cid:63) RR . (9)If µ < µ is controlled by two factors. The first one is the sign of (cid:101) and the second one is thesign of the Pontryagin invariant. Taking the explicit form of (cid:63) RR in axial symmetry (7) it is easyto show that eq. (5) is invariant under the change of the sign of (cid:101) combined with the change ofthe angle θ to π − θ . Thus, the solutions with different (cid:101) are just mirrored one with respect to theother that gives us the freedom to consider only the case of (cid:101) = III. RESULTSA. Scalar field coupling and domain of existence
Even though the adopted method can be very powerful and accurate in determining the fi-nal scalarized rotating black hole equilibrium solutions, the decoupling limit does no allow torigorously prove the existence of solutions. In GB gravity for example the solutions can quicklydisappear after the bifurcation point because the condition for the regularity of the scalar field andthe metric functions at the black hole horizon is violated, and this behavior is strongly dependenton the choice of coupling function and the scalar field potential [30–35]. This violation can not bereproduced in our case if we assume that the background spacetime geometry remains the Kerrone. Luckily, such violation of the regularity conditions is not observed for the non-perturbativedCS black holes with linear coupling [15] (moreover extremal black holes exist in dCS gravity[36]). This gives us the confidence that the domain of existence of the scalarized solutions willprobably span from the bifurcation point all the way to the extremal limit and we will adopt thisassumption in the calculations below, i.e. we will examine the scalarization of both slowly rotatingKerr black holes and near-extremal solutions.In our studies we will consider more conservative and “secure” choices of coupling functionsbased on the studies of scalarized black holes in GB theory. In GB theory the simplest possiblecase that leads to scalarization f ( ϕ ) = ϕ causes instabilities that can be cured either by adding anadditional quartic term in the function f ( ϕ ) [31, 37] or by considering a self-interaction scalar fieldpotential [32]. Still the domain of existence, the stability, etc. is strongly dependent on the weightof these stabilizing terms, i.e. on the values of the associated additional parameters. Much saferand easier to handle numerically in GB gravity are the cases when the coupling has an exponentialform, where the solutions are often stable in the whole domain of existence and almost all of thecases of rotating black holes in GB theory were calculated for such couplings [34, 39, 40] (with theexception of [35]). That is why, similar to the studies in GB theory, we will adopt the followingtwo choices of the coupling function [26, 30] f I ( ϕ ) = β (cid:16) − e − βϕ (cid:17) , (10) f II ( ϕ ) = β (cid:18) − ( βϕ ) (cid:19) , (11)where β is a parameter. For scalarized GB black holes the value of β practically controls the degreeof scalarization – the increase of β suppresses the scalar field and leads to solutions close to theKerr one [30]. This is also fulfilled in dCS gravity and according to our numerical experiments β =
12 leads to relatively weak scalar field where the decoupling limit is a good approximationeven away from the bifurcation point. That is why we will adopt this value in the followingstudies.The domain of existence of scalarized solutions in dCS gravity, that is practically independenton the particular form of the coupling function as long as the condition d fd ϕ ( ) = M as a normal-ization parameter. Thus the family of solutions is described by two parameters (assuming that β is fixed) – the normalized black hole angular momentum a / M and the constant associated withcoupling function λ / M . Thus, for a fixed λ / M there exists a threshold a / M (a bifurcation point)above which the Kerr solution loses stability and gives rise to a new class of scalarized solutions.This threshold a / M decreases with the increase of λ / M , i.e. for larger coupling constant slowerrotating black holes can scalarize. In the limit when λ / M → a / M →
0, i.e. there is no lower limit on a / M for the development of scalarization [23, 24]. Inaddition, the growth time of the scalar field, i.e. the characteristic time required for the scalarfield to develop from a small perturbation, tends to infinity at the bifurcation point and quicklydecreases as the angular momentum or the coupling parameter is increased. Below we will firstdiscuss the dynamics of scalarization followed by an investigation of the scalarized equilibriumblack hole properties. B. Dynamics of the scalarization
The time evolution of the scalar field equations is performed with the numerical code devel-oped in [26] for the case of GB gravity with the necessary adjustments to handle the Pontryaginscalar. The boundary conditions we impose are the standard ones – the scalar field should havethe form of an ingoing wave at the black hole horizon and an outgoing wave at infinity. Not considering of course the region where the hyperbolic character of the field equations is lost [38]. - 5 - 4 - 3 - 2 - 1 (cid:1) H t / M C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:1) / M = 0 . 9 , a / M = 0 . 4 5 - 3 - 2 - 1 (cid:1) H t / M C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:1) / M = 0 . 9 , a / M = 0 . 7 - 2 - 1 (cid:1) H t / M C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:1) / M = 0 . 9 , a / M = 0 . 9 9 9 FIG. 1. Time evolution of the scalar field on the horizon at θ = π /4 for λ / M = a / M = a / M = a / M = β = -20 -15 -10 -5 (cid:2) (cid:1) x / (cid:1) C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:3) / M = 0 . 9 , a / M = 0 . 4 5 t / M = 0 t / M = 1 1 2 t / M = 2 2 4 t / M = 3 3 6 t / M = 4 4 8 -2 0 -15 -10 -5 (cid:2) (cid:1) x / (cid:1) C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:3) / M = 0 . 9 , a / M = 0 . 7 t / M = 0 t / M = 3 7 t / M = 7 4 t / M = 1 1 1 t / M = 1 4 8 -20 -15 -10 -5 (cid:2) (cid:1) x/M C o u p l i n g I I , (cid:2) (cid:1) = 1 2 (cid:3) / M = 0 . 9 , a / M = 0 . 9 9 9 t / M = 0 t / M = 3 7 t / M = 7 4 t / M = 1 1 1 t / M = 1 4 8 FIG. 2. Snapshots of the two-dimensional scalar field profile at different times during the evolution for thesame values of the parameters as in Fig. 1. The last profile in time is always chosen to represent the finalstate after an equilibrium is reached.
The time evolution of the scalar field at the horizon is depicted in Fig. 1 for three differentvalues of the normalized angular momentum, keeping fixed λ / M = λ / M thebifurcation happens at a / M = a / M = a / M = a / M = a / M = a / M = ( r , θ ) domain. Close to the extremal limit ( a / M = C. Properties of the equilibrium black hole solutions
Even though the evolution discussed above can slightly change depending on the initial data,we have explicitly tested that the final equilibrium state is independent. A sample of solutionsfor a variety of parameters is plotted in Fig. 3. The most prominent feature one can notice is thatthey do not have a symmetry with respect to the plane θ = π /2. This behavior is quite differentcompared for most of the example to dCS black holes with linear coupling [15], the GB black holes(with or without scalarization) [34, 41], or the rotating black holes in Einstein-Maxwell-dilatongravity [42] and is connected with the symmetries of the Chern-Simmons term appearing in thescalar field equation (5) if one considers coupling functions of the form (10) and (11).The qualitative behavior of the scalar field for different θ can be understood if one examinesmore closely the criterion for scalarization µ < µ is defined by eq. (9). One shouldkeep in mind that our discussion is specifically about the case of (cid:101) >
0. The Pontryagin invariant -2 0 -1 5 -1 0 -5 (cid:1) (cid:1) x / (cid:1) C o u p l i n g I , (cid:1) = 1 2 (cid:2) / M = 0 . 4 a / M = 0 . 8 -2 0 -1 5 -1 0 -5 (cid:1) (cid:1) x / (cid:1) C o u p l i n g I , (cid:1) = 1 2 (cid:2) / M = 0 . 4 a / M = 0 . 8 5 -2 0 -1 5 -1 0 -5 (cid:1) (cid:1) x / (cid:1) C o u p l i n g I , (cid:1) = 1 2 (cid:2) / M = 0 . 4 a / M = 0 . 9 0 -2 0 -1 5 -1 0 -5 (cid:1) (cid:1) x / (cid:1) C o u p l i n g I , (cid:1) = 1 2 (cid:2) / M = 0 . 4 a / M = 0 . 9 5 -2 0 -1 5 -1 0 -5 (cid:1) (cid:1) x / (cid:1) C o u p l i n g I , (cid:1) = 1 2 (cid:2) / M = 0 . 4 a / M = 0 . 9 9 FIG. 3. Some representative equilibrium scalar field profiles for fixed λ / M = β =
12. Models with several different normalized angular momenta are plotted showing thechange of the scalar field from models close to the bifurcation point all the way to the extremal limit. (cid:2) (cid:1) = 1 2 (cid:1) / M a / M
00 . 0 60 . 0 80 . 10 . 1 20 . 1 4 D / (cid:1) S t a b l e K e r r B H s
S calarized d C S B H s (cid:2) (cid:1) = 1 2 (cid:1) / M a / M
00 . 0 20 . 0 2 80 . 0 3 60 . 0 4 40 . 0 5 2 D / (cid:1) S t a b l e K e r r B H s
S calarized d C S B H s
FIG. 4. Contour plot of the normalized scalar charge D / M as a function of the black hole angular momen-tum a / M and coupling parameter λ / M . The results are for the both coupling functions (10) (left panel)and (11) (right panel) with β = drops very rapidly as we go away from the black hole and that is why we will concentrate onlyon the region in the immediate vicinity to the black hole horizon. Clearly, the condition for thepresence of tachyonic instability µ < θ dependent. The sign of the Potryagin invariant iscontroller by the cos θ term in eq. (7) for low and moderate values of the angular momentum. Inthis case µ < ≤ θ < π /2. Clearly, once the condition for scalarization is fulfilledfor some θ a nontrivial scalar field will develop for the whole black hole, but as the numericalsimulations demonstrate, it is much stronger close to the θ = a / M (more precisely for a / M > ( r − a cos θ ) term in eq. (7) can change sign as well allowing for stronger scalarizationclose to the θ = π rotational axis. Indeed, as we see in the figure, a second maximum of the scalarfield forms at θ = π for very large rotational rates.Another quantity of interest is the scalar charge D . It is defined through the asymptotic of thescalar field at infinity. For the considered coupling functions, the leading order asymptotic has theform ϕ ∼ D / r . This is very different from the case of Chern-Symons gravity with linear couplingwhere the scalar charge is zero [15]. A contour plot of D as a function of the angular momen-tum a / M and the parameter λ / M is given in Fig. 4. We have depicted the cases for both couplingfunctions. Close to the bifurcation point the accuracy of the calculations is somewhat reduced andthus the contour plots are not given in greater detail. The reason is purely computational. Closeto the bifurcation point the growth time of the modes strongly increases and thus a much longerevolution time is required to see the development of accurate scalar field asymptotic. In addition,determining the scalar charge needs very good precision of the solution at large distances and anyundesired reflected signal from infinity will spoil the procedure. Even though we have imposedoutgoing wave boundary conditions, due to numerical inaccuracies there is always such unde-sired reflected signal that eventually spoils the asymptotic behavior of the solution. Our way ofcompletely “eliminating” the reflected signal is the simplest and most straightforward one – wepush the outer boundary to very large values and put an upper limit on the computational time so A more sophisticated method to eliminate the ingoing wave from infinity was implemented in [23]. a / M or λ / M the scalar charge increases monotonically. We havelimited the figure up to λ / M = λ / M the scalar charge quickly increases reachinglarge values for which the backreaction of the scalar field on the spacetime geometry can not besafely neglected anymore. An important property we can observe in Fig. 4 is that the behaviorof the scalar charge is qualitatively the same for both coupling functions and only the range of D / M changes. Moreover, our studies show that the magnitude of D scales with β , i.e. for thesame mass models the increase of β leads to a decrease of D and vise versa. This suggests that theobservations we have made in the present paper are relatively generic, more or less independentof the particular form of the coupling (as long as it leads to scalarization of course), and willremain approximately valid even if we consider the backreaction of the scalar field on the blackhole metric. IV. CONCLUSION
In the present paper, we have studied non-perturbatively the dynamics of black hole scalar-ization in dynamical Chern-Symons gravity neglecting the backreaction of the scalar field on thespacetime metric. This approach has proven to give good results as long as the scalar field is keptsmall enough. Thus, we have limited ourselves to coupling functions and regions of the parame-ter space that fulfill this criterion. Apart from the study of the process of scalar field development,we have paid special attention to the properties of the newly formed equilibrium scalarized dCSblack holes since such solutions were not obtained in the literature until now.The results show that close to the bifurcation point, i.e. for lower angular momentum of theblack hole, the initial exponential increase of the scalar field is more or less monotonic until it set-tles to an equilibrium state. Close to the extremal limit, though, the growth time of the scalar fieldis very large and as a result oscillations with relatively large amplitude around the equilibrium areobserved. Thus the evolution towards this equilibrium state is not a monotonic one. Such behav-ior will have an influence on the metric perturbations and thus the gravitational wave emission,producing an interesting observational signature. In order to study this effect thoroughly, though,one has to take into account also the backreaction of the scalar field on the spacetime metric.The resulting equilibrium solutions at late times are also studied in detail. The scalar fieldprofile is substantially different compared to most of the other solutions in dCS gravity (withlinear coupling) or GB theory (including the case of scalarization). First of all, because of thesymmetries of the Chern-Simons term in the scalar field equation, the Z symmetry of the scalarfield is broken. In addition, for low and moderate rotational rates, the condition for scalarizationis fulfilled only in one of the black hole “hemispheres” leading to a strong maximum of the scalarfield at the corresponding rotational axis, and minimum at the other. Close to the extremal limitthis behavior changes because the condition for the presence of tachyonic instability becomesmore complicated, and thus two maxima (with different magnitude) appear at the two rotationalaxes. This behavior is very different than most of the beyond-Kerr black holes and can havevarious astrophysical implications. For example, the breaking of the symmetry with respect to the1equator will have a clear signature for thick accretion discs, in the quasinormal mode spectrum,etc.We have studied as well the scalar charge for a large domain of the parameters and for twocoupling functions. The results show that the scalar charge gets stronger as we go to larger rota-tional rates and larger coupling parameter λ . We have limited our study to moderate values of λ ,more precisely we focused on λ / M < λ / M . The qualitative behavior of the results is the same for thecoupling function we have considered suggesting that the main observations we have made inthe paper would remain correct even if one considers coupling functions that lead to a strongerscalar field beyond the decoupling limit approximation. ACKNOWLEDGEMENTS
We would like to thank C. Herdeiro and E. Radu for reading the manuscript and the usefulsuggestions. DD acknowledges financial support via an Emmy Noether Research Group fundedby the German Research Foundation (DFG) under grant no. DO 1771/1-1. SY would like to thankthe University of T ¨ubingen for the financial support. The partial support by the Bulgarian NSFGrant KP-06-H28/7 and the Networking support by the COST Actions CA16104 and CA16214are also gratefully acknowledged. [1] R. Jackiw and S. Y. Pi, Phys. Rev. D , 104012 (2003), arXiv:gr-qc/0308071.[2] V. Taveras and N. Yunes, Phys. Rev. D , 064070 (2008), arXiv:0807.2652 [gr-qc].[3] S. Mercuri, Phys. Rev. Lett. , 081302 (2009), arXiv:0902.2764 [gr-qc].[4] A. Ashtekar, A. P. Balachandran, and S. Jo, Int. J. Mod. Phys. A , 1493 (1989).[5] S. Weinberg, Phys. Rev. D , 123541 (2008), arXiv:0804.4291 [hep-th].[6] S. Alexander and N. Yunes, Phys. Rept. , 1 (2009), arXiv:0907.2562 [hep-th].[7] E. Berti et al. , Class. Quant. Grav. , 243001 (2015), arXiv:1501.07274 [gr-qc].[8] N. Yunes and F. Pretorius, Phys. Rev. D , 084043 (2009), arXiv:0902.4669 [gr-qc].[9] K. Konno, T. Matsuyama, and S. Tanda, Prog. Theor. Phys. , 561 (2009), arXiv:0902.4767 [gr-qc].[10] M. Cambiaso and L. F. Urrutia, Phys. Rev. D , 101502 (2010), arXiv:1010.4526 [gr-qc].[11] K. Yagi, N. Yunes, and T. Tanaka, Phys. Rev. D , 044037 (2012), [Erratum: Phys.Rev.D 89, 049902(2014)], arXiv:1206.6130 [gr-qc].[12] L. C. Stein, Phys. Rev. D , 044061 (2014), arXiv:1407.2350 [gr-qc].[13] K. Konno and R. Takahashi, Phys. Rev. D , 064011 (2014), arXiv:1406.0957 [gr-qc].[14] R. McNees, L. C. Stein, and N. Yunes, Class. Quant. Grav. , 235013 (2016), arXiv:1512.05453 [gr-qc].[15] T. Delsate, C. Herdeiro, and E. Radu, Phys. Lett. B , 8 (2018), arXiv:1806.06700 [gr-qc].[16] P. V. P. Cunha, C. A. R. Herdeiro, and E. Radu, Phys. Rev. D , 104060 (2018), arXiv:1808.06692 [gr-qc].[17] M. Okounkova, L. C. Stein, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D , 104026 (2019),arXiv:1906.08789 [gr-qc].[18] I. Z. Stefanov, S. S. Yazadjiev, and M. D. Todorov, Mod. Phys. Lett. A , 2915 (2008), arXiv:0708.4141[gr-qc].[19] D. D. Doneva, S. S. Yazadjiev, K. D. Kokkotas, and I. Z. Stefanov, Phys. Rev. D , 064030 (2010),arXiv:1007.1767 [gr-qc]. [20] D. D. Doneva and S. S. Yazadjiev, Phys. Rev. Lett. , 131103 (2018), arXiv:1711.01187 [gr-qc].[21] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti, Phys. Rev. Lett. , 131104 (2018),arXiv:1711.02080 [gr-qc].[22] G. Antoniou, A. Bakopoulos, and P. Kanti, Phys. Rev. Lett. , 131102 (2018), arXiv:1711.03390 [hep-th].[23] Y.-X. Gao, Y. Huang, and D.-J. Liu, Phys. Rev. D , 044020 (2019), arXiv:1808.01433 [gr-qc].[24] Y. S. Myung and D.-C. Zou, Phys. Lett. B , 136081 (2021), arXiv:2012.02375 [gr-qc].[25] Y. Brihaye, C. Herdeiro, and E. Radu, Phys. Lett. B , 295 (2019), arXiv:1810.09560 [gr-qc].[26] D. D. Doneva and S. S. Yazadjiev, (2021), arXiv:2101.03514 [gr-qc].[27] J. L. Ripley and F. Pretorius, Class. Quant. Grav. , 134001 (2019), arXiv:1903.07543 [gr-qc].[28] H. O. Silva, H. Witek, M. Elley, and N. Yunes, (2020), arXiv:2012.10436 [gr-qc].[29] H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, Phys. Rev. D , 064035 (2019), arXiv:1810.05177[gr-qc].[30] D. D. Doneva, S. Kiorpelidi, P. G. Nedkova, E. Papantonopoulos, and S. S. Yazadjiev, Phys. Rev. D ,104056 (2018), arXiv:1809.00844 [gr-qc].[31] H. O. Silva, C. F. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti, Phys. Rev. D , 064011(2019), arXiv:1812.05590 [gr-qc].[32] C. F. B. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou, Phys. Rev. D , 104041(2019), arXiv:1903.06784 [gr-qc].[33] D. D. Doneva, K. V. Staykov, and S. S. Yazadjiev, Phys. Rev. D , 104045 (2019), arXiv:1903.08119[gr-qc].[34] P. V. Cunha, C. A. Herdeiro, and E. Radu, Phys. Rev. Lett. , 011101 (2019), arXiv:1904.09997 [gr-qc].[35] L. G. Collodel, B. Kleihaus, J. Kunz, and E. Berti, Class. Quant. Grav. , 075018 (2020),arXiv:1912.05382 [gr-qc].[36] B. Chen and L. C. Stein, Phys. Rev. D , 084012 (2018), arXiv:1802.02159 [gr-qc].[37] M. Minamitsuji and T. Ikeda, Phys. Rev. D , 044017 (2019), arXiv:1812.03551 [gr-qc].[38] J. L. Bl´azquez-Salcedo, D. D. Doneva, J. Kunz, and S. S. Yazadjiev, Phys. Rev. D , 084011 (2018),arXiv:1805.05755 [gr-qc].[39] C. A. Herdeiro, E. Radu, H. O. Silva, T. P. Sotiriou, and N. Yunes, (2020), arXiv:2009.03904 [gr-qc].[40] E. Berti, L. G. Collodel, B. Kleihaus, and J. Kunz, (2020), arXiv:2009.03905 [gr-qc].[41] B. Kleihaus, J. Kunz, S. Mojica, and E. Radu, Phys. Rev. D , 044047 (2016), arXiv:1511.05513 [gr-qc].[42] C. A. Herdeiro, E. Radu, N. Sanchis-Gual, and J. A. Font, Phys. Rev. Lett.121