Measuring the ringdown scalar polarization of gravitational waves in Einstein scalar Gauss-Bonnet gravity

We model the scalar waves produced during the ringdown stage of binary black hole coalescence in Einstein scalar Gauss-Bonnet (EsGB) gravity, using numerical relativity simulations of the theory in the decoupling limit. Through a conformal coupling of the scalar field to the metric in the matter-field action, we show that the gravitational waves in this theory can have a scalar polarization. We model the scalar quasi-normal modes of the ringdown signal in EsGB gravity, and quantify the extent to which current and future gravitational wave detectors could observe the spectrum of scalar radiation emitted during the ringdown phase of binary black hole coalescence. We find that within the limits of the theory's coupling parameters set by current theoretical and observational constraints, the scalar ringdown signal from black hole remnants in the $10^1 - 10^3 \, M_{\odot}$ mass range is expected to be well below the detectability threshold with the current network of gravitational-wave detectors (LIGO-Virgo-KAGRA), but is potentially measurable with next-generation detectors such as the Einstein Telescope.


I. INTRODUCTION
General relativity (GR) has passed all observational tests to date [1,2] and, as of 2015, observations of gravitational wave signals in LIGO [3] and Virgo [4] data have been consistent with its predictions [5][6][7][8][9][10]. Nevertheless, there has been an increasing amount of work that tests the predictions of non-GR theories of gravity (modified theories of gravity) with gravitational wave measurements of binary black hole (BH) and neutron star (NS) mergers [11][12][13][14][15]. These model-dependent tests of GR have allowed for remarkably strong constraints to be placed on a wide array of modified gravity theories.
Here we focus on the gravitational wave implications for black hole binaries in Einstein-scalar-Gauss-Bonnet (EsGB) gravity. EsGB gravity can be motivated from effective-field-theory-styled arguments [16,17], and has been shown to appear in the low-energy limits of heterotic string theories [18][19][20]. What makes EsGB interesting from a phenomenological point of view is that black holes in the theory can have scalar hair [21,22] (for a review, see [23]). Because of this, binary black hole (BBH) systems can emit scalar radiation [24]. From a theoretical point of view, EsGB gravity is a Horndeski theory [25] and, as such, has second order equations of motion, circumventing the Ostrogradsky instability. In addition, the theory has been shown to have a well-posed initial value problem in the weak coupling limit [17] which allows for numerical relativity simulations of the full theory [26][27][28][29]. 1 In this paper, we focus on the leading order effects of EsGB on the gravitational waveform emitted during the ringdown stage of BBH mergers. It is expected that the strongest constraints on most modified theories of gravity will come from the inspiral, due to the cumulative dephasing of the signal from the GR predicted waveform, over a large number of observed cycles [33,34]. Nevertheless, correctly modeling the ringdown is crucial to understanding the properties of the remnant black hole formed from binary coalescence, and to performing consistency tests with theory-specific estimates of the final black hole properties inferred from the inspiral and merger [35]. While the dephasing of the inspiral can be a common feature of many different theories of modified gravity and the theory-specific details can be difficult to model at high accuracy (see [13,36] for perturbative approaches in EsGB), a spectroscopic analysis of the ringdown signal may reveal characteristic features that will help narrow down the origin of an observed deviation from GR.
Working in the decoupling limit of EsGB gravity, we extract first-order tensorial and scalar waveforms and focus on the ringdown portion of the signal. We study the quasinormal mode (QNM) spectrum of the resulting scalar and tensor waveforms and numerically compute the amplitudes of the excited modes. Unlike most previous studies on EsGB gravity, we study the impact of the radiated scalar field on the scalar polarization measured by the gravitational wave detectors, through the presence of a nonminimal scalar field coupling to the metric tensor in the matter field action (without such a coupling, gravitational waves would not have a scalar polarization component in EsGB gravity [37]). While there has been work on model-independent tests of the polarization of gravitational waveforms [38][39][40][41], we present the first model-dependent analysis of the feasibility of measuring a nontensor polarization from black hole binaries with current and future gravitational wave detectors (we note though that Ref. [42] considers model-dependent tests with stochastic gravitational wave backgrounds). Ultimately, we find the scalar polarization-given current constraints on the EsGB gravity and weak-field tests of general relativity-is unobservable with current gravitational wave detectors, although there is a higher chance that it could be measured with next-generation detectors.
The paper is organized as follows. In Sec. II A we present the theoretical framework of EsGB gravity in the Einstein frame. In Sec. II B, we review how to extract the polarization content of gravitational waves (GWs), and show how the conformal scalar field coupling can give rise to a scalar polarization in gravitational waves. Next, we review and justify the application of the decoupling limit for solving the equations of motion of EsGB gravity in Sec. II C. For reference, we present a summary of current observational and theoretical constraints on the couplings of the theory we consider in Sec. II D. In Secs. III A-III C we overview our numerical set-up and the diagnostics we used to extract physical quantities from our simulations. In Sec. III D, we present the scalar waveforms measured from our numerical relativity (NR) simulations. We analyse the QNM spectrum of our resulting scalar waveforms in Sec. IV and compute their signal-to-noise ratio (SNR) in Sec. V A to assess the measurability of the scalar signal. We then move on to a full Bayesian data analysis of our simulated signals by performing NR injections of the scalar waveform into simulated current and future-generation detector noise using a tailored ringdown pipeline in Sec. VI. We finish with concluding remarks in Sec. VII.
Our notation and conventions are as follows. The metric signature is (− + ++) and we use M to denote the Arnowitt-Deser-Misner (ADM) mass. We use lower case Greek indices to index spacetime components (indexed 0, 1, 2, 3, with 0 being the timelike index) and lower case Latin indices to index spatial components. The Riemann tensor is R α µβν = ∂ β Γ α µν − · · · .

II. THEORETICAL BACKGROUND
A. Formulation of four-derivative scalar-tensor (4∂ST) gravity The general action for Einstein-scalar-Gauss-Bonnet gravity is where V (φ) is the scalar field potential,ᾱ(φ), β(φ) are arbitrary coupling functions of the scalar field, X ≡ − 1 2 (∇φ) 2 , S M is the action functional for the matter fields (which we denote schematically with Ψ), A 2 (φ) is the conformal coupling of matter to the metric, and R GB is the Gauss-Bonnet scalar If A ̸ = 0, the scalar field couples directly to matter. This coupling generally leads to violations of the weakequivalence principle (barring the presence of screening mechanisms [43]), as the effective gravitational field felt by matter fields will depend on the value of the scalar field φ. Weak field tests of gravity have strongly constrained violations of the weak-equivalence principle, which can be translated to constraints on the strength of the coupling A [1]. In geometric units, the coupling parameter α GB has dimensions length squared, L 2 . It is then natural to expect strong constraints on the scalar Gauss-Bonnet coupling will come from regions of high spacetime curvature [33,34,44]. In this work we examine a subset of theories for which V =ᾱ = 0. We assume the beyond-GR corrections give only relatively small deviations to the background GR solution. Therefore, we impose that |φ|, |α GB |, |β| ≪ 1 and to leading order we have where we have Taylor expanded β and A in φ, set A(0) = 1, A ′ (0) ≡ a 0 /2, and β 0 ≡ α GB β ′ (0). We have discarded the term β(0) as the Gauss-Bonnet scalar is a total derivative in four spacetime dimensions, so it does not contribute to the equations of motion.

