Numerical binary black hole collisions in dynamical Chern-Simons gravity

TL;DR

This paper presents the first numerical relativity simulations of binary black hole collisions in dynamical Chern-Simons gravity, revealing unique modifications to the quasi-normal mode spectrum that could aid gravitational wave tests of gravity.

Contribution

It provides the first numerical waveforms for black hole collisions in a higher-curvature gravity theory beyond general relativity, focusing on dynamical Chern-Simons gravity.

Findings

Modified quasi-normal mode frequencies depend on remnant spin.

Modifications are not degenerate with Kerr black hole frequencies.

Results can help test general relativity with gravitational waves.

Abstract

We produce the first numerical relativity binary black hole gravitational waveforms in a higher-curvature theory beyond general relativity. In particular, we study head-on collisions of binary black holes in order-reduced dynamical Chern-Simons gravity. This is a precursor to producing beyond-general-relativity waveforms for inspiraling binary black hole systems that are useful for gravitational wave detection. Head-on collisions are interesting in their own right, however, as they cleanly probe the quasi-normal mode spectrum of the final black hole. We thus compute the leading-order dynamical Chern-Simons modifications to the complex frequencies of the post-merger gravitational radiation. We consider equal-mass systems, with equal spins oriented along the axis of collision, resulting in remnant black holes with spin. We find modifications to the complex frequencies of the quasi-normal mode spectrum that behave as a power law with the spin of the remnant, and that are not degenerate with the frequencies associated with a Kerr black hole of any mass and spin. We discuss these results in the context of testing general relativity with gravitational wave observations.

Figures

Tables

Table 1. Table 1: Fitted QNM parameters for each head-on collision configuration considered in this study. All configurations have mass ratio q = 1 𝑞 1 q=1 and final background Christodolou mass M f = 0.9896 subscript 𝑀 f 0.9896 M_{\mathrm{f}}=0.9896 . The first column corresponds to the (equal) initial spins of the background black holes, which are oriented in the same direction along the axis of collision (cf. Fig.  1 ). The second column corresponds to the dimensionless spin χ f subscript 𝜒 f \chi_{\mathrm{f}} of the final background black hole. The third, fourth, and fifth columns correspond to the leading-order dCS correction to the QNM frequencies of the ( 2 , 0 , n ) 2 0 𝑛 (2,0,n) modes, ω ( 2 , 0 ) ( 2 ) subscript superscript 𝜔 2 2 0 \omega^{(2)}_{(2,0)} (multiplied by the final background mass, and with the dCS coupling scaled out), for the n = 0 , 1 , 2 𝑛 0 1 2 n=0,1,2 overtones. The sixth, seventh, and eighth column similarly correspond to τ ( 2 , 0 , n ) ( 2 ) subscript superscript 𝜏 2 2 0 𝑛 \tau^{(2)}_{(2,0,n)} , the leading-order dCS correction to the QNM damping times (divided by the final background mass and with the dCS coupling scaled out) for the n = 0 , 1 , 2 𝑛 0 1 2 n=0,1,2 overtones. We detail our sign convention for τ ( 2 ) superscript 𝜏 2 \tau^{(2)} in Sec.  III.4 . We provide a maximum allowed value of ( ℓ / G ​ M ) 4 superscript ℓ 𝐺 𝑀 4 (\ell/GM)^{4} for each configuration (cf. Sec.  IV.2 ) in the last column. In order to be physically meaningful, the dCS QNM parameters must be multiplied by this factor. We have checked that adding max ​ ( ( ℓ / G ​ M ) 4 ) ​ τ ( 2 , 0 , n ) ( 2 ) max superscript ℓ 𝐺 𝑀 4 subscript superscript 𝜏 2 2 0 𝑛 \mathrm{max}((\ell/GM)^{4})\tau^{(2)}_{(2,0,n)} to the GR background solution τ ( 2 , 0 , n ) ( 0 ) subscript superscript 𝜏 0 2 0 𝑛 \tau^{(0)}_{(2,0,n)} does not change the sign (meaning that all modes are decaying).

${\chi }_{1,2}$	${\chi }_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(2,0,0)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(2,0,1)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(2,0,2)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(2,0,0)}^{\left(2\right)}/M_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(2,0,1)}^{\left(2\right)}/M_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(2,0,2)}^{\left(2\right)}/M_f$	max ${\left(\ell /GM\right)}^4$
0.1	0.05106	$6.\left(1\right)\times {10}^{-3}$	$-1.0\left(1\right)\times {10}^{-1}$	$-2.1\left(3\right)\times {10}^0$	$-4.4\left(2\right)\times {10}^0$	$8.\left(1\right)\times {10}^0$	$3.1\left(4\right)\times {10}^1$	$1.26\times {10}^{-1}$
0.2	0.1021	$4.1\left(2\right)\times {10}^{-2}$	$-6.9\left(5\right)\times {10}^{-1}$	$-1.3\left(5\right)\times {10}^1$	$-2.8\left(1\right)\times {10}^1$	$5.4\left(2\right)\times {10}^1$	$1.9\left(1\right)\times {10}^2$	$3.01\times {10}^{-2}$
0.3	0.1532	$1.1\left(5\right)\times {10}^{-1}$	$-1.9\left(5\right)\times {10}^0$	$-3.5\left(2\right)\times {10}^1$	$-7.5\left(7\right)\times {10}^1$	$1.3\left(1\right)\times {10}^2$	$5.1\left(1\right)\times {10}^2$	$1.18\times {10}^{-2}$
0.4	0.2042	$2.6\left(5\right)\times {10}^{-1}$	$-6.1\left(2\right)\times {10}^0$	$-7.0\left(2\right)\times {10}^1$	$-1.4\left(4\right)\times {10}^2$	$3.0\left(1\right)\times {10}^2$	$1.0\left(5\right)\times {10}^3$	$5.68\times {10}^{-3}$
0.5	0.2553	$4.9\left(6\right)\times {10}^{-1}$	$-1.1\left(6\right)\times {10}^1$	$-1.2\left(5\right)\times {10}^2$	$-2.7\left(1\right)\times {10}^2$	$5.6\left(1\right)\times {10}^2$	$1.9\left(1\right)\times {10}^3$	$2.97\times {10}^{-3}$
0.6	0.3062	$8.9\left(2\right)\times {10}^{-1}$	$-2.1\left(1\right)\times {10}^1$	$-2.2\left(1\right)\times {10}^2$	$-4.8\left(3\right)\times {10}^2$	$9.7\left(2\right)\times {10}^2$	$3.4\left(2\right)\times {10}^3$	$1.61\times {10}^{-3}$
0.7	0.3574	$1.5\left(2\right)\times {10}^0$	$-4.1\left(1\right)\times {10}^1$	$-3.8\left(1\right)\times {10}^2$	$-8.2\left(4\right)\times {10}^2$	$1.6\left(2\right)\times {10}^3$	$5.7\left(3\right)\times {10}^3$	$8.79\times {10}^{-4}$

