Perturbed black holes in Einstein-dilaton-Gauss-Bonnet gravity: Stability, ringdown, and gravitational-wave emission

Gravitational waves emitted by distorted black holes---such as those arising from the coalescence of two neutron stars or black holes---carry not only information about the corresponding spacetime but also about the underlying theory of gravity. Although general relativity remains the simplest, most elegant and viable theory of gravitation, there are generic and robust arguments indicating that it is not the ultimate description of the gravitational universe. Here we focus on a particularly appealing extension of general relativity, which corrects Einstein's theory through the addition of terms which are second order in curvature: the topological Gauss-Bonnet invariant coupled to a dilaton. We study gravitational-wave emission from black holes in this theory, and {\bf(i)} find strong evidence that black holes are linearly (mode) stable against both axial and polar perturbations; {\bf(ii)} discuss how the quasinormal modes of black holes can be excited during collisions involving black holes, and finally {\bf(iii)} show that future ringdown detections with large signal-to-noise ratio would improve current constraints on the coupling parameter of the theory.


I. INTRODUCTION
The historical detection of gravitational waves (GWs) by the LIGO/Virgo Collaboration has marked the beginning of a new era in astrophysics and the birth of GW astronomy [1]. The next generation of detectors will routinely observe the coalescence of compact objects, such as black holes (BHs) and neutron stars. These observations will probe, for the first time, the highly dynamical regime of strong-field gravity and may provide the answer to long-standing issues [2][3][4]. Is cosmic censorship preserved in violent gravitational interactions? Do GW observations carry incontrovertible evidence for the event horizon of BHs? Can we pinpoint, in gravitational waveforms, the signature of the light ring or of ergosurfaces?
Simultaneously, the entire coalescence process can be used to constrain gravity theories in novel ways [4][5][6]. That general relativity (GR) is not the ultimate theory of gravity is a possibility that should be entertained in the light of several observations (such as those related to the dark-matter and the dark-energy problems), and of the difficulty to reconcile GR with quantum field theory [5]. Although such an extension of GR is unknown-and a robust spacetime parametrization in strong-field gravity is lacking-GW observations will help us to exclude or to strongly constrain wide classes of alternative theories. The inspiral stage, for example, when the two objects are far apart, can teach us about possible extra radiation channels [4,[6][7][8], while the final ringdown stage-when the end-product is relaxing to its final state-provides for remarkable tests of GR, through the measurement of the characteristic quasinormal modes (QNMs) [9]. In GR, as well as in essentially any relativistic theory of gravity, BHs are extremely simple objects described by only a handful of parameters. Accordingly, their QNMs are completely characterized by only a few parameters as well. For example, Kerr BHs in GR are characterized by their mass and angular momentum, and so are their QNMs. In a nutshell, measurement of one single QNM (i.e, a ringing frequency and a decay time scale [9,10]) allows for a determination of the BH mass and angular momentum. The measurement of a second QNM tests GR [10][11][12][13]. In the context of modified theories of gravity, a second QNM can be used to measure possible extra coupling parameters, as was shown recently for a theory with an extra vector degree of freedom [8].
Some of the most viable and appealing modifications of gravity are those obtained via the inclusion of extra scalar fields-such as scalar-tensor theories of gravity-or of highercurvature terms in the action, or both. Higher-order gravity is generically motivated by UV corrections, which also arise naturally in some low-energy truncations of string theories. The paradigmatic case, and the one we focus on here, is Einsteindilaton-Gauss-Bonnet (EDGB) gravity, described by the action [5,14] where is the Gauss-Bonnet topological term, S m represents the matter sector and we use (throughout this work) units for which the Newton's constant and the speed of light are unity, G = c = 1. Current best constraints on the coupling constant α are √ α < 10 km [5,15] 1 .
We will close two important gaps in the literature concerning BHs in this theory: we first find strong evidence that static EDGB BHs are linearly stable, and then compute gravitational waveforms from plunging particles, which can be argued to be an indicator of how BH collisions proceed in this theory. Finally, we discuss the constraints on the coupling parameter α from current and future ringdown observations.

II. FRAMEWORK
The equations of motion obtained by extremizing (1) with respect to the metric and dilaton field are given by [14] where G ab = R ab − 1 2 g ab R is the Einstein tensor, T ab is the matter stress-energy tensor, and K ab = (g ac g bd + g ad g bc ) idjk ∇ l R cl jk ∂ i e φ , where abcd is the contravariant Levi-Civita tensor,R ab cd = abij R ijcd .

A. BH solutions in EDGB gravity
BHs in EDGB gravity are scalar-vacuum solutions of the above equations, which were first constructed analytically in spherical symmetry, in the small-coupling regime [16], where M is the BH mass. In spherical symmetry, the line element reads where dΩ 2 is the standard unit 2-sphere line element. Both functions A(r) and B(r) and the scalar field φ can be expanded in powers of the coupling parameter ζ, and the corresponding solutions can be found by solving the field equations (3)-(4) perturbatively (see also Ref. [17]). Further details are given in Appendix A. Nonperturbative solutions were investigated numerically in Ref. [14] for static geometries and in Ref. [18] for slowly rotating BHs to first order in the spin. It was shown that static BH solutions exist only up to a maximum value of ζ, namely [18] 0 ≤ ζ 0.691 .
Because ζ is strictly less than unity, higher-order perturbative expansions [19] are accurate almost in the entire parameter space.
Slowly rotating solutions were described numerically in Ref. [18], and analytically in the small-coupling regime in Refs. [20,21], and were recently extended to higher order in the coupling and in the spin parameter [19]. In the latter case, the line element (as well as the scalar field) can be expanded in a complete basis of orthogonal functions according to their symmetry properties [19].
Numerical solutions describing rotating BHs for arbitrary coupling and spin were found in Ref. [22] and have been recently thoroughly discussed in Ref. [23].
The linearized mode stability of spherically symmetric BHs against radial fluctuations was studied in Ref. [24], whereas axial gravitational perturbations were studied in Ref. [18]. In the polar sector, the linear stability of EDGB BHs was analyzed in Ref. [25], focusing in the particular regimes where perturbations are dominantly gravitational or scalar, as well using a high-frequency analysis of the perturbations. In addition, axial quasinormal modes of neutron stars in EDGB theory were studied in Ref. [26].
Due to the cumbersome field equations, there is currently no result concerning the stability of nonrotating EDGB BHs for the most relevant gravitational polar sector, and there is no stability analysis for rotating solutions. In this work we partly fill this gap by performing a full linear (mode) stability analysis of static EDGB BHs.

B. Perturbed BHs in EDGB
We are interested in understanding how BHs in EDGB theory respond to small perturbations, as those induced by a fluctuation of the metric or of the dilaton or by a small external perturbing object. We focus our attention in both dilaton fluctuations in vacuum and those induced by a small pointlike particle plunging into the BH. The first case will describe the latetime behavior of perturbed EDGB BHs, which also dominates the ringdown signal from a distorted BH formed in a coalescence. Pointlike particles, on the other hand, are a good proxy for small BHs or neutron stars falling into massive BHs, but are also known to provide reasonably accurate estimates even for equal-mass BH collisions [27][28][29].
Pointlike particles of mass µ M are modeled by [30][31][32] with dτ = dλ −g abẋ aẋb . We will consider pointlike objects with a scalar charge which is entirely due to the Gauss-Bonnet coupling. In other words, we will investigate BH spacetimes of which the scalar charge arises purely from the Gauss-Bonnet term. The theory we consider contemplates no other couplings to matter, and therefore pointlike particles follow geodesics. This need not be the case generically, and nontrivial couplings to matter can be envisioned [31,33]. These couplings will certainly influence the motion of particles and the radiation in collisions but will not affect the intrinsic ringdown properties of the spacetime. We consider a spherically symmetric EDGB BH distorted by either the pointlike particle or through some fluctuation in the metric or scalar field. At the linearized level, the full geometry is described by where ε 1 is a bookkeeping parameter, g (0) ab is described by (7) and φ 0 (r) is the corresponding background scalar. As background solution, we consider both a perturbative solution [19] up to O(ζ 6 ), and a numerical solution for arbitrary values of ζ.
The fluctuations h ab and δφ are functions of (t, r, θ, ϕ). Einstein's equations can be further simplified by Fourier transforming these quantities and by expanding them in (tensor and scalar) spherical harmonics, e.g., where Y lm are the standard spherical harmonics and ω is a Fourier frequency. The metric perturbations can be decomposed in terms of tensorial spherical harmonics [30,34]. By using this decomposition, perturbations naturally split into two sectors according to their parity (either axial or polar). Axial perturbations of EDGB BHs are simpler, because they are decoupled from the scalar-field perturbations [18]. Here, we shall consider the two sectors of the perturbations.
In the Regge-Wheeler gauge [34], the polar sector of the metric perturbations is given by while the axial sector reads In the above definitions, the metric perturbations depend only on t and r. Note that we have already specialized the spacetime to axial symmetry: for the cases handled here, one can always rotate the coordinate axis such that the spacetime is axially symmetric. We shall Fourierdecompose these perturbation functions as, e.g., X(t, r) = (2π) −1/2 dωX(ω, r)e −iωt . Likewise, the stress-energy tensor can be decomposed into spherical harmonics [30,35]. In the radial plunging case considered here, the particle only disturbs the spacetime in the polar sector, and its stress-energy tensor can be written as where the functions A lm and A lm , like the dilaton and metric perturbations, can be Fourier decomposed, such that they depend only on r and ω. The explicit form of these functions is given in Appendix B.
The linearized dynamics is governed by a coupled system which can be obtained by expanding the field equations (3)-(4) up to first order in the perturbation functions h ab and δφ. The equations take the schematic form (no sum on j) where j = (p, a) for the polar and axial sectors, respectively, Ψ j is a column vector, the components of which are Ψ p ≡ (H 1 , K, φ 1 , φ 1 ) for the polar sector and Ψ a ≡ (h 1 , h 0 ) for the axial sector, respectively. The matrices V j describe the coupling among the perturbations and depend on the background fields. The vectors S j represent the source terms associated to the point-particle stress-energy tensor. The explicit forms of V j and S j are derived in Appendix B, where we also show that the functions H 2 and H 0 can be completely determined in terms of Ψ p . For the axial sector, the source terms vanish identically and we can write the coupled firstorder equations into a second-order Schrödinger-like equation with an effective potential.

C. Computation of the QNMs
The late-time gravitational signal from a perturbed BH is dominated by a sum of exponentially damped sinusoids, the QNMs, which correspond to the characteristic vibration modes of the spacetime [9,36]. The QNMs are solutions of the sourceless wave equations (16) (i.e., S j = 0), along with proper boundary conditions-purely outgoing waves at infinity and ingoing waves at the horizon. The latter correspond to the following behavior of the wave function at the boundaries where r h is the horizon radius in these coordinates and the "tortoise" coordinate r * is defined by To compute the QNMs, we use a direct-integration method. We construct a square matrix X (four and two dimensional in the polar and axial case, respectively) of which the columns are independent solutions of Eq. (16). This matrix can be constructed by a certain combination of the solution which is regular at the horizon and the one which is regular at infinity (cf. Refs. [37,38] for details). For the polar sector, we note that the boundary conditions defining the QNM eigenvalue problem depend on two parameters, related to the amplitudes of the scalar and gravitational perturbations. Two of the columns of X can be constructed by integrating the equations from the horizon outward. Likewise, two other solutions can be constructed by integrating the equations from infinity inward. In general, the solutions integrated from the horizon are linearly independent from the ones integrated from infinity, unless ω is the QNM frequency. In other words, the QNM frequencies are obtained by imposing where r m is an arbitrary matching point of the order of the horizon radius. The same procedure can be done for the axial modes, with the difference that the boundary conditions apply only to the gravitational amplitude and therefore the problem is technically less involved.
We have computed the lowest QNMs of static BHs using both the full numerical background and the perturbative analytical solution given by Eqs. (A6)-(A8), up to N = 4, i.e. up to O(ζ 6 ). We checked the numerical stability of the QNM frequencies against changes in the values of the numerical horizon r h , numerical infinity and the matching radius r m (typically, we use r m ∼ 4r h ). One of the advantages of the above procedure is that it can be applied to both the perturbative solution and the full numerical one.

D. Plunging particles
To obtain the metric and dilaton perturbations due to a particle plunging into a BH, we need to solve the inhomogeneous system (16). Two common methods used in the literature to solve this problem are a direct integration and a Green's function approach (see, e.g., Ref. [8]).
The Green's function method relies on the fundamental matrix X, as constructed by using the homogeneous solutions as discussed in the previous section. The formal solution of (16) can be written as [39] where β is a constant vector to be determined by imposing the proper boundary conditions. Thus, once the fundamental matrix X of the homogeneous problem is computed, its convolution with the source term S in (20) yields the solution of the inhomogeneous problem. However, in many problems, the source term might converge slowly at the boundaries or even diverge, as in our case in which the source term diverges at the event horizon. These problems can be avoided by performing a suitable nontrivial transformations of the perturbation functions [40,41], or with a careful choice of Green's function [42]. A different scheme which avoids this problem consists in integrating directly the full inhomogeneous system, by imposing the proper boundary conditions for the full solution. First, we expand the perturbation functions near the horizon as where p is a constant chosen such that the above expansion satisfies the inhomogeneous equations near the horizon and Ψ j,H is the (ingoing) solution of the homogeneous equation. The above expansion is solved iteratively near the event horizon for the coefficients ψ j,k up to k = N , and they generically depend on constants (say ψ j,0 ) which are related to the amplitude of the fields at the horizon. The amplitudes at the horizon are then used as shooting parameters; i.e. we chose them such that the numerical solution satisfies the proper boundary conditions also at infinity. For arbitrary amplitudes at the horizon, the numerical solution far from the BH is a combination of ingoing and outgoing waves, i.e., and the required solution is obtained by setting the amplitude of the ingoing waves to zero. In the EDGB case, the amplitudes are related to the gravitational and dilaton perturbations, and therefore the problem is a two-parameter shooting problem for the amplitudes at the horizon. Note that for radial plunging, since the source terms for the axial sector are zero, the particle only induces perturbations in the polar sector.
With the numerical solution at hand, one can compute the gravitational and scalar energy spectra at infinity. This is achieved through the effective stress-energy tensor for the gravitational perturbations (i.e., the Isaacson tensor) and the stress-energy tensor of the dilaton field [43]. As shown in Ref. [44], the Isaacson tensor is the same in GR and in a class of theories including EDGB gravity. The gravitational flux, for a given multipole l, can be written as [30] where K(r) is the polar perturbation given by Eq. (13). 2 The dilaton flux reads where φ 1 (r) is given by the dilaton perturbation in Eq. (12). As shown in Appendix B, since the particle does not have a direct coupling to the dilaton field, dilaton radiation only exists due to the gravitational perturbations. In this sense, the gravitational perturbations work as a "source" for the dilaton radiation. Since in the dilaton equation they are proportional either to the parameter ζ or the background scalar field (and derivative), it is natural to expect that the dilaton radiation scales with ζ 2 . Additionally, for l = 0 and l = 1 the perturbations of EDGB BHs can be written as a single second-order ordinary differential equation, which represents the dilaton perturbation (see Appendix B).

III. STABILITY AND QNMS OF EDGB BLACK HOLES
As previously discussed, the QNMs govern the late-time behavior of any small fluctuation away from axisymmetry of the BH. Since the time dependence is of the form e −iωt , the absence of a QNM frequency with the positive imaginary part in the spectrum implies that all fluctuations decay exponentially with time. Thus, a criterion for linearized mode stability of the spacetime is that all its QNM frequencies have a negative imaginary part [9]. In addition, the late-time dynamics is controlled by the fundamental QNM, i.e., the mode with the smallest imaginary component (equivalently, with the longest decay time).
The Schwarzschild spacetime is stable, and its fundamental QNM frequency, ω S = M ω S R + iM ω S I , for the l = 2 mode reads [9,46,47] M ω S ≈ 0.3737 − i 0.08896 gravitational for the gravitational and scalar fundamental modes, respectively, where M is the BH mass. In EDGB gravity, because the coupling ζ is smaller than unity, we expect that the fundamental QNM frequencies 3 are only slightly different from Eqs. (25) and (26). In other words, we expect EDGB BHs to be stable for sufficiently small coupling. We will now study if these BHs are stable throughout all values of ζ, and also quantify the deviation from the corresponding Schwarzschild value.
We have computed the QNMs of EDGB BHs with two independent codes, both in a small-ζ expansion and using the full numerical background. In the small-coupling limit, where the expressions can be expanded in powers of ζ (see Appendix A), the real and imaginary parts of the QNM frequencies can be written as where ω S R and ω S I are, respectively, the real and imaginary parts of the modes of Schwarzschild BH with the same mass M . The coefficients R j and I j are obtained by fitting the numerical data with the above expression and depend on l and on the nature of the mode.
A. QNMs, light ring, and geodesic correspondence Computing the QNMs of spacetimes known only numerically can be challenging. For spherically symmetric spacetimes, a WKB-type analysis has shown that, in the eikonal (l 1) limit, the QNMs can be obtained using only properties of the light ring (which defines the radius of the photon sphere of the BH). This "null geodesic correspondence" [48][49][50] is useful as it requires only manipulation of background quantities which are easy to obtain, and provides a clear physical insight into the QNMs of BHs: they correspond to waves trapped near the peak of the potential barrier for null particles (i.e., within the photon sphere), slowly leaking out on a time scale given by the geodesic instability time scale.
The geodesic correspondence only works, formally, in the l 1 regime, but can be used even at low l using an appropriate calibration. In fact, this approach has proven to provide reliable results also for low-l modes for a variety of BH spacetimes [49], including Kerr-Newman BHs [8].
By extending the analysis of Refs. [48][49][50], Ref. [8] recently showed that the complex QNMs of a stationary and axisymmetric BH in the eikonal limit can be written as where n is the overtone number, is the orbital frequency at the light ring on the orbital plane, and is the Lyapunov coefficient evaluated at the light-ring location on the equatorial plane. In the above expression, a prime denotes radial derivative, whereas E and L are the (conserved) specific energy and angular momentum of the geodesic, and V is the effective radial potential. The expression (28) is valid for l = m 1 modes. A more involved result for the QNMs of a Kerr BH with generic l 1 and |m| ≤ l is derived in Ref. [50].
The geodesic correspondence has been formally proven for Kerr BHs and for a variety of spacetimes, but it has never been checked for BH solutions in modified gravity. This is particularly interesting in light of the breaking of the isospectrality of the axial and polar QNMs of BHs in EDGB theory, as we discuss in the next section. It is important at this stage to point out that the geodesic approximation must fail to capture some of the features of the full problem. It is well adapted, in principle, to describe the effects of rotation, but cannot take into account (at least not blindly) the presence of extra degrees of freedom like scalars in EDGB. Thus, using this correspondence to extract the QNMs of BHs in modified gravity should be done carefully.
As previously discussed, the background solution and the metric of a spinning BH in EDGB theory is known analytically up to O(χ 5 , ζ 7 ) [19] and numerically for any value of χ and ζ [22,23], where χ is the dimensionless angular momentum parameter, We will use this result in Sec. V, to estimate the modes of spinning EDGB BHs.

B. Axial modes
The axial sector of gravitational perturbations is decoupled from the scalar-field perturbations, and hence is simpler to study. Using the direct-integration procedure described in Sec. II C, we have computed the axial QNMs both in a small-ζ expansion and in the full numerical background. Our results are summarized in Fig. 1 and in Table I. In Fig. 1, we show the behavior of the axial l = 2 fundamental mode as a function of the coupling ζ, normalized by the corresponding Schwarzschild quantity. In most of the range of ζ, the behavior of the modes is smooth and given by a corresponding small deformation of the Schwarzschild QNMs. The only exception occurs close to the critical value of the coupling constant ζ ≈ 0.691, where the QNMs have a very sensitive dependence on ζ. For small values of ζ, say for ζ 0.4, analytically expanded backgrounds [up to O(ζ 6 )] yield QNMs which are in very good agreement with the full numerical solution. This provides a nontrivial check for both our (independent) codes. Figure 1 also shows the result of the geodesic algorithm described in Sec. III A. For the axial modes, which are decoupled from the scalar perturbations, the geodesic predictions are in good agreement with the results of the full numerical solution. Although not shown in Fig. 1, the agreement is better for higher multipoles, as expected (cf . Table I).
Finally, Table I shows the results of the polynomial fit (27) to the axial QNM of EDGB BHs, which is specially accurate for small ζ. By analyzing the perturbation equations, it is easy to show that R 1 = I 1 = 0 in the expansion (27). The full numerical results are available online [47]. In Table I, we also give the coefficients obtained by the geodesic algorithm, which help to quantify the accuracy of the geodesic approximation; for instance at ζ = 0.5 the difference in percentage between the real and the imaginary parts of the (notnormalized) frequency, relative to the numerical result, is, respectively, less than 3% and less than 8% for l = 2, and improves for l > 2. Similar deviations are obtained in the GR limit, ζ = 0.
As explained in Appendix B, similarly to what happens in GR, there are no axial QNMs for l = 0 and l = 1.

C. Polar modes
Unlike the axial sector, the polar gravitational sector of the metric perturbations couples to the scalar-field perturbations. The system of ODEs is more complex and finding the QNM frequencies is therefore more challenging. Even for arbitrarily small values of ζ the QNMs contain two families: (i) gravitational-led modes, which reduce to the gravitational QNMs of Schwarzschild BHs in the ζ → 0 limit, and (ii) scalar-led modes, which reduce to the QNMs of a test scalar field on a Schwarzschild metric when ζ → 0 (see Ref. [51] for a similar situation in another theory). The extent to which each of these modes is excited in actual physical setups is discussed in the next section.
Our results for the quadrupole modes (l = 2) are summarized in Fig. 2 and Table II, where we show the fundamental gravitational-led and scalar-led QNMs. The deviations from the GR case are larger than in the axial case. This is probably due to the extra coupling between gravitational and scalar degrees of freedom. From the fits given in Table II, it is interesting to note that the leading-order correction to the l = 2 gravitational-led mode has an opposite sign compared to the axial mode. In particular, since the geodesic correspondence predicts only one type of modes and the latter are in good agreement with the axial modes, we find that the behavior of the polar modes is not captured by the geodesic correspondence; therefore, we do not plot the geodesics results in the figure. This qualitative difference is expected, because the axial potential resembles the geodesic potential at large l, whereas in the polar sector the coupling to scalar pertur-TABLE I. Numerical value of the coefficients Rj and Ij for the expansions in the small-coupling limit, cf. (27) for the axial QNMs. The geodesic coefficients are computed from the exact analytical solution for small ζ limit, while the QNM frequencies coefficients are obtained through a polynomial fit with the data.  Due to the coupling between the dilaton and gravitational perturbations, there are also nontrivial l = 0, 1 scalar-led modes for EDGB BHs. These reduce to their respective scalar modes in the Schwarzschild spacetime when ζ → 0. The results at finite coupling follow the trend of higher multipoles.

D. Mode stability
From the above results, it is clear that the fundamental QNMs of an EDGB BH change at most by a few percent relative to the Schwarzschild case. As a consequence, these modes are stable for any value of ζ in the domain of existence of static EDGB BHs. We have investigated this issue also for higher multipoles (l ≥ 2) and our numerical search has found no unstable modes in the entire parameter space. This strongly indicates that static EDGB BHs are linearly mode stable, just like Schwarzschild BHs.

IV. RADIAL PLUNGE
In this section, by using the procedures depicted in Sec. II D, we discuss the gravitational and dilaton radiation FIG. 4. Dilaton quadrupolar flux for radial plunges into an EDGB BH with different boosts for ζ = 0.1. As expected, the scalar flux starts to present an exponential suppression for frequencies roughly larger than the dilaton mode. In the inset, we plot the ratio between the fluxes for ζ = 0.1 and ζ = 0.05, which confirms that the scalar flux scales dominantly with ζ 2 .
emitted by a particle plunging radially into an EDGB BH. For practical reasons, this computation was done within the small-ζ approach for the background, and the equations are expanded up to order O(ζ 6 ). As discussed in the previous section, this higher-order perturbative approach gives precise results even for relatively large values of ζ. In Fig. 3 we show the gravitational flux, by considering different initial boosts for the particle (see Appendix A). The deviations from the GR case are very small, of order of ∼ 1%, at least for ζ 0.1. Moreover, we note that all additional terms appearing in the metric perturbation equations (see Appendix B)-sources included-are at least of order ∼ ζ 2 , and therefore the corrections to the flux are proportional to ∼ ζ 2 . Therefore, in the small-coupling limit, the gravitational flux can be written as where dE S g /dω is the corresponding Schwarzschild flux. Although the overall gravitational corrections are small, the dilaton perturbations can also be radiated by the plunging particles, similarly to the case of a neutral particle plunging into a charged BH [8,52,53]. This is due to the coupling between the gravitational and dilaton perturbations, cf. Eq. (B7). We show the dilaton flux for ζ = 0.1 in Fig. 4. Note that the flux displays a cutoff roughly at ω ∼ ω φ R , where ω φ R is the scalarled QNM. Also in this case the dilaton flux scales dominantly as ζ 2 . Additionally, because the source terms in the dilaton field are only due to the gravitational perturbations, the dilaton radiation is also dominantly quadrupolar (see Ref. [52] for a similar setup with perturbations in Reissner-Nordström BHs).
Although the total radiated energy due to the dilaton radiation scales as ζ 2 , it can still be considerably high, depending on the radiating source. For instance, for the source GW150914 the luminosity due to the GWs was dE g /dt ≈ 3.6 × 10 56 erg/s [1]. Therefore, even if the scalar radiation is small compared to the gravitational one, it can still have a considerable value. The implication of a burst of dilaton radiation depends on how the environment (plasma, surrounding stars, etc.) interacts with the dilaton field.

V. CONSTRAINTS ON THE EDGB COUPLING FROM RINGDOWN OBSERVATIONS
The results of the previous section show that plunges of point particles do not excite considerably the scalar-led QNMs. This suggests that such modes might be only mildly excited during the coalescence of two BHs of equal mass and that the main signature of EDGB theory would be a shift of the ringdown frequencies, the latter being governed by the fundamental gravitational-led modes. Thus, the use of ringdown measurements in EDGB theory to estimate the magnitude of the coupling ζ would rely only on the deviation of the gravitational-led modes from their GR counterpart. These deviations are parametrized in terms of the fit (27) of which the coefficients are given in Tables I and II for the l = 2, 3 fundamental axial and polar modes, respectively.
The results of the previous sections refer to nonspinning BHs, whereas the end product of the coalescence is a spinning compact object. Thus, ringdown tests require the knowledge of the first dominant modes of a spinning BH as a function of the spin χ and of the coupling constant of the modified theory of gravity. Computing the QNMs of generic spinning BHs in modified gravity is a very challenging task 4 , which has witnessed some developments only recently (cf. Ref. [37] for an overview). Nonetheless, most of the results are obtained within a perturbative expansion valid for χ 1 and quickly become intractable at the higher perturbative order. The latter is required to extrapolate the perturbative result up to χ ≈ 0.7, which is roughly the spin of the final BH measured in the two coalescence events detected by aLIGO to date [1,54].
To overcome this limitation and estimate how rotation affects the QNMs of an EDGB BH, we rely on the geodesic correspondence and on our knowledge of the metric of a spin-ning BH in EDGB theory. As previously mentioned, the latter is known analytically up to O(χ 5 , ζ 7 ) [19] and numerically for any value of χ and ζ [22,23]. In the previous section we have checked that the geodesic correspondence works reasonably well for axial modes, so in this section we will focus on the latter. However, we may argue that the order of magnitude of our estimates should be correct also for the more relevant polar modes.
The detection of two ringdown modes is necessary to estimate the mass, spin and coupling ζ. Here we adopt the same Fisher-matrix technique presented in Ref. [8]. Let us consider first the axial modes, which are, generically, excited during the merger phase. The measurement of the two most dominant modes (which we take to be l = m = 2, 3, excited with a relative amplitude of 3 : 1 [11]) gives us sufficient information to extract the mass and spin of the BH, as well as the coupling parameter of the theory. To do this, we use the geodesic correspondence, Eq. (34). Assuming that both modes are excited to detectable amplitude (this can be quantified using the methods of Ref. [11]), then a Fisher-matrix computation allows us to estimate the uncertainties in the parameters determining the ringdown [8]: M, χ, ζ. This then allows us to constrain the magnitude of the coupling parameter ζ, where ρ is the signal-to-noise ratio in the ringdown waveform and the O(χ 4 ) and O(χ 5 ) terms are negligible. The numerator in the expression above is a (mildly) decreasing function of the spin and ranges from ≈ 4 to ≈ 3.5 in the region 0 ≤ χ ≈ 0.8. This analysis can be extended to polar modes. In fact, one might even argue that corrections to polar modes are higher, since the polar QNMs of nonrotating BHs are more affected than the axial, cf. Tables I and II. However, there is little evidence that they follow the geodesic correspondence, but we expect that the order of magnitude change in the modes remains the same. Therefore, it is reasonable to believe that, even when polar modes are excited, a mode analysis of the ringing signal can set constraints on ζ of the order of those given by Eq. (35).
As a reference, the final BH spin of GW150914 was χ ≈ 0.67 [1] and ρ ≈ 7.7 in the ringdown part [6]. This [assuming, as a first approximation, that Eq. (35) also holds for polar modes] yields roughly ζ 1.3, which is weaker than the theoretical bound [18] for the finite-ζ solution, ζ 0.691, and therefore meaningless. This constraint is nevertheless comparable with that derived from the orbital decay rate of lowmass x-ray binaries [15] and it is slightly larger than the projected bound achievable in the near future from the measurements of quasiperiodic oscillations in the spectrum of accreting BHs [55], although the latter might be affected by astrophysical systematics. Our estimate suggests that ρ 25 in the ringdown waveform is needed to obtain an upper limit which is more stringent than the theoretical bound using GW ringdown detections. Of course, in order to set these upper limits, an accurate computation of both polar and axial modes of rotating EDGB BHs will be needed.

VI. DISCUSSION AND CONCLUSIONS
EDGB gravity is a simple and viable higher-curvature correction to GR which predicts BH solutions with scalar charge. Since strong-curvature corrections are suppressed at large distance, it is natural to expect that the most stringent constraints on this theory come from the strong-curvature, highly dynamical regime as the one involved in a BH coalescence. Furthermore, compact stars in this theory possess only a very small scalar charge [56,57] and therefore EDGB gravity evades the stringent constraints on the dipole radiation coming from current binary-pulsar systems [5].
The estimate (35) translates to the following upper bound on the dimensionful EDGB coupling 5 , where the prefactor changes by less than 10% depending on the final BH spin. This result is in agreement with the simple estimates derived in Ref. [58]. As a consequence, our analysis also confirms that in most cases modified-gravity effects can be distinguished from environmental effects [58]. Future GW detectors will greatly increase the signal-tonoise ratio, a large value of which is necessary to perform ringdown tests of the Kerr metric [59]. The signal-to-noise ratio of a ringdown waveform scales approximately (among its dependence on other quantities not shown here) as ρ ∼ M 3/2 /S n (f ) 1/2 [10], where M is the final BH mass and S n (f ) is the detector noise power spectral density at a given frequency f . The best sensitivity of the future Voyager [60] and Einstein Telescope [61] detectors will be, respectively, roughly a factor of 10 and a factor of 100 better than in the first aLIGO observing run at the same optimal frequency f ∼ 10 2 Hz [59]. Thus, the Einstein Telescope with an optimal design can achieve a signal-to-noise ratio of roughly ρ ≈ 100 for the ringdown signal of a GW150914-like event. From Eqs. (36) and (35), this would translate into the bound km and ζ 0.4. As expected, lighter BHs would provide a significantly more stringent constraint on α, although their ringdown frequency might not fall into the optimal frequency range for ground-based detectors. Due to the small exponent of ρ in Eq. (36), even an increase of ρ of 1 order of magnitude will not provide a significantly more stringent constraint on the EDGB coupling. A stronger constraint may be set if future observations detect a light BH with a very large signal-to-noise ratio.
Given this scenario, electromagnetic observations of accreting BHs (like the one discussed in Ref. [55]) might provide more stringent constraints in the future, although the latter are affected by astrophysical systematics that are absent in the ringdown case.
Our estimates in the case of spinning BHs rely on the geodesic analogy for QNMs, which we verified only for axial modes in the static case and for Kerr BHs with any spin [8]. It would be interesting to compute the modes of slowly rotating EDGB BHs (e.g. by adapting the methods discussed in Ref. [37]) and to check the geodesic approximation in the spinning case. This computation will be required to place precise constraints on the EDGB coupling through future detections of BH ringing with high signal-to-noise ratio.
Another interesting extension of our work concerns the scalar waves emitted during the coalescence. Although the luminosity in scalar waves is significant, this radiation may 5 Since one of the parameters of our Fisher-matrix analysis is ζ = α/M 2 , propagation of errors implies a relative uncertainty δα/α = δζ/ζ + 2δM/M . However, in the large-ρ limit the term δM/M is negligible because it scales as 1/ρ, compared to the 1/ √ ρ behavior of the error on ζ [cf. Eq. (35)]. Therefore, in this limit δα δζM 2 .
be possibly detected only if the dilaton is coupled to matter. Such coupling is presumably small and would not give rise to any effects in the detectors. Nonetheless, if the dilaton-matter coupling is non-negligible, the scalar radiation might be investigated through the same techniques developed to study the scalar emission in scalar-tensor theories, e.g. by using a network of ground-based detectors [62]. the metric functions and the scalar field behave as The coefficients a j , b j and φ j are constants, but they are not free parameters. In fact they depend only on the coupling parameters of the theory and horizon radius of the BH (related to the total mass M ) through some complicated algebraic relations [14].
In order to study the quasinormal modes of the full geometry, we build numerically the solutions of the background metric. We do so by integrating the system of ordinary differential equations for A(r), B(r) and φ 0 (r), from the horizon up to infinity. We perform this integration in a compactified coordinate x = 1 − r h /r, with x ∈ [0, 1]. With a suitable reparametrization of the functions we can impose as boundary conditions the correct behavior at infinity (x = 1) and at the horizon (x = 0). The integration is performed with the package COLSYS [63], and typically background solutions are generated with 1000-10000 points and required relative precisions of the functions smaller than 10 −6 .
Another approach, which allows us to obtain analytical expressions for the background metric and dilaton field is the small-coupling limit. In the small-coupling limit, in which ζ ≡ α/M 2 1, the equations can be greatly simplified, and the background field can be computed analytically. Considering the expansion for the background fields and expanding the background equations in ζ, one has that each order j +2 for the metric and j +1 for the scalar field can be solved consistently, imposing regularity at the horizon r = r h and that the scalar field goes to zero at infinity. Typically, we consider solutions up to N = 4, verifying that the results converge in the small-ζ region.
To obtain the system as (16) we do the following: we use Eqs. (B3)-(B5) to eliminate H 0 , H 0 , and H 2 in terms of K, H 1 , φ 1 , and their first derivatives, as well as source terms, substituting them in Eqs. (B1), (B2) and (B7). In this way we can write which is the extended form of Eq. (16) for j = p. Obviously, V 3k = 0 for k = 4, V 34 = −1, and S 3 = 0. Due to the complexity of the components of V p , we shall not show them here, but they can be seen in the supplementary Mathematica notebook [47]. The boundary conditions for the perturbations can be found with the aid of the expansions of the background metric and dilaton field. Note that, for QNMs, due to the natural divergence at spatial infinity, one must consider a highorder expansion of the perturbations at infinity.
For l = 0, 1, simpler gauge choices can be chosen [30]. First, for l = 1 we notice that we can pick a gauge in which K vanishes identically. Therefore, we can use the first two equations in (B8) to eliminate H 1 and H 1 in favor of φ 1 , φ 1 . For l = 0, we can choose a gauge in which both K and H 1 vanish and again the equations are reduced to a second-order equation for the dilaton perturbation. Note that this is possible because some harmonics in the expansion are identically zero in the l = 0 and l = 1 cases.

Axial sector
The equations for the axial sector are much simpler. The fundamental equations can be written as where h 0 , h 1 are the axial parity perturbation functions defined in Sec. II B. The above differential equations are already in the form of Eq. (16). As mentioned in the main text, the particle does not induce perturbation in the axial sector, and hence there are no source terms in the above equations.
For l = 0, the tensorial harmonics multiplying the functions h 0 and h 1 vanish identically, and therefore the axial perturbation vanish identically.
For l = 1, the analysis of the perturbations follows in a very similar manner as the one in GR [30,35]. One of the harmonics vanishes identically, and one is left with only two equations, namely (B10) and (B11). We can exploit the gauge freedom to set either h 1 or h 0 to zero. The remaining equations, which are asymptotically GR, have a nonradiative behavior at infinity and contribute only to give an infinitesimal angular momentum to the BH.

Source terms
As mentioned in the main text, the particle stress-energy tensor can be expanded in spherical harmonics. The procedure to obtain the coefficients is outlined in Refs. [30,35,64]. In the frequency domain, the source functions for a particle falling radially into the BH are given by A lm = l + 1 2 where γ(≥ 1) is the specific energy of the particle (boost parameter) and T is the time trajectory of the particle, as a function of the radial coordinate. The function T can be obtained by solving the differential equation We can solve this equation together with the perturbed equations, imposing that T = 0 at the numerical horizon, without loss of generality. We note that in the small-ζ approximation for the metric Eq. (B15) can be solved analytically, expanding the right-hand side in powers of ζ.