B. Polarization content
In alternative metric theories of gravity, there may exist up to six polarizations associated with additional degrees of freedom: tensorial (plus and cross), vector (x and y) and scalar (breathing and longitudinal) [45]. Because of their symmetries, the breathing and longitudinal modes are fully degenerate to networks of quadrupolar antennas [46], so we will refer to them jointly as the scalar polarization. The scalar polarization state is generally a mixture of the transverse breathing and longitudinal polarizations and is excited by a massive scalar field. In the massless limit, the longitudinal polarization disappears, while the breathing one persists [47]. Since the scalar field we consider is massless, we only study the breathing polarization in this work.
The six polarization states are encoded in six "electric" components of the Riemann tensor, R i0j0 (i, j spatial), which can be written in terms of the following Newman-Penrose scalars [45] In Fig. 1 we illustrate the form of the different polarization states. Perturbing the metric and scalar field, the metricĝ µν = (1 + a 0 φ)g µν that couples to the matter fields (such as a gravitational wave detector) in the wave zone isĝ where η µν is the Minkowski metric. The tensorial perturbation h µν gives the plus and cross polarizations, while the coupling constant a 0 , along with a nonzero value of φ, gives a new scalar breathing-mode polarization to gravitational waves 2 [42,48].

C. Perturbative equations of motion
Unless otherwise stated, in this section we work in geometric units G = c = 1. Our discussion roughly follows that of [49,50]. We start with the full equations of motion derived from action (3) where 2 We note that if a 0 = 0, then gravitational waves have no scalar polarization in EsGB gravity [37].
We expand the metric and the scalar field as where ϵ is a book-keeping parameter for our perturbative expansion. By considering the equations of motion order by order in ϵ, one finds that up to O(ϵ) [50] φ (0) = 0, We note that the zeroth order equations of motion are the Einstein field equations minimally coupled to a massless scalar field. The choice of φ (0) = 0 in the zeroth order solution (15) is motivated by the fact that asymptotically flat black holes cannot carry scalar hair [51], which if initially present would be radiated away at late times, and that a cosmological value of ϕ for EsGB gravity has been constrained by measurements of the speed of gravitational waves [12,52]. Using Eqs. (15) and (16), the equations of motion to O(ϵ) reduce to We note that there is no backreaction of the scalar field on the background geometry, however at O(ϵ 2 ) the metric will be corrected by the backreaction of the scalar field φ (0) on the spacetime [50,53,54]. We reintroduce dimensions by choosing the characteristic length scale of the system to be L, and introduce a dimensionless coordinatex = Lx. By redefining x →x, we acquire a factor of L −2 in front of the wave operator and a factor of L −4 in front of the Gauss-Bonnet scalar, leading toˆ We then may rescale the scalar field so thatˆ This implies that a solution with a given β 0 /L 2 parametrizes a family of solutions with different L and β 0 . Further, Eq. (20) naturally introduces the choice for the dimensionless parameter ϵ, which we use to parametrize the weak coupling of our solutions. Since we are interested in binary black hole systems and the corrections become larger for smaller masses, we choose L = m 1 , The classification of the six polarizations encoded in the four Newman-Penrose scalars. The shapes represent displacement that each mode induces on a sphere of test particles. The wave is assumed to propagate in the z-direction indicated by the icons/arrows in the right corners.
which is the mass of the smaller black hole in our binary black hole configuration. We then define By combining Eqs. (20) and (8), we note that the scalar GW amplitude is then controlled by the couplings of our theory via a 0 β 0 /m 2 1 . The dimensionless EsGB parameter ϵ plays an important role in EsGB gravity. Depending on the coupling β 0 chosen in the theory, a stationary BH solution of mass m 1 cannot exist above a certain threshold max β 0 /m 2 1 = ϵ thr > ϵ. Such threshold is controlled by the regularity of the scalar field on the horizon and the finite radius singularity being hidden behind it for sufficiently small β 0 . If the coupling becomes too large, a naked singularity emerges from within the horizon in solutions to the full theory. For instance, in the spherical-symmetric case of shift-symmetric EsGB, ϵ thr ∼ 0.3 [22]. Hyperbolicity constraints on the theory suggest a somewhat more stringent bound of ϵ thr ∼ 0.2 [55]. In Fig. 2, we outline schematically the range of validity of our EsGB theory in the parameter space of a 0 × β 0 /m 2 1 , taking into consideration the observational constraints on a 0 from Table I and theoretical constraints on β 0 /m 2 1 . The contour lines therefore represent possible scalar GW amplitude values.

D. Current observational constraints
Here we discuss current experimental tests and observational bounds on the parameters a 0 and β 0 . The EsGB coupling β 0 has been most strongly constrained through the analysis of inspiral and merger of binary systems of compact objects using GW data released by the LIGO-Virgo-KAGRA (LVK) collaboration [56]. Weakfield tests of GR also place a constraint on the coupling, but do not constrain it as significantly [57], as the Gauss-Bonnet coupling does not lead to much scalarization of stellar solutions [24]. This being said, as we are considering a 0 ̸ = 0, stars can scalarize in the theory we consider. In this work, we use the first constraint on a 0 , obtained from the Viking relativity experiment on verification of signal retardation by solar gravity [58] 3 , a 2 0 ≲ 10 −3 . Table I references the studies constraining the EsGB coupling 4 β 0 and the conformal coupling a 0 . We additionally summarize theoretical constraints on the dimensionless coupling β 0 /m 2 1 as discussed in Sec. II C.

A. Numerical relativity code and evolution equations
Here we describe our numerical set-up to solve Eqs. (17) and (18). For our simulations, we use a modified version of GRChombo [70][71][72], which is a publicly available finite difference numerical relativity code built on the Chombo [73] adaptive mesh refinement libraries. The 3 We note however that the most stringent constraint of a 2 0 ≲ 10 −5 is known to be from the Cassini measurements of the Shapiro time delay [1,59]. Here we use a weaker constraint from the Viking experiment to see how well in principle the scalar polarization could be measured by present GW detectors. 4 If one assumes that both objects in GW190814 are BHs, then an even more stringent constraint on the EsGB coupling can be found [60]. However, in this work we take a conservative approach, by quoting the constraints from GW190814 that assume the secondary object is a neutron star.    metric and the scalar field are evolved with the method of lines using fourth order finite difference stencils and Runge Kutta time integration. We use the covariant and conformal Z4 (CCZ4) formulation [74] with the moving puncture gauge [75,76] for evolving the Einstein equations (17). The evolution equations for the Einstein field equations in the CCZ4 formulation can be found in Sec. III F of [77], where we replace κ 1 → κ 1 /α, in order to stably evolve BHs, and choose the constraint damping parameters to be κ 1 = 0.1, κ 2 = 0 and κ 3 = 1.
The evolution equations for the scalar field (18) read , α is the lapse function, β i is the shift vector,γ ij ≡ χγ ij is the conformally rescaled metric constituent of the conformal factor χ = det(γ ij ) −1/3 and physical spatial metric γ ij . We rewrite the Gauss-Bonnet scalar in the gravito-electric, E ij , and gravito-magnetic, B ij , counterparts [50] where in CCZ4 variables they take the form of [72] Here R ij is the 3D Ricci tensor, K ij is the extrinsic curvature, Θ = −n µ Z µ is the projection of the CCZ4 vector Z µ onto the timelike unit normal n µ , [..] TF denotes the tracefree part and D i is the covariant derivative compatible with the physical metric γ ij . In the case of Θ = 0, we recover the Baumgarte-Shapiro-Shibata-Nakamura-Oohara-Kojima formulation [78][79][80].