Table 2. Table 2: Same as Table  1 , but for the ( 4 , 0 , n ) 4 0 𝑛 (4,0,n) modes.

${\chi }_{1,2}$	${\chi }_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(4,0,0)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(4,0,1)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\omega }_{(4,0,2)}^{\left(2\right)}{M}_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(4,0,0)}^{\left(2\right)}/M_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(4,0,1)}^{\left(2\right)}/M_f$	${\left(\ell /GM\right)}^{-4}$ ${\tau }_{(4,0,2)}^{\left(2\right)}/M_f$	max ${\left(\ell /GM\right)}^4$
0.1	0.05106	$-1.2\left(5\right)\times {10}^{-2}$	$4.\left(2\right)\times {10}^{-2}$	$8.5\left(1\right)\times {10}^{-2}$	$1.\left(3\right)\times {10}^0$	$-7.\left(1\right)\times {10}^{-1}$	$-6.\left(3\right)\times {10}^{-1}$	$1.26\times {10}^{-1}$
0.2	0.1021	$-6.\left(2\right)\times {10}^{-2}$	$2.\left(2\right)\times {10}^{-1}$	$4.8\left(8\right)\times {10}^{-1}$	$1.\left(1\right)\times {10}^1$	$-5.\left(1\right)\times {10}^0$	$-4.2\left(7\right)\times {10}^0$	$3.01\times {10}^{-2}$
0.3	0.1532	$-1.5\left(4\right)\times {10}^{-1}$	$5.\left(5\right)\times {10}^{-1}$	$1.2\left(2\right)\times {10}^0$	$3.\left(4\right)\times {10}^1$	$-1.2\left(4\right)\times {10}^1$	$-1.0\left(3\right)\times {10}^1$	$1.18\times {10}^{-2}$
0.4	0.2042	$-3.\left(1\right)\times {10}^{-1}$	$1.\left(1\right)\times {10}^0$	$2.5\left(6\right)\times {10}^0$	$8.\left(9\right)\times {10}^1$	$-2.8\left(1\right)\times {10}^1$	$-2.2\left(6\right)\times {10}^1$	$5.68\times {10}^{-3}$
0.5	0.2553	$-5.\left(2\right)\times {10}^{-1}$	$1.\left(2\right)\times {10}^0$	$4.\left(1\right)\times {10}^0$	$1.\left(1\right)\times {10}^2$	$-5.6\left(2\right)\times {10}^1$	$-4.\left(1\right)\times {10}^1$	$2.97\times {10}^{-3}$
0.6	0.3062	$-9.\left(6\right)\times {10}^{-1}$	$2.\left(3\right)\times {10}^0$	$7.\left(2\right)\times {10}^0$	$3.\left(3\right)\times {10}^2$	$-1.0\left(1\right)\times {10}^2$	$-8.\left(1\right)\times {10}^1$	$1.61\times {10}^{-3}$
0.7	0.3574	$-1.4\left(7\right)\times {10}^0$	$2.\left(8\right)\times {10}^0$	$1.2\left(5\right)\times {10}^1$	$8.\left(5\right)\times {10}^2$	$-1.8\left(1\right)\times {10}^2$	$-1.6\left(3\right)\times {10}^2$	$8.79\times {10}^{-4}$

Equations

$$S\equiv \int d^{4}x\sqrt{-g}\left(\frac{m_{\text{pl}}^2}{2}R-\frac{1}{2}(\partial \vartheta {)}^2-\frac{m_{\text{pl}}}{8}{\ell }^2\vartheta {}^{\ast }RR\right).$$

$${}^{\ast }RR\equiv {}^{\ast }{R}^{abcd}{R}_{abcd},$$

$$\square \vartheta \equiv {\nabla }_a{\nabla }^a\vartheta =\frac{m_{\text{pl}}}{8}{\ell }^2{}^{\ast }RR,$$

$$m_{\text{pl}}^2{G}_{ab}+m_{\text{pl}}{\ell }^2{C}_{ab}=T_{ab}^{\vartheta },$$

