Numerical binary black hole collisions in dynamical Chern-Simons gravity

I Introduction

At some length scale, Einstein’s theory of general relativity (GR) must break down and be reconciled with quantum mechanics in a beyond-GR theory of gravity. Binary black hole (BBH) mergers probe the strong-field, non-linear regime of gravity, and gravitational waves from these systems could thus contain signatures of such a theory. Current and future gravitational wave detectors have the power to test GR Berti et al. 2015, and BBH observations from LIGO and Virgo have given a roughly 96% agreement with GR Abbott et al. 2016; Abbott et al. 2017.

These tests of GR, however, are presently null-hypothesis and parametrized tests Yunes et al. 2016; Abbott et al. 2016, which use gravitational waveforms produced in GR with numerical relativity. An open problem is the simulation of BBH systems through full inspiral, merger, and ringdown in beyond-GR theories. Waveform predictions from such simulations would allow us to perform model-dependent tests, and to parametrize the behavior at merger in beyond-GR theories.

In this study, we consider dynamical Chern-Simons (dCS) gravity, a beyond-GR effective field theory that adds a scalar field coupled to spacetime curvature to the Einstein-Hilbert action, and has origins in string theory, loop quantum gravity, and inflation Alexander and Yunes 2009; Green and Schwarz 1984; Taveras and Yunes 2008; Mercuri and Taveras 2009; Weinberg 2008. Computing the evolution of a binary system requires first specifying suitable initial conditions. Because the well-posedness of the initial value problem in full dCS gravity is unknown Delsate et al. 2015, we work instead in a well-posed order-reduction scheme, in which we perturb the metric and scalar field around a GR background Okounkova et al. 2017. The leading-order modification to the spacetime metric, and hence gravitational radiation, occurs at second order, which is precisely the order we consider in this study, building on our previous work Okounkova et al. 2017; Okounkova et al. 2018; Okounkova et al. 2019.

While our ultimate goal is to produce full inspiral-merger-ringdown waveforms relevant for astrophysical BBH systems, in this study we consider the leading-order dCS corrections to binary black hole head-on collisions. Such configurations, while less astrophysically relevant than orbiting binaries, serve as a proof of principle for our method of producing BBH waveforms in a beyond-GR theory Okounkova et al. 2019, and are fast and efficient to run. Head-on collisions also contain interesting science in their own right, as they cleanly probe the quasi-normal mode (QNM) spectrum of the post-merger gravitational radiation Anninos et al. 1993; Anninos et al. 1995; Baker et al. 2000; Sperhake et al. 2005. In this study, we thus produce the first BBH waveforms in a higher-curvature beyond-GR theory, and probe the leading-order dCS modification to the QNM spectrum of a head-on BBH collision.

This paper is organized as follows. We give an overview of our methods in Sec. II, and refer the reader to previous papers, Okounkova et al. 2019 and Okounkova et al. 2018, as well as Appendices B and C, for technical details. We discuss fitting perturbed quasi-normal modes in Sec. III. We present and discuss our results, including quasi-normal mode fits, in Sec. IV. We discuss the implications of this study on testing GR in Sec. V. We conclude in Sec. VI.

We set $G=c=1$ throughout. Quantities are given in terms of units of $M$, the sum of the Christodoulou masses of the background black holes at a given relaxation time Boyle and Mroue 2009. Latin letters in the beginning of the alphabet $\{a,b,c,d\ldots\}$ denote 4-dimensional spacetime indices, while Latin letters in the middle of the alphabet $\{i,j,k,l,\ldots\}$ denote 3-dimensional spatial indices (present in the appendices). $g_{ab}$ refers to the spacetime metric with connection $\Gamma^{a}{}_{bc}$, while $\gamma_{ij}$ (used in the appendices) refers to the spatial metric from a 3+1 decomposition with corresponding timelike unit normal one-form $n_{a}$ (cf. Baumgarte and Shapiro 2010 for a review of the 3+1 ADM formalism).

II Methods

Full details about order-reduced dynamical Chern-Simons gravity and our methods to simulate black hole spacetimes in this theory are given in Okounkova et al. 2019; Okounkova et al. 2018; Okounkova et al. 2017. Here we only briefly summarize.

The first term is the Einstein-Hilbert action of GR, with the Planck mass denoted by $m_{\textrm{\tiny{pl}}}$. The second term in the action is a kinetic term for the (axionic) scalar field. The third term, meanwhile, couples $\vartheta$ to spacetime curvature via the parity-odd Pontryagin density,

The equations of motion for $\vartheta$ and $g_{ab}$ have the form

and $T_{ab}^{\vartheta}$ is the stress energy tensor for a canonical, massless Klein-Gordon field

Because of $C_{ab}$ in Eq. (4), the equation of motion is different from that of a metric in GR sourced by a scalar field.

Each order in $\varepsilon$ leads to an equation of motion with the same principal part as GR, and therefore is known to be well-posed at each order. Order $\varepsilon^{0}$ gives the Einstein field equations of general relativity for $g_{ab}^{(0)}$, the background GR metric, minimally coupled to a massless scalar $\vartheta^{(0)}$, which we can consistently treat as “frozen out” and thus set to zero. The scalar field is unfrozen at order $\varepsilon^{1}$ (cf. Okounkova et al. 2017), and it takes the form of a sourced wave equation