B. Grid setup and initial data
For the results presented in this paper we set up a computational domain of size 1024M . We use a grid spacing of ∆x = 2M in the outermost refinement level with eight additional refinement levels and a Courant factor of 1/4. In this setup we roughly have 64 points covering the smallest BH in the most asymmetric run (mass ratio 1:2). For boundary conditions, we use Sommerfeld boundary conditions and take advantage of the bitant symmetry to evolve only half of the grid. To validate our results, we estimated the discretization error for the amplitudes of Ψ 4 and φ to be ≲ 2% and for phase of Ψ 4 to be ≲ 0.15 radians. We further estimate the finite radius extraction error to be ≲ 2% and thus a total error budget of ≲ 4%. See Figs. 14-16 of Appendix A 1 for a quantitative illustration of the convergence of the code.
For our initial data, we consider two nonspinning black holes with initial masses m 1 and m 2 (we assume m 2 > m 1 ), initial linear momenta P 1 and P 2 and initial separation D. We prescribe quasicircular initial data for a black hole binary by using puncture initial data [81] of Bowen-York [82] type provided by the spectral initial data solver TwoPunctures [83,84].
We have investigated the effects of different initial separations on the eccentricity and found that D = 8M (resulting in four orbits) gave us an eccentricity estimator of no more than 0.012 for all of the configurations considered here (we used Eq. (20) of [85] to calculate the eccentricity).
The initial data for the scalar field is chosen to be φ = 0 and ∂ t φ = 0. The authors in [50] found that using an initially vanishing scalar or superposing two hairy solutions gave identical results after sufficiently long evolution (see Fig.4 of [50]). Moreover, this choice of initial data solves the exact constraint equations of EsGB gravity [26,29]. We therefore do not investigate any other forms of initial data for the scalar field. We further pick three configurations for black holes of different mass ratios q = m 1 /m 2 , (q ≤ 1) summarized in Table II.

C. Extraction of physical quantities
To order O(ϵ), the metric sector of our spacetime is determined by the Einstein equations of GR (see Eq. (17)), so that the gravitational waveform will be no different from GR. To calculate gravitational radiation we use the Newman-Penrose formalism [86,87], and determine the outgoing gravitational radiation by computing complex scalar Ψ 4 . By interpolating the real and imaginary parts of Ψ 4 onto a sphere of fixed radius R ext = 110M , we decompose them in terms of spin-weighted spherical harmonics Y s lm with s = −2 [88,89], where dΩ = sin θdθdϕ. By similar decomposition into spherical harmonics Y lm and integrating over the sphere, we construct the scalar waveform In what follows, we use the subscript "rad" to denote quantities calculated from the start of our simulation, which includes spurious "junk" radiation. Otherwise, we consider quantities from t = R ext + 50M onwards. The mass and spin of the final BH are estimated from balance arguments [90] where E GW rad denotes the total radiated energy of GWs, J z rad is the radiated angular momentum in the z direction (by symmetry) and L = P y D. In the above balance arguments the radiated energy and angular momentum  II: Two puncture initial data used for three binary BH configurations considered in this paper. M fin and j fin denote the final mass and final spin of the remnant BHs respectively, which were calculated from balance arguments given in Eqs. (29) and (30). We estimate 0.01% error in the final mass and 0.05% error in the final spin when compared to their Richardson-extrapolated values. In Appendix A 2 we present convergence plots for the radiated energy and angular momentum.
are determined through [91,92] where e r is the unit radial vector of a sphere and J is the angular momentum operator for spin weight s = −2 Apart from tensor energy, we expect some of the scalar field energy to be emitted too. Since the scalar field is of order φ ∼ ϵ, the radiated energy must be of order E φ ∼ ϵ 2 . We note that if we were to continue the order reduction scheme for our solution to O(ϵ 2 ), there would be an additional radiated tensor energy sourced by O(ϵ 2 ) corrections. At least in post-Newtonian (PN) theory this term is suppressed compared to the scalar radiation, so we will ignore it as it is done in [49,50]. The energy flux of the scalar field through a two-sphere of radius R is given via where n a is the timelike unit normal and T ab representing the right hand side of the modified Einstein's equations in Eq. (9). Assuming that a spherical surface is located in the radiation zone of an asymptotically flat spacetime and φ behaves as outgoing radiation in the radiation zone, to leading order we find that T (φ) ab = 1 2 T ab . Then, to leading order Eq. (35) reads where the dot denotes time differentiation. We therefore integrate Eq. (36) to find the total radiated scalar energy to leading order.

D. Numerical results: Scalar waveforms
In this section we present our numerical results. The initial data for the simulations are listed in Table II. We find good agreement with numerical results of [67], which made use of the same decoupling approximation as we do. We show scalar waveforms for l ≥ 1 modes in Fig. 3 and for l = 0 mode in Fig. 4, where we align the waveforms so that the amplitude peaks of the scalar field coincide att = 0 and choose β 0 /m 2 1 for convenience. For equal-mass binaries we have only even l-modes present in the waveform, as equal-mass nonspinning binaries are symmetric under parity transformations (⃗ r → −⃗ r). Since Y m l=odd (ϑ, φ) is antisymmetric and the EsGB coupling is invariant under parity, there must be no emission from odd l-modes for equal-mass binaries. On the other hand, in the unequal-mass case, odd l-modes are present and the dipole emission dominates all other higher modes of the scalar waveform for EsGB gravity (it enters at −1PN) order) [24,36,93]. In Fig. 4, we see that the l = 0 mode looks unlike the other modes and is roughly constant before merger time (this is consistent with leading order PN theory [24,67]). After merger it settles to almost the same value for all the mass ratios, since the remnant BHs formed have quite similar final masses and spins, and thus the same amount of remnant scalar hair. Overall, from Figs. 3 and 4, at fixed β 0 /m 2 1 , the scalar field radiation becomes stronger for more extreme unequal-mass binaries, especially in the inspiral and merger portions of the signal. End state BHs with smaller spin (i.e. BHs formed from more extreme mass ratios) have faster damping times than BHs with higher spin. As a result, the scalar waves in the ringdown are more pronounced at late times for black holes produced from more symmetric-mass-ratio (q ∼ 1) progenitors. To illustrate this point, we align the (22) modes in the amplitude peak of the scalar field for different mass ratios in Fig. 5, and zoom-in on the ringdown, where the amplitude for more symmetric-mass-ratio binaries becomes larger than for more asymmetric systems later in time.
Finally, the scalar energy flux, as computed from For unequal mass cases, the scalar field amplitude is more pronounced. The waveforms have been aligned so that the amplitude peaks of the scalar field coincide att = 0, and we set β 0 /m 2 1 = 1 for convenience (recall that since we are using an order-reduction scheme, there is no restriction on the value of β 0 /m 2 1 in our evolution). In our analysis, we choose a physically allowed value of β 0 /m 2 and rescale the scalar waveforms using Eq. (20). Eq. (36), is shown in Fig. 6. We see that it is dominated by the contribution from l = 0 mode, which is larger and q = 1/2 waveforms become smaller in amplitude than the q = 1 waveform.
for more unequal-mass-ratio binaries, as expected from PN theory [24,67]. This is unlike the tensorial gravitational energy flux, where the strongest flux comes from the equal-mass binary. This can be explained by the fact that to leading order the total radiated gravitational energy E GW rad ∼ η 2 , where η = q/(1 + q) 2 is the symmetric mass ratio [94]. We note that the scalar energy flux measured here using Eq. (36) is the second-order effect and so is significantly suppressed-at merger, it accounts for much less than 1% of the gravitational energy flux.