$$C_{ab}\equiv {ϵ}_{cde(a}{\nabla }^d{R}_{b)}{}^c{\nabla }^e\vartheta +{}^{\ast }{R}^c{}_{(ab)}{}^d{\nabla }_c{\nabla }_d\vartheta ,$$

$$T_{ab}^{\vartheta }$$

$$g_{ab}$$

$$\vartheta$$

$${\square }^{(0)}{\vartheta }^{(1)}=\frac{m_{\text{pl}}}{8}{\ell }^2{}^{\ast }R{R}^{(0)},$$

$$m_{\text{pl}}^2{G}_{ab}^{(0)}[h_{ab}^{(2)}]=-m_{\text{pl}}{\ell }^2{C}_{ab}^{(1)}+T_{ab}^{(\vartheta (1))},$$

$$T_{ab}^{(\vartheta (1))}$$

$$C_{ab}^{(1)}$$

$$+{}^{\ast }{R}^c{}_{(ab)}{}^d{}^{(0)}{\nabla }_c{}^{(0)}{\nabla }_d{}^{(0)}{\vartheta }^{(1)}.$$

$$h_{ab}^{(2)}\equiv \frac{(\ell /GM{)}^4}{8}\Delta g_{ab},{\vartheta }^{(1)}\equiv \frac{m_{\text{pl}}}{8}(\ell /GM{)}^2\Delta \vartheta$$

$${\square }^{(0)}\Delta \vartheta ={}^{\ast }R{R}^{(0)}.$$

$$G_{ab}^{(0)}[\Delta g_{ab}]=-C_{ab}^{(1)}[\Delta \vartheta ]+\frac{1}{8}{T}_{ab}^{(1)}[\Delta \vartheta ].$$

$${\Psi }_4={\Psi }_4^{(0)}+\sum _{k=1}^{\infty }{\epsilon }^k{\Psi }_4^{(k)}.$$

$${\Psi }_4^{(2)}=\frac{(\ell /GM{)}^4}{8}\Delta {\Psi }_4,$$

$$r{\Psi }_4(t)=\sum _{l,m,n}{\tilde{A}}_{(l,m,n)}{e}^{-i{\tilde{\omega }}_{(l,m,n)}t}{}_{-2}{Y}_{(l,m,n)}^{a\omega }.$$

$$\tilde{\omega }=\omega -i/\tau .$$

$$\tilde{A}=A{e}^{i\theta }.$$

$$r{\Psi }_4$$

$$r{\Psi }_4^{(0)}$$

$$+i{A}^{(0)}\sin (-{\omega }^{(0)}t+{\theta }^{(0)})e^{-t/{\tau }^{(0)}}.$$

$$\omega$$

$$A$$

$$\mathrm{Re}(r{\Psi }_4^{(2)})$$

$$-{\theta }^{(2)}{A}^{(0)}\sin (-{\omega }^{(0)}t+{\theta }^{(0)})e^{-t/{\tau }^{(0)}}$$

$$+t{\omega }^{(2)}{A}^{(0)}\sin (-{\omega }^{(0)}t+{\theta }^{(0)})e^{-t/{\tau }^{(0)}}$$

$$+t\frac{{\tau }^{(2)}}{({\tau }^{(0)}{)}^2}{A}^{(0)}\cos (-{\omega }^{(0)}t+{\theta }^{(0)})e^{-t/{\tau }^{(0)}},$$

$$s$$

$$c$$

$$\mathrm{Re}(r{\Psi }_4^{(2)})$$

$$\displaystyle\quad+tA^{(0)}\left(\omega^{(2)}s+\frac{\tau^{(2)}}{(\tau^{(0)})^{2}}c\right)\Big{]}\,.$$

$$\mathrm{Im}(r{\Psi }_4^{(2)})$$

Full text

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, 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, 2018, 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, 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.

I.1 Roadmap and conventions

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

II.1 Order-reduced dynamical Chern-Simons gravity

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, 2018, 2017). Here we only briefly summarize.

The full dCS action takes the form

[TABLE]

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,

[TABLE]

where $\,{}^{*}\!R^{abcd}=\frac{1}{2}\epsilon^{abef}R_{ef}{}^{cd}$ is the dual of the Riemann tensor, and $\epsilon_{abcd}\equiv\sqrt{-g}[abcd]$ is the fully antisymmetric Levi-Civita tensor Misner et al. (1973). This interaction is governed by a coupling constant $\ell$, which has dimensions of length, and physically represents the length scale below which dCS corrections become relevant.

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

[TABLE]

and

[TABLE]

where

[TABLE]

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

[TABLE]

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.

$C_{ab}$, as given in Eq. (5), contains third derivatives of the metric, and it is thus unknown whether dCS has a well-posed initial value formulation Delsate et al. (2015). But as is typical in the modern treatment of beyond-GR theories of gravity, we assume that dCS is a low-energy effective field theory (EFT) of some UV-complete theory that is well-posed. Therefore we work instead in well-posed order-reduced dCS, in which we perturb the metric and scalar field about an arbitrary GR spacetime, and obtain perturbed equations of motion. In particular, we introduce a dimensionless formal order-counting parameter $\varepsilon$ which keeps track of powers of $\ell^{2}$ (this formal order-counting parameter can later be set to one). We then write

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

II.1.1 Scaled variables

Because the evolution equations at some order $\varepsilon^{k}$ are homogeneous in $\ell^{2k}$, we can scale out the $\ell$ dependence by defining new variables

[TABLE]

(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

[TABLE]

Eq. (10) similarly becomes

[TABLE]

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)}$.