where $\square^{(0)}$ is the d’Alembertian operator of the background and $\,{}^{*}\!RR^{(0)}$ is the Pontryagin density of the background.

Because $\vartheta^{(0)}$ vanishes, there is no correction to the metric at order $\varepsilon^{1}$. The leading-order dCS correction to the spacetime metric, which will produce the leading-order dCS correction to the gravitational radiation, occurs at order $\varepsilon^{2}$ (cf. Okounkova et al. 2017), and takes the linear form

where $G_{ab}^{(0)}$ is the linearized Einstein field equation operator of the background, and

where $\nabla_{a}{}^{(0)}$ denotes the covariant derivative associated with $g_{ab}^{(0)}$. Meanwhile,

To produce beyond-GR gravitational waveforms, our goal is thus to simultaneously evolve fully nonlinear vacuum Einstein equations for $g_{ab}^{(0)}$, Eq. (9) for $\vartheta^{(1)}$, and Eq. (10) for $h_{ab}^{(2)}$, to obtain the leading-order dCS correction to the spacetime metric and corresponding gravitational radiation.

(recall from Sec. I.1 that $M$ is the sum of the Christodoulou masses of the background black holes at a given relaxation time Boyle and Mroue 2009). With these substitutions, Eq. (9) becomes

where $T_{ab}^{(1)}[\Delta\vartheta]$ refers to the Klein-Gordon stress-energy tensor in Eq. (11) computed from $\Delta\vartheta$ instead of $\vartheta^{(1)}$, and $C_{ab}^{(1)}[\Delta\vartheta]$ similarly refers to the $C$-tensor in Eq. (12) computed with $\Delta\vartheta$ instead of $\vartheta^{(1)}$.

II.2 Evolution

To evolve the first-order dCS metric perturbation, we evolve three systems of equations simultaneously: one for the GR background BBH spacetime, one for the scalar field $\Delta\vartheta$ [cf. Eq. (14)] sourced by the background curvature, and one for the metric perturbation $\Delta g_{ab}$ [cf. Eq. (15)], sourced by the background curvature and $\Delta\vartheta$. We evolve all variables concurrently, on the same computational domain.

All variables are evolved using the Spectral Einstein Code SpE, a pseudo-spectral code. The GR BBH background is evolved using a well-posed generalized harmonic formalism, with details given in Lindblom et al. 2006; Scheel et al. 2009; Szilagyi et al. 2009; Hemberger et al. 2013. The first-order scalar field is evolved using the formalism detailed in Okounkova et al. 2017. Finally, the metric perturbation is evolved using the formalism given in Okounkova et al. 2019, a well-posed perturbed analogue of the generalized harmonic formalism. When evolving the metric perturbation, we have the freedom to choose a perturbed gauge, which we choose to be a harmonic gauge. We give details on perturbed gauge choices in Appendix B. We use the boundary conditions detailed in Cook and Pfeiffer 2004; Rinne et al. 2007; Okounkova et al. 2017; Okounkova et al. 2019.

We use the standard computational domain used for BBH simulations with the Spectral Einstein Code SpE (such as used in Boyle et al. 2019). The computational domain initially has two excision regions (one for each black hole), and the post-merger grid has one excision region (for the final black hole) (cf. Hemberger et al. 2013 for mode details). The outer boundary is chosen to be $\sim 700\,M$. We use adaptive mesh refinement (as detailed in Szilagyi et al. 2009), with the background GR variables governing the behavior of the mesh refinement. This is justified, as high gradients in the background will source higher gradients in both the scalar field and the metric perturbation. For all of the evolved variables, in spherical subdomains we filter the top four tensor spherical harmonics, while we use an exponential Chebyshev filter in the radial direction Szilagyi et al. 2009. We similarly filter the variables in subdomains with other topologies according the the prescriptions in Szilagyi et al. 2009. For the constraint damping parameters (cf. Lindblom et al. 2006; Okounkova et al. 2019), we choose the standard values for BBH simulations.

Because the code is pseudo-spectral, we expect roughly exponential convergence with numerical resolution in all of the evolved variables. We specify numerical resolution by choosing adaptive mesh refinement tolerances Hemberger et al. 2013; Szilagyi et al. 2009; because mesh refinement is based on thresholds, this means that in practice errors decay roughly but not rigorously exponentially—see Boyle et al. 2019 for further discussion. In Okounkova et al. 2019, we performed detailed tests of the metric perturbation system, showing exponential convergence of evolved variables. We will quote all physical extracted quantities (cf. Tables 1 and 2) with error bars given by comparing the highest two numerical resolutions.

II.3 Initial data

To perform an evolution, we must generate initial data for the background (metric) fields, the scalar field, and the metric perturbation. The background initial data for a BBH system are given by a constraint-satisfying superposition of black hole metrics in Kerr-Schild coordinates Lovelace 2009; Ossokine et al. 2015. The scalar field initial data are given by a superposition of slow-rotation solutions Okounkova et al. 2017; Yunes and Pretorius 2009; Yagi et al. 2012. The constraint-satisfying initial data for $\Delta g_{ab}$ are generated using the methods outlined in Okounkova et al. 2018. For head-on collisions, we start with a separation of $25\,M$, assuming that the contributions to the gravitational radiation and energy flux from times $t\lesssim 25\,M$ are negligible.