IV. MODELING THE SCALAR RINGDOWN WAVEFORM
When exploring the GW phenomenology of astrophysical systems in modified gravity, the ultimate aim is to obtain an accurate waveform model for the signal (here the scalar ringdown) and to be in a position to observe the effects in question and infer information on the underlying theory from the observed GW data. In this section we model the scalar ringdown based on the NR data of Sec. III D. The gravitational waves emitted by the remnant BH formed by a merger, known as the ringdown signal, can be approximated as a superposition of damped sinusoids [95]. In GR, each mode of the tensorial ringdown signals can then be written as where t 0 is the start time of the ringdown, J = (lmn) is the mode number, H J and p J are real amplitude and phase of mode J at t = t 0 , (which depend on the binary configuration and dynamics near merger), and ω J are the complex frequencies of the N most dominant QNMs. In a similar fashion, the scalar ringdown signal can be written as, where A J now denote scalar real amplitudes. For a mode J corresponding to a given set of QNM indices (lmn), we write the real and imaginary parts as The integers (lm) describe the angular properties of the emission, while n is the overtone number [96]. Modes with n = 0 are referred to as the fundamental modes, while with n ≥ 1 as overtones, which (at least in GR) have lower frequencies and greater damping times than their corresponding fundamental modes. We have denoted the total number of quasinormal modes with which we model the ringdown signal by N . The no-hair theorem of GR states that the damping times, τ lmn , and frequencies, ω lmn , for gravitational wave perturbations around astrophysical (noncharged) BHs in GR, are uniquely determined by the mass and spin of the final black hole. The excitation coefficients also depend on the initial conditions of the perturbation [96], in our case the progenitor parameters. For modified theories of gravity, the QNM spectrum of BH spacetimes generally depends on the theory in question, and so observations of the ringdown signal can provide a way to discriminate between GR and possible modified theories of gravity. Unlike GR [97], for alternative theories of gravity there is typically no set of separable master equations for linearized perturbations of black holes. Instead, the calculation of QNMs has relied on either the assumption of spherical symmetry [98,99], the slow spin approxi- mation [100,101], 5 or on numerical fits of time-domain waveforms computed from NR simulations [50,103]. In this work we focus on the last method by performing numerical fits to the scalar and tensor waveforms.
In our analysis, we take a similar approach of modeling the scalar ringdown as a superposition of damped sinusoids given by Eq. (38), except our QNM frequencies no longer correspond to just scalar QNMs of GR. We note that the operator acting on φ in its equation of motion consists of the background GR metric in the decoupling approximation we use (see Eq. (18)). Therefore, the homogeneous solution of the beyond-GR scalar perturbation equation around a BH spacetime should contain the same scalar QNMs as the corresponding Kerr BH in GR, ω Kerr,scalar lmn . However, our equation is nonhomogeneous as we also have the Gauss-Bonnet scalar evaluated on the background GR metric, acting as a source term. Particular solutions to the scalar equation introduce additional frequencies into the ringdown QNM spectrum [50,103]. The spectrum of these "driven" modes by the source term for l ≥ 2 has been shown to coincide with the quasinormal modes for tensor perturbations for Kerr, ω Kerr,grav lmn ; while for l = 1 the driven mode has been fitted empirically [50].
Can our scalar ringdown waveforms be modeled by the gravitational and/or scalar QNMs of a Kerr black hole in GR as above findings suggest? To address this question we perform one-mode (N = 1) and two-mode (N = 2) 5 However, there has been recent work on extending the Teukolsky formalism to modified gravity theories to compute their quasinormal modes for fast-rotating black hole solutions [102].
fits, whilst keeping the QNMs fixed to their values in GR. These fixed values employed in our analysis can be found in Table VI  and ω Kerr,grav lmn simultaneously. We note that for the l = 1 ringdown modes fitted with the N = 2 mode fit we take a conservative approach and only use a fundamental scalar QNM and its first overtone: Equations (40) and (41) fit to the real parts of the scalar multipoles, as their imaginary parts will have the same amplitudes and phases. We make use of all scalar mode waveforms presented in Sec. III D, and therefore perform the fitting on the (lm) = (11), (22), (33), (44) modes.
We sample the parameter space using emcee [104]. Dropping the (lmn) indices in Eqs. (40) and (41) and keeping just the mode number N ∈ [1, 2], we use uniform priors in our inference, namely log A i ∈ [−3, 1] and p i ∈ [0, 2π] for i ∈ [1, N ]. We compute the frequencies ω lmn with the qnm package [105] using our final spin calculations from Table II. There is an ongoing debate around when the optimal time is to begin fitting quasinormal ringdown (see e.g. [106]), as different start times can give different answers for the mode fits, and have also caused some controversy regarding the confidence of a subdominant QNM observation in GW data [107]. While overtones are typically quick to decay, they can play an important role as we approach the time of merger [108][109][110][111][112][113]. To avoid the intricacies of having to account for overtones, as well as nonlinear effects in the early ringdown [114,115], we take a conservative approach and start the ringdown analysis t 0 = 16M after the time of the peak amplitude of the scalar mode 6 . By starting too early we risk including nonlinear effects from the merger and by starting too late the SNR of the signal may significantly be lower. In our analysis, the damped sinusoid model is most accurate as an approximation in the late ringdown and so we sacrifice on the strength of the signal. The end time of our analysis is determined by the point in the late ringdown when the signal becomes very faint and the numerical noise becomes more significant. The point when this happens varies on the mode-by-mode basis, with higher modes being prone to numerical noise earlier than others. However, the typical length of our ringdown varies in the range of ∼ 15 − 30M .
In the optimization process, we use a least-squares likelihood in the time domain, with a flat noise spectral density where d(t i ) is the time sequence of numerical data, σ 2 is the variance and θ = {A i , p i }. We estimate the variance by taking our total error budget percentage (i.e. 4%) from the peak of the amplitude for each mode. As such, the error we account for through the variance parameter varies on a mode by mode basis. This choice for error estimation in our analysis assumes that the total NR error is not correlated in time. In reality, this is not likely to be an accurate assumption to make as in certain regions of the waveform we may in fact be overestimating or underestimating the allowed deviation of the data from the true values. Ideally, our likelihood function should contain some sort of correction to account for such correlation of the data with time. However, this is out of the scope of this work and we leave this for future study. We have verified the robustness of our choice of variance in the likelihood function of Eq. (42) by gradually decreasing its value and assessing its effect on the estimation of the modes' amplitudes. As expected, with decreasing variance, the amplitudes converge to some fixed values with smaller error bars, representing an ideal scenario where our NR error is very small. These idealized values provide the relative error we make in our fitting procedure when estimating amplitudes using our actual total NR error budget percentage to calculate fixed variance. Figure 7 shows this convergence of A i for i ∈ [1, 2] of the (22) mode of the BBH-12 configuration, where the fixed variance of σ 2 = 0.003 used in the fitting procedure results in median amplitudes lying within the error bars of the smallest variance.
In Table III we summarize the results of our fits. We note that we fixed the dimensionless coupling β 0 /m 2 1 to unity for convenience there. Therefore, when choosing a different value for the coupling, appropriate rescaling to the amplitudes has to be done according to Eq. (20). We find that the ratio of amplitudes of scalar Kerr QNMs to gravitational Kerr QNMs varies on a mode-by-mode basis for l > 1 modes. In the case of the l = 1 mode fit, the amplitude of the overtone n = 1 dominates the amplitude of its fundamental mode as also found in the analysis of gravitational ringdown of Ref. [109]. Unless indicated by dots in Table III, we find that for certain modes a N = 1 mode fit is able to describe the data accurately enough. This happens in the cases where the amplitude of one of the modes in the N = 2 mode fit is dominant, or when the frequencies of the N = 2 mode fit lie close to each other.
We additionally fit the amplitudes of our gravitational waveforms in Table IV, where we use (lm) = (11), (22), (33), (44), (55), (66) modes in the fitting procedure. We find that the N = 1 mode fit of Eq. (37) is sufficient to describe the data. Since the Einstein equations to leading order remain unchanged (see Eq. (17)), we have fixed the frequency of the N = 1 mode fit to the fundamental gravitational Kerr QNM of GR.
In summary, the leading order QNMs of the scalar and tensor components of the ringdown waveform in EsGB gravity are accurately described by QNMs of GR. The next correction in ϵ to the scalar QNM spectrum would depend on the nature of the coupling to the Gauss-Bonnet term, while the correction to gravitational QNM spectrum would happen at O(ϵ 2 ). In the case of subleading corrections to the scalar QNMs, for EdGB gravity with exponential coupling, β(φ) = e φ , the next correction will come from second order in ϵ, while for shiftsymmetric Gauss-Bonnet gravity with coupling β(φ) = φ, at the third order in ϵ [54]. Given the smallness of ϵ, we expect any kind of deviation from GR in the scalar QNM spectrum in Gauss-Bonnet gravity to be very small.