We thus need to solve Eqs. (14) and (15) only once for each BBH background configuration, and then multiply our results for $\Delta g_{ab}$ and $\Delta\vartheta$ by appropriate powers of $\ell/GM$ and factors of 8 afterward.

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, 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

[TABLE]

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.

In practice, we write [cf. Eq. (13)]

[TABLE]

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.

III.1 Quasi-normal modes in general relativity

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

[TABLE]

Here, $l$ and $m$ label the spin-weight $-2$ spheroidal harmonic under consideration, while $n$ refers to the overtone, ordered by largest damping time. The quantities $\tilde{A}$ and $\tilde{\omega}$ are the complex amplitude and frequency of the $(l,m,n)$ mode under consideration. For simplicity, we will henceforth omit the $(l,m,n)$ indices on $\tilde{\omega}$ and $\tilde{A}$, and consider each mode separately. We can write $\tilde{\omega}$ in terms of a real frequency, $\omega$, and a damping time, $\tau$, to give

[TABLE]

Let us similarly write

[TABLE]

where $A\equiv|\tilde{A}|$ is the norm of $\tilde{A}$, and $\theta$ is the complex phase of $\tilde{A}$. Then we obtain, for a single mode,

[TABLE]

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:

[TABLE]

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

[TABLE]

Recall that the leading-order correction to the gravitational radiation is $r\Psi_{4}^{(2)}$, which is linear in $h^{(2)}_{ab}$ and has a coupling factor of $(\ell/GM)^{4}$. Thus, the leading-order correction to $\omega$ will be $\omega^{(2)}$, as computed from a linearization of Eq. (18), with a coupling factor of $(\ell/GM)^{4}$. Each mode of $r\Psi_{4}^{(2)}$, the leading-order dCS correction to the gravitational radiation, will thus be parametrized by and linear in $\{\omega^{(2)},\tau^{(2)},A^{(2)},\theta{}^{(2)}\}$.

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

[TABLE]

which can be more compactly written as

[TABLE]

The imaginary part is similarly modified as

[TABLE]

We thus see both an amplitude modification to the background QNM spectrum from the $A^{(2)}$ and $\theta^{(2)}$ terms, and a modification linear in time from the $\omega^{(2)}$ and $\tau^{(2)}$ terms. We fit precisely the functional form in Eqs. (28) and (29) to the $r\Psi_{4}^{(2)}$ obtained from the simulation. Our fit determines four free parameters for each mode: $\{A^{(2)},\theta^{(2)},\omega^{(2)},\tau^{(2)}\}$; the other free parameters $\{A^{(0)},\theta^{(0)}\}$ in Eqs. (28) and  (29) are determined by the fit to $r\Psi_{4}^{(0)}$ using Eq. (22). Note that this is different from simply fitting a damped sinusoid to $r\Psi_{4}^{(2)}$. The presence of a term that behaves as $\sim te^{-i\tilde{\omega}t}$ beside the ordinary damped sinusoids is a sign of secular breakdown, which is discussed in more detail in Sec. IV.2.2.

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

Because the simulations (cf. Sec. II.1.1) are independent of the coupling parameter $\ell/GM$, the resulting waveforms for $\vartheta^{(1)}$ and $r\Psi_{4}^{(2)}$ have the coupling scaled out. We will thus report our results as

[TABLE]

and so on.

Much of the QNM literature reports $\tilde{\omega}$ in terms of its real and imaginary parts, $\tilde{\omega}=\mathrm{Re}(\omega)+i\mathrm{Im}(\omega)$, without invoking a damping time $\tau$. We can transform our results for $\tau$ into $\mathrm{Im}(\omega)$ as

[TABLE]

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

[TABLE]

III.5 Mass and spin definitions

Since $\tau^{(0)}$ and $\omega^{(0)}$ are (by the no-hair theorem) dependent only on $M$ and $\chi$, the mass and dimensionless spin of the final black hole in GR, we should similarly expect $\omega^{(2)}$ and $\tau^{(2)}$ to be dependent on some final mass and final spin. In the full dCS theory, we expect the mass and spin of a dCS black hole to be modified with respect to those of a GR black hole (recall that Kerr is not a solution of dCS Yunes and Pretorius (2009)). The formulae used to compute mass and spin, because they are derived using properties of GR (cf. Baumgarte and Shapiro (2010)), may themselves be modified in the full dCS theory. If we had access to the full theory, we could parametrize the QNM spectra in terms of $\chi_{\mathrm{dCS}}$ and $M_{\mathrm{dCS}}$, as well as $\ell/GM$. Since we are working in an order-reduction scheme, we can instead linearize the formulae used to compute the spin and mass of the final background black hole, and compute the corrected mass and spin. In this study, however, we choose to parametrize the QNM spectra in terms of the Christodoulou mass and dimensionless spin of the final background black hole, which we will call $M_{\mathrm{f}}$ and $\chi_{\mathrm{f}}$.

III.6 Fitting window

When fitting $r\Psi_{4}^{(0)}$ and $r\Psi_{4}^{(2)}$, we must be careful about the time window of the post-merger waveform used for the fit. For $r\Psi_{4}^{(0)}$, if we choose a starting time $t_{\mathrm{start}}$ too close to merger, then our assumption that each mode can be fit by a function of the form in Eq. (18) breaks down Bhagwat et al. (2018). The later we choose $t_{\mathrm{start}}$, the less data (not contaminated by numerical noise) are available to perform the fit. However, in Giesler et al. (2019), the authors found that, when including overtones, the post-merger spectrum could be fit with QNMs as early as the peak of the gravitational waveform. We similarly fit enough overtones so that we can faithfully choose $t_{\mathrm{start}}$ to be the peak of $r\Psi_{4}^{(0)}$.111Reference Giesler et al. (2019) uses the peak of the gravitational waveform strain as the start of the fit, while we use the peak of $r\Psi_{4}^{(0)}$.

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