In this study, we will consider axisymmetric configurations where the background spins of the black holes are oriented along $\hat{x}$, the axis along which they are colliding. Moreover, we will choose configurations where the two spins have the same orientation along the axis of collision so that the system has a reflection symmetry for $x\to-x$ (recall that spin is a pseudo-vector). We illustrate this configuration in Fig. 1. We consider equal mass, equal spin configurations, with dimensionless spins $\chi$ between $0.1$ and $0.8$, in steps of $0.1$. Kerr with $\chi\neq 0$ is not a solution of dCS, and hence the initial configurations will have a non-zero dCS metric perturbation Yunes and Pretorius 2009. However, Schwarzschild is a solution of the theory, and hence we do not consider $\chi=0.0$, as there will be no metric perturbation in that case.

As a check, we also consider the opposite configuration to Fig. 1, where the spins have opposite orientations. For the equal mass, equal spin systems considered in this study, the final remnant in this case (for all spins) is a Schwarzschild black hole. As Schwarzschild is a solution of dCS, there is zero (to within numerical error) final dCS metric perturbation or scalar field in the spacetime.

II.4 Wave Extraction

In the order reduction scheme, $\Psi_{4}$, the Newman-Penrose scalar measuring the outgoing gravitational radiation is expanded about a GR solution as

If we substitute the expanded metric given in Eq. (7) into the expression for $\Psi_{4}$ (cf. Appendix C), we can match the terms order-by-order. $\Psi_{4}^{(1)}$, the first-order correction, will have pieces linear in $h_{ab}^{(1)}$. Recall, however, that $h_{ab}^{(1)}=0$, so $\Psi_{4}^{(1)}$ vanishes. $\Psi_{4}^{(2)}$, the second-order correction, will have pieces quadratic in $h_{ab}^{(1)}$, which will similarly vanish, and pieces linear in $h_{ab}^{(2)}$. Thus, the leading-order correction to the gravitational radiation will be linear in the leading-order correction to the spacetime metric.

and we compute $\Delta\Psi_{4}$ using the methods detailed in Appendix C.

Throughout the evolution, we extract $\Psi_{4}^{(0)}$ and $\Delta\Psi$ on a set of topologically spherical shells using the methods given in Taylor et al. 2013. We similarly extract the scalar field $\Delta\vartheta$ radiation on these spherical shells (cf. Okounkova et al. 2017). $\Psi_{4}^{(0)}$ and $\Delta\Psi_{4}$ are then fit to a power series in $1/r$ (where $r$ is the radius of the spherical shell) and extrapolated to infinity using the methods given in Taylor et al. 2013; Boyle and Mroue 2009. We report all of the quantities as $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$.

III Perturbations to quasi-normal modes

Once we have obtained $r\Psi_{4}^{(0)}$, the background gravitational radiation, and $r\Psi_{4}^{(2)}$, the leading order dynamical Chern-Simons deformation to the gravitational radiation, we can analyze the quasi-normal mode spectrum. As discussed in Sec. I, head-on BBH collisions cleanly probe the quasi-normal mode (QNM) spectrum of the post-merger spacetime. We are thus most interested in fitting for the QNM spectrum of $r\Psi_{4}^{(0)}$, and the leading-order deformation to this spectrum in $r\Psi_{4}^{(2)}$. A more technical/abstract derivation can be found in Appendix A.

A GR QNM waveform takes the form of a superposition of damped sinusoids

Since the GR background gravitational radiation is composed of QNMs, we can use the form above to fit for $\Psi_{4}^{(0)}$ for each mode:

The quantities $\omega^{(0)}$ and $\tau^{(0)}$ are known from perturbation theory for each $(l,m,n)$ Stein 2019. Our fit thus determines two free parameters for each mode: $A^{(0)}$, and $\theta^{(0)}$.

III.2 Perturbed quasi-normal modes

Let us now consider how to fit $r\Psi_{4}^{(2)}$ after the merger. Note that all fitting is performed in the time domain (cf. London et al. 2014, Giesler et al. 2019). The QNM frequency, damping time, and amplitude will all be corrected from the background values as

Let us focus on the real part of Eq. (22). Computing the leading-order perturbation to this expression gives us the form

The imaginary part is similarly modified as

III.3 Predictions for particular and homogeneous solutions

The metric perturbation $h_{ab}^{(2)}$ satisfies a linear inhomogeneous differential equation. Its general solution will be a linear combination of a homogeneous and particular solution. Shortly after merger, the source driving $h_{ab}^{(2)}$ is constructed from both the dCS scalar field $\vartheta^{(1)}$ and the nonstationary background spacetime. At very late times, when $\vartheta^{(1)}$ settles down to a stationary configuration, the source for $h_{ab}^{(2)}$ will also be stationary, sourcing just the stationary deformation $h_{ab}^{\textrm{Def}}$ away from Kerr, plus any remaining homogeneous solution. That homogeneous solution coincides with the GR homogeneous solution, and thus has the same frequency and decay time as QNMs in GR (see Appendix A for a more rigorous derivation).