A. Measurability of the scalar ringdown
Having formulated a model for the scalar polarization component of the EsGB ringdown waveform, we now make a first attempt to assess the measurability of its presence in compact binary coalescence observations by the current and future networks of GW detectors. Returning to Eq. (8), we recall that the second term predicts an additional scalar contribution, whose strength is controlled by the magnitude of the conformal coupling a 0 defined in Sec. II B and the amplitude of the scalar field φ, which may be rescaled using the dimensionless parameter β 0 /m 2 1 according to Eq. (20). Therefore, the choice of the rescaling applied to the amplitude of the scalar field will affect the strength of the scalar polarization, and our choices of the re-scaling will be detailed in the following sections. However, before making predictions on the strength of the scalar signal, it is necessary to make the conversion to physical units. The scalar polarization now takes the form of where M fin is the final BH mass and D L is its luminosity distance from the Earth. We note that φ is the sum of all (lm) modes and is dependent on the inclination angle ι, which is defined as the angle between orbital angular momentum and the line-of-sight In our analysis we set ι = 60 • . As we discussed in Sec. II B, each polarization has a distinct geometrical im-print on a GW detector. We therefore begin with projecting the scalar strain, h S , onto the detector D located at Here F S is the response, or antenna pattern, of a detector D to scalar polarization. It can be given in terms of the polar angle θ and the azimuthal angle ϕ in the detector frame (i.e. θ and ϕ measured with respect to detector arms along the x and y axes) in the following form Note the absence of dependence on the polarization angle ψ that defines the orientation of the source projected on the celestial sphere; this is a direct consequence of the spin-0 nature of the scalar field, in contrast with the tensorial GWs that behave as a spin-2 field. The response to scalar polarization depends on the polar and azimuthal   (40) and (41). As in Fig. 3, we set β 0 /m 2 1 = 1 for convenience. The dots indicate that a specific fit could not describe the data accurately. In all cases with an N = 2 mode fit, A 1 represents the amplitude of the scalar Kerr QNM. For (11) mode, A 2 is the amplitude of the scalar Kerr overtone, while for all other modes, A 2 represents the amplitude of gravitational Kerr QNM. The upper and lower limits on the estimated parameters lie within 90% confidence interval and the central value is the median.  angles at which the GW comes in with respect to the detector frame. Figure 8 illustrates the change of the antenna pattern for various angles. We position the source at an optimal sky location for one of the detectors (LIGO Hanford in particular), which is straight down one of its arms, i.e. θ = π 2 , ϕ = 0. We then compute the squared SNR of the scalar signal defined via whereh D is the signal in the Fourier domain and S n (f ) is the noise power spectral density (PSD) of the detector, which we set to the estimated noise curve of i. the current network of LIGO (Hanford and Livingston) and Virgo detectors at design sensitivity [116,117] and ii. the Einstein Telescope (ET) [118] in the ET-D configuration [119]. We list luminosity distance and the range of masses chosen for each of the detectors in Table V.  In the computation of SNR from Eq. (47), the signal is weighted by the noise PSD curve in the frequency domain. Having a frequency of the scalar waveform close to the minima of the PSD would therefore increase the integrand of SNR for a fixed amplitude of the waveform. However, for varying mass, the amplitude is not fixed but scales according to Eq. (43). Further, the amplitude of the scalar signal is affected by the mass ratio: as suggested by our NR simulations from Sec. III D, amplitude increases for more asymmetric binaries.
We now move on to estimate the actual signal strength for two different scenarios with respect to the values of coupling parameters in EsGB.