IV.1 Waveforms

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}$.222Note 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.

We show the dominant modes of $r\Psi_{4}^{(0)}$ for a representational case with $\chi=0.1\hat{x}$ in Fig. 2. We similarly show the dominant modes of $r\Psi_{4}^{(2)}$ for this configuration in Fig. 3. Recall that the physical gravitational radiation includes a coupling factor $(\ell/GM)^{4}$, which is scaled out in the numerical computation, and thus we report the waveforms as $(\ell/GM)^{-4}r\Psi_{4}^{(2)}$. This configuration is axisymmetric about the $\hat{x}$ axis, and thus we expect only $m=0$ modes to be excited. Since the spins have the same orientation, there is reflection symmetry about the $\hat{y}-\hat{z}$ plane, so we expect only the $l=\mathrm{even}$ modes to be excited.

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

IV.2.1 Instantaneous regime of validity

As discussed in Sec. II.1.1, the leading-order scalar field $\Delta\vartheta$ and metric perturbation $\Delta g_{ab}$ as computed from the code are independent of the coupling constant $\ell/GM$. In order to make our results physically meaningful, we must multiply the leading-order scalar field by $(\ell/GM)^{2}$ and the leading-order metric correction by $(\ell/GM)^{4}$. Similarly, we must multiply the computed leading-order correction to the gravitational radiation, $\Delta\Psi$, by a factor of $(\ell/GM)^{4}$.

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:

[TABLE]

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

[TABLE]

and hence

[TABLE]

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

We show the regime of validity for a $\vec{\chi}=0.7\hat{x}$ head-on collision in Fig. 5. $\ell/GM$ takes its smallest allowed value in the strong-field region, outside the apparent horizon of each black hole. We see that closer to merger, where there is power in the metric perturbation, the maximal allowed value of $\ell/GM$ decreases. After merger, the maximal allowed value of $\ell/GM$ increases as the dCS metric perturbation partially radiates away, and the final constant value is governed by the strength of the dCS metric perturbation around the final black hole.

We show the behavior of the minimum allowed value of $\ell/GM$, over the entire simulation, as a function of final dimensionless spin $\chi_{\mathrm{f}}$ in Fig. 6. The regime of validity decreases with spin, as the magnitude of $\Delta g_{ab}$ increases with spin. This scaling serves as a proxy for the allowed values of $\ell/GM$ when considering gravitational waveforms.

IV.2.2 Secular regime of validity

In addition to an instantaneous regime of validity, our computations will also have a secular regime of validity. This arises from our perturbative scheme and is discussed in detail in Okounkova et al. (2017). In order for the perturbation scheme to remain valid, we want the perturbation that is going as $\sim te^{-i\tilde{\omega}^{(0)}t}$ in Eqs. (28) and (29) to remain small relative the background part $\sim e^{-i\tilde{\omega}^{(0)}t}$ in Eq. (22). This gives us the condition

[TABLE]

where $|\tilde{\omega}^{(2)}|$ denotes the norm of the leading-order correction to the complex frequency [cf. Eqs. (19) and (33)]. This condition bounds how long the perturbation is valid by the secular time $t_{\mathrm{sec}}$, roughly

[TABLE]

The (physical) value of $\ell/GM$ thus determines the time window over which the perturbative scheme is appropriate. Since this value is the same for all modes, we must take the most conservative estimate, minimizing the RHS of Eq. (38) over all modes.

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.

We show representative fits in Fig. 7. We find that for head-on collisions, we can most successfully fit each mode from the peak of the waveform using three overtones. We set $t_{\mathrm{end}}=25\,M$ after the peak. We take a closer look at this fit in Fig. 8, where we give an illustration of the linear-in-time behavior of $r\Psi_{4}^{(2)}$ predicted by Eqs. (28) and (29). We divide the terms of the form $t\cos()$ and $t\sin()$ by $\cos()$ and $\sin()$ respectively, finding that the result is indeed linear with time.

For each of the configurations, the final Christodolou mass of the background spacetime is $M_{\mathrm{f}}=0.9896$. However, $\chi_{\mathrm{f}}$, the final background spin, varies with configuration, and we include these values in Tables 1 and  2.

We plot the values of $\tau^{(2)}_{(2,0,0)}(\ell/GM)^{-4}$, the leading-order dCS correction to the damping time of the least-damped $(2,0)$ mode of the gravitational radiation, and $\omega^{(2)}_{(2,0,0)}(\ell/GM)^{-4}$, the leading-order dCS correction to the frequency, as functions of $\chi_{\mathrm{f}}$ in Figs. 9 and 10. We see that $\tau^{(2)}(\ell/GM)^{-4}$ and $\omega^{(2)}(\ell/GM)^{-4}$ behave as a power law with spin. This behavior can be expected by considering analytical results in dCS theory. In the slow-rotation approximation, the horizon (and hence the light ring) is modified at quadratic order in spin Yagi et al. (2012), while containing no modifications at first order in spin Yunes and Pretorius (2009). We additionally plot $\tau^{(2)}(\ell/GM)^{-4}$ and $\omega^{(2)}(\ell/GM)^{-4}$ for the $(4,0,0)$ mode in Figs. 11 and 12. Again, we see that these quantities behave as a power law with spin.

