Black holes and binary mergers in scalar Gauss-Bonnet gravity: scalar field dynamics

We study the nonlinear dynamics of black holes that carry scalar hair and binaries composed of such black holes. The scalar hair is due to a linear or exponential coupling between the scalar and the Gauss--Bonnet invariant. We work perturbatively in the coupling constant of that interaction but nonperturbatively in the fields. We first consider the dynamical formation of hair for isolated black holes of arbitrary spin and determine the final state. This also allows us to compute for the first time the scalar quasinormal modes of rotating black holes in the presence of this coupling. We then study the evolution of nonspinning black-hole binaries with various mass ratios and produce the first scalar waveform for a coalescence. An estimate of the energy loss in scalar radiation and the effect this has on orbital dynamics and the phase of the GWs (entering at quadratic order in the coupling) show that GW detections can set the most stringent constraint to date on theories that exhibit a coupling between a scalar field and the Gauss--Bonnet invariant.

The history of alternative theories of gravity is almost as old as that of general relativity (GR) itself (see, e.g., [1,2]). For more than a century, each astrophysical revolution and the corresponding observational opportunity led to a new milestone test of gravity. Einstein's theory has remained consistent with observations (although dark matter and dark energy may be considered as indication to the contrary), while several modified theories of gravity have been strongly constrained or ruled out [1][2][3].
The recent gravitational-wave (GW) revolution provides yet another opportunity to test gravity in a new regime: the highly-dynamical, strong-curvature regime probed by black holes (BHs), compact objects, and binaries thereof [4][5][6]. While the recent GW events are all consistent with GR [7,8], the constraints one can extract on alternative theories are rather weak [7], due to the lack of complete waveforms that correspond to binary evolution and mergers in these theories. Obtaining such waveforms is necessary to go beyond performing null tests of GR. Using theory-specific waveforms could constrain the corresponding theory to unprecedented levels, or uncover new effects, using data that is already available. Moreover, our ability to probe the highly nonlinear regime of gravity will improve further when LIGO/Virgo detections will perform routinely at design sensitivity, and when future instruments such as third-generation ground-based interferometers [9,10] and the future space mission LISA [11] will become operational.
Motivation for testing gravity are manifold [1][2][3]6], but arguably the most pressing one is of fundamental nature: finding an underlying, consistent description of quantum gravity is still the "holy grail" in modern physics. GR itself fails at this task -e.g., it is non renormalizableand is believed to be the leading-order manifestation of a more fundamental (possibly quantum) theory. Remarkably, it has been shown [12] that including terms quadratic in the curvature in the gravitational action can render the theory renormalizable (such terms also arise in the low-energy limit of string theories [13]).
Uncovering a deviation from the predictions of GR would provide the first experimental insight into this fundamental theory. There is no indisputable argument suggesting that new physics should make an appearance in BH binaries -the curvature involved is not high enough by fundamental physics standards. Nonetheless, it is significantly higher than the curvature scale tested by any other observation and experiment, and it is particularly hard to argue against fully exploring data from a new regime.
As we will see in more detail later, the field equations of sGB gravity, though second order, contain highly nonlinear quadratic terms of those second derivatives. Similar PDE structures, e.g. in hydrodynamics, are known to lead to shocks. Hence, this raises reasonable concerns about the predictivity of the theory in the strongfield regime. Indeed, it has been shown that sGB combined with a generalized harmonic gauge is not wellposed [39,40] although this is a gauge dependent result and, hence, not conclusive. A potential cure was suggested in [41,42] and will be adressed in future work.
These problems only arise if action (2) is taken at face value. In this article we take a different route and treat the coupling between Φ and R GB as the leading-order term of a low-energy expansion and α GB , or more correctly the dimensionless ratio ≡ α GB /l 2 (where l is some reference length) as the control parameter of this expansion. Within this framework, the theory is known to O ( ) only, solutions that are not smooth in the limit α GB → 0 should be considered spurious, and hence one can solve the equations perturbatively in the coupling. This effective-field theory inspired approach has been popular in the literature when obtaining stationary solutions [24,26,29] and recently it has also been used in the context of dynamical evolution [33,34,43].
In what follows we will study the dynamics of isolated BHs and BH binaries for sGB gravity, working perturbatively and up to linear order in the coupling, and for the choices f = Φ/8, e Φ /8. As we will argue below, within the perturbative treatment and at this order, these two choices are actually identical. Moreover, they are also equivalent (modulo rescaling of the coupling α GB ) to any other choice of the function f (Φ), as long as f (0) = 0. 1 We summarize relations between different conventions for the coupling constant that are common in the literature in Appendix A. For instance, the values of |α GB | in the notation of [22] and of this article differ from those of [7,24,25] by a factor of 4 4 √ π.
Dimensional analysis suggests that there are terms other than f (Φ)R GB that could appear at O ( ) in an effective action (such as those including derivatives of the scalar field) which we are implicitly neglecting. However, the fact that ΦR GB is the only shift-symmetric term that leads to BH hair suggests strongly that including O( ) terms that respect this symmetry will lead to secondary corrections only. Indeed, we will demonstrate below that, within our setup, such corrections are O( 2 ) at least. Hence, we do consider our results as rather generic, at least at the qualitative level. When shift symmetry is abandoned though, the scalar can acquire a mass which can lead to effects that our analysis does not capture.
Our results on isolated BHs will provide new insights on the quasi-normal ringing of spinning BHs. Perturbations of spherically symmetric and static BHs in EdGB gravity have been studied in [44,45], where it has been shown that gravitational and scalar perturbations are coupled 2 and the quasinormal-mode (QNM) spectrum contains two branches of modes which have been respectively called gravitational-led and scalar-led, according to which sector they reduce to when the coupling constant vanishes. Toy models obtained with point particles plunging into hairy BHs in EdGB gravity suggest that both gravitational-led and scalar-led QNMs are present in the post-merger GW ringdown signal [44,45].
In the case of rotation, which is so far unexplored at perturbative level, it is reasonable to expect the same qualitative structure for the QNM spectrum. An explicit perturbative computation is extremely challenging due to the lack of (an extended version of) the Teukolsky formalism for BHs in EdGB gravity. Current estimates are based on the assumption that, in the eikonal limit, the light-ring modes are related to the QNMs even in modified gravity theories [44,47], but in fact the coupling of the gravitational perturbations with the scalar field breaks this analogy [44]. Our nonperturbative analysis circumvents the aforementioned issues.
Our numerical simulations of BH binaries give the first complete waveforms for scalar radiation, including the ringdown, within the class of theories we are studying. Moreover, we estimate the impact of the scalar radiation and the corresponding loss of energy on the binary evolution and its imprint on the phase of GWs. We use this estimate to argue that future GW detections can place stronger bounds on the theory than the known ones.
Most of the current constraints on the coupling constant of sGB gravity theories, α GB , have been derived for EdGB gravity and shift-symmetric sGB gravity. In this case, observations on the orbital decay of low-mass X-ray binaries lead to |α GB | 10 km [25,48] 3 . Slightly more stringent constraints could be set from the measurements of quasiperiodic oscillations in the X-ray emission from accreting BHs, although the latter might be affected by large systematics [49,50]. However, in several sGB theories a stronger bound can be obtained by a theoretical constraint. In these theories, a stationary BH solution with mass M only exists if where we have fixed the reference length l to 2M , and M is the Arnowitt-Deser-Misner (ADM) mass which coincides with the BH mass in the case of an isolated BH, and with the binary total mass in the case of BH binaries. The threshold max ∼ O(1) depends on the specific sGB theory [22,32,[35][36][37]. For spherically symmetric BHs, 4 max 0.619 in EdGB gravity [46], 4 max ∼ 0.3 in shift-symmetric sGB gravity [32]. When this bound is reached, a curvature singularity emerges from within the horizon [32], and the solution does not describe a BH anymore. In the case of rotating BHs, the bound becomes stronger (at least, in the case of EdGB gravity, where it has been shown [30] that max decreases as the spin increases, and vanishes at extremality). Therefore, the mere existence of a BH of mass M implies that We remark that the above constraints do not apply to the class of theories found in [35][36][37], which predict the existence of both Kerr BHs and "scalarized" BHs. It should be noted that here we are not requiring the scalar field to have any cosmological significance, so we will not discuss cosmological constraints in any detail. It is also worth emphasising the following: the coincident detection of GWs and gamma rays from the binary neutron star merger GW170817 have constrained the speed of GWs to extremely high accuracy [52,53]. This has in turn been used to place stringent bounds on a class of theories that included sGB gravity [54][55][56][57][58][59] but these bounds rely on the assumption that the scalar field is to account for dark energy [60]. Since we are making no such assumption here, these bounds are inapplicable. Note that in asymptotically flat spacetimes, the speed of GWs approaches unity asymptotically (see e.g. [61]).
In the present paper we report new GW-based constraints on the GB coupling from fully nonlinear simulations covering the inspiral, merger and ringdown. In particular, we compare estimates of the expected GW dephasing against that of current GW detections and the forecast for future third-generation detectors. For a BH binary like GW151226 with mass ratio ∼ 1/2 and total mass ∼ 20M , we find |α GB | 2.7km -comparable to the theoretical constraints and about one order of magnitude stronger than those based on the inspiral only [7,25].

A. Action and field equations
Varying the action (2) with respect to the scalar field Φ and metric g ab yields their field equations where f ≡ df /dΦ. G ab ≡ (4) R ab − 1/2g ab (4) R is the Einstein tensor, the canonical part of energy-momentum tensor for the scalar field is and the modification due to the GB term is [22,46] where * R ab cd = abef (4) R ef cd is the dual Riemann tensor, and abcd is the totally anti-symmetric Levi-Civita pseudo-tensor.
B. Perturbative treatment in the coupling

Preliminaries
Since we want to use a perturbative treatment in the coupling, we will assume that 1 and formally expand any tensor X as In particular, the spacetime metric and the scalar field are expanded as It should be stressed that this is not a weak-field expansion. We raise indices of all tensorial quantities with g (0)ab , e.g. we define and likewise for all other tensors. Since is dimensionless, any tensor X (k) has the same dimensions as the background tensor X (0) -for instance, the scalar field perturbations Φ (k) are dimensionless, as is the background scalar field, and so on.

Field equations
Applying the perturbative treatment to the field equations (4) yields, to O ( ), ab refer to the k-th order correction to the corresponding quantity. The crucial feature of the equations above is that higher-curvature corrections at any given order always enter only as source terms computed from the metric and the scalar field at lower order. Hence, at any given order the system of partial differential equations can be made well-posed by an appropriate gauge choice or reformulation [62,63].

Zero-th order
The zero-th order in the perturbative expansion, equivalent to taking the limit → 0 of Eqs. (4), leads to Einstein's equations minimally coupled to a massless scalar field, Eq. (10a). It has been shown that this system can be cast into a well-posed initial value formulation [62], which is a necessary condition for numerical stability.
Stationary, asymptotically flat BHs cannot carry hair if they satisfy Eqs. (10a) [15]. That is, they would be solutions of vacuum Einstein's equations and any nontrivial initial scalar configuration would be shedded away. One also expects that the scalar field would not be excited in binaries composed of such BHs, as there is no scalar charge to begin with and the equations are linear in the scalar. This suggests that, (at least) at late times, the solution to the zero-th order equation should be of the form where g GR ab is a solution of the vacuum Einstein equations. However, there is a subtlety in this argument. As is evident from Eqs. (10b), Φ (0) effectively sources the first-order (and subsequent order) equations. Hence, a nontrivial initial Φ (0) configuration could in principle leave some imprint on the evolution. Though our expectation is that this effect would be rather small, we have not explored this in any detail. Instead we focus on the late-time behaviour and we enforce Φ (0) = 0. This choice will affect the form of the first-order equations.

First order
Using the solution (11) the first-order field equations (10b) reduce to where f (0) = 1/8 for both EdGB and shift-symmetric sGB gravity, (0) and R GB are, respectively, the d'Alembertian and Gauss-Bonnet invariant evaluated from the background metric g ab . Hence, the metric itself is not deformed and it is safe to set h (1) ab = 0. As indicated in Eq. (12) the scalar field Φ (1) is sourced by the curvature of the background spacetime and, therefore, develops a nontrivial profile. Then, the solution at O( ) is where Φ (1) can be solved for analytically in certain approximations discussed below, or numerically in the general case. Since Eqs. (12) are the Einstein-scalar field equations sourced by tensors computed from (g GR ab , 0), they can be cast into a well-posed initial value formulation.
We remark that under the assumption Φ (0) = 0 and for vanishing ordinary matter the scalar field is O( ) and the Ricci tensor is O( 2 ), hence the GB invariant is equivalent to the Kretschmann scalar up to O( 4 ) terms. Any other term which we have neglected in the action would be of the same order or higher in the perturbation expansion (with the exception of the parity-violating Pontryagin density abcd R ab ef R f ecd leading to Chern-Simons gravity [64][65][66]). Therefore, within our perturbative approach and excluding parity violation, sGB gravity provides the most general theory with higher-order curvature corrections.

Second order
Although we shall solve the field equations up to O ( ) only, in order to assess the validity of our perturbative approach and to estimate backreaction effects on the system's dynamics we inspect the field equations at order O 2 . Since Φ = O( ), its corrections to the metric appear at second order, as we now show. Energy-momentum tensor: When Φ (0) = 0, one gets whereas the first nonvanishing contribution to the scalar stress-energy tensor reads Gauss-Bonnet correction and invariant: Both quantities enter the field equations at a given order (k) only as source terms, i.e., computed from lower-order terms. Up to the linear level considered in (10) and inserting solution (11) we have (neglecting for simplicity the superscript (4) in front of the curvature tensor and its contractions) i.e., the Gauss-Bonnet invariant depends only on the curvature of the background GR spacetime. Instead, the GB correction does contribute to the energy-momentum content of the system at second order (which we will use to estimate the GW dephasing), given by Field equations: The field equations at order O 2 read We remark that the second-order equations are different for EdGB gravity (f (1) = Φ (1) /8) and for shift-symmetric sGB gravity (f (1) = 0). The right-hand side of (18a) defines an effective energy-momentum tensor where T (2) ab and G (1) ab are given in Eqs. (15) and (17), respectively.
We remark that, if one wishes to compute the first nonvanishing corrections to the metric components (including GW emission), it is sufficient to solve the modified Einstein equations (18a) with the linear-order scalar field as an input, whereas the quadratic correction to the scalar field Φ (2) does not affect the metric to leading order.

Summary
In the following, we set f (0) = 1 8 . This choice corresponds to the coupling functions describing, respectively, EdGB gravity or shiftsymmetric sGB gravity. However, we should stress that, within our perturbative scheme, our results are far more general. Indeed, based on the previous discussion, the choice of f (0) completely determines the form of the field equation for the scalar field up to order O( ), cf. Eqns. (12), and also the form of Einstein's equations at O( 2 ), cf. Eq. (18a). Moreover, the precise value of f (0) can be absorbed in α GB . Hence, one does not need to specify f (Φ) any further to fully determine the evolution to O( ) and to estimate how the scalar emission affects the gravitational waveform at O( 2 ), as we will do here. Note however that theories for which f (0) = 0 are not covered by our analysis. It has been found recently that such theories can yield interesting phenomena such as BH scalarization [35][36][37], so binary evolution in these theories deserves further consideration.
Summarizing, the set of field equations (up to O( )) that we evolve numerically is As discussed in Sec. III, we evolve the scalar field simultaneously with the background, i.e. GR spacetime, that is set up either as a single rotating BH or a BH binary. Unless needed for clarity, we will drop the superscripts (0) , (1) in the following, i.e. Φ ≡ Φ (1) and g ab ≡ g

A. Spacetime split
To write the field equations (21) as a time evolution problem we perform a spacetime decomposition. Specifically, we foliate the background spacetime (M, g ab ) into a set of spatial hypersurfaces (Σ t , γ ij ) labelled by a time parameter t and with 3-metric γ ab = g ab + n a n b , where n a denotes the timelike unit vector normal to the hypersurface and is normalized such that n a n a = −1. The spatial metric defines a projection operator with γ a b n b = 0 by construction. The line element takes the form where α and β i are, respectively, the lapse function and shift vector. We denote the covariant derivative and Ricci tensor associated to the 3-metric γ ij as D i and R ij , respectively. The extrinsic curvature is defined as where L n is the Lie derivative along n a .

B. Background spacetimes
We consider two types of background spacetimes: (i) Rotating BHs that will allow us to study the dynamical formation of hairy BHs and their quasi-normal ringdown in the time domain. This also allows us to benchmark our code at late times against the analytic solutions summarized in Appendix B 1. In this case we consider an isolated BH with total mass M and dimensionless spin χ ≡ a/M = J/M 2 . (ii) BH binaries that will enable us to explore, for the first time, the scalar excitation in sGB gravity induced by the strong-field dynamics of a coalescing compact binary in the background, and estimate its potential effect on the GW emission in this theory. In this case we consider a nonspinning binary with total mass M , mass ratio q and initial separation d/M . We briefly describe each of the settings below.

Isolated rotating black holes
We are interested in tracking the formation of scalar hair around a rotating BH. As previously discussed, to linear order in the background spacetime is a solution of vacuum Einstein's equations and, hence, the unique stationary solution is the Kerr metric with mass M and spin J = aM . In Boyer-Lindquist coordinates (t, r, θ, φ) its line element reads where the metric functions are defined as Σ =r 2 + a 2 cos 2 θ , (26b) and the inner and outer horizons are located at In light of the numerical evolutions of the scalar field dynamics it is convenient to introduce a quasi-isotropic radial coordinate R [67,68], such that In the remainder of this section we use r and r ± purely as shorthand notation. In these new coordinates the outer horizon is located at whereas the inner horizon is not part of the domain R ∈ (0, ∞). In contrast to the definition of Refs. [69,70], the transformation (28) yields a finite (nonzero) horizon radius in the extremal limit lim a→M R + = M/4, thus allowing us to set up highly spinning BH backgrounds. After applying the coordinate transformation (28) and performing the spacetime split (23), the 3-metric and gauge functions are given by Here we have analytically continued the lapse function, and its positive (negative) sign corresponds to the exterior (interior) region. The nonvanishing components of the extrinsic curvature (24) read The Gauss-Bonnet invariant that sources the scalar field dynamics reduces to the Kretschmann scalar. For a single Kerr BH it is In practice, our (spatial) numerical domain is described in terms of Cartesian coordinates X i = {x, y, z}. Their relation to spherical coordinates X a = {R, θ, ϕ} is given by x =R sin θ cos ϕ , y = R sin θ sin ϕ , z = R cos θ . (33) Applying the coordinate transformation (33), the spatial line element can be written explicitly as where η ij is the flat space metric and we introduced The extrinsic curvature and shift vector transform according to where Λ i a = ∂X i /∂X a is the Jordan matrix, and K K ab and β a K are given in Eqs. (31) and (30b), respectively.

Black hole binaries
Even within GR, the (near merger) two-body dynamics of BHs with comparable masses has to be solved numerically. The techniques are by now standard and regularly employed in GW source modelling, so we only give a brief summary of the specific ingredients that we use and refer the interested reader to textbooks, e.g. [71][72][73]. Specifically, we follow the ADM-York approach [74,75]. Applying the spacetime split discussed in Sec. III A we can rewrite Einstein's equations -a system of coupled partial differential equations of mixed character -as a set of elliptic-type constraint equations and a set of hyperbolictype evolution equations.
The Hamiltonian and momentum constraints in vacuum GR are To provide initial data (γ ij , K ij )| t=0 describing quasicircular BH binaries, we employ the Bowen-York construction [76]. The time development of the 3-metric and extrinsic curvature is determined by the evolution equations where d t = ∂ t − L β and L β is the Lie derivative along the shift vector. To obtain a strongly hyperbolic and, hence, well-posed initial value formulation of Einstein's equations, we employ the W -version of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [77,78] whose variables are where γ ≡ det (γ ij ) is the determinant of the physical 3-metric, A ij = K ij − 1 3 γ ij K is the tracefree part of the extrinsic curvature, and the last relation for the conformal connection functionΓ i holds becauseγ = 1 by construction. The resulting evolution equations are given explicitly, e.g., in Refs. [79,80] and we complement them with the moving puncture gauge [81,82] where we typically choose β Γ = 0.75 and η β = 1.

C. Gauss-Bonnet invariant in 3+1 form
Next, we express the Gauss-Bonnet invariant (1) in terms of spatial quantities. A particular convenient reformulation is that in terms of the electric and magnetic parts of the Weyl tensor defined as where * W abcd denotes the dual Weyl tensor. By construction the electric and magnetic parts of the Weyl tensor are symmetric, tracefree and spatial, i.e., E ij = E (ij) , γ ij E ij = 0, E ab n a = 0, and likewise for B ij . Comparing Eqs. (41) to the spacetime decomposition of the Weyl tensor and of the GR field equations (38) we find is the tracefree part of the spatial Ricci tensor. Assuming a Ricci-flat background spacetime, i.e. (4) R ab = 0, the Gauss-Bonnet invariant (1) simplifies to This equation, evaluated on the background metric, is the source term in Eq. (21).

D. Scalar-field evolution equations
To simulate the scalar field's dynamics we rewrite its field equation (21) as a set of time evolution equations. To this end, we introduce the scalar field momentum in analogy to the extrinsic curvature. Then, the time evolution is determined by where the Gauss-Bonnet invariant is calculated from Eq. (43).

E. Scalar-field initial data
We will focus on different types of scalar-field initial data that appear most relevant. Initial data 1 (ID1): The first set is trivial initial data This will allow us to verify the formation of nontrivial scalar hair around rotating BHs or around a BH binary solely sourced by the spacetime curvature. Initial data 2 (ID2): To investigate perturbations around isolated BHs we initialize the scalar field and its momenetum as a condensate with a Gaussian profile centered around R 0 with width σ and amplitude A. Specifically, we set where Σ lm is a superposition of spherical harmonics, typically set to Initial data 3 (ID3): These initial conditions represent data for multiple (hairy) BHs. For simplicity, we neglect any linear or angular momenta of the BHs. Since the scalar field equation (21) is linear, we can superpose the static solution of [19]; see (B1) with χ = 0. Then, for N BHs, we have where the field associated to the (a)-th BH with ADM mass m (a) and at position R (a) is in the same quasi-isotropic coordinates (t, R, θ, ϕ) employed to construct the background (GR) initial data. Finally, we note that our initial data for the scalar fields are strictly valid only when the BHs are all at rest. When considering a BH binary, we are neglecting the initial velocity of each BH, which could be taken into account by boosting each of the scalar-field profiles. Neglecting this boost introduces some spurious initial-data effect that is neglible as the scalar field adjusts to its actual configuration during the evolution.

F. Analysing the data
To analyse and interpret the numerical data we extract a number of observables from our simulations each of which we summarize in the following. Waveforms: To calculate the gravitational radiation produced in our background spacetime we employ the Teukolsky formalism [83,84]. In this spinor-inspired approach one defines a null tetrad and a set of complex scalars that contain information about the radiative degrees of freedom. They are constructed from contractions of the Weyl tensor with the tetrad vectors. With an appropriate choice of the tetrad one of these complex scalars, Ψ 4 , encodes the outgoing gravitational radiation. For details of the construction see, e.g., Refs. [68,71,85]. In practice, we measure the Newman-Penrose scalar on spheres of fixed extraction radius R ex and decompose it into multipoles where −2 Y lm are s = −2 spin-weighted spherical harmonics.
In a similar fashion we extract scalar radiation: we interpolate the scalar field Φ onto a sphere of radius R ex and perform a multipole decomposition using spherical harmonics Y lm . In particular, we measure Energy and momentum fluxes: In addition to the waveforms, the energy and momentum fluxes provide crucial insight into the phenomenology of the system and allow us to estimate the order-of-magnitude of metric deformations and radiation at second order without actually evolving it. First, let us recap the energy and momentum fluxes of the gravitational radiation in GR [63,71]. These are fluxes present in the background spacetime and given by where = − (sin θ cos ϕ, sin θ sin ϕ, cos θ). Furthermore, we consider the energy and momentum fluxes of the scalar field, which are associated to the scalar stress-energy tensor T ab and to the Gauss-Bonnet correction G ab . For a generic energy-momentum tensor T ab they can be defined as [63,71,84] where j R and S iR are the radial energy-momentum flux and stress tensor, generically computed from Since our code is implemented explicitly in Cartesian coordinates, we transform these quantities to spherical coordinates where Λ i a = dX i dX a is defined by transformation (33). Since Φ = O( ), the leading-order components of the scalar energy-momentum tensor and of the Gauss-Bonnet correction [i.e., the source terms of modified Einstein's equations (18a)] are O( 2 ). Using the effective stress-energy tensor defined in Eq. (19) we can write the (leading order) energy flux carried by the scalar field as fluxes as and likewise for the momentum flux. Here the energy and momentum fluxes for the scalar field associated to the canonical energy-momentum tensor, indicated by (Φ) , and those associated to the Gauss-Bonnet correction, indicated by GB , are obtained from Eqs. (53) using the relevant flux densities and spatial stress tensors. Using (54) and replacing T ab with T ab (of the scalar) defined in Eq. (15) we find 4 If, instead, we replace T ab with G ab defined in Eq. (17) and insert (20) we get In practice, we take the following steps: (i) we compute j i and S ij during the numerical evolution from Eqs. (57) and (58); (ii) since the numerical code is in Cartesian coordinates we perform a coordinate transformation (55); (iii) we use the radial fluxes to calculate the energy and momentum fluxes through Eqs. (53).
We remark that at O( 2 ) the energy flux contains a contribution from gravitational radiation coming from "mixed" terms ∼ ... h ab denotes the background contribution. We do not explicitly compute h (2) ab , since it would require to numerically evolve Einstein's equations at O( 2 ). However, it is reasonable to assume that this contribution is at most comparable to the energy flux carried by the scalar field. In the inspiral phase, it is expected to be sub-dominant with respect to the scalar field flux, since it is of higher post-Newtonian (PN) order [43]. Therefore, Eq. (56) provides a reliable estimate of the energy flux at O( 2 ).

A. Code description
To simulate BHs in sGB gravity we implemented the field equations (37), (38) and (45) in Canuda [79,86]. Canuda 5 is a novel numerical relativity library that is compatible with the open-source Einstein Toolkit [87][88][89], and capable to simulate BHs in extensions of GR; see, e.g. [68,79,90]. The Einstein Toolkit itself is a community code originally designed to solve the two-body problem in GR. It is based on the Cactus computational toolkit [91,92] and the Carpet boxes-in-boxes adaptive mesh refinement package [93].
In practice, we evolve the background BH spacetime and the scalar field simultaneously, as described in Sec. III. We set up BH initial data using either the TwoPunctures spectral code for binaries [94] or fix the background as a single rotating BH as described in Sec. III B 1. To evolve BH binaries we typically employ Canuda-Lean (an upgraded version of [95]), although our implementation can be combined with other Einstein Toolkit evolution thorns such as McLachlan [96]. We extract apparent or isolated horizon properties of the BHs with AHFinderDirect [97,98] and Quasilocalmeasures [99]. For the wave extraction we typically use our own implementation of the Weyl scalars or the built-in Einstein Toolkit version thereof.
The core thorns 6 developed for the purpose of the present paper are already publicly available [86], and consist of an initial data thorn implementing the prescription of Sec. III E and an evolution thorn implementing Eqs. (45) in the W -version of the BSSN formulation and analysis capabilities as described in Sec. III F.

B. Hair formation around rotating black holes
Following up on previous studies [33,34] we explore the formation of scalar hair around rotating BHs. Stationary solutions have been obtained analytically or generated numerically (see Appendix B 1). Their key feature of the most general family of such solutions is that it is generically singular on the horizon. Imposing that the solution is regular across the horizon selects a sub-family of solutions that is uniquely characterised by the mass and the spin of the black hole. The corresponding scalar configuration is unique, i.e. there are no independent scalar charges [22,32]. This kind of nontrivial scalar configuration around BHs is also known as scalar hair of the second kind. Though it seems natural to discard the singular configurations (as taking the stationary limit of a PDE can lead to dynamically spurious singular solution), it is crucial to demonstrate explicitly that the regular solutions are indeed the endpoints of dynamical processes, such as collapse and mergers. Setup: To explore this process we performed a set of time domain simulations for different initial configurations of the scalar field around rotating BHs with total mass M = 1 and various dimensionless spins, χ ∈ [0, 0.99]. We initialized the scalar field either as vanishing or as a Gaussian composed of the dipole or quadrupole mode located at R 0 = 10M with width σ = 1M and amplitude A = 1M denoted, respectively, as Initial data 1 or 2 in Sec. III E. The different configurations are summarized in Table I.
Our numerical domain with outer boundaries at 256M consisted of nine refinement levels centered around the BH. On the outermost level we typically set the gridspacing to dx = 2.0M . Scalar evolution: In Fig. 1 we illustrate the formation of scalar hair exemplarily around a highly rotating BH with dimensionless spin χ = 0.99. Specifically, we present the scalar field profile Φ along the x-axis (equatorial plane) and z-axis (direction aligned with BH spin) at different instances in time. To compare the numerically obtained field after an evolution time of t = 300M to the analytic solution in a large radius expansion (see Eq. (B2)), we transform the latter to quasi-isotropic coordinates that coincide with the numerical ones for R 5; see Fig. A1 of [34]. We find excellent agreement between the numerical and analytic solution at large distances.
Near the horizon, the analytical solution is less accurate mostly in the polar direction, whereas on the equatorial plane it provides a good approximation to the exact numerical result also close to the horizon. We further quantify that statement in Fig. 2 where we benchmark the numerical scalar field solution after t = 300M against the analytic one, Eq. (B2). We focus on highly spinning BH backgrounds and find excellent agreement for sufficiently large distance where the analytic approximation is valid.  Waveforms: As the system settles down to a hairy BH in sGB gravity, the scalar field sheds away all non axisymmetric modes. This is illustrated in Fig. 3, where we show the l = m = 1, 2 modes of the scalar, measured at R ex = 60M , around a BH with dimensionless spin χ = 0.7 and χ = 0.99. The latter case exhibits a decay pattern that is consistent with the quasi-normal ringdown of a free scalar around a Kerr BH [100], i.e., it is determined by the dominant (scalar) QNM frequency. In contrast, the χ = 0.7 case shows a more complex, modulated ringdown that indicates the presence of multiple modes with comparable excitation frequency and amplitude. Inspecting Eq. (21) this is not surprising; the scalar field is sourced by the Gauss-Bonnet invariant. Hence, we expect the superposition of two types of modes 7 : (i) a scalar mode corresponding to the decay of a free scalar field around a Kerr BH; (ii) a mode driven by the background curvature. We refer to the former as "scalar-led" and the latter as "gravitational-led" modes. The signal's specific morphology then depends on the amplitude with which each of those modes are excited. To identify the composition of the ringdown signal we perform a twomode fit where we suppress multipole indices (lm), we indicate the dominant and subdominant mode by numeral subscripts, their real and imaginary parts by subscripts "R" or "I", and the relative amplitude between the modes is δA 2 = A 2 /A. We summarize the results in Table II. In particular, we provide estimates of the dominant ringdown frequency ω 1 , the real part of the secondary mode ω 2 , and their relative amplitude. In practice, we cannot accurately estimate the decay rate ω 2I of the secondary mode. Additionally, we compare our time-domain estimates to frequency domain calculations of scalar QNMs [100] and find agreement within 6%. If the l = 2 multipole is present in the scalar initial data and, hence, in the ringdown signal, the secondary, gravitational-led mode's frequency coincides with that of a l = m = 2 gravitational perturbation. Although there is no gravitational analogue of the dipole mode, we find good agreement with a back-of-the-envelope estimate ω G 11 ∼ ω G 22 /2. Both findings support our expectation of the presence of a gravitational-led mode in the scalar field emission.
There are, however, some caveats in performing these fits: the early response is followed by a (in most cases) relatively short ringdown signal that transitions to the latetime tail (not shown in the plots). Combined with the uncertainty regarding the end of the direct response and the starting point of the actual ringdown -a hindrance that is not even completely resolved within GR [102,103] -this results in the uncertainties quoted above.
We expect that the occurrence of a secondary, gravitational-led mode to be a generic feature for all spins, but in some cases (indicated by a dash in Table II) with excitation factors that could be too small to  We denote the dominant and the (real part of the) first subdominant mode, estimated from the time domain data, as ω 1,lm and ω 2R,lm , respectively. Their estimated relative amplitude is δA2 = A2/A. A dash indicates that the ringdown is dominated by a single mode and we could not accurately estimate a secondary one. In the last column we provide the scalar QNM frequency ω f,S lm of a Kerr BH calculated in the frequency domain [100]. Where a secondary mode is present it agrees well with the corresponding gravitational mode (in case of l = m = 2) or with ∼ ω f, be extracted from our fit. This seems to be the case for small spin and near extremality. Interestingly, the maximum excitation factor of the secondary mode seems to occur for the case χ ≈ 0.7, which is also phenomenologically relevant since it is approximately the final spin of a BH remnant from the coalescence of two slowly-spinning BHs. In this case our analysis predicts the emission of scalar radiation at two dominant frequencies, the scalar and the gravitational fundamental QNMs of the corresponding Kerr BH. We also expect that this effect occurs at higher order. In particular, the corrections to the gravitational ringdown (not computed here since they are of O( 2 )) should contain modes of the O( ) scalar field, in particular the l = m = 1 and l = m = 2 scalar QNMs a Kerr BH. Computing the excitation factors of the latter requires to solve the field equations at O( 2 ) and is left for future work.

C. Black hole binaries
Setup: We focus on nonrotating BHs because they are the simplest BH solutions in sGB gravity that develop scalar hair. The BHs initial separation for the results presented below is d = 10M , though we have considered other initial distances as well. The system's total mass M = m 1 + m 2 = 1 and the mass ratio varies between q = m 1 /m 2 = 1, 1/2, 1/4. We summarize details of the inital BHs' parameters and the final state in Table III. This includes the dimensionless spin χ f of the final BH computed from [85] where A AH and C e denote the area and equatorial circumference of the apparent horizon, its Christodoulou mass where M irr = A AH /16π is the irreducible mass, and the energy E GW /M radiated in GWs. These results are in III. Parameters for nonspinning BH binaries with total mass M = m1 + m2 = 1, mass ratio q = m1/m2 and symmetric mass ratio η = m1m2/M 2 . mi denote the physical BH masses. Their initial separation is d = x1 − x2 = 10M (situated along the x-axis). We denote their orbital and tangential momenta Pr = P2x = −P1x and P ⊥ = P1y = −P2y. We also denote the dimensionless spin χ f and mass M f of the remnant BH, and the energy E GW ∞ /M radiated in form of GWs and the second-order We considered scalar initial data ID1 and ID3 of Sec. III E, i.e. either starting with an initially vanishing scalar or superposing two hairy solutions. Except for an early transition or built-up period both setups yield equivalent results as we illustrate in Fig. 4 for the case of BHs with a mass ratio q = 1/2. Therefore, in the following we present results only for initially hairy BHs. To simulate these systems we set up a numerical grid with outer boundary at 256M consisting of 8 refinement levels. For the results presented here we used a grid spacing of dx= 1.6M in the outermost refinement level. To validate our results we performed a convergence analysis of run BBH q1 using additional grid spacings of TABLE IV. Postmerger ringdown frequencies of the scalar field obtained from a two-mode fit (60). We list the scalar-led mode frequency M f ω S lm and the real part of the gravitationalled mode M f ω G lm , rescaled by the final BH mass M f . We also denote the relative amplitude δAG = AG/AS. A dash indicates that the ringdown is dominated by a single mode and we could not accurately estimate a secondary one.
Waveforms: We present the background gravitational waveform and the O( ) scalar waveform for binaries with mass ratio q = 1 and q = 1/2, 1/4 in Figs. 6 and 7, respectively. All presented waveforms are shifted in time such thatt/M = (t − t merger − R ex )/M = 0 indicates the maximum in the dominant gravitational mode as measure for the time of merger. The waveforms exhibit the typical morphology: a sinusoid with increasing frequency that is driven by the orbital motion of the BHs, the highly nonlinear merger followed by the exponentially damped ringdown. During the inspiral we compare the numerical results to the analytical expressions obtained at leading PN order in [24] (see Appendix B 2). We remark that these expressions depend on the time-dependent orbital frequency Ω(t), which, at this PN order, can not be obtained with good approximation. Therefore, we extract the orbital frequency from the numerical data (measuring, at each half-cycle, the wavelength of the gravitational waveform). Within this approach, which is similar to that used in [43], the comparison between PN and numerical results concerns the amplitudes of the scalar waveforms, while their phases agree by definition.
Interestingly, while the scalar signal for l ≥ 2 is qualitatively similar to the gravitational waveform and displays the classical chirp, the dipole is qualitatively different. As shown in Fig. 7, the frequency of the dipole mode grows as expected during the merger, but the amplitude remains almost constant. This is a strong-field behavior that is not captured by the PN approximation. A potential explanation of this behaviour is that the scalar configuration ceases to be dominantly dipolar before the merger, i.e. the dynamical evolution of the scalar is more complex and it involves additional oscillations and reconfiguration. Our simulations, and in particular the time evolution of the scalar distribution, do seem to be consistent with this explanation, though limitations in resolution do not allow us to make a conclusive statement.
In the post-merger phase the background approaches a stationary spinning BH, so we expect to observe the same multiple ringing discussed in the previous section for an isolated BH. This is confirmed in the insets of Figs. 6 and 7 and by the postmerger ringdown frequencies extracted from the scalar waveform using the twomode fit (60) and presented in Table IV. Note that, in contrast to the single BH case, the background is now a perturbed BH plus gravitational radiation, both of which modify the source term of the scalar field. In particular, gravitational radiation seems to cause an enhancement of the gravitational-led quadrupole modes, which dominate over the scalar-led one in some configurations (see, e.g., the l = m = 2 case for q = 1, 1/2 in Table IV).
Finally, the scalar field monopole for the same values of the mass ratio is shown in Fig. 8. The pre-merger amplitude is larger for smaller values of the mass ratio q, while the final amplitude is approximately independent of q. This behaviour can be understood noting that the post-merger amplitude is, as a first approximation, Φ α GB /(2M r) (see App. B 1 a), where (by construction) the total mass M is the same in all cases, i.e., independent of the mass ratio. Instead, when the two BHs are well separated, the scalar field amplitude is for sufficiently large radii encompassing the entire binary. Therefore, the ratio between the pre-merger and the postmerger amplitude is expected to be determined by the (inverse of the) symmetric mass ratio η. In particular we have 1/η = 4, 4.5 and 6.25 for mass ratios q = 1, 1/2 and 1/4. These values are in agreement with Fig. 8. Energy and momentum fluxes: Next, we investigate the energy radiated in gravitational and scalar waves. We compute their energy fluxes using (52)-(53) with (56), i.e. accounting for both the canonical scalar's and Gauss-Bonnet contributions to the energy flux. We furthermore estimate the total radiated energy by integrating Eqs. (52) and (53) in time, measuring it at different extraction radii and performing the extrapolation and likewise for the second-order flux transported by the scalar field. We estimate the extrapolation error by comparing to E/M = E ∞ + B/R ex + C/R 2 ex and find ∆E GW ∞ /E GW ∞ 0.8% and ∆E 7%. The results are summarized in Table III. In Fig. 9 we present the fluxes for all three configurations. The background, i.e., GW flux (black solid lines) follows the common pattern: it increases monotonically in amplitude as the BHs circle around each other for the last few orbits, culminates in a peak during their merger, and decays exponentially as the newly born BH rings down to a Kerr BH. We also show the second-order energy flux carried by the scalar waves (56) (blue dashed lines) together with the canonical scalar field energy fluẋ E (Φ) (red dot-dashed lines), rescaled by the appropriate power 2 of the expansion parameter. Exemplarily, we set = 0.01. We observe that the second-order scalar flux is dominated by the contribution of the scalar stress-energy tensor.
The morphology of the signal is determined by the orbital dynamics and monotonically increases during the inspiral of the background BHs. The canonical scalar field flux also exhibits a peak during the merger that is predominantly determined by the monopole mode, as is illustrated by the green dot-dashed lines in Fig. 9. This is because the system changes rapidly from two nonrotating BHs, each with its own scalar hair determined by (63) to a single rotating BH with a new scalar configuration of this form but with larger mass and, hence, smaller scalar charge.
The ratio between scalar and gravitational radiation dramatically increases as the mass ratio decreases. This is because the scalar charge is determined by the smallest mass scale in the system yielding the largest curvatures, and then undergoes a transition to the final BH mass. So while this characteristic scale changes at most by a factor of two in the equal-mass case, it can be vastly different as we decrease the mass ratio. Whether and how this trend continues for higher mass ratios is beyond the scope of the paper and will be presented in a more detailed parameter study elsewhere.  9. (Color online) Gravitational (black solid lines) and scalar second-order (blue dashed lines) energy fluxes measured at Rex = 100M for mass ratios q = 1 (left), q = 1/2 (middle) and q = 1/4 (right). They have been shifted in time such thatt = 0 coincides with the time of merger. We rescaled the scalar field flux by the appropriate powers of the expansion parameter, and set = 0.01. We also show the canonical scalar field flux (red dot-dashed lines) and its monopole contribution (green dotted lines) that exhibits a peak during the merger when the scalar field adjusts to the final, single BH solution.

V. RANGE OF VALIDITY AND OBSERVATIONAL CONSTRAINTS
Our results allow us to put new constraints on the Gauss-Bonnet coupling with previous and upcoming GW detections. However, before doing so we need to quantify the validity of the low-energy perturbative expansion presented in Sec. II B up to first order. Therefore, we consider the "instantaneous" range of validity as well as integrated, secular effects. In the former case we demand that the scalar energy flux at a given instant t be much smaller than the GW flux, whereas in the latter case we require that the dephasing due to scalar emission accumulated during the inspiral be smaller than the GW phase. Although we have not explicitly evolved the second-order scheme we can estimate deviations at this order from the source terms in Eqs. (18a) as they only depend on background or first-order quantities that we obtained in our numerical simulations.

A. Instantaneous range of validity
We start by investigating the instantaneous range of validity, that is, we check for which couplings the perturbative expansion (8) remains applicable at every timestep in our simulation. To this end we compare the secondorder energy flux (carried by the scalar waves) with the background GW flux. Then, a necessary condition for the perturbative expansion to apply iṡ In Fig. 10 we present the (instantaneous) bounds on the dimensionless coupling obtained from We observe that the tightest constraints come from binaries with small mass ratio, because in such case the first of the binary components has a smaller mass than in the equal-mass case (recall that the total mass is fixed to unity) and therefore yields a larger dimensionless coupling ∼ α GB /m 2 1 . This is consistent with the fact that sGB gravity is a strong-curvature correction to GR, so the strongest effects come from small BHs for which the near-horizon curvature is large.
At the merger nonlinear effects dominate, the scalar field transitions to a new configuration (see Fig. 8) while the final BH forms, and we observe a burst of scalar energy flux. This is indicated by the dip aroundt = 0 in Fig. 10, where the allowed value of the dimensionless coupling drops by about an order of magnitude.

B. Secular effects and dynamical range of validity
While the instantaneous range of validity is a first check, we find it more instructive to explore the influence of the scalar field on the binary's evolution. As we have seen in Sec. IV the (background) spacetime dynamics source scalar radiation that needs to be accounted for in the full energy budget. Because a BH binary in sGB gravity emits not only gravitational but also scalar waves, the inspiral is accelerated as compared to the same system in GR. This yields an increase in the orbital angular velocity of the binary and causes a dephasing ∆φ in the gravitational waveform when compared to GR. Since we adopted a perturbative approach up to first order only, we cannot compute this dephasing -a second-order effect -directly, but we can provide qualitative estimates based on our numerical results.
In accordance with this approach we expand the orbital phase φ and orbital frequency Ω =φ as in Eq. (7). Recalling that the first-order contribution to the metric vanishes [see Eq. (13)] we note that there is also no firstorder correction to the orbital phase and frequency, and the first nonvanishing shifts appear at second order in the perturbative expansion, i.e., Ω Ω (0) + 1 where φ (0) and Ω (0) = |d×ḋ|/d 2 are the orbital phase and frequency of the BH binary calculated from the distance d between the puncture positions, and ∆φ = 2 φ (2) /2, ∆Ω = 2 Ω (2) /2 are the Gauss-Bonnet corrections to the phase and frequency. The validity of the perturbative expansion requires φ (0) ∆φ, where Therefore, the threshold is We evaluate φ (i) (t) througḣ whereĖ (0) =Ė GW andĖ (2) are the GW flux in GR and the second-order energy flux carried by the scalar waves, respectively. While these are computed from the numerical data using Eqs. (52a) and (56), we can only estimate the last term dE (2) /dΩ in Eq. (71) using the PN approximation. As shown in [27], the scalar interaction binding energy of a compact binary is where s, t are the (leading) multipole numbers of the scalar field emission from the two BHs, and µ 1 , µ 2 are their multipole charges (see App. A for the different notations for charges). Note that, while in the case of Chern-Simons gravity studied in [43] the leading-order contribution to the binding energy is the dipole-dipole interaction (only present when the BHs are rotating), in the case of sGB gravity the leading-order binding energy contribution is the monopole-monopole interaction E M M , which does not depend on the BH spins. Therefore, s = t = 0 in Eq. (72), and where µ i = α GB /(2m i ) = 2 M 2 /m i is the scalar charge (see [24] and App. B 1 a) of the i-th body, and, at leading PN order, Ω = (M/r 3 ) 1/2 8 . Therefore, and Now we have all the ingredients to compute the condition (69) for the expansion to be consistent also during the binary's evolution. We illustrate it in Fig. 10 as function of time (shifted by the time of merger). As in the previous case, the bounds on the validity of the perturbative approach up to first order become more stringent with decreasing mass ratio and near the merger, as one might expect in this highly nonlinear regime.

C. Observational bounds
While an accurate computation of the observational constraints on the Gauss-Bonnet coupling would require computing the O( 2 ) corrections to the gravitational waveform (which is a higher-order effect than the ones computed in this work), we can estimate an observational constraint based on the assumption that no GW dephasing ∆φ GW = 2∆φ = 2 φ (2) + O( 3 ) induced by some non-GR extension is observed by a given (present or future) GW detector and, hence, it must be at least below the detector's GW phase uncertainty ∆φ det . This is a conservative estimate since: (i) smaller effects could be constrained by comparing directly waveform models in GR and in modified theory; and (ii) the contribution of the energy flux at O( 2 ) due to gravitational radiation, which we are neglecting, would further increase the dephasing. Computing the full energy flux to O( 2 ) is an interesting problem which we leave for future work.
We remark, however, that in this estimate we are not taking into account O( 2 ) shift in the physical masses of BHs in sGB gravity. For instance, in EdGB gravity the physical mass of an isolated, static BH measured from the asymptotic limit of the metric is M ∼ M (1 + 49α 2 GB /(20480πM 4 )) [21,22,105]. Since the orbital energy −Gm 1 m 2 /(2r) and the monopole scalar binding energy (73) have similar expressions, the mass shift and the scalar monopole energy may have degenerate effects 9 . On the other hand, we have verified that neglecting the contribution of the monopole binding energy dE (2) dΩ in Eq. (71) would modify the dephasing by 20%. Moreover, the dissipative and conservative contributions to the dephasing (i.e. the two terms in Eq. (71)) have different dependence on the parameters of the binary system. Thus, unless fine-tuning cancelations occur, our computations provide a correct order-of-magnitude estimate of the dephasing due to sGB gravity.
Absence of non-GR dephasing appears to be the case in previous LIGO detections [8,106], which can therefore be used to put actual upper bounds. Using the same assumption of a negative search, we can also forecast observational bounds for third-generation detectors such as the Einstein Telescope or Cosmic Explorer and for the space-based LISA mission. This can be translated into the bound where we again evaluate φ (2) as discussed in the previous section. Specifically, we consider the reference systematic phase errors (i) ∆φ LIGO 0.1 for LIGO's O2 run [106]; (ii) ∆φ 3G 0.01 for third generation ground-based detectors [107]; and (iii) ∆φ LISA 0.01 for LISA [11,108].
While LIGO allows us to constrain the dimensionless Gauss-Bonnet coupling to be O 10 −3 , thirdgeneration ground-based detectors and the space-based LISA mission will provide bounds that are about one order of magnitude more stringent. Note, however, that since α GB ∼ M 2 the bounds become much weaker for heavier BHs, so supermassive BHs are not good probes of higher-curvature corrections to GR.
This can be seen in Fig. 11, where we present bounds on the dimensionful GB coupling that depends on the mass of the system. We choose a binary with M = 20M for ground-based detectors and M = 10 5 M for LISA. Furthermore, we summarize the constraints on both the dimensionless and dimensionful coupling constants in Table V. For ground-based detectors, we consider total binary masses of M = 20M or M = 60M , corresponding to the lightest and most massive BH binary detected so far, whereas we consider M = 10 5 M for space-based detectors.
Using this phase information we can put the most stringent observational constraints on the Gauss-Bonnet coupling to date. We choose to consider a binary with mass ratio q = 1/2 and M ∼ 20M because its characteristics strongly resemble X GW151226 [109]. Hence, our results, within the caveats of our approach, place the constraint 10 √ α GB 2.7 km . 10 We recall (see Sec. I) that the low mass -ray binary constraint |α YYSP | 1.9km [25] and GW constraint |α YYSP | 5.1km [7] from GW151226 translate to |α GB | 10km and |α GB | 27km, respectively, in our notation; see Eq. (A1).
This bound can improve with detections of less massive systems or binaries with a lower mass ratio. As indicated in Table V and Fig. 11, upcoming third-generation ground-based detectors will have the potential to place even more stringent constraints on Gauss-Bonnet-type modifications to GR.
Finally, in Fig. 12 we show how the absence of a dephasing from LIGO's O2 run and from third generation ground-based detectors would constrain the dimensionful coupling parameter α GB for a range of possible source masses. Let us emphasize that these observational bounds are already stronger than previous ones and can cover the entire range of total masses up to M ∼ 100M for mass ratios q = 1/4.

VI. SUMMARY
In this paper we have considered spinning isolated BHs and binary BH systems in sGB gravity, described by the action (2). We have worked perturbatively in the coupling constant α GB but nonperturbatively in the fields. We have solved the field equations to first order in α GB and computed the right-hand side of the modified Einstein equations to second order. This allowed us to simulate, for the first time, the dynamics of the scalar field in the highly dynamical strong field regime probed by binary BHs, to obtain complete scalar waveforms in this context, and to check the range of validity of our scheme.
We first investigated the formation of scalar hair around a Kerr BH with arbitrary spin and studied the scalar QNM ringing. The latter contains at least two dominant modes, corresponding to the scalar and gravitational QNMs of a Kerr BH. We expect the same result to hold qualitatively at second order as well. That is, we expect that the post-merger ringdown waveform from a BH binary coalescence in sGB gravity will contain both gravitational-led and scalar-led QNMs.
We then investigated the emission of scalar waves from the coalescence of nonspinning BH binaries with various mass ratios. While the scalar radiation generally displays the typical chirp signal, the dipole mode (radiated only for unequal mass binaries) displays a more peculiar behavior and no chirping. This suggests that the scalar field  (76) and, exemplarily, considering a GW151226-type source. We present estimates for LIGO's O2 run (black lines) and forecasts for third generation ground-based detectors (orange lines) and LISA (red lines). We consider all simulated mass ratios q = 1 (solid lines), q = 1/2 (dashed lines) and q = 1/4 (dashed-dotted lines). For comparison we also show the so far most stringent bound coming from the existence of light BHs (green dotted line).  (76) and estimated near the merger, for a range of source masses M = m1 q+1 q . We start with the smaller BH's mass at m1 = 3M . We present estimates for LIGO's O2 run (black lines) and forecasts for third generation ground-based detectors (orange lines) for all simulated mass ratios q = 1 (solid lines), q = 1/2 (dashed lines) and q = 1/4 (dashed-dotted lines). For comparison we also show the so far most stringent bound coming from the existence of light BHs (green dotted line). exhibits interesting dynamics in the pre-merger phase.
In all of our simulations, the axisymmetric components of the scalar field approach the profile of the stationary hairy BH. This BH remnant is characterized by the mass and spin only (primary hair) and all other multipole moments, including the scalar charge (secondary hair), can be written in terms of those. Any deviation from this stationary multipolar structure is radiated away during the merger and ringdown, just like in GR, and there is a unique scalar configuration that acts an the endpoint of dynamical evolution for a BH of given mass and spin. This is a significant generalization of similar results obtained previously [33,34] for spherically symmetric spacetimes.
We calculated the scalar energy flux emitted in our binary simulations. This flux leads to modifications of the binary evolutions with respect to GR that are imprinted on the standard tensor modes as dephasing. The effect is of second order in the coupling so we could not compute it directly. However, the calculation of the flux allowed us to estimate this effect and derive some qualitative observational bounds on the sGB coupling constant. Requiring that the dephasing due to the GB corrections is smaller than LIGO's phase sensitivity, our results imply the preliminary constraint √ α GB 2.7km on the sGB coupling constant for a GW151226-type source, i.e., a BH binary with mass ratio q = 1/2 and total mass M = 20M . The detection of the GW151226 event does enforce this (approximate) bound, which may be viewed as the strongest constraint to date on the coupling constant of sGB gravity.
Our results indicate that systems with smaller mass ratios and smaller total mass can put even more stringent bounds. Future third-generation detectors will improve such a constraint by at least a factor ∼ 3; taking into account the fact that 3G will detect several sources, possibly with a larger distribution of mass ratios, the improvement may reach one order of magnitude. On the other hand, projective bounds on α GB based on dephasing from the future space-based LISA detector are unlikely be competitive: even though LISA will provide the most stringent constraints on the dimensionless coupling constant, = α GB /(4M 2 ), it will probe larger total masses/ curvatures and this will significantly weaken the constraints on α GB itself.
It should be stressed that dephasing due to dipolar emission is not the only way to constrain α GB . In fact, the dephasing constraint relies only on energy loss. One expects to obtain significantly stronger bounds by attempting to fit complete inspiral-merger-ringdown waveforms (when available) to specific events or even by stacking events.
As we have argued in detail, within a perturbative approach in the coupling constant (but not in the fields), the leading-order effects of a scalar field on BHs and their binaries are driven by the coupling of the scalar to the Gauss-Bonnet invariant in the absence of parity violations. The coupling to the Pontryagin density comes at the same order if parity invariance is broken. Hence, our results, potentially combined with those of Ref. [43] in Chern-Simons gravity provide the complete scalar dynamics, to linear order in the coupling constant.
Work on second-order effects and on the metric backreaction is ongoing. Future work will also focus on a more detailed analysis of the parameter space for BH binaries, on the mode excitation for the remnant BHs with various spin values, and on a more rigorous analysis on the detectability of sGB corrections in the GW signal from BH coalescences.

ACKNOWLEDGMENTS
We thank K. Yagi for sharing his notes used in Sec. B 1 b and for useful discussions. We also thank M. Okounkova and L. Stein for useful discussions and comments. H.W. acknowledges financial support provided by the European Union's H2020 research and innovation program under Marie Sklodowska-Curie grant agreement no. BHstabNL-655360, by the Royal Society University Research Fellowship UF160547 and the Royal Society Research Grant RGF\R1\180073. She also acknowledges partial support by grants FPA-2016-76005-C2-2-P, AGAUR SGR-2017-754. P.P. acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement no. DarkGRA-757480. T.P.S. acknowledges partial support from the STFC Consolidated Grant No. ST/P000703/1. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. The authors would like to acknowledge networking support by the COST Action CA16104. We thankfully acknowledge the computer resources at Marenostrum IV, Finis Terrae II and LaPalma and the technical support provided by the Barcelona Supercomputing Center via the PRACE grant Tier-0 PPFPWG, and via the BSC/RES grants AECT-2017-2-0011, AECT-2017-3-0009 and AECT-2018-1-0014. The code developed for this paper is based on the Einstein Toolkit [89] and public at [86]. The xTensor package for mathematica [110,111] has been used. The contributions of m = 0 and m = 1 agree with those given in Ref. [24], whereas we explicitly present also the other multipoles that are relevant in our case. As expected, in the equal-mass case only the even multipoles are nonvanishing.
From Eq. (51), it is easy to check that the leading-order contribution to the radiative mode Φ ll comes only from Φ m=l , i.e.
Finally, the orbital parameters needs to be evolved adiabatically, i.e. b → b(t) and ω → ω(t), where ω(t) is extracted from the evolution at zeroth order (i.e., the GR coalescence).