Meanwhile, at earlier times just after merger, the oscillating $\vartheta^{(1)}$ will generate a source term for $h_{ab}^{(2)}$ that oscillates at the scalar field’s frequency, and decays at the rate of the scalar field’s decay. Thus at early times just after merger, there can be a substantial particular solution with a different frequency and decay time than the late-time behavior. We observe this behavior in Sec. IV.4.

III.4 Scaling

Similarly, given $\tau^{(2)}$, we can perturb the above expression to give

III.5 Mass and spin definitions

III.6 Fitting window

III.7 Practical considerations

To perform these linear fits, we use a least-squares method Jones et al. 01. We fit a sum of overtones to each mode $(l,m)$. We shift $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$ to align at the peak of $r\Psi_{4}^{(0)}$ for each mode. We compute errors in our estimates of the parameters by considering the fitted values for a medium numerical resolution simulation and a high numerical resolution simulation (for the same initial configuration).

IV Results

During each simulation, we extract $r\Psi_{4}^{(0)}$, the Newman-Penrose scalar measuring the outgoing gravitational radiation of the background spacetime, decomposed into spin-weight $-2$ spherical harmonics labelled by $(l,m)$. Similarly, we extract and decompose $r\Psi_{4}^{(2)}$, the leading-order dCS correction to the gravitational radiation. Since the computational domain is of finite extent, both quantities are extrapolated to infinity. We additionally extract $\vartheta^{(1)}$, the scalar field, decomposed into spherical harmonics. In all cases, the spherical harmonics’ azimuthal axes are oriented along the collision axis of the black holes, which we will call $\hat{x}$. Note that we decompose into spin-weighted spherical harmonics Taylor et al. 2013, not spheroidal harmonics (which do not form a basis), and ignore spherical-spheroidal mode mixing.

Finally, we plot the dominant modes of $\vartheta^{(1)}$, the leading-order dCS scalar field for this configuration, in Fig. 4. Because the scalar field around each black hole takes the form of a dipole oriented around $\hat{x}$ (cf. Yagi et al. 2012), and the spins are pointing in the same direction (cf. Fig. 1), we expect power only in the odd $l$ modes. Because of the axisymmetry of the configuration, we expect only the $m=0$ modes to be excited. We see the $(1,0)$ mode asymptotes to a value that corresponds to the remnant dipolar profile of the scalar field on the final black hole.

IV.2 Regime of validity

Recall, however, that the order-reduction scheme is perturbative. The modifications to the spacetime must actually form a convergent perturbation series around GR. We thus require that $g_{ab}$, the background metric, have a larger magnitude than $h_{ab}^{(2)}$ at each point in the spacetime:

for some tolerance $C$. This gives an instantaneous regime of validity. Following Eq. (15), we can compute

In practice, the ratio is taken point-wise on the computational domain. We choose $C=0.1$ as a rough tolerance.

IV.2.2 Secular regime of validity

IV.3 Quasi-normal mode fits

We perform the quasi-normal mode fits detailed in Sec. III to $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$. We fit three overtones to each $(l,m)$ mode. For each mode of $r\Psi_{4}^{(0)}$, we use the perturbation theory results for the corresponding $\omega^{(0)}$, the GR QNM frequency, and $\tau^{(0)}$, the GR damping time Stein 2019, and fit for the the QNM amplitudes (cf. Eq. (22)). From $r\Psi_{4}^{(2)}$, we extract $\omega^{(2)}$, the leading-order dCS correction to the QNM frequency, and $\tau^{(2)}$, the leading-order correction to the QNM damping time, as well as the leading-order corrections to the QNM amplitudes (cf. Eqs. (28) and (29)).

We tabulate all of our fit results in Tables 1 and 2. We quote errors on each of the quantities by comparing the results from the two highest numerical resolutions of the NR simulation.

However, the dCS-modified QNMs should always be exponentially decaying, meaning that we require

Let us consider the sources of error in these computations. The resolution of the simulation is the dominant source of error (for example, varying the fitting window as detailed in Sec. III.6 does not significantly change the results). In each of Figs. 9 and 11, as well as the tabulated values in Tables. 1, and 2, the error bars on a fitted quantity $Q$ are computed by comparing the value of $Q$ for two simulations with different numerical resolutions (cf. Sec. IV.3). The error bars on the fits for $\tau^{(2)}$, and $\omega^{(2)}$ increase with $l$, being lowest for the $(2,0)$ mode, and highest for the $(6,0)$ mode. Higher modes are more difficult to resolve numerically SpE; Scheel et al. 2009, and thus it takes higher resolution for the error bars on the $(6,0)$ mode to decrease to those on the $(2,0)$ mode at lower resolution. The errors also increase with the spin of the system. This is because it is more difficult to resolve higher spin systems numerically Lovelace et al. 2008; Lovelace et al. 2011.

IV.4 Particular and homogeneous solutions