Let us think about the signs of $\tau^{(2)}_{(l,m,n)}(\ell/GM)^{-4}$. In the GR case, given the definition of $\tau^{(0)}$ in Eq. (19), we always expect $\tau^{(0)}_{(l,m,n)}>0$ in order for the modes to be exponentially damped (as opposed to exponentially growing). The dCS modifications to the damping time can either be positive or negative. In Fig. 9, for example, we see that for the dominant mode, $(\ell/GM)^{-4}\tau^{(2)}_{(2,0,0)}<0$, meaning that the damping time is decreased. We can physically think about this as the dCS scalar field removing additional energy from the system, leading to faster damping. We see in Fig. 11, however, that $(\ell/GM)^{-4}\tau^{(2)}_{(4,0,0)}>0$.

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

[TABLE]

Note that the dCS coupling $(\ell/GM)^{4}$ enters into the above expression. Thus, if we multiply our numerical results for $(\ell/GM)^{-4}\tau^{(2)}_{(l,m,n)}$ (no matter what signs they carry) by the maximum allowed value of $(\ell/GM)^{4}$ to obtain $\tau^{(2)}_{(l,m,n)}$, Eq. (39) must hold. We have checked that this is indeed the case (cf. Tables 1 and  2).

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, 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.

Let us first suppose that $\ell=0$, meaning that there is no modification from GR. In GR, assuming the no-hair theorem is true, the frequency and damping time for each mode of ringdown should be parametrized purely by the mass, $M_{\mathrm{f}}$, and spin, $\chi_{\mathrm{f}}$, of the final black hole. Given two observed modes, we have access to four quantities — two frequencies and two damping times. We can check that the fitted $\omega_{(l,m,n)}$ and $\tau_{(l,m,n)}$ are consistent with the predicted GR values for $M_{\mathrm{f}}$ and $\chi_{\mathrm{f}}$ Berti et al. (2015); Gossan et al. (2012); Abbott et al. (2016); Kamaretsos et al. (2012); Yunes et al. (2016).

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.1 Checking non-degeneracy: projected

Let us consider whether it is possible in principle to measure dCS corrections $\omega_{(l,m,n)}^{(2)}$ and $\tau_{(l,m,n)}^{(2)}$ from a given waveform signal, or instead whether these quantities are degenerate with the GR values $\omega_{(l,m,n)}^{(0)}$ and $\tau_{(l,m,n)}^{(0)}$ corresponding to a different remnant mass and spin. Consider the 4-dimensional parameter space $\mathbb{P}$ of $\{M\omega^{(0)},\tau^{(0)}/M,M\omega^{(2)},\tau^{(2)}/M\}$ for a given mode. GR solutions exist purely in the 2-dimensional submanifold $\mathbb{S}_{\mathrm{GR}}$ specified by $\{\omega^{(2)}=0,\tau^{(2)}=0\}$, with coordinates $\{M_{\mathrm{f}},\chi_{\mathrm{f}}\}$ on the manifold. Suppose at some point $(M_{1},\chi_{1})$ on $\mathbb{S}_{\mathrm{GR}}$, we introduce a dCS deviation with some coupling $(\ell/GM)^{4}$. In other words, we will have

[TABLE]

where we have explicitly written out the dependence on the coupling constant with $\Delta\omega\equiv(\ell/GM)^{-4}\omega^{(2)}$ and $\Delta\tau\equiv(\ell/GM)^{-4}\tau^{(2)}$, which can be read off of Figs. 9–12. If dCS modifications and GR are degenerate, then this modification will move purely within $\mathbb{S}_{\mathrm{GR}}$. However, if dCS modifications and GR are non-degenerate, then the new point will be off $\mathbb{S}_{\mathrm{GR}}$ in $\mathbb{P}$ and the dCS modifications will form a 3-dimensional submanifold of $\mathbb{P}$, $\mathbb{S}_{\mathrm{dCS}}$, with coordinates $\{M_{\mathrm{f}},\chi_{\mathrm{f}},\ell/GM\}$.

Let us now consider this statement in the context of our numerical results. For simplicity, let us first consider holding $M_{\mathrm{f}}$ constant in the comparisons. In Fig. 14, we plot values of $M_{\mathrm{f}}\omega_{(2,0,0)}$ and $\tau_{(2,0,0)}/M_{\mathrm{f}}$ for various spins. We similarly plot $M_{\mathrm{f}}\omega_{(4,0,0)}$ and $\tau_{(4,0,0)}/M_{\mathrm{f}}$ in Fig. 15. For $\ell=0$, we can use perturbation theory to compute the values of $\omega=\omega^{(0)}$ and $\tau=\tau^{(0)}$ in GR. Holding $M_{\mathrm{f}}$ fixed, the GR solutions form a curve $L_{\mathrm{GR}}$ in the $\omega-\tau$ plane, parametrized by $\chi_{\mathrm{f}}$.

Now let us introduce $\ell\neq 0$. For each simulation that we have performed, with a given $\chi_{\mathrm{f}}$ (recall all of the $M_{\mathrm{f}}$ are equal), we compute $\omega(\chi_{\mathrm{f}})$ and $\tau(\chi_{\mathrm{f}})$ via Eqs. (40) and (41) using our results for $\Delta\omega=(\ell/GM)^{-4}\omega^{(2)}$ and $\Delta\tau=(\ell/GM)^{-4}\tau^{(2)}$. This computation requires specifying a value of $(\ell/GM)$. If we vary $(\ell/GM)$ over an allowed range (cf. Sec. IV.2), for each $\chi_{\mathrm{f}}$ we obtain a curve $L_{\mathrm{dCS}}(\chi_{\mathrm{f}})$ in the $M_{\mathrm{f}}\omega-\tau/M_{\mathrm{f}}$ plane parametrized by $(\ell/GM)$.