Maximally allowed couplings
Our prior for the range allowed for the couplings a 0 and β 0 /m 2 1 , is controlled by the theoretical and observational constraints described in Sec. II D. In this section we make the conservative choice of max(β 0 /m 2 1 ) = ϵ thr ≈ 0.2 and a max 0 = √ 10 −3 ≈ 0.0316 to be maximally allowed values for our couplings (i.e. the blue region of Fig. 2 is where the theory can be said to be no longer predictive) and investigate in more detail the dependence of the SNR of the scalar ringdown signal on the total mass M fin and mass ratio q of the progenitors. Following common practice [120,121], we assume SNR ≥ 8 is required for a detection. Figure 9 shows the network SNR for the current LIGO (Livingston and Hanford) and Virgo detector network, calculated via quadrature, ρ network = j ρ 2 j , where ρ j represents the SNR in a single detector. Here the values of numerically calculated SNR span three horizontal lines of constant q = {1/2, 2/3, 1} and are interpolated on the parameter space of M fin × q. The typical scalar signal is weaker than its tensorial counterpart by at least 3 orders of magnitude, and the choice of maximally allowed couplings results in increased SNR for larger masses. The largest network SNR of 1.8 is produced from the smallest mass ratio of q = 1/2 and the final BH mass of M fin = 1000M ⊙ . We note, however, that our choice of t 0 = 16M to be the ringdown start time reduces the strength of the signal. As seen in Sec. IV, the amplitudes for highly asymmetric mass ratios become significantly damped at later times. Therefore, one may expect larger SNR, if an earlier ringdown start time was chosen.
The strength of the signal significantly improves in the ET, which we illustrate in Fig. 10. Here we further disentangle contributions of each mode to the total SNR by plotting SNR for the strongest individual modes, i.e. (lm) = (11), (22), (33). We observe that some cross terms between modes may contribute destructively in the integrand of SNR (see Eq. (47)) and as such, the total SNR is not simply a sum in quadrature of the individual modes' contributions. We find that the strongest total SNR calculated from all modes is recovered from the most extreme mass ratio considered in our simulations, q = 1/2. This is as expected since its modes' amplitudes are significantly larger than for milder mass ratios. Assuming luminosity distance of D L = 1Gpc for BBH-12 configuration, the total SNR of the scalar ringdown is weak for BHs with masses M fin ≲ 50M ⊙ and falls beyond the detectability threshold. However, taking parameters of GW151226 event with q ∼ 1/2, D L = 440 Mpc and M fin ∼ 20.8M ⊙ [122], we find the SNR in ET to be around 7.15. Interestingly, the SNR of individual modes has a more complicated structure and there are several competing factors contributing to their strength. First, the final BH mass determines where the QNM frequencies lie in relation to the high-sensitivity region of the detector's noise curve. Second, the hierarchy and relative amplitudes of the scalar modes are determined by the mass ratio of the progenitor (e.g. with odd-m modes vanishing for symmetric binaries). Third, higher modes are less long-lived for more extreme mass ratios, meaning that at later ringdown start times these modes become significantly damped. Finally, the inclination angle contributes an additional l-dependent geometric factor to the relative mode amplitudes as seen by the observer. As such, in our simulated analysis, the SNR of the (22) mode is the strongest for q = 2/3 and (11) and (33) modes are the strongest for the q = 1/2 one.
2. Imposing the most pessimistic constraint, √ β0 = 1.18 km So far we explored the strength of the signal in the optimistic scenario, where the strength of the EsGB dimensionless parameter is set at the limit of the theoretical bound and appears to be strong enough for certain binary configurations. Constraints on the EsGB coupling β 0 coming from data analyses of astrophysical observations (see Table I) further limit the range of the coupling and therefore are expected to significantly lower the strength of the scalar signal. We now follow the analysis of [56] on GW signals from the inspiral of BBH and BH-NS binaries and set the value of the EsGB coupling to their quoted 90% combined upper bound, i.e. √ β 0 = 1.18 km. This result is based on an analytically approximated inspiral model, where the leading-order correction to the waveform at −1PN (i.e. δΨ ∝ (β 2 0 /M 4 )v −7 ) is inferred and mapped to a bound on β 0 . Higher-order corrections to an incomplete 2PN level can vary the resulting bounds by more than 10%, depending on the binary parameters, therefore this value should not be considered as a robust upper bound at 90% confidence. Nevertheless, we expect that the upcoming fourth observing run (O4) of the LVK collaboration [117] will probe this value range with much higher confidence.
In this scenario the SNR of the scalar ringdown becomes significantly suppressed, which we demonstrate in Fig. 11 on a logarithmic scale. Unlike the previous case, where we bounded a 0 and β 0 /m 2 1 , here, smaller values of the final BH mass and more asymmetric mass ratios yield a stronger scalar polarization signal, as expected. In particular, the largest SNR of 0.42 is observed for a binary of mass ratio q = 1/2 with the smallest final mass of M fin = 10M ⊙ . Unless the source is observed at very small luminosity distance of D L ∼ 52.5 Mpc, the scalar ringdown signal of this binary configuration will not be detectable even with third generation detectors. Our predictions are less pessimistic for lower-mass binaries at even higher mass ratios, however, in this work we do not consider BH-NS binaries or BHs with masses lower than 5M ⊙ .

VI. BAYESIAN INFERENCE ON THE SCALAR RINGDOWN SIGNAL
We now perform full Bayesian inference on a hypothetical scalar signal present in the bands of the LVK detector network and ET. We model the data stream in each detector as where h(t) is the signal as described in Sec. IV and n(t) is a realization of detector noise. Although we use continuum-appropriate notation, in practice the above are realized as discrete time series with a sampling rate of 4096 Hz, which is sufficient for our purposes in this work.