There is interesting behavior later on in the $r\Psi_{4}^{(2)}$ waveforms. As we can see from e.g. Fig. 3, there is a change in the overall slope that occurs in $r\Psi_{4}^{(2)}$ around $40\,M$ after the peak time. This change in slope is convergent with resolution, and is present with and without adaptive mesh refinement. Later in the waveform, after this change in slope, both $r\Psi_{4}^{(2)}$ and $r\Psi_{4}^{(0)}$ are well-described by damped sinusoids, and have the same decay time and frequency. In other words, at late times $r\Psi_{4}^{(2)}$ has the same QNM spectrum as $r\Psi_{4}^{(0)}$ a GR QNM on a Kerr spacetime (we know from previous work that the resulting GR spacetime of the numerical simulations with this code is Kerr Owen 2009; Bhagwat et al. 2018). This suggests that $r\Psi_{4}^{(2)}$ switches from being dominantly driven by the dCS scalar field entering the source term, to being source-free, as suggested in Sec. III.3 and discussed further in Appendix A.4. In other words, the early post-merger dCS waveform correction is dominated by a particular solution of Eq. (10), whereas later it is dominated by a homogeneous solution.

We illustrate this behavior schematically in Fig. 13. We consider the slopes of the logarithms of $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$, which is equivalent to finding a decay time for each. Note that this is not the same as the perturbed fits for $r\Psi_{4}^{(2)}$ given in Eqs. (28) and (29), which we use to extract $\tau^{(2)}$ and $\omega^{(2)}$. At late times, the decay times of $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$ are the same. In other words, the late-time leading order dCS modification to the gravitational radiation is the same as a GR QNM on Kerr. This behavior is consistent across all spins considered in this study.

We can corroborate this interpretation by looking at the scalar field in the strong field region, whose dynamics drive the radiative part of $h_{ab}^{(2)}$. As the scalar field settles down, it is no longer the dominant source driving $h_{ab}^{(2)}$, and the metric perturbation is dominantly driven by the Kerr background. However, making this interpretation more precise would be tricky: we must keep in mind that mapping between the strong-field region and a gravitational waveform at infinity requires utmost care (cf. Bhagwat et al. 2018).

V Implications for testing general relativity

Let us now discuss this work in the context of testing GR with gravitational wave observations. Suppose that we were to observe a post-merger gravitational wave, given by $r\Psi_{4}$. We can check the consistency of this waveform with GR, and with dCS.

If the values are not consistent with GR, meaning that there is a shift away from the predicted GR frequencies and damping times, we can check whether the fitted $\omega_{(l,m,n)}$ and $\tau_{(l,m,n)}$ agree with the leading-order dCS-corrected frequencies and damping times computed using the methods in this study.

V.2 Checking non-degeneracy: full case

Non-degeneracy in this context means that the dimension of the span of $\{\phi_{*}\partial/\partial\chi,\phi_{*}\partial/\partial M,\phi_{*}\partial/\partial\varepsilon^{2}\}$ is 3. This can be checked by looking at the rank of the $2k\times 3$ dimensional (Jacobian) matrix

VI Conclusion

In this study, we have produced the first beyond-GR BBH gravitational waveforms in full numerical relativity for a higher-curvature theory. We have considered head-on collisions of BBHs in dynamical Chern-Simons gravity. While these are not likely to be astrophysically relevant configurations, they serve as a proof of principle of our ability to produce beyond-GR waveforms Okounkova et al. 2019. Future work in this program thus involves adding initial orbital angular momentum to the system and producing beyond-GR gravitational waveforms for inspiraling systems. We have previously evolved the leading order dCS scalar field for an inspiraling BBH background Okounkova et al. 2017, and can use our (fully-general) methods given in Okounkova et al. 2018 and Okounkova et al. 2019 to produce initial data for and evolve an inspiraling BBH system.

We have also studied modifications to the post-merger BBH head-on collision QNM spectra. We found that at leading order, the damping time of each QNM receives a modification that increases with the spin of the final black hole as a power law. The frequency of each QNM receives a similar modification. These modifications are not degenerate with GR.

When performing inspiraling BBH simulations, we can repeat the analysis outlined in this paper to learn about the dCS modification to the QNM spectrum of an astrophysically relevant system. These results can then be applied to beyond-GR tests of BBH ringdowns Abbott et al. 2016; Yunes et al. 2016. In particular, the investigation of modification to overtones is useful for an analysis of the form detailed in Isi et al. 2019. Note that the GR-dCS non-degeneracy results found in this paper assume infinite signal to noise ratio. Thus, future work also includes checking degeneracy in the presence of gravitational wave detector noise. Inspiraling simulations will also allow us to perform even more powerful tests of GR using full inspiral-merger-ringdown waveforms, thus taking advantage of the entire gravitational wave signal.

Acknowledgements

This work was supported in part by the Sherman Fairchild Foundation, and NSF grants PHY-1708212 and PHY-1708213 at Caltech and PHY-1606654 at Cornell. Computations were performed using the Spectral Einstein Code SpE. All computations were performed on the Wheeler cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by Caltech.

Appendix A Formalism for QNMs beyond GR

In this Appendix we give an abstract formalism for QNM modeling in theories beyond GR. For simplicity we will present perturbative expansions with leading power $\varepsilon^{1}$, but the generalization to the behavior of dCS (where the leading metric correction is at $\mathcal{O}(\varepsilon^{2})$) is straightforward.