If the dCS corrections to the quasi-normal mode spectrum were degenerate with GR, then $L_{\mathrm{dCS}}(\chi_{\mathrm{f}})$ would lie purely along $L_{\mathrm{GR}}$. In other words, the resulting QNM spectrum for $\chi_{\mathrm{f}}$ would be degenerate with that of GR for some other spin $\chi^{\prime}$. However, we see in Figs. 14 and 15 that in all cases $L_{\mathrm{dCS}}(\chi_{\mathrm{f}})$ does not lie purely along $L_{\mathrm{GR}}$, meaning that the QNM spectra are non-degenerate. This in turn means that dCS modifications to QNM spectra can in principle be observed (in the limit of infinite signal-to-noise ratio). Note that in this analysis, we have held $M_{\mathrm{f}}$ fixed, since all of the head-on collision simulations in this paper have the same final black hole mass.

V.2 Checking non-degeneracy: full case

We can perform a more rigorous analysis, checking for full degeneracy, rather than the simpler check that holds $M_{\mathrm{f}}$ fixed. Let us think about the 3-dimensional space $\mathbb{F}$ with coordinates $\{\chi,M,\varepsilon^{2}\}$ (where $\varepsilon$ is our dCS order-reduction parameter). Suppose we observe $k$ QNMs, which gives us $2k$ quantities ($\omega$ and $\tau$ for each mode). Let $\mathbb{Q}$ be the $2k$-dimensional space with these coordinates.

Let us consider the map $\phi:\mathbb{F}\to\mathbb{Q}$, which maps each set of parameters $\{\chi,M,\varepsilon^{2}\}$ to the QNM values. The image $\phi(\mathbb{F})$ will form a 3-dimensional submanifold of $\mathbb{Q}$, and the tangent space of the image will be spanned by the pushforwards of $\{\partial/\partial\chi,\partial/\partial M,\partial/\partial\varepsilon^{2}\}$. That is, $\{\phi_{*}\partial/\partial\chi,\phi_{*}\partial/\partial M,\phi_{*}\partial/\partial\varepsilon^{2}\}$.

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

[TABLE]

Let us consider how to evaluate this matrix, working near $\varepsilon^{2}=0$ for each $\chi_{\mathrm{f}}$ and $M_{\mathrm{f}}$ for which we have performed a head-on collision. Suppose we are considering some mode with QNM frequency $\omega_{lmn}$ and damping time $\tau_{lmn}$. This will give us two rows in the matrix $J$.

Let us first compute

[TABLE]

This can be done by computing the values of $\omega_{lmn}M_{\mathrm{f}}$ and $\tau_{lmn}/M_{\mathrm{f}}$ from perturbation theory Stein (2019), varying only $\chi$ around $\chi_{\mathrm{f}}$, and then taking a numerical derivative. We work with a step-size of $10^{-10}$, which is the default precision of Stein (2019).

Now let us compute the $\partial/\partial M$ column. For fixed $\chi_{\mathrm{f}}$ and $\ell^{4}=0$, the dependence on $M$ is

[TABLE]

We evaluate these expressions at $M=M_{\mathrm{f}}$.

Finally, for the last column, for fixed $\chi_{\mathrm{f}},M_{\mathrm{f}}$, we have

[TABLE]

where $(\ell/GM)^{-4}M_{f}\omega_{lmn}^{(2)}$ and $(\ell/GM)^{-4}\tau_{lmn}^{(2)}/M_{f}$ are the quantities we compute from our numerical fits.

We put the matrix $J$ together with these values. Note that the $\omega$ rows all have a factor of $1/M_{\mathrm{f}}$, while the $\tau$ rows have a factor of $M_{\mathrm{f}}$. Since each row is homogeneous in a power of $M_{\mathrm{f}}$, we can divide through without changing the rank of the matrix. We evaluate the rank of this matrix using a singular-value decomposition (SVD) Jones et al. (01). For all values of $\{\chi_{\mathrm{f}},M_{\mathrm{f}}\}$ in our head-on collisions dataset, we find that the rank of $J$ is 3. The lowest singular value is $10^{-2}-10^{-1}$, while the condition numbers (the 2-norm, computed from the SVD) are of order $10^{3}$. Therefore, we conclude that dCS corrections to ringdowns have a distinct observational signature, as they are not degenerate with changing GR parameters of final mass and spin.

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.

A.1 QNMs in GR

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

[TABLE]

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,

[TABLE]

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, 1973), which we will denote as

[TABLE]

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

[TABLE]

Here $\tilde{A}_{lmn}$ is a complex mode amplitude with complex frequency $\tilde{\omega}_{lmn}$. Here $l,m$ label the angular harmonics, and $n\geq 0$ labels the QNMs with the same $lm$ in terms of increasing damping time, related to the imaginary part of the complex frequency (henceforth we suppress the $lmn$ labels for brevity). The damping time $\tau>0$ is found via

[TABLE]

The functions $e^{im\phi}\,{}_{-2}S_{lm}^{a\tilde{\omega}}(\theta)$ are spin $-2$-weighted spheroidal harmonics. However our numerical simulations compute $\Psi_{4}^{(0)}$ and $\Psi_{4}^{(2)}$ decomposed into spin-weighted spherical harmonics, where the spheroidals are a deformation of the sphericals. The spherical-spheroidal mixing coefficients can be computed Hughes (2000); Cook and Zalutskiy (2014); Berti et al. (2015); Stein (2019), and the amount of mixing is small at low spins and frequencies, so we do not model the mixing in this paper.

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