A. Analysis setup
Before delving into the details of our analysis, let us first describe the simplifications performed here. First, we assume that the GR part of the signal is well modeled and has been successfully reconstructed and removed from the data, without being particularly interested in the details of that part of the analysis. In practice this means that we have reduced the problem of performing an apples-to-apples comparison between GR and EsGB to the problem of searching for the presence of our EsGB scalar-ringdown signal in the residual data. This, in turn, makes use of the decoupling-limit assumption, that is, GW generation is driven by GR dynamics and therefore the tensorial ringdown signals in GR and EsGB are identical. Furthermore, we factorize the analysis of multiple QNMs into individual analyses for each (lm) pair, here performed using N = 1 or N = 2 modes for each pair that can capture the presence of a "gravity-led" mode (or even an unusually strong overtone). This simplification assumes that the frequency content of e.g. the l = m = 1 and l = m = 2 modes that are studied here can be separated into disjoint regions in the frequency domain, as was done in recent searches for subdominant modes [107]. Although this latter assumption may not always hold, prior information on the final BH parameters, and hence on the expected scalar spectrum, will always be available from the analysis of the much louder tensorial signal ("+" and "×" polarizations), thus allowing for the scalar analysis setup to be adapted accordingly. Details related to the above considerations will not be discussed here, but are left to be explored in future work.
We thus define our two competing hypotheses as follows: 1. H 0 is our null hypothesis, corresponding to the data consisting of pure Gaussian noise, according to the PSD of each detector in our network.
2. H 1 is our scalar-ringdown signal hypothesis, corresponding to the presence of a scalar signal in addition to noise in our data, where the scalar model is described in Sec. IV and is parametrized by the real QNM amplitudes and phases (A lm k , p lm k ), while the complex mode frequencies are fixed to their expected values (see details in Eqs. (40) and (41)).
Our first task is then to infer whether our hypothetical scalar signal could be identified as being present in the data-according to (48), where h(t) = h S (t; ⃗ θ) is consistent with our scalar signal model-or the data is consistent with pure noise i.e. h(t) = 0. For scalar signals that are detected with confidence against the noise hypothesis, our analysis will also estimate the parameters ⃗ θ of our scalar ringdown model as posterior probability distributions on the model parameter space p( ⃗ θ|d, H 1 ). To accomplish this, we inject our NR scalar waveforms h S (t) into data streams of simulated Gaussian noise and perform time-domain ringdown analysis using the pyRing pipeline [123][124][125], which relies on a Python implementation of the nested sampling algorithm [126] called CPNest [127]. The built-in calculation of the evidences P (d|H i ) in nested sampling, by statistically integrating each likelihood function p(d| ⃗ θ, H i ) over the respective parameter space Σ i , enables model selection between the signal and noise hypotheses, by means of the Bayes factor, Since H 0 here is the noise hypothesis, its parameter space is the empty set and the integral in the denomina-tor reduces to a product of noise likelihoods p(d|H 0 ) ∝ exp{− 1 one for each detector, where C ij is the two-point autocovariance matrix, which characterizes the detector noise in the time domain as a discrete stochastic process that is assumed to be Gaussian and stationary [124]. Similarly, the pointwise likelihood for H 1 is computed using the same expression, where we replace d → d − h S (t; ⃗ θ). That is, the signal hypothesis states that the residual is pure noise.
We adopt two approaches in the injection of the scalar signal. First, we perform a zero-noise injection by setting n(t) = 0 in (48), in order to isolate possible biases due to realization-specific effects of the noise in our analysis and disentangle them from potential systematic biases inherent to our model. Note, however, that since the likelihood is weighted by the noise PSDs of our GW detectors, the overall effect of the noise on the uncertainty of our measurements is still incorporated in our method. Second, we inject the scalar waveform into simulated noise and assess its effects on the parameter estimation. We have followed the standards of the LIGO Scientific Collaboration Algorithm Library (LAL) [128] to perform NR scalar waveform injections, for which we have used a modified version of LALsuite [129] to accommodate the construction of the scalar waveform.
We present results for the injection of the BBH-12 binary configuration, which gives us the strongest scalar wave signal. We place the source along one of the arms of the detector, as described in Sec. V A, at a luminos-ity distance of D L = 100 Mpc, and choose a final mass of M fin = 500M ⊙ . We choose to inject the two loudest modes (lm) = (11) and (lm) = (22) and further take maximally allowed values for the theory couplings β 0 /m 2 1 and a 0 as described in Sec. V A, i.e. max(β 0 /m 2 1 ) = 0.2 and a max 0 = √ 10 −3 ≈ 0.0316. The SNR value for the configuration involving the (11) mode is roughly 700 and for the (22) is around 100. We use Table III of Sec. IV to determine which recovery model of Eqs. (40) and (41) to use for each of our injections. As such, for the (11) mode, we use the N = 2 fit with frequencies fixed to the fundamental and first overtone scalar Kerr QNMs, while for the (22) mode we use the N = 1 mode fit with the frequency fixed to the fundamental gravitational Kerr QNM.

B. Results
Using the two hypotheses, H 1 and H 2 , as defined in Sec. VI A, in Fig. 12, we present the results of the recovery of the (11) and (22) modes with zero-noise and colored Gaussian noise injections in ET. The results suggest that the scalar-ringdown signal hypothesis H 1 is favored for both injections, i.e. the data contains a scalar-polarized GW present at the expected frequencies within the detector's sensitivity band. The Bayes factors for the zero-noise and Gaussian noise injections are ln B 1 0 ∼ 16771 and ln B 1 0 ∼ 16127, and ln B 1 0 ∼ 47 and ln B 1 0 ∼ 15 for the (11) and (22) modes, respectively. Finally, in Fig. 13, we demonstrate the effect of Gaussian noise on the estimation of the modes' amplitudes and phases. Overall, with the source parameters of our choice we conclude that Gaussian noise does not compromise our ability to measure scalar polarization, neither does it introduce significant biases in recovering the modes.
As we have seen in Sec. V A, for a fixed set of coupling parameters a 0 and β 0 /m 2 1 , larger total masses or more extreme mass ratios are needed to achieve higher SNRs at a given distance. From the observational point of view, the astrophysical population statistics for BBH systems allow for a rather faint probability of either placing stringent bounds on EsGB or making a direct detection, based on the scalar ringdown signal of a single BBH source, even with a 3G/XG detector. Having multiple detectors in the network and combining information from multiple BBH events across a wide mass range will improve our overall sensitivity to the consistent presence of such signals. As the rate of BBH events with SNR>100 in the nextgeneration network is expected to reach O(10 2 − 10 3 ) per year [130], joint information from multiple events, even if they may be weak, will statistically improve the confidence of detecting a scalar signal as well as the precision of the inference by roughly a factor of ∼ √ n, n being the number of detected scalar signals [131]. While the presence of scalar polarization in the ringdown spectrum will be a smoking gun (unlike the inspiral dephasing), indicating the existence of a fundamental scalar field, more care will be needed in disentangling this weak effect from noise artifacts and other systematics.

VII. CONCLUSION
Through the presence of a conformal coupling to matter, black hole binaries in EsGB gravity emit gravitational waves that have a scalar polarization. The family of theories we consider in this work is controlled by two coupling parameters: the Gauss-Bonnet coupling β 0 and the conformal coupling a 0 , which both impact the strength of the scalar-polarized gravitational wave. We have studied the extent to which these couplings could be constrained through gravitational wave observations of the scalar polarization during binary black hole ringdown. By exploring the parameter space in the range of validity of our couplings, and calculating respective signal-to-noise ratios, we conclude that the scalar polarization has a low chance of being detected. There is a stronger possibility of detection with next generation detectors like the Einstein Telescope, targeting larger masses or more asymmetric binary progenitors, given the observed amplification in the scalar polarization amplitude for those systems from our NR simulations. Furthermore, future observations with LISA [132] may possibly allow for more precise measurements of ringdown from extreme mass-ratio inspirals, although the much higher mass range that LISA will probe may not be favorable for sourcing sufficiently strong scalar modes in theories like the ones studied here. The SNR is greatly dependent on the value of the dimensionless EsGB parameter, β 0 /m 2 1 . Even at the largest allowed theoretical value for the dimensionless coupling (β 0 /m 2 1 ∼ 0.2), one would need to have total mass of roughly ≳ 50M ⊙  and a fairly asymmetric mass ratio (q ∼ 1/2) in order for there to be a detection with a future generation GW detector. For smaller values of β 0 /m 2 1 , the strength of the signal drops beyond the detectability threshold for any inspiral scenario we considered, even for the Einstein Telescope. Given the present observational constraint of √ β 0 ≲ 1.18km, and that M ⊙ ∼ 1.5km in geometrized units, detecting a scalar polarization in a 50M ⊙ event requires the smaller mass black hole to have a mass of ∼ 2.5M ⊙ in order for a competitive constraint on β 0 to be possible from a ringdown measurement with the Einstein Telescope.
The most stringent constraint on most modified theories of gravity comes from the inspiral, through measurements of the phase of the gravitational wave signal [14,56]. This being said, a measurement of the inspiral phase does not completely determine the nature of a potential deviation from GR, as the dephasing could come from the emission of energy in other energy channels beyond (for example) that of a gravitational scalar field. Measuring the ringdown signal provides a complementary test to the measurement of the inspiral phase.
Apart from performing a targeted test on the GW data for signs of our EsGB model, the method that we propose here can be applied to a more general setting. EsGB with a conformal scalar coupling in its decoupling limit [133,134] is a prime example of a well-motivated theory, in which the spin-0 QNMs of the Kerr metric are excited and observed via radiation of scalar GWs. Since backreaction can be neglected in the limit of a small coupling β 0 /m 2 1 , part of the scalar quasinormal mode spectrum is identical to that of the spectrum of a (massless) scalar about a Kerr black hole [50,96]. More generally, any alternative theory featuring an additional scalar field that couples to the matter-field metric will also radiate spin-0 QNM frequencies, provided that the theory respects the "no-backreaction" property (i.e. to leading order, the GR sector determines the background metric and is driving the dynamics). In short, the general strategy amounts to using the information from the observed tensorial QNM frequencies to measure the remnant BH parameters, fixing the spin-0 QNM frequencies, and subsequently searching for the presence of any residual scalar GW signal on those frequencies, while being agnostic about the amplitude of each mode. The mode amplitudes are sensitive to the details of a given theory, since the field equations will determine how strongly each mode is sourced from the initial perturbation. Therefore, following the detection of a scalar signal, measurements of the relative mode amplitudes, combined with measurements on progenitor parameters, will narrow down the specifics of the theory and its couplings. A similar approach could also cover theories featuring a gravitational vector field. A thorough analysis of this strategy is beyond the scope of this article.
There are several more directions for future work. We have only considered the impact of the scalar polarization on the ringdown; it would be interesting to see if stronger constraints could be placed on β 0 and a 0 from polarization measurements of the inspiral waveform. We have solved for the equations of motion of EsGB gravity in the decoupling limit. While this approximation has been shown to be accurate for modeling the scalar waveform in EsGB gravity, at least during the inspiral phase of evolution (and ignoring the effects of dephasing) [50], it would be interesting to compare it to full solutions to EsGB gravity [26,27,54]. It will also be important to extend our study to spinning progenitors and see whether the scalar mode amplitudes' dependence on spins follows a functional form similar to that of the tensorial modes [135,136]. Finally, while we have argued that EsGB gravity with a conformal scalar coupling to the metric may give the largest scalar polarization signal for black hole binaries among scalar-tensor gravity theories, it would be interesting to numerically model other modified theories of gravity featuring scalar or vectorpolarized gravitational waves and devise similar targeted tests.

ACKNOWLEDGMENTS
We thank Ulrich Sperhake and Daniela Doneva for useful conversations on aspects of this project, Miren Radia for the help with GRChombo, Danny Laghi and Gregorio Carullo for getting started with pyRing and Maxence Corman for reviewing our preprint. We also thank the entire GRChombo collaboration for their support and code development work.
T In this section we follow the methodology of [72]. We start with our results on the convergence of the code and the discretization error due to finite resolution for simulation configuration BBH-11. We use ∆x high = 2M , ∆x med = 2.67M and ∆x low = 3.2M on the coarsest level of high, medium and low resolution runs, respectively. We note that the results presented in the main body were obtained from evolving the higher resolution configuration. The order of convergence is estimated for the amplitude and phase of Ψ 4 : In Fig. 14, where we have already multiplied the difference between medium and high resolutions by an appropriate factor, we plot the convergence order for the amplitude and phase of Ψ 4 . The convergence factor of order n for a given quantity τ is calculated by its difference at low/medium and medium/high resolutions: In the continuum limit the truncation error typically approaches zero as a power of the discretization parameter ∆x, that is for a given approximation error of order n, lim ∆x→0 τ ∆x = ∆x n × const.
(A3) Therefore, in the continuum limit the convergence factor reduces to For the amplitude of Ψ 4 we find the convergence order to be between second and third order in the inspiral, which then increases to between third and fourth in the late inspiral and merger. For the phase of Ψ 4 we find convergence consistent with fourth order dropping to third around merger. We perform a similar analysis for the amplitude of the scalar field amplitude and find the order of convergence fluctuating between second and third as indicated in the left panel of Fig. 15.
To estimate the discretization error for the amplitude e A , we measure the relative percentage error between the finest resolution result and Richardson extrapolation, while for the phase error e ϕ we simply give it as the difference, Finally, to estimate the error due to finite radius extraction we compute a given radiated quantity at several finite extraction radii and extrapolate to infinity by fitting a polynomial in 1/r. We use first order extrapolation of polynomial order n = 1: We then compute the corresponding errors for the amplitude and phase as given in Eqs. (A5) and (A6), where Richardson extrapolated quantities are now replaced by the extrapolated quantities prescribed by Eq. (A7). Additionally, we present convergence results for BBH-12 configuration, where we now use ∆x high = 2M , ∆x med = 2.29M and ∆x low = 2.67M on the coarsest level of high, medium and low resolution runs respectively. We present our results in Fig. 16 and the right panel of Fig. 15 for the gravitational (22) mode of Ψ 4 and the scalar field amplitude respectively. We find similar results to BBH-11 configuration. We observe between second and third orders of convergence in the amplitude and between third and fourth orders of convergence for the phase. The convergence order for the scalar field amplitude is second order in the inspiral and fourth order in the merger and ringdown.

GW energy and angular momentum flux convergence
Finally, we present convergence plots for BBH-11 of GW energy flux and GW angular momentum flux computed with the use of Eqs. (31) and (32), respectively. We employ the same set of resolutions as in Sec. A 1 and compute the total radiated energy and the total angular momentum flux, as described in Sec. III C of the main text.
From Fig. 17, we conclude that both the energy flux and angular momentum flux have an order of convergence between third and fourth. By using third order Richardson extrapolation, we estimate the error for final mass and spin to be of 0.01% and 0.05%, respectively.

Appendix B: Code Validation
In this section we address the tests on the implemented Gauss-Bonnet term in GRChombo. First of all, we performed a self-convergence test by evolving a scalar field on a Schwarzschild background. We use ∆x high = 1.347M , ∆x med = 1.6M and ∆x low = 2M on the coarsest level of high, medium and low resolution runs respectively. Fig. 18 demonstrates our results and overall we find second order of convergence.
Finally, we compare our numeric result for a scalar field evolved on a Schwarzschild background to the analytic solution in the Painlevè-Gullstrand (PG) coordinates. By starting with the scalar field equation we impose a regularity condition of ∂ r ϕ at the horizon r = 2M and the falloff of the scalar field at infinity, lim r→∞ ϕ = 0. We find an analytic solution given in terms of PG coordinates by (B2) We then calculate the L1-norm for the analytic and numeric solutions and normalize it over the total volume. Figure 19 shows good agreement between solutions after t ∼ 100M , when the scalar has grown and settled to a constant value. In Sec. IV and in particular Eqs. (40) and (41), we utilize QNMs of a Kerr black hole in GR to perform numerical fits to our data. In Table VI we present the values of the frequencies (ω lmn ) and damping times (τ lmn ), for BBH-11. The difference for medium/high resolutions is rescaled according to third and forth orders of convergence. Right: same as left but the diagnostic quantity is the angular momentum flux in the z direction extracted at R ext = 110M for BBH-11.
of the gravitational quasinormal modes of the remnant black hole. These modes were computed using the qnm package [105].