Under the final state conjecture Penrose 2002; Chrusciel et al. 2012; Klainerman 2002, the result of a merger of two Kerr black holes will uniquely be a perturbed Kerr BH, with the perturbation decaying with time. Therefore, the waveform after merger is typically modeled using linear black hole perturbation theory. That is, we treat the post-merger metric as

where $g_{ab}^{K}$ is the Kerr metric, and $\xi$ is a small formal order-counting parameter. The full metric satisfies the nonlinear Einstein equations up to order $\xi^{2}$, yielding the linear partial differential equation (PDE) that $h_{ab}$ satisfies,

where $G^{(1)}[\cdot]$ is the Einstein operator linearized about the Kerr background. Notice that QNMs are homogeneous solutions to this linear PDE.

In practice, the metric perturbation equations (50) of GR are intractable, whereas curvature perturbations for the Weyl scalar $\Psi_{0}$ and $\Psi_{4}$ can be decoupled, yielding the Teukolsky equation Teukolsky 1972; Teukolsky 1973, which we will denote as

In the Teukolsky formalism, QNMs are still homogeneous solutions. As Teukolsky showed, this partial differential equation is amenable to separation of variables. The most general homogeneous solution, at large $r$, is a linear combination

A.2 QNMs beyond GR

Now suppose we are interested in some (unknown) beyond-GR theory that is UV complete. We assume that there is a parameter $\varepsilon$ such that as $\varepsilon\to 0$, this theory recovers GR. Therefore we can expand around $\varepsilon=0$ to control the calculation. If we keep some leading number of terms, this theory will coincide with an order-reduced EFT, for example with order-reduced dCS as presented in Sec. II.1. There are now two small parameters: the size of the GW perturbations $\xi$, and the amount of deformation away from GR, $\varepsilon$. The metric and any additional new degrees of freedom will now be a bivariate expansion in $(\varepsilon,\xi)$. We will demonstrate this abstractly.

Suppose the nonlinear, UV-complete EOMs can be written as

where $\boldsymbol{u}$ is a vector of all the field variables [for example, any UV completion of dCS must include at least $\boldsymbol{u}=(g_{ab},\vartheta,\ldots)$]. The new degrees of freedom should be frozen out in the $\varepsilon\to 0$ limit.

We assume these equations have a family of nonlinear solutions for stationary, axisymmetric BHs $\boldsymbol{u}^{{\textrm{\tiny{BH}}}}(\varepsilon,M,a,\ldots)$. As in GR, we conjecture that the final state of a merger of two BHs will be a unique perturbation to a single BH of this family. Therefore, we can perform linear perturbation theory about one solution, positing

and similarly for all other fields in $\boldsymbol{u}(\varepsilon)=\boldsymbol{u}^{{\textrm{\tiny{BH}}}}(\varepsilon)+\xi\boldsymbol{u}^{(1)}(\varepsilon)+\mathcal{O}(\xi^{2})$. Linearizing Eq. (54) about $\boldsymbol{u}^{{\textrm{\tiny{BH}}}}$ would derive the linear EOMs that parallel GR’s metric perturbation equations (50), with some linear PDE

The QNMs in this beyond-GR theory need not be diagonal in field space $(g_{ab},\vartheta,\ldots)$. For example, mixed QNM modes are present in dCS, as discussed in Molina et al. 2010, and in Einstein-dilaton-Gauss-Bonnet gravity, as claimed in Blázquez-Salcedo et al. 2016; Blázquez-Salcedo et al. 2017; Pani and Cardoso 2009; Witek et al. 2019 (though see below for further discussion). However, since all new degrees of freedom beyond the metric should freeze out in the limit $\varepsilon\to 0$, the mode structure should become diagonal in field space in this limit. Regardless of the diagonal basis, in the $r\to\infty$ limit, $r\Psi_{4}$ will be a linear combination of the form

A.3 Perturbative treatment of QNMs beyond GR

Now turn to the bivariate expansion of Eq. (55),

Here $h^{{\textrm{\tiny{Def}}}}$ is the deformation of the full BH metric away from the Kerr metric, $g^{{\textrm{\tiny{BH}}}}=g^{K}+\varepsilon h^{{\textrm{\tiny{Def}}}}+\mathcal{O}(\varepsilon^{2})$, and similarly the QNMs can be expanded into their GR parts $h^{Q\,(0)}$ and deformation $h^{Q\,(1)}$. Doing so requires expanding all quantities in powers of $\varepsilon$,

If we plug the leading order corrections from Eqs. (59) and (60) into the ansatz Eq. (57), we will arrive at the form

The quantity denoted by “mode mixing” is proportional to $\frac{d}{d\varepsilon}S_{\lambda}(\varepsilon,\theta)$. Since we are already ignoring the difference between spherical and spheroidal harmonics in GR, we also ignore this extra mode mixing term in this manuscript.

A.4 Particular, homogeneous modes, and modes that are mixed in field space

Now including the terms at $\mathcal{O}(\xi)$, we see