[TABLE]

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

[TABLE]

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

[TABLE]

At the level of the “full” equations, QNMs are once again homogeneous solutions to this linear PDE. Since the background is axisymmetric, the perturbations will again decompose into azimuthal modes proportional to $e^{im\phi}$, labeled by $m$. Since the background is stationary, the perturbations will also decompose into modes proportional to $e^{-i\tilde{\omega}t}$, labeled by complex frequencies $\tilde{\omega}$, with $\text{Im}[\tilde{\omega}]<0$ if we assume the nonlinear BH solutions are stable.

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, 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

[TABLE]

where $\lambda$ labels the modes.

A.3 Perturbative treatment of QNMs beyond GR

We can now ask how these quantities are deformed away from their GR values. Notice that in Eq. (57), the frequencies $\tilde{\omega}(\varepsilon)$ and the amplitude of each mode $\tilde{A}_{\lambda}(\varepsilon)$ should depend continuously on $\varepsilon$. Some modes $\lambda$ will reduce to GR modes $lmn$ in the limit as $\varepsilon\to 0$, and their frequencies and amplitudes will reduce to the appropriate GR quantities. Additional modes (e.g. the “scalar-led” modes of Blázquez-Salcedo et al. (2016, 2017)) will have amplitudes that vanish at $\varepsilon=0$.

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

[TABLE]

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$,

[TABLE]

and similarly for $S_{lmn}(\varepsilon,\theta)$. The quantities with superscript ${(0)}$ are the GR expression appearing in Eq. (52). The real and imaginary parts of $\tilde{\omega}^{(k)}$ follow the same expansion, but we also use $\tau\equiv-1/\text{Im}[\tilde{\omega}]$. Similarly expanding $\tau=\sum_{k}\varepsilon^{k}\tau^{(k)}$ and using the chain rule gives

[TABLE]

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

[TABLE]

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.

Let us make a few notes about the functional form of $r\Psi_{4}^{(1)}$ in Eq. (64). First, although the “full” QNM frequency $\tilde{\omega}(\varepsilon)$ at finite $\varepsilon$ is shifted relative to $\tilde{\omega}^{(0)}$, everything in Eq. (64) is proportional to $e^{-i\tilde{\omega}^{(0)}u}$, where $u\sim t-r$ is retarded time. Second, Eq. (64) is not just a sum of damped sinusoids, because of the term going as $\sim\tilde{\omega}^{(1)}\tilde{A}^{(0)}ue^{-i\tilde{\omega}^{(0)}u}$. Since $\text{Im}[\tilde{\omega}^{(0)}]<0$, this term remains bounded, but it grows in importance relative to the term $\tilde{A}^{(1)}e^{-i\tilde{\omega}^{(0)}u}$. This is a typical symptom of a secular breakdown of perturbation theory Bender and Orszag (1978); Kevorkian and Cole (1996).

A.4 Particular, homogeneous modes, and

modes that are mixed in field space

How does the new term $\sim ue^{-i\tilde{\omega}^{(0)}u}$ arise, which differs in form from both the GR and “exact” beyond-GR behaviors? To understand we have to look at the perturbative treatment of the equations of motion. We can model the metric sector of $\mathcal{F}$ as some deformation of Einstein’s equations,

[TABLE]

For example in dCS, $\varepsilon H_{ab}=m_{\textrm{\tiny{pl}}}\ell^{2}C_{ab}-T_{ab}^{\vartheta}$ which is divergence-free on shell. Inserting Eq. (58) into Eq. (65) and expanding will give order-reduced equations for the deformations. For example, setting $\xi=0$ reduces to the BH background with no QNMs. The deformation $h^{{\textrm{\tiny{Def}}}}$ satisfies the inhomogeneous equation

[TABLE]

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

[TABLE]

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.

This equation is different from the linearized EOM in GR (50) or in the “exact” theory (56): those are both homogeneous, whereas Eq. (67) is inhomogeneous (it has a source term on the RHS). The solution will therefore have two parts: a homogeneous solution, and a particular solution. The linear differential operator $G^{(1)}_{ab}[h_{cd}^{Q\,(1)}]$ appearing on the LHS of Eq. (67) is the same operator, expanded about the same Kerr background, as in the GR Eq. (50). Therefore the homogeneous solutions are the same. This accounts for the $e^{-i\tilde{\omega}^{(0)}u}$ terms in Eq. (64). The $ue^{-i\tilde{\omega}^{(0)}u}$ terms are particular solutions. Their secular behavior is due to the source term having support at the poles of the Green’s function.

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.333Working 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)),

[TABLE]

where $r_{*}$ is the usual Schwarzschild “tortoise” radial coordinate, $\vartheta$ is a spherical harmonic mode of the full dCS scalar field, $\Psi$ is a spherical harmonic mode of the metric, and the $V_{ij}$ potentials are functions dependent on $r$ and the dCS coupling parameter. Since $\vartheta$ and $\Psi$ are coupled in this 2-dimensional linear ODE problem, we find two types of QNMs (scalar-led and gravitational-led, which become diagonal in field space in the limit $\ell\to 0$).

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

[TABLE]

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

[TABLE]

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,

[TABLE]

must be satisfied.

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

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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

Perturbing Eq. (85), we obtain

[TABLE]

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

[TABLE]

and

[TABLE]

We can then perturb the projection operators,

[TABLE]

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}$.