At linear order in $\xi$, $H_{ab}^{(1)}[h^{Q\,(0)}_{cd}]$ (expanded about $g^{K}$) will generate a source term for $h^{Q\,(1)}$. In our numerical implementation we have the full GR metric solution, not an expansion in powers of $\xi$, so there can also be nonlinear (QNM)2 and higher terms appearing in the source.

In Witek et al. 2019, the authors simulated the leading-order scalar field behavior on a binary black hole background in Einstein-dilaton-Gauss-Bonnet gravity, and found that the solution for the scalar field during ringdown contained two parts, similar to how we observed two pieces in Sec. IV.4. The authors referred to these as “scalar-led” and “gravitational-led” modes which they suggest are due to mixing in field space.

This nomenclature was introduced in Blázquez-Salcedo et al. 2016 but was seen earlier in dCS in Molina et al. Molina et al. 2010. In Molina et al. 2010, the authors investigated QNMs of Schwarzschild black holes in full dCS gravity. For zero spin, the system is well-posed, and thus can be solved in the full theory, without working in an order-reduction or other perturbative scheme. Working on a Schwarzschild background within the order-reduction scheme also faithfully reproduces $\vartheta^{(1)}=0$ and $g_{ab}^{(2)}=0$. The radial parts of the scalar and gravitational QNMs for each mode are governed by a set of fully coupled ordinary differential equations (ODEs) of the form (cf. Eqs. (2.8) and (2.9) in Molina et al. 2010),

However, we and Witek et al. 2019 are both working in a perturbative scheme. In the perturbation scheme, the leading-order dCS metric perturbation $h_{ab}^{(2)}$ does not back-react onto the leading scalar field $\vartheta^{(1)}$. If one applies the perturbation scheme to Eq. (68) on the Schwarzschild background, the ODEs take the form

This matrix is triangular, so the solution for $\vartheta^{(1)}$ can be found independently of $\Psi^{(2)}$. Meanwhile, the QNMs of $\vartheta^{(1)}$ enter into the source term for $\Psi^{(2)}$.

Thus, the presence of the “scalar-led” mode (now on a nonlinear, perturbed Kerr background) seen by Witek et al. 2019 is not surprising, because the homogeneous solution of the beyond-GR scalar perturbation equation will contain the same scalar QNMs (homogeneous solutions) as in GR. Meanwhile, the “gravitational-led” mode in $\vartheta^{(1)}$ is surprising. However, looking back at Eq. (67) suggests its origin (keeping in mind that everything in order-reduced dCS and EDGB is pushed up by a power of $\varepsilon$, so we should think of Eq. (67) of suggesting the form of the scalar equation for $\vartheta^{(1)}$). The nonlinear, perturbed Kerr background enters the source term on the RHS of (67), generating a source that oscillates at the frequency of the background (GR) gravitational waves. This sources a particular solution at this frequency.

Therefore we conjecture that the “gravitational-led” mode appearing the scalar field in Witek et al. 2019 was actually a particular solution, though it is fair to still consider it as part of the QNM spectrum. We can more generally conjecture that when the perturbative approach is applied to field-space mixed QNM modes, they will appear as combinations of homogeneous and particular solutions of the linear equations. Investigating this conjecture is beyond the scope of this work.

Appendix B Choosing a perturbed gauge

Throughout this appendix, as well as Appendix C, we use the notation developed in Okounkova et al. 2019, and standard 3+1 ADM decomposition notation Baumgarte and Shapiro 2010. Recall that $g_{ab}$ refers to the 4-dimensional spacetime metric, while $\gamma_{ij}$ refers to the 3-dimensional spatial metric. $\Delta Q$ is the leading-order perturbation to quantity $Q$.

The generalized harmonic evolution for the background follows from the equation

where $\Gamma_{a}\equiv g^{bc}\Gamma_{bca}$, and $H_{a}$ is known as the gauge source function (cf. Lindblom et al. 2006 for more details). Throughout the evolution, the gauge constraint,

When generating initial data for $g_{ab}$ and $\partial_{t}g_{ab}$, we are free to choose $\partial_{t}\alpha$ and $\partial_{t}\beta^{i}$, the initial time derivatives of the lapse and shift. These quantities appear in $\Gamma_{a}$, so choosing them is equivalent to choosing initial values of $H_{a}$, via Eq. (71). For example, for initial data in equilibrium, we can set $\partial_{t}\alpha=0$ and $\partial_{t}\beta^{i}=0$, and set $H_{a}$ to initially satisfy Eq. (71). Alternatively, we can choose to work in a certain gauge, such as harmonic gauge with $H_{a}=0$, and set $\partial_{t}\alpha$ and $\partial_{t}\beta^{i}$ to satisfy Eq. (71).

As the evolution progresses, we can either leave $H_{a}$ fixed, or continuously ‘roll’ it into a different gauge, with the restriction that it contains only up to first derivatives of $g_{ab}$ to ensure well-posedness. In practice, for BBH in GR, we work in a damped harmonic gauge, with $H_{a}$ specified using the methods given in Szilagyi et al. 2009.

The perturbed generalized harmonic evolution takes a similar form as Eq. (70), with

where $\Delta\Gamma_{a}$ is the first-order perturbation to $\Gamma_{a}$, and $\Delta H_{a}$ is a perturbed gauge source function. Similar to Eq. (71), we have a perturbed gauge constraint,

At the start of the evolution, we similarly have the freedom to choose $\Delta H_{a}$, provided that it contains no higher than first derivatives of $\Delta g_{ab}$, and satisfies the perturbed gauge constraint Eq. (73). When solving for perturbed initial data (cf. Okounkova et al. 2018), we similarly have the freedom to choose $\partial_{t}\Delta\alpha$ and $\partial_{t}\Delta\beta^{i}$, the time derivatives of the perturbed lapse and shift. An easy choice, for example, is to work in a perturbed harmonic gauge,

Let us now work out how to set $\partial_{t}\Delta\alpha$ and $\partial_{t}\Delta\beta^{i}$ in order to satisfy Eq. (73) for some desired perturbed gauge source function $\Delta H_{a}$. Let us first consider the unperturbed case, setting $\partial_{t}\alpha$ and $\partial_{t}\beta^{i}$ for some gauge source function $H_{a}$. We will work with the $\kappa_{abc}$ variable, which is the fundamental variable encoding the spatial and time derivatives of the metric (cf. Okounkova et al. 2019) as

where $n^{c}$ denotes the timelike unit normal vector. We can use our freedom to set $\partial_{t}\beta^{i}$ and $\partial_{t}\alpha$ to modify $\kappa_{abc}$ to satisfy $\Gamma_{a}=-H_{a}$ as

where $\Gamma_{ijk}$ is the spatial Christoffel symbol of the first kind, and

where $\kappa_{00j}$ in the above expression is given by Eq. (77). We can then use this modified $\kappa_{abc}$ to compute $\Gamma_{a}$ and ensure that Eq. (71) is satisfied for $H_{a}=H_{a}$.

Perturbing Eqs. (77) and (78), we can get an expression for a modified $\Delta\kappa_{abc}$ to satisfy Eq. (73) for some desired perturbed gauge source function $\Delta H_{a}$. We thus obtain

Note that this computation also uses the gauge source function of the background, $H_{a}$. Assuming that the background is in a satisfactory gauge, we set $H_{a}$ to the initial background gauge source function. All of the perturbed quantities in Eqs. (79) and (80) are given in Okounkova et al. 2019.

In this study, we choose to work in a perturbed harmonic gauge, with $\Delta H_{a}=0$.

Appendix C Computing perturbed gravitational radiation

The outgoing gravitational radiation of a spacetime is encoded in the Newman-Penrose scalar $\Psi_{4}$. In order to compute the leading-order correction to the binary black hole background radiation due to the metric perturbation $\Delta g_{ab}$, we need to compute $\Delta\Psi_{4}$, the leading-order correction to $\Psi_{4}$.

$\Psi_{4}$, a scalar, is computed on a topologically spherical surface from a rank-two tensor $U_{ij}$, contracted with a tetrad (in our case, a coordinate tetrad that converges to a quasi-Kinnersley tetrad at large radii). $U_{ij}$ on a surface with normal vector $\hat{n}^{i}$ takes the form

where $E_{ij}$ is the electric Weyl tensor, $B_{ij}$ is the magnetic Weyl tensor, $\epsilon_{ijk}$ is the (spatial) Levi-Civita tensor, and the projection operators are given by

Here, the vector $\hat{n}^{i}$ and the one form $\hat{n}_{i}$ are normalized using $N\equiv\sqrt{\gamma^{ij}n_{i}n_{j}}$ with $n^{i}=\gamma^{ij}n_{j}$.

In order to perturb $\Psi_{4}$, let us write the electric and magnetic Weyl tensors in Eq. (81) in terms of the extrinsic curvature $K_{ij}$,

where $R_{ij}$ is the spatial Ricci tensor and $D_{i}$ is the spatial covariant derivative associated with $\gamma_{ij}$.

All of the perturbed quantities $\Delta g^{ij},\Delta K_{ij},\Delta(D_{k}K_{ij})$, and $\Delta R_{ij}$ are given in terms of the perturbation to the spatial metric, $\Delta\gamma_{ij}=\Delta g_{ij}$, its spatial derivative $\partial_{k}\Delta\gamma_{ij}=\partial_{k}\Delta g_{ab}$, and its time derivative, $\partial_{t}\Delta\gamma_{ij}=\partial_{t}\Delta g_{ij}$ in Okounkova et al. 2018. Note that since we use a first-order scheme, we have access to $\Delta g_{ab}$, $\partial_{c}\Delta g_{ab}$ throughout the evolution (cf. Okounkova et al. 2019).

Let us now work through the perturbations to the normal vectors and projection operators. Because we want the perturbation to the gravitational radiation to be extracted on the same surface as the background gravitational radiation, we will hold the unnormalized one-form to the surface, $n_{i}$, fixed. That is, $\Delta n_{i}=0$. From this, we can then compute

We can then perturb the projection operators,

where $\Delta\gamma^{i}{}_{j}=\Delta\gamma^{ik}\gamma_{kj}+\gamma^{ik}\Delta\gamma_{kj}$

Once we obtain $\Delta U_{mn}$, we use the same tetrad to generate $\Delta\Psi_{4}$ from $\Delta U_{ij}$ as we do for $\Psi_{4}$.

References