Numerical binary black hole mergers in dynamical Chern-Simons: I. Scalar field

TL;DR

This paper develops a numerical scheme to simulate binary black hole mergers in dynamical Chern-Simons gravity, revealing new phenomenology and estimating bounds on the theory's parameters from gravitational wave observations.

Contribution

It introduces a well-posed perturbation scheme for beyond-GR theories and applies it to simulate BBH mergers in dynamical Chern-Simons gravity, uncovering new merger phenomenology.

Findings

Good agreement with analytic predictions at early times

Discovery of a burst of dipole radiation during merger

Estimated bounds on dCS length scale from LIGO data

Abstract

Testing general relativity in the non-linear, dynamical, strong-field regime of gravity is one of the major goals of gravitational wave astrophysics. Performing precision tests of general relativity (GR) requires numerical inspiral, merger, and ringdown waveforms for binary black hole (BBH) systems in theories beyond GR. Currently, GR and scalar-tensor gravity are the only theories amenable to numerical simulations. In this article, we present a well-posed perturbation scheme for numerically integrating beyond-GR theories that have a continuous limit to GR. We demonstrate this scheme by simulating BBH mergers in dynamical Chern-Simons gravity (dCS), to linear order in the perturbation parameter. We present mode waveforms and energy fluxes of the dCS pseudoscalar field from our numerical simulations. We find good agreement with analytic predictions at early times, including the absence of pseudoscalar dipole radiation. We discover new phenomenology only accessible through numerics: a burst of dipole radiation during merger. We also quantify the self-consistency of the perturbation scheme. Finally, we estimate bounds that GR-consistent LIGO detections could place on the new dCS length scale, approximately $\ell \lesssim \mathcal{O}(10)~\mathrm{km}$.

Figures

Tables

Table 1. Table 1: Parameters of numerical runs. Each run was performed at low, medium, and high resolutions. We give the mass ratio m 1 / m 2 subscript 𝑚 1 subscript 𝑚 2 m_{1}/m_{2} where the subscripts label the black holes. All of the spins are aligned in the z 𝑧 z -direction, so we give the z ^ ^ 𝑧 \hat{z} component of the dimensionless spin vector χ → A subscript → 𝜒 𝐴 \vec{\chi}_{A} for each black hole. The initial orbital frequency is Ω 0 subscript Ω 0 \Omega_{0} . Initial orbital parameters were chosen so that the eccentricity was below 5 × 10 − 4 5 superscript 10 4 5\times 10^{-4} . The time simulated to merger is t Merger subscript 𝑡 Merger t_{\mathrm{Merger}} , and the amount of ringdown simulated thereafter is t RD subscript 𝑡 RD t_{\mathrm{RD}} , both in units of G ​ M 𝐺 𝑀 GM . The final mass of the remnant black hole is m Final subscript 𝑚 Final m_{\mathrm{Final}} , in units of M 𝑀 M . The remnant spins are in the z 𝑧 z -direction, and thus we give the z ^ ^ 𝑧 \hat{z} component χ Final subscript 𝜒 Final \chi_{\mathrm{Final}} of the dimensionless spin.

Name	$\frac{m_1}{m_2}$	${\chi }_1$	${\chi }_2$	${\Omega }_0\left(GM\right)$	$\frac{t_{\mathrm{Merger}}}{GM}$	$\frac{t_{\mathrm{RD}}}{GM}$	$\frac{m_{\mathrm{Final}}}{M}$	${\chi }_{\mathrm{Final}}$
Spin 0.3	3.0	0.30	0.30	0.0163	5841	764	0.96	0.68
Spin 0.1	3.0	0.10	0.10	0.0164	5452	817	0.97	0.59
Spin 0.0	3.0	0.00	0.00	0.0190	3457	697	0.97	0.54

Equations

$$I=\int d^{4}x\sqrt{-g}\left[L_{\text{EH}}+L_{\vartheta }+L_{\text{int}}+L_{\text{mat}}+\dots \right],$$

$$L_{\text{EH}}$$

$$L_{\text{int}}$$

$${}^{\ast }RR\equiv {}^{\ast }{R}^{abcd}{R}_{abcd}=\frac{1}{2}{ϵ}^{abef}{R}_{ef}{}^{cd}{R}_{abcd},$$

$${\kappa }^{\text{YSYT}}$$

$${\kappa }^{\text{AY}}$$

$${\kappa }^{\text{MSY}}$$

$$I_{\epsilon }=\int d^{4}x\sqrt{-g}\left[L_{\text{EH}}+L_{\vartheta }+\epsilon L_{\text{int}}+L_{\text{mat}}+\dots \right],$$

$$\square \vartheta =\epsilon \frac{m_{\text{pl}}}{8}{\ell }^2{}^{\ast }RR,$$

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

$$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 }={\nabla }_a\vartheta {\nabla }_b\vartheta -\frac{1}{2}{g}_{ab}{\nabla }_c\vartheta {\nabla }^c\vartheta .$$

$$g_{ab}$$

$$\vartheta$$

$$m_{\text{pl}}^2{G}_{ab}[g^{(0)}]$$

$${\square }^{(0)}{\vartheta }^{(0)}$$

$$G_{ab}[g^{(0)}]=0,$$

$$(g^{(0)},{\vartheta }^{(0)})=(g^{\text{GR}},0),$$

$$m_{\text{pl}}^2{G}_{ab}^{(1)}[h^{(1)};g^{(0)}]$$

$${\square }^{(0)}{\vartheta }^{(1)}+{\square }^{(1)}{\vartheta }^{(0)}$$

$$C_{ab}^{(0)}[{\vartheta }^{(0)}=0]$$

$$m_{\text{pl}}^2{G}_{ab}^{(1)}[h^{(1)};g^{(0)}]$$

$${\square }^{(0)}{\vartheta }^{(1)}$$

$$m_{\text{pl}}^2{G}_{ab}^{(1)}[h^{(2)}]$$

$${\square }^{(0)}{\vartheta }^{(2)}$$

$$\mathcal{O}({\epsilon }^0):$$

$$\mathcal{O}({\epsilon }^1):$$

$$\mathcal{O}({\epsilon }^2):$$

$$(GM{)}^{-2}{\square }^{(0)}{\vartheta }^{(1)}=\frac{m_{\text{pl}}}{8}{\ell }^2(GM{)}^{-4}[{}^{\ast }RR{]}^{(0)}.$$

$${\vartheta }^{(1)}=\frac{m_{\text{pl}}}{8}{\left(\frac{\ell }{GM}\right)}^2\Psi .$$

$${\square }^{(0)}\Psi =[{}^{\ast }RR{]}^{(0)}.$$

$$h_{ab}^{(2)}\equiv {\left(\frac{\ell }{GM}\right)}^4{{\rm Y}}_{ab}.$$

$${\nabla }^2{\rm Y}+L.O.T.\sim (\nabla \Psi {)}^2+(\nabla \Psi )(\nabla R)+({\nabla }^2\Psi )R,$$

$${\vartheta }_A^{(1)}=\frac{{\mu }_A^i{n}_A^i}{R_A^2},$$

Full text

I Introduction

General relativity has been observationally and experimentally tested for almost a century, and has been found consistent with all precision tests to date Will (2014). But no matter how well a theory has been tested, it may be invalidated at any time when pushed to a new regime. Indeed, there are many theoretical reasons to believe that general relativity (GR) cannot be the ultimate description of gravity, from non-renormalizability to the black hole information problem.

Moreover, from the empirical standpoint, all precision tests of GR to date have been in the slow-motion, weak-curvature regime. With the Laser Interferometer Gravitational Wave Observatory (LIGO) now detecting the coalescence of compact binary systems Abbott et al. (2016a, b, 2017), we finally have direct access to the non-linear, dynamical, strong-field regime of gravity. This is an arena where GR lacks precision tests, and it may give clues to a theory beyond GR. The LIGO collaboration has already used the detections of GW150914, GW151226, and GW170104 to perform some tests of GR Abbott et al. (2016c, 2017), but these are not yet very precise: a model-independent test gives 96% agreement with GR.

Both black hole (BH) and neutron star (NS) binaries probe the strong-field regime. However, NSs have the added complication that the equation of state of dense nuclear matter is presently unknown. Until more is known about the equation of state, we must rely on binary black holes (BBHs) for precision tests of GR. Yunes, Yagi, and Pretorius argued Yunes et al. (2016) that the lack of understanding of BBH merger in beyond-GR theories severely limits the ability to constrain gravitational physics using GW150914 and GW151226. Thus, to perform tests of GR with BBHs, we require inspiral, merger, and ringdown waveform predictions for these systems, which can only come from numerical simulations.

To date, BBH simulations have only been performed in GR and scalar-tensor gravity Berti et al. (2015) (note that BBHs in massless scalar-tensor gravity will be identical to GR, under ordinary initial and boundary conditions). There are a huge number of beyond-GR theories Berti et al. (2015), and for the vast majority of them, there is no knowledge of whether there is a well-posed initial value formulation, a necessity for numerical simulations. Indeed, there is evidence that dynamical Chern-Simons gravity, the beyond-GR theory we use here as an example, lacks a well-posed initial value formulation Delsate et al. (2015).

Our goal is to numerically integrate BBH inspiral, merger, and ringdown in theories beyond GR that are viable but that do not necessarily have a well-posed initial value problem. This goal is relevant even for those only interested in parametric, model-independent tests, because there is presently no theory guidance for late-inspiral and merger waveforms in theories beyond GR.

We are only interested in theories that are sufficiently “close” to GR: for a theory to be viable, it has to be able to pass all the tests that GR has passed. This motivates an effective field-theory (EFT) approach. We assume that there is a high-energy theory whose low-energy limit gives GR plus “small” corrections. The effective theory of GR with corrections does not need to capture arbitrarily short-distance physics. Such a theory is valid up to some cutoff, and modes shorter than this distance scale are said to be outside of the regime of validity of the EFT. The EFT only needs to be well-posed for the modes within the regime of validity. This can be accomplished with perturbation theory.

We present a perturbation scheme for numerically integrating beyond-GR theories that limit to GR. For such a theory, we perturb it about GR in powers of the small coupling parameter. We collect equations of motion at each order in the coupling, creating a tower of equations, with each level inheriting the same principal part as the background GR system. The well-posedness of the initial value problem in GR Wald (1984) thus ensures the well-posedness of this framework, even if the “full” underlying theory may not have a well-posed initial value formulation.

In this study, we apply our perturbation framework to BBH mergers in dynamical Chern-Simons gravity (dCS) Alexander and Yunes (2009), to linear order in perturbation theory. This theory involves a pseudoscalar field coupled to the parity-odd Pontryagin curvature invariant with a small coupling parameter, and at linear order gives a scalar field evolving on a GR BBH background.

There are a number of theoretical motivations for considering dynamical Chern-Simons. The dCS interaction arises when cancelling gravitational anomalies in chiral theories in curved spacetime Delbourgo and Salam (1972); Eguchi and Freund (1976); Alvarez-Gaume and Witten (1984), including the famous Green-Schwarz anomaly cancellation in string theory Green and Schwarz (1984) when compactified to four dimensions Alexander and Yunes (2009); Polchinski (1998a, b). DCS also arises in loop quantum gravity when the Barbero-Immirzi parameter is allowed to be a spacetime field Taveras and Yunes (2008); Mercuri and Taveras (2009). From an EFT standpoint, dCS is the lowest-mass-dimension correction that has a parity-odd interaction. All other EFTs at the same mass dimension have parity-even interactions, so the phenomenology of dCS is distinct Yagi et al. (2012a). The dCS interaction was also included in Weinberg’s EFT of inflation Weinberg (2008).

From a practical standpoint, there are already a large number of dCS results in the literature that we can compare against Yunes and Pretorius (2009); Yagi et al. (2012a, b, c, 2013); Konno and Takahashi (2014); Stein (2014a), including post-Newtonian (PN) calculations for the BBH inspiral. One of the more important results is that scalar dipole radiation is highly suppressed in dCS during the inspiral Yagi et al. (2012a). Dipole radiation is present in scalar-tensor theory and Einstein-dilaton-Gauss-Bonnet (EdGB), and enters with two fewer powers of the orbital velocity (i.e. 1 PN order earlier) than the leading quadrupole radiation of GR. This leads to gross modifications of the inspiral, but dCS avoids this problem because the dipole is suppressed. As a result, the perturbative treatment of dCS will be valid for a longer period of inspiral than scalar-tensor or EdGB.

The paper is organized as follows. Section II covers the analytical and numerical formalisms. More specifically, in Sec. II.1 we introduce dynamical Chern-Simons, and in Sec. II.2 we present the perturbation scheme, which is valid for any theory with a continuous limit to GR. We discuss the numerical scheme in Sec. II.3 (some numerical details are in the Appendix). We present the results of numerically implementing this formalism in dCS on three different binary mergers in Sec. III. Section III.1 reviews some previously-known analytic phenomenology of the BBH inspiral problem in dCS. Section III.2 presents the waveform results, and III.3 presents the energy fluxes, both including comparison to PN. In Sec. III.4 we use the numerical results to assess the validity of the perturbation scheme. In Sec. III.5 we use the numerical results to estimate the detectability of dCS and the bounds that could be placed by LIGO detections. We conclude and discuss in Sec. IV, and lay out plans for future work.

II Formalism

Throughout this paper, we set $c=1$ and $\hbar=1$ so that $[M]=[L]^{-1}$. Since there will be more than one length scale, we explicitly include factors of the reduced Planck mass $m_{\textrm{\tiny{pl}}}^{-2}=8\pi G$ and the “bare” gravitational length $GM$, though quantities in our code are non-dimensionalized with $GM=1$. Latin letters in the middle of alphabet $\{i,j,k,l,m,n\}$ are (3-dimensional) spatial indices, while Latin letters in the beginning of the alphabet $\{a,b,c,d\}$ refer to (4-dimensional) spacetime indices. We follow the sign conventions of Wald (1984), and $g_{ab}$ refers to the 4-dimensional spacetime metric, with signature $(-{}+{}+{}+)$, and with $\nabla$ its Levi-Civita connection.

II.1 Action and equations of motion

The method we present in this paper applies to a large number of beyond-GR theories that have a continuous limit to GR, but for concreteness we focus on dCS. We start with the four-dimensional action

[TABLE]

where the omitted terms $(\ldots)$ are above the cutoff of our EFT treatment. Here $g$ without indices is the determinant of the metric, $L_{{\textrm{\tiny{EH}}}}$ is the Einstein-Hilbert Lagrangian, $L_{\vartheta}$ is the Lagrangian of a minimally coupled (pseudo-)scalar field $\vartheta$ (also referred to in the literature as the axion), $L_{{\textrm{\tiny{int}}}}$ is a beyond-GR interaction between $\vartheta$ and curvature terms, and $L_{{\textrm{\tiny{mat}}}}$ is the Lagrangian for ordinary matter. In this paper, we are considering a binary black hole (BBH) merger in dCS, so we ignore $L_{{\textrm{\tiny{mat}}}}$.

Explicitly, these action terms are given by

[TABLE]

Here the Ricci scalar of $g_{ab}$ is $R$. With our unit system, $[g]=[L]^{0}$, coordinates carry dimensions of length, $[x]=[L]^{1}$, and note that the scalar field $\vartheta$ has been canonically normalized, $[\vartheta]=[L]^{-1}$. We have omitted any potential $V(\vartheta)$, so $\vartheta$ is massless and long-ranged, as appropriate for a “gravitational” degree of freedom. In the interaction Lagrangian $L_{{\textrm{\tiny{int}}}}$, the scalar field $\vartheta$ is coupled to the 4-dimensional Pontryagin density (also known as the Chern-Pontryagin density) ${}^{*}\!RR$,

[TABLE]

where $\epsilon_{abcd}$ is the fully antisymmetric Levi-Civita tensor.

The coupling strength of this interaction is governed by the new parameter $\ell$ with dimensions of length. This parameter takes on specific values if this EFT arises from the low-energy limit of certain string theories Green and Schwarz (1984) or to cancel gravitational anomalies Alvarez-Gaume and Witten (1984); Polchinski (1998a, b). However, here we simply take it as a “small” coupling parameter. In the limit that $\ell\to 0$, we recover general relativity with a massless, minimally coupled scalar field.

The coupling parameter conventions vary throughout the literature. To enable comparisons, we express the couplings of a number of works in terms of our conventions. To put Yagi et al. Yagi et al. (2012a) into our conventions, use

[TABLE]

To convert Alexander and Yunes Alexander and Yunes (2009) into our conventions,

[TABLE]

To compare with McNees et al. McNees et al. (2016), use

[TABLE]

The conventions of Stein Stein (2014a) agree with ours (except for an inconsequential sign change in the definition of ${}^{*}\!RR$, which is compensated for by an additional sign everywhere ${}^{*}\!RR$ appears).

Below we will perform an expansion in powers of $\ell^{2}$. To simplify matters, we insert a dimensionless formal order-counting parameter $\varepsilon$ that will keep track of powers of $\ell^{2}$. Expanding in a dimensionless parameter ensures that field quantities at different orders have the same length dimension.

Specifically, we replace the action in Eq. (1) with

[TABLE]

a one-parameter family of actions parameterized by $\varepsilon$. Formally, we recover the action in Eq. (1) when $\varepsilon=1$.

Varying the action Eq. (7) with respect to the scalar field, we have the sourced wave equation

[TABLE]

where $\square=\nabla_{a}\nabla^{a}$ is the d’Alembertian operator. Varying with respect to the metric gives the corrected Einstein field equations,

[TABLE]

where $G_{ab}$ is the Einstein tensor of $g_{ab}$, and the tensor $C_{ab}$ includes first and second derivatives of $\vartheta$, and second and third derivatives of the metric,

[TABLE]

Since we are focusing on BBH mergers, $T^{{\textrm{\tiny{mat}}}}_{ab}=0$. The scalar field’s stress-energy tensor $T^{\vartheta}_{ab}$ is given by the expression for a canonical, massless Klein-Gordon field,

[TABLE]

From here forward we will drop the superscript $\vartheta$.

The “full” system of equations for dCS is thus the pair of Eqs. (8) and (9).

II.2 Perturbation scheme

Because $C_{ab}$ in Eq. (9) contains third derivatives of the metric, the “full” system of equations for dCS likely lacks a well-posed initial value formulation Delsate et al. (2015). In the language of particle physics, this is equivalent to the appearance of ghost modes above a certain momentum scale Dyda et al. (2012).

From the EFT point of view, though, the ghost modes and ill-posedness are nothing more than the breakdown of the regime of validity of the theory, which should be valid for long wavelength modes in the decoupling limit $\ell\to 0$. To excise the ghost modes and arrive at a well-posed initial value formulation, we expand about $\varepsilon=0$, which is simply GR coupled to a massless minimally-coupled scalar field and certainly has a well-posed initial value problem Wald (1984). As a result, all higher orders in $\varepsilon$ will inherit the well-posedness of the zeroth-order theory, by inheriting the principal parts of the differential equations.

We begin this order-reduction scheme by expanding the metric and scalar field in power series in $\varepsilon$,111Note that this is not a Taylor series, since there is no factor of $1/k!$ in the $k$th term. These factors must be tracked if using standard perturbation theory, e.g. with the xPert package xAc ; Brizuela et al. (2009).

[TABLE]

Note that since $\varepsilon$ is dimensionless, each $\vartheta^{(k)}$ has the same units as $\vartheta$, and similarly for $h_{ab}^{(k)}$. This expansion is now inserted into the field equations, which are likewise expanded in powers of $\varepsilon$, and we collect orders homogeneous in $\varepsilon^{k}$, as below. This results in a “tower” of systems of equations that must be solved at progressively increasing orders in $\varepsilon$. This scheme is quite general, and should apply to any theory that has a continuous limit to GR.

II.2.1 Order $\varepsilon^{0}$

Zeroth order comes from taking $\varepsilon\to 0$, which simply gives the system of GR coupled to a massless, minimally coupled scalar field,

[TABLE]

where $G_{ab}[g^{(0)}]$ is the Einstein tensor of the background metric $g^{(0)}$, $\square^{(0)}$ is the associated d’Alembert operator, and $T^{(0)}$ is the stress-energy of $\vartheta^{(0)}$. This system certainly has a well-posed initial value problem.

Because of the explicit presence of $\varepsilon$ in front of $L_{{\textrm{\tiny{int}}}}$ in the action [Eq. (7)], $C_{ab}$ does not appear in the metric equation (13a), and the Pontryagin source does not appear on the right-hand side of the scalar equation (13b). These terms have been pushed to one order higher and will appear below.

On general grounds, we expect that any initially non-vanishing scalar field will radiate away within a few dynamical times. Similarly, if we start with a $\vartheta^{(0)}=0$ initial condition and impose purely outgoing boundary conditions, $\vartheta^{(0)}$ will remain zero throughout the entire simulation. Therefore, rather than simulating a vanishingly small $\vartheta^{(0)}$, we simply analytically assume that $\vartheta^{(0)}=0$.

Therefore, at order $\mathcal{O}(\varepsilon^{0})$, the system will simply be

[TABLE]

and the solution will be

[TABLE]

where $g^{{\textrm{\tiny{GR}}}}$ is a GR solution to the BBH inspiral-merger-ringdown problem.

II.2.2 Order $\varepsilon^{1}$

Continuing to linear order in $\varepsilon$, we find the system

[TABLE]

As noted above, the explicit presence of $\varepsilon$ in the action (7) and equations of motion [(8) and (9)] lead to $C^{(0)}$ and $[{}^{*}\!RR]^{(0)}$ appearing in these $\varepsilon^{1}$ equations strictly as source terms. By construction, the principal part of this differential system is the same as the principal part of the $\mathcal{O}(\varepsilon^{0})$ system, and thus it inherits its well-posedness property. This is true at all higher orders in perturbation theory.

Here, $G^{(1)}[h^{(1)};g^{(0)}]$ is the linearized Einstein operator, built with the covariant derivative $\nabla^{(0)}$ compatible with $g^{(0)}$, acting on the metric deformation $h^{(1)}$. The d’Alembert operator receives the correction $\square^{(1)}$, which depends on the metric deformation $h^{(1)}$. The quantity $C^{(0)}_{ab}$ is the same as the definition given in Eq. (10), evaluated on the background quantities $(g^{(0)},\vartheta^{(0)})$. Similarly, $[{}^{*}\!RR]^{(0)}$ is the Pontryagin density evaluated on the background spacetime metric $g^{(0)}$. Finally, $T^{(1)}_{ab}$ is the first-order perturbation to the stress-energy tensor; since $T_{ab}$ is quadratic in $\vartheta$, $T^{(1)}_{ab}$ has pieces both linear and quadratic in $\vartheta^{(0)}$ (the quadratic-in-$\vartheta^{(0)}$ pieces are linear in $h^{(1)}$).

The crucial property at this order is that both $C^{(0)}$ and $T^{(1)}$ are built from pieces linear and quadratic in $\vartheta^{(0)}$. At order $\mathcal{O}(\varepsilon^{0})$, we found that $\vartheta^{(0)}=0$. Therefore, when evaluated on the $\mathcal{O}(\varepsilon^{0})$ solution [Eq. (15)], these both vanish,

[TABLE]

Therefore, at order $\mathcal{O}(\varepsilon^{1})$ in perturbation theory, evaluating on the background solution, we have the system

[TABLE]

In the metric perturbation equation (18a), starting with $h^{(1)}=0$ initial conditions and imposing purely outgoing boundary conditions will enforce $h^{(1)}=0$ throughout the entire simulation. Similarly, we can argue that small perturbations of $h^{(1)}$ would radiate away on a few dynamical times, since there is no potential to confine the metric perturbations. Once again, rather than simulating a vanishingly small field, we will just analytically assume that $h^{(1)}=0$. Therefore, at order $\mathcal{O}(\varepsilon^{1})$, there is no metric deformation, and the system is only Eq. (18b), driven by the background system (14) which generates the source term $[{}^{*}\!RR]^{(0)}$.

II.2.3 Order $\varepsilon^{2}$

This perturbation scheme can be extended to any order desired. Although this paper reports only on work extending through $\mathcal{O}(\varepsilon^{1})$, we sketch the derivation of $\mathcal{O}(\varepsilon^{2})$, since that is the lowest order where a metric deformation is sourced.

Schematically, the system at $\mathcal{O}(\varepsilon^{2})$, after accounting for the vanishing of $\vartheta^{(0)}$ and $h^{(1)}$, is

[TABLE]

The operator $C^{(1)}[\vartheta^{(1)}]$ is linear in its argument, and $T^{(2)}[\vartheta^{(1)},\vartheta^{(1)}]$ is linear in each slot. Various other combinations have vanished. In (19a), vanishing source terms were quadratic in $h^{(1)}$ or built from the product of $h^{(1)}\times\vartheta^{(1)}$. In (19b), $\ell^{2}[{}^{*}\!RR]^{(1)}$ is proportional to $h^{(1)}$ and thus vanishes, as do terms such as $\square^{(1)}\vartheta^{(1)}$ (linear in $h^{(1)}$) and $\square^{(2)}\vartheta^{(0)}$ (linear in $\vartheta^{(0)}$).

We leave detailed discussion of order $\mathcal{O}(\varepsilon^{2})$ to future work Okounkova and Stein (2017).

II.2.4 Summary and scaling

Let us briefly summarize the perturbative order-reduction scheme and discuss the scaling of different orders. The system at orders $\varepsilon^{0}$ and $\varepsilon^{1}$ is

[TABLE]

where $T^{{\textrm{\tiny{eff}}}}_{ab}$ may be determined from the right hand side of Eq. (19a).

Zeroth order (20a) is just vacuum GR, which has no intrinsic scale. As is very common in numerical relativity simulations, the coordinates used in the simulation are dimensionless and in units of the total ADM mass, $X^{a}=x^{a}/(GM)$. This means that $\nabla$ may be non-dimensionalized by pulling out a factor of $(GM)^{-1}$, Riemann may be non-dimensionalized by pulling out a factor of $(GM)^{-2}$, etc.

Meanwhile, the new length scale and coupling parameter $\ell$ enters at first order. If we non-dimensionalize the derivative operator and curvature tensors in Eq. (20b), we will find

[TABLE]

We therefore define the dimensionless scalar field $\Psi$ via

[TABLE]

Then $\Psi$ will satisfy

[TABLE]

Thus the analytic dependence of $\vartheta^{(1)}$ on $(\ell/GM)$ has been extracted. The solution $\Psi$ can later be scaled to reconstruct $\vartheta^{(1)}$ for any allowable value of $(\ell/GM)$.

All of the results that we present will be given in terms of powers of the dimensionless coupling $(\ell/GM)$. We will also compare to known post-Newtonian results Yagi et al. (2013), that were presented in terms of $\alpha_{4}^{{\textrm{\tiny{YSYT}}}}$. To perform the comparison, we use the conversion given in Eq. (4).

Finally, though we do not address $\mathcal{O}(\varepsilon^{2})$ simulations in this paper, we should still study how $h^{(2)}$ scales with $\ell$ and $(GM)$. Since the perturbative scheme preserves the units of length of fields, $[h^{(k)}]=[g]=[L]^{0}$ is already dimensionless; however, it still depends on $(\ell/GM)$ in a specific way. When we move to units in which we measure lengths and times in units of $(GM)$, we find it is appropriate to define a scaled metric deformation $\Upsilon$ via

[TABLE]

Then this dimensionless quantity $\Upsilon$ will satisfy an equation that is schematically

[TABLE]

where L.O.T. stands for lower order terms, and all derivatives and curvatures are $\mathcal{O}(\varepsilon^{0})$ dimensionless quantities.

II.3 Numerical scheme

For the order $\varepsilon^{1}$ part of the order reduction scheme, our overall goal is to solve Eq. (23) on a dynamical background metric. We co-evolve the metric and the scalar field, where Eq. (23) is driven by Eq. (20a). The whole system is simulated using the Spectral Einstein Code (SpEC) SpE , which uses the generalized harmonic formulation of general relativity in a first-order, constraint-damping system Lindblom et al. (2006) in order to ensure well-posedness and hence numerical stability. We have added a scalar field module that is similarly a first-order, constraint-damping system, following Holst et al. (2004), as outlined in App. A.

The code uses pseudospectral methods on an adaptively-refined grid Lovelace et al. (2011); Szilágyi (2014), and thus numerical convergence with resolution of both the metric variables and the scalar field is exponential. We demonstrate the numerical convergence of the scalar field in App. A.

The initial data for the binary black hole background is a superposition of two Kerr-Schild black holes with a Gaussian roll-off of the conformal factor around each black hole Lovelace et al. (2008). The initial data for the scalar field is similarly given by a superposition of approximate dCS solutions around isolated black holes, and is given in more detail in Sec. III.2.

The metric equations are evolved in a damped harmonic gauge Szilágyi et al. (2009); Lindblom and Szilágyi (2009), with excision boundaries just inside the apparent horizons Hemberger et al. (2013); Scheel et al. (2015), and minimally-reflective, constraint-preserving boundary conditions on the outer boundary Rinne et al. (2007). The scalar field system, meanwhile, uses purely outgoing boundary conditions modified to reduce the influx of constraint violations into the computational domain Holst et al. (2004).

The Pontryagin density source term ${}^{*}\!RR$ is computed throughout the simulation in a 3+1 split from the available spatial quantities as outlined in App. B.

III Results

III.1 Background: Phenomenology of binary black hole inspirals in dCS

To give the proper context for our numerical results, we first review the previously-known phenomenology relevant to this problem. Analytical and numerical results are known for isolated black holes in the decoupling limit, and analytical results are known for the binary black hole problem in the decoupling limit and at slow velocities ($v/c\ll 1$).

Any spherically-symmetric metric will have vanishing Pontryagin density.222This is straightforward to verify with a computer algebra system, using the canonical form for a spherically symmetric metric, $ds^{2}=-e^{2\alpha(t,r)}dt^{2}+e^{2\beta(t,r)}dr^{2}+r^{2}d\Omega^{2}$. Since it is true in this coordinate system, it is true in general. This is also proven in App. A of Grumiller et al. (2002) following a tensorial approach. Finally, one can appeal to a symmetry argument. If the metric is invariant under an $O(3)$ isometry, then the curvature tensor and ${}^{*}\!RR$, being tensorial objects built only from $g$, must also be invariant under this symmetry. Therefore ${}^{*}\!RR$ must be a constant on each 2-sphere. The group $O(3)$ also contains the reflection symmetry, sending points to their antipodes. The metric is invariant under this reflection, but ${}^{*}\!RR$ must flip sign, as it is a pseudo-scalar. But then we must have ${}^{*}\!RR=-{}^{*}\!RR$, so ${}^{*}\!RR=0$.

Thus the Schwarzschild solution with vanishing scalar field is already a solution to the “full” dCS system. An isolated spinning black hole in dCS, however, is not given by the Kerr solution of GR Campbell et al. (1991); Yunes and Pretorius (2009); Konno et al. (2009); Yagi et al. (2012b); the scalar field is sourced, and the metric acquires corrections. Analytical results for the leading-order, small-coupling corrections to the Kerr metric have been found in the slow-rotation approximation ($a\ll M$) Yunes and Pretorius (2009); Konno et al. (2009); Yagi et al. (2012b); Maselli et al. (2017). Additionally, numerical results have been found for the scalar field for general rotation Konno and Takahashi (2014); Stein (2014b). The leading-order correction to Kerr is dipolar scalar hair, while the scalar monopole vanishes. This vanishing scalar monopole means that scalar dipole radiation is heavily suppressed in dCS. At a large radius away from an isolated black hole labeled by $A$, the dipolar scalar field goes as

[TABLE]

where $R_{A}$ is the distance from black hole $A$, $n_{A}^{i}$ is the spatial unit vector pointing away from BH $A$, and $\mu_{A}^{i}$ is the scalar dipole moment of the BH. This scalar dipole moment is given by Yagi et al. (2012a)

[TABLE]

where $\chi_{A}^{i}$ is the dimensionless spin vector of black hole $A$, $\chi_{A}^{i}=J_{A}^{i}/GM_{A}^{2}$ (this factor of $G$ in the denominator arises from our usage of natural units, where angular momentum is dimensionless, $[J]=[L]^{0}$, in units of $\hbar$).

The dCS binary inspiral problem in the post-Newtonian regime ($v\ll c$) was first treated by Yagi et al. Yagi et al. (2012a). When two spinning BHs with scalar dipole hair are placed in proximity with each other, the hair is responsible for a number of effects. First, there is a correction to the binding energy due to the dipole-dipole interaction. Second, as the BHs orbit each other, the net quadrupole of the binary system has a time derivative on the orbital timescale. The binary’s combined dipole moment is also time-varying, but only on the spin-precession timescale, so it is heavily suppressed. Thus in the far zone of the binary, the scalar field exhibits predominantly quadrupole and higher radiation, and no $l=0$ monopole radiation.

The dominant far-zone multipole moments for the scalar field have $|m|=l-1$ with $l\geq 2$ and the $l=1$ modes radiate on the spin-precession timescale. To make comparing to PN simpler, we are simulating aligned-spin systems, so the $l=1$ mode will in fact be non-radiative at early times. Yagi et al. Yagi et al. (2012a) gave expressions for the scalar field $\vartheta^{(1)}$ due to spinning and non-spinning binaries, presented in terms of symmetric tracefree (STF) tensors. In most numerical relativity work, however, we decompose fields into spherical harmonics,

[TABLE]

Using Blanchet and Damour (1986), we convert the STF expressions from Yagi et al. (2012a); Stein and Yagi (2014) into spherical harmonics at extraction radius $R$ for a spin-aligned binary, when the post-Newtonian approximation is valid (the early inspiral), giving

[TABLE]

Here $\phi=\phi(t)$ is the orbital phase, $\omega=\omega(t)=\dot{\phi}$ is the orbital frequency, $m_{A}$ is the mass of each black hole, $M=m_{1}+m_{2}$ is the total mass,333In PN literature, $m$ is often used as the total mass. We use $M$ here in order to be consistent with numerical relativity literature. and $\mu_{A}$ is the $z$ component (the only component since this calculation is for a spin-aligned binary) of the scalar dipole moment from Eq. (27). Note that the $(1,0)$ mode is time-independent (and hence non-radiative), since we are focusing on spin-aligned systems.

The behavior of the scalar field during the late inspiral and merger was previously unknown, and is part of the motivation for the present numerical study.

III.2 Scalar field waveforms

We performed three numerical simulations in this formalism, each at low, medium, and high numerical resolutions, with parameters given by Table 1. We chose three values for the BHs’ dimensionless spins of $0.0$, $0.1$, and $0.3$, to qualitatively see the effect of spin on the physics, and to allow for comparison with analytical calculations. While SpEC can simulate very high spins Scheel et al. (2015), the analytics we compare against use the small-spin expansion and stop at linear order in spin. Therefore the $\mathcal{O}(\chi^{2})$ errors should be at most $\sim 30$% of the $\mathcal{O}(\chi)$ effects we compare against. Similarly, while modeling spin precession is possible Ossokine et al. (2015), it is not the focus of this study, and thus we have eliminated this complication by aligning all of the spins with the orbital angular momentum.

As mentioned in Sec. III.1, the scalar field around an isolated, slowly spinning black hole in dCS is approximately a dipole. We use this analytic approximation as the basis for our initial data, as mentioned in Sec. II.3. The initial scalar field is a superposition of two slow-rotation dipole solutions (since all of the dimensionless spins are $\leq 0.3$), one around each black hole. We apply a boost to account for the initial velocity of each black hole. As our scalar field evolution system is first-order (see App. A), we also initialize the variables corresponding to the spatial and time derivatives of $\Psi$ to the analytical derivatives of the approximate dipole solution. For the non-spinning simulation, we set the initial value of $\Psi$ and its derivatives to zero.

We plot mode-decomposed waveforms extracted from the highest resolution simulations of the three simulations in Figs. 1, 2, and 3. Each figure shows the $(l=2,m=2)$ mode of the Newman-Penrose quantity $\Psi_{4}$ decomposed into spin-weight $-2$ spherical harmonics, and the dominant $(l,m=l-1)$ modes of the scalar $\vartheta^{(1)}$ for $l=1,2,3$, along with the PN comparisons from Eq. (29).

We immediately see that at early times, there is good qualitative agreement between the numerical waveforms and the PN predictions, with the $(l=2,m=1)$ mode dominating, as expected. In the PN formulas of Eq. (29), we used the instantaneous coordinate orbital frequency and phase calculated from the black hole trajectories for $\omega$ and $\phi$. Since the starting phase is arbitrary, we perform a phase alignment (by eye) between the numerical results and the PN waveforms.

As expected, because the spins are not precessing, there is no dipole radiation at early times. The offset away from zero seen in the $(l=1,m=0)$ panel of Fig. 1 is a real physical effect: it is due to the combined dipole moments of the two individual black holes and their orbital angular momentum. After merger, the $l=1$ moment settles down to a new non-zero value (below the resolution of this figure) determined by the spin of the final black hole, again via Eq. (27). In between, there is a burst of scalar dipole radiation. This is a newly discovered phenomenon that could not have been computed with analytic post-Newtonian calculations. Scalar monopole radiation, meanwhile, is consistent with zero within the numerical errors of the simulation.

III.3 Energy fluxes

Having solved for the scalar field $\vartheta^{(1)}$, we can evaluate physical quantities such as its stress-energy tensor, Eq. (11). From $T_{ab}^{(\vartheta)}$, we can compute the energy flux through some 2-sphere $S^{2}_{R}$ at coordinate radius $R$ via

[TABLE]

Here $n^{a}$ is the timelike unit normal to the spatial slice, and $dS^{b}$ is the proper area element of $S^{2}_{R}$, i.e. $dS^{b}=N^{b}\sqrt{\gamma}dA$, where $N^{b}$ is the spacelike unit normal to $S^{2}_{R}$, $\gamma$ is the determinant of the induced 2-metric, and $dA$ is the coordinate area element.

Like the metric and scalar field, we similarly expand $T_{ab}^{(\vartheta)}$ and $\dot{E}^{(\vartheta)}$ in powers of $\varepsilon$,

[TABLE]

where each $\dot{E}^{(\vartheta,k)}$ includes the appropriate orders of both the scalar field and metric. Since $\vartheta^{(0)}=0$ and $T_{ab}^{(\vartheta)}$ is quadratic in $\vartheta$, we have $T_{ab}^{(\vartheta,0)}=T_{ab}^{(\vartheta,1)}=0$, and similarly $\dot{E}^{(\vartheta,0)}=\dot{E}^{(\vartheta,1)}=0$. The first non-vanishing order is $T_{ab}^{(\vartheta,2)}$, which is given by

[TABLE]

Using the results of the simulations, we compute $T_{ab}n^{a}$, interpolate it onto surfaces of fixed coordinate radius $R$, compute $T_{ai}n^{a}N^{i}$ by contracting with the normal, and perform spectral integration with the induced area element to obtain $\dot{E}^{(\vartheta,2)}$. That is, we compute

[TABLE]

We also compute the energy flux at order $(\ell/GM)^{0}$, which for vanishing $\vartheta^{(0)}$ consists purely of the background gravitational energy flux, as (c.f. Ruiz et al. (2008))

[TABLE]

where numerically we set the lower bound of the time integral to the start of the simulation, assuming there was comparatively little radiation before the start.

We plot the numerical values of $\dot{E}^{(\vartheta,2)}(R)$ and $\dot{E}^{(0)}(R)$ in Fig. 4, keeping (spin-weighted) spherical harmonics up through $l=8$. We check for the convergence of the flux quantities with increasing extraction radius, and present the results at $R=300~{}GM$, which agree with the results at $R=200~{}GM$.

In Fig. 4 we also plot a post-Newtonian approximation to $\dot{E}^{(\vartheta,2)}$. This is computed using the far-zone PN solution for $\vartheta^{(1)}$ from Yagi et al. (2012a), which only includes the $l=2$ quadrupole radiation. We impose circular orbits and aligned spins, convert to our conventions via Eq. (4), and re-insert the appropriate factors of $G$. The result for at least one non-zero spin is

[TABLE]

and for two non-spinning black holes,

[TABLE]

In these expressions, $\chi_{A}$ is the dimensionless spin of black hole $A$, $\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio, and $\delta m=m_{1}-m_{2}$ is the mass difference.

Although the gravitational flux at order $(\ell/GM)^{0}$ is by far the largest energy flux, the scalar field flux at order $(\ell/GM)^{4}$ sharply increases before merger. The spin contributions are dominant, as the scalar flux for the spin-0 simulation is comparatively small until the merger, when nonlinearities become very important. At early times, our fully numerical results qualitatively agree with the PN results of Yagi et al. (2012a), validating our and their calculations. We expect the $\mathcal{O}(1)$ ratio between PN and full numerics in Fig. 4 stems from the PN expressions (35), (36) only including $l=2$, whereas our numerics include all modes up through $l=8$.

III.4 Regime of validity

Since this method is perturbative, we expect that it breaks down—becomes invalid—at some point. There are two types of breakdown. First, at every instant of time, there is the question of whether the series converges. We expect that the series should only converge when $\ell\ll GM$, and we assess this in Sec. III.4.1. Second, over much longer times, there will be a secular drift between the perturbative solution and the “true” solution, so that the two solutions become out of phase. We assess the dephasing below in Sec. III.4.2.

III.4.1 Instantaneous validity

The perturbative scheme is valid pointwise at every instant in time if the series for the metric (12a) and scalar (12b) are convergent. Roughly, we can assess this by comparing the magnitudes of successive terms in the series. As shown in Sec. II.2, up through order $\varepsilon^{2}$, the metric and scalar are expanded as

[TABLE]

Thus we cannot assess the convergence of $\vartheta$ without going to $\mathcal{O}(\varepsilon^{3})$, but at $\mathcal{O}(\varepsilon^{2})$ we can compare the magnitudes of $g_{ab}^{(0)}$ and $h_{ab}^{(2)}$. A rough condition for convergence is that

[TABLE]

where $\left\|\cdot\right\|$ is an $L^{2}$ norm.

The magnitude of $h_{ab}^{(2)}$ depends on the strength of the coupling parameter $\ell$, as discussed in Sec. II.2.4, via $h_{ab}^{(2)}=(\ell/GM)^{4}\Upsilon_{ab}$, where $\Upsilon_{ab}$ is independent of $\ell$. Thus we translate Eq. (38) into a condition on the maximum allowed value of $\ell/GM$,

[TABLE]

where $C$ is some factor of order unity, and on the right-hand side, the ratio is evaluated pointwise, and then the minimum is taken over the domain outside of the apparent horizons, at each coordinate time. At values of $\ell/GM$ larger than this estimate, we expect the perturbative approach fails to converge somewhere in the spacetime.

In these order $\varepsilon^{1}$ simulations, we have not simulated $\Upsilon_{ab}$. We can, however, make scaling estimates from its schematic equation of motion, Eq. (25). The source term $m_{\textrm{\tiny{pl}}}\ell^{2}C_{ab}^{(1)}$ should be of the same order of magnitude as $T_{ab}^{(2)}$ (which we do compute in the simulations), so to within an order of magnitude, we estimate

[TABLE]

Here $L$ is a characteristic curvature length scale, and $T_{ab}[\Psi]$ is shorthand for the “stress-energy” $T_{ab}[\Psi]=\nabla_{a}\Psi\nabla_{b}\Psi-\frac{1}{2}g_{ab}(\nabla\Psi)^{2}$. Therefore, we estimate the allowed value for $\ell/GM$ as

[TABLE]

We plot this estimate in Fig. 5 for each of the spin configurations considered in this study. During inspiral, the curvature is highest around the smaller black hole, so we let $L=\min(Gm_{1},Gm_{2})$. After merger, we let $L=Gm_{\mathrm{Final}}$ (see Table 1 for values).

We can compare our estimates for the regime of validity $|\ell/GM|_{\max}$ to those computed in Stein Stein (2014a). Stein computed $|\ell/Gm|_{\max}$ of a stationary, isolated black hole as a function of $\chi$ of the body, using methods that are independent of ours. At late times, we find direct agreement, at the 5% level, by setting $C=(32)^{1/4}\approx 2.38$. At early times, after including a factor of $M/m_{2}$ to convert from $|\ell/GM|$ to $|\ell/Gm_{2}|$, we again find agreement. At early times, the low-spin simulation has a very large regime of validity, because the Pontryagin density is small, and hence Chern-Simons effects are also small. However, approaching the time of merger, orbital motion and nonlinearities source enough energy density in the scalar field to restrict the regime of validity of $|\ell/GM|$ to order unity.

III.4.2 Secular validity (dephasing)

The true physical system at $\varepsilon>0$ radiates energy more quickly than the GR-only ($\varepsilon=0$) solution that we are using as the background for perturbation theory. As a result, the true solution will inspiral more quickly, so the orbital phase will have a secularly growing deviation away from the background. A post-Newtonian scaling estimate (see below) says that the standard solution will break down over a secular timescale of order $T_{\mathrm{sec}}\sim T^{\textrm{GR}}_{\mathrm{RR}}(\ell/GM)^{-2}v^{-2}$, where $T^{\textrm{GR}}_{\mathrm{RR}}$ is the radiation-reaction timescale in GR. This scaling $(\ell/GM)^{-2}$ is characteristic of singular perturbation theory Bender and Orszag (1978); Chen et al. (1996); Galley and Rothstein (2017).

If the length of a detected gravitational waveform is long compared to the secular breakdown time, then we will need a method to extend the secular regime of validity of the calculation—for example, multiple-scale analysis (MSA) Bender and Orszag (1978) or the dynamical renormalization group Chen et al. (1996); Galley and Rothstein (2017). We save this issue for future work. Here, we will estimate the dephasing time (secular breakdown time).

Let us focus on quasi-circular, adiabatic inspirals. Similarly to the scalar field and metric variables in Eqs. (12a) and (12b), we can expand the accumulated orbital phase $\phi(t)$ and the orbital frequency $\omega(t)=\dot{\phi}(t)$ of the binary in powers of $\varepsilon$,

[TABLE]

where $\phi^{(0)}$ corresponds to the phase of the binary in pure GR, and $\phi^{(1)}$ contains the dCS corrections at order $\varepsilon^{1}$ and so on. Since the metric deformation at $\mathcal{O}(\varepsilon^{1})$ vanishes, the phase correction at $\mathcal{O}(\varepsilon^{1})$ also vanishes, $\phi^{(1)}=0=\omega^{(1)}$. The first non-vanishing orbital phase correction is

[TABLE]

We can use $\Delta\phi$ to assess the secular regime of validity, and in Sec. III.5 we will also use it to assess the detectability of dynamical Chern-Simons.

We do not have $\Delta\phi$ directly from the simulation, as we do not evolve the $\varepsilon^{2}$ system. However, we can estimate it from previously-known analytical results combined with numerical quantities available during the simulation.

Consider the local-in-time expansion of the orbital phase correction $\Delta\phi$ around any ‘alignment time’ $t_{0}$,

[TABLE]

If our simulation had started at reference time $t_{0}$, then we would have $\Delta\phi(t_{0})=0$. The linear piece $(t-t_{0})\omega^{(2)}(t_{0})$ corresponds to a perturbative, instantaneous frequency shift, which is completely degenerate with a renormalization of the physical mass $M(\varepsilon)$ in terms of the ‘bare’ mass $M(\varepsilon=0)$. Therefore, the constant and linear pieces of this expansion are not observable.

However, the curvature $\frac{1}{2}(t-t_{0})^{2}d\omega^{(2)}/dt|_{t=t_{0}}$ cannot be redefined or scaled away. Therefore, within a sufficiently short window of time around the alignment time $t_{0}$, the deformation to the orbital phase is given by

[TABLE]

We use this to define the perturbative secular time $T_{\mathrm{sec}}(t_{0})$ at any instant $t_{0}$ via

[TABLE]

roughly the time to dephase by order one radian.

Thus we need to estimate $d\omega^{(2)}/dt$ from our simulation. Under the assumption of quasi-circular, adiabatic orbits, there is a one-to-one correspondence between the orbital frequency $\omega$ and orbital energy $E$. In other words, there exists a function of one variable, $E(\omega)$ or $\omega(E)$. Therefore from the chain rule we can find the time derivative

[TABLE]

This depends on the conservative sector through the frequency-dependence of orbital energy, $dE/d\omega$, and on the dissipative sector through the radiated power, $dE/dt$. Just as with the frequency, we expand the orbital energy in powers of $\varepsilon$,

[TABLE]

We can then use this to expand Eq. (51) in powers of $\varepsilon$. The $\mathcal{O}(\varepsilon^{2})$ piece is given by

[TABLE]

The prefactor $d\omega^{(0)}/dt$ is simply the background (GR) evolution of the orbital frequency. The first term in square brackets in Eq. (53) comes from the dissipative sector of the dynamics, since it depends on the radiated power $dE^{(2)}/dt$. The second term, meanwhile, comes from the conservative sector, as it depends on the correction to the orbital energy $E^{(2)}(\omega)$. Both of the factors in square brackets scale as $(\ell/GM)^{4}v^{4}$ Yagi et al. (2012a, 2013) for BBHs with spin. Plugging this scaling into Eq. (50) recovers $T_{\mathrm{sec}}\sim T^{\textrm{GR}}_{\mathrm{RR}}(\ell/GM)^{-2}v^{-2}$.

We find it useful to rewrite $dE^{(0)}/d\omega$ in the second term using the chain rule (51) to give

[TABLE]

Now we can discuss how to evaluate these factors from our numerical simulation and previously-known analytical results. The background energy flux $dE^{(0)}/dt$ comes from the numerical simulation via Eq. (34). We also have the background frequency evolution $d\omega^{(0)}/dt$ from the numerical simulation, via a time derivative of the coordinate orbital frequency.

The two $\mathcal{O}(\varepsilon^{2})$ quantities require approximations. In the dissipative sector, there are two contributions to $dE^{(2)}/dt$: the first from scalar radiation, and the second from gravitational radiation. We expect these to be the same order of magnitude. Since we do not have access to the gravitational radiation, we approximate that to within an order of magnitude,

[TABLE]

where $\dot{E}^{(\vartheta,2)}$ was given in Eq. (33). This is further justified during the inspiral, where the $\mathcal{O}(\varepsilon^{2})$ dissipative correction due to gravitational waves is higher-PN than the scalar radiation Yagi et al. (2012a).

In the conservative sector, we can approximate $E^{(2)}(\omega)$ from a post-Newtonian calculation Yagi et al. (2013); Stein and Yagi (2014). The (PN-approximate) correction to the orbital energy $E^{(2)}$ also has two pieces: the scalar binding energy and the metric-deformation binding energy. Again we are going to make an approximation and ignore the metric deformation piece, approximating

[TABLE]

where $E^{(\vartheta)}_{\mathrm{DD}}$ is the scalar dipole-dipole interaction. After accounting for a missing minus sign in Yagi et al. (2013); Stein and Yagi (2014), this is given by

[TABLE]

where again $\mu_{A}^{i}$ is the scalar dipole moment given in Eq. (27). In our case the spins are in the $\hat{z}$ direction, so the $(\mu_{A}\cdot n_{12})$ term vanishes. To leading PN order, we use the Kepler relation $\omega^{2}=GM/r_{12}^{3}$ and obtain

[TABLE]

where $\mu_{A}$ now refers to the $\hat{z}$ component. For $\omega$ we again use the coordinate orbital frequency from the simulation.

To summarize this calculation: we are approximating the secular breakdown time $T_{\mathrm{sec}}$ [Eq. (50)] by assuming a quasi-circular, adiabatic inspiral, and thus we compute $d\omega^{(2)}/dt$, Eq. (54). We approximate the dissipation $\dot{E}^{(2)}$ from only the scalar flux, Eq. (55). We approximate the conservative correction $E^{(2)}(\omega)$ from the post-Newtonian scalar dipole-dipole interaction, Eq. (56).

In Fig. 6 we plot $(\ell/GM)^{2}T_{\mathrm{sec}}(t_{0})$, the time to secularly dephase by about $\sim 1$ radian, around various alignment times $t_{0}$. We have checked that at early times, this numerical estimate agrees with an analytic PN estimate. As expected, $T_{\mathrm{sec}}$ is parametrically longer than the GR radiation-reaction time. The time window for secular validity shrinks approaching merger, but does not vanish.

The value of $T_{\mathrm{sec}}$, and hence secular regime of validity, is smallest near merger. For the spin 0.3 simulation, just before merger, we find the time to dephase by about 1 radian from the GR background is $T_{\mathrm{sec}}\sim 15~{}GM(\ell/GM)^{-2}$. If Advanced LIGO detects a gravitational waveform of length, say, $200~{}GM$, then a perturbative calculation without MSA/renormalization would be valid for $(\ell/GM)\lesssim 1/4$. For longer waveforms or larger values of $(\ell/GM)$, MSA or renormalization would be required. However, larger values of $(\ell/GM)$ will be very close to the limit on the instantaneous regime of validity, Fig. 5.

III.5 Detectability and bounds estimates

We now turn to the issue of how well Advanced LIGO/Virgo would be able to detect or bound the effects of dynamical Chern-Simons gravity from observations of a binary black hole merger. As we do not yet have metric waveforms [that arise at $\mathcal{O}(\varepsilon^{2})$], we make order-of-magnitude projections of detectability and bounds from the dephasing estimates in the previous section.

Suppose that LIGO detects a gravitational waveform similar to one of those we have simulated, with approximately 5 cycles of inspiral in band before merger—similar to GW150914 Abbott et al. (2016a), with a total mass approximately $M\approx 60M_{\odot}$. Such a detection would come with errors due to noise and calibration uncertainty; let us define the overall waveform phase uncertainty $\sigma_{\phi}$. Let us further assume that the dCS corrections to the full waveform are not degenerate with redefining ‘bare’ binary parameters. Upon detection there are two distinct possibilities: (i) the detected waveform is consistent with GR predictions; or (ii) the detection is inconsistent with any point in the GR parameter space.

In the case of consistency, we would be able to place bounds on the size of $\ell$. Crudely, we would be able to say

[TABLE]

where the factor of two comes from the gravitational wave being at twice the orbital frequency. If we have consistency with GR, then the quadratic approximation for $\Delta\phi$ in Eq. (48) holds.

We plot the quadratic approximation to the orbital phase difference (relative to GR) in Fig. 7. By taking the maximum value of $\Delta\phi$ over the length of the waveform, and taking into account the scaling with $(\ell/GM)^{4}$, we can derive a projected bound on $\ell$. For example, from the spin 0.3 simulation and $M\approx 60M_{\odot}$, we would find

[TABLE]

and from the spin 0.1 simulation,

[TABLE]

The spin 0.0 simulation would only give $(\ell/GM)\lesssim 1.4(\sigma_{\phi}/0.1)^{1/4}$. Such a bound would be past the instantaneous regime of validity limit during merger for this simulation (see Fig. 5). It is not internally self-consistent to use this perturbative result to claim a constraint on the regime past perturbative validity, so conservatively, no statement can be made. The higher spin simulations do not suffer from this problem.

These bounds forecasts can immediately be turned around into detectability forecasts. We can forecast that dynamical Chern-Simons corrections would be detectable in a $M\approx 60M_{\odot}$ binary with parameters consistent with our spin 0.3 simulation if $\ell\gtrsim 11$ km, and similarly for the spin 0.1 simulation if $\ell\gtrsim 18$ km.

We can draw three simple lessons on detectability and bounds from these results. First, better phase sensitivity (smaller $\sigma_{\phi}$) is an obvious way to improve the odds of detectability, or place stronger bounds. This comes from improved detector sensitivity, but also from higher signal-to-noise ratio (SNR) events. Second, at fixed phase sensitivity, lower-mass events would be better than higher mass events, to a point. Lower mass events obviously have smaller $GM$, but they also spend more time in band, and thus have more time for dephasing. There is a tradeoff, though, because lower mass events are quieter, and also because most of the dephasing comes right before merger—so the mass must be high enough for merger to be in band. Finally, we can easily see that higher spin systems would lead to stronger constraints or a better chance of detecting dCS effects.

Let us compare our projected bounds to those appearing previously in the literature. Ali-Haïmoud and Chen Ali-Haïmoud and Chen (2011) used solar system data from Gravity Probe B and the LAGEOS satellites to constrain the characteristic length scale to $\ell\lesssim 10^{8}~{}\mathrm{km}$. Yagi, Yunes and Tanaka Yagi et al. (2012b) found a similar bound from table-top experiments. This is comparable to the curvature radius in the solar system.

Yunes and Pretorius Yunes and Pretorius (2009) applied a precession calculation from the extreme mass-ratio limit to PSR J0737–3039 to estimate a constraint of $\ell\lesssim 10^{4}~{}\mathrm{km}$. However, this calculation missed some effects (such as the scalar binding energy), and the mass ratio of PSR J0737–3039 is very close to 1. Moreover, the curvature radius at the surface of one of the NSs in this system should be order $\sim 10$ km, which means there is room between $10-10^{4}$ km where $\ell$ could be large compared to the curvature length, and thus the calculation would not be internally self-consistent. Yagi, Stein, Yunes, and Tanaka Yagi et al. (2013) performed a more careful analysis, using post-Newtonian theory for binary NS systems. They concluded that even PSR J0737–3039, with its high orbital velocity and exquisite timing, would not be able to yield a constraint on dCS for the foreseeable future, and that gravitational wave measurements would be the best hope.

Yagi, Yunes, and Tanaka Yagi et al. (2012c) used post-Newtonian calculations to project the level of constraints that might arise from second and third generation GW detectors. If next-generation detectors such as Einstein Telescope Punturo et al. (2010) were to observe binary black hole inspirals consistent with GR, then YYT project a bound of $\ell\lesssim\mathcal{O}(10-100)~{}\mathrm{km}$. Second-generation ground-based detectors could place a similar constraint. The only caveat here is that YYT use post-Newtonian estimates, stopping at the ISCO frequency, for systems that would be seen not only in the inspiral, but also in the merger and ringdown, where PN is invalid. The additional SNR contributed by merger and ringdown will likely improve constraints.

Stein and Yagi Stein and Yagi (2014) projected a number of constraints on $\ell$ based on both pericenter precession in pulsar binaries and gravitational wave measurements. For a LIGO detection of a $(10+11)M_{\odot}$ BBH inspiral, consistent with GR, at an SNR of 30, they projected a bound on the order of $\ell\lesssim 10~{}\mathrm{km}$. Note that this is the same order of magnitude as the projected bound we estimate here.

Finally, Stein Stein (2014a) projected a bound based on the observations of the black hole candidate GRO J1655–40. Assuming observations were consistent with GR, Stein projected a constraint of $\ell\lesssim 22~{}\mathrm{km}$. However, such a constraint would require (for example) accretion disk modeling in the presence of the dCS correction, which has not been simulated.

IV Discussion and future work

In this study, we have performed the first fully nonlinear inspiral, merger, and ringdown numerical simulations of a binary black hole system in dynamical Chern-Simons gravity. These are the first BBH simulations in a theory besides general relativity and standard scalar-tensor gravity. BBH in scalar-tensor is identical to that in GR, unless one imposes an external scalar field gradient Healy et al. (2012); Berti et al. (2013). Therefore these are also the first numerical simulations in a theory where the BBH dynamics differ from GR under ordinary initial and boundary conditions.

The “full” equations of motion for dCS, and many other corrections to GR, probably lack a well-posed initial value formulation. This is not an obstacle if the corrections are treated as being a low-energy effective field theory. In Sec. II we formulated a perturbation scheme which guarantees a well-posed initial value problem. We stress that this scheme is applicable not just to dCS, but also any deformation of general relativity which has a continuous limit to GR.

We performed fully nonlinear numerical simulations through order $\mathcal{O}(\varepsilon^{1})$ in the perturbation scheme. We simulated binaries with mass ratio $q=3$ and aligned spins with equal dimensionless spin parameters $\chi_{1}=\chi_{2}$, taking on three values, $\chi=0.0,0.1,0.3$. The background ($\varepsilon^{0}$) metric radiation and perturbative ($\varepsilon^{1}$) scalar radiation waveforms are presented in Sec. III.2. We found good agreement with PN waveform predictions Yagi et al. (2012a, 2013) during the early inspiral.

We have also discovered new phenomenology in dCS. In agreement with PN predictions, dCS does not suffer from dipole radiation during the early inspiral. However, during merger, there is a burst of dipole radiation. This phenomenon can only be studied with full numerical simulations.

We extracted energy fluxes in Sec. III.3, finding good agreement with PN at early times. We found that the scalar field’s $\mathcal{O}(\varepsilon^{2})$ energy flux during the inspiral was approximately $10^{-6}(\ell/GM)^{4}$ times smaller than the corresponding $\mathcal{O}(\varepsilon^{0})$ (GR) energy flux for the highest spin simulation, rising to a $10^{-3}(\ell/GM)^{4}$ fraction of GR during merger. This energy flux enters into our detectability estimate.

Since we use a perturbative scheme, it is important to understand where perturbation theory breaks down. In Sec. III.4.1 we estimated the maximum values of $\ell/GM$ for the perturbation theory to be convergent at each time during the simulation. During the inspiral and ringdown, the regime of validity agrees with estimates from Stein (2014a). The tightest bound on the instantaneous regime of validity comes during merger, and is comparable for spinning and non-spinning black hole mergers, close to $\ell/GM\lesssim 1$.

The additional radiation in the scalar field $\vartheta^{(1)}$ leads to a secular drift in orbital phase between the “true” orbital dynamics and the GR background from which we perturb. Therefore even if perturbation theory is instantaneously under control, the perturbative solution will dephase after a sufficiently long time. We numerically estimated this dephasing time in Sec. III.4.2, and it agrees with post-Newtonian scaling at early times. At times approaching merger, the dephasing time becomes shorter, but remains nonzero.

This dephasing calculation served as the basis for estimating detectability and predicting bounds that LIGO would be able to place on $\ell$, in Sec. III.5. For $q=3$, $M\approx 60M_{\odot}$, and aligned dimensionless spins of $\chi_{1}=\chi_{2}=0.3$, we estimated that a GR-consistent detection would yield a bound of

[TABLE]

where $\sigma_{\phi}$ is LIGO’s statistical phase uncertainty on the detected waveform, which depends on the SNR of the detection. Conversely, an $\ell$ above this value would be detectable by LIGO. Lower spins lead to poorer detectability and/or bounds. Better bounds come from three places: (i) improved phase sensitivity (higher SNR), (ii) lower mass events (while keeping merger in band), and (iii) higher spin systems.

IV.1 Future work

The natural next step in this program is to continue to the order $\varepsilon^{2}$ system, as outlined in Sec. II.2.3. This is the lowest order where gravitational radiation is modified, and would involve solving for $h_{ab}^{(2)}$, which is sourced by $g_{ab}^{(0)}$ and $\vartheta^{(1)}$.

With the solution for the deformation to the metric $h_{ab}^{(2)}$, we will be able to directly compare dCS predictions against LIGO data. We will also have a more complete assessment of the convergence of the perturbation scheme.

Comparing dCS predictions against LIGO data will yield the first direct bounds on the theory from the strong-field, dynamical regime of gravity. To do so will involve extending GR parameter estimation Veitch et al. (2015) with one additional parameter, $\ell$, which will be simultaneously inferred or constrained from the data.

A complete analysis would involve thorough exploration of the 7-dimensional parameter space of quasicircular BBHs (mass ratio and two spin vectors; the $\ell$ dependence is analytic in the perturbative approach). For example, in this work, we have focused on aligned-spin binaries in order to simplify comparisons with analytic predictions. The scalar energy flux in the case of misaligned binaries may be an order of magnitude larger than in the spin aligned case (see Yagi et al. (2012a) and the erratum). Building a surrogate waveform model Blackman et al. (2015, 2017) would simultaneously allow for an efficient exploration of parameter space and efficient parameter estimation/constraints with LIGO data.

The standard perturbation theory approach we used here will be sufficient if we find that the dephasing time is long compared to LIGO signals. However, if we need to extend the secular regime of validity, some form of multiple-scale analysis Bender and Orszag (1978) or dynamical renormalization group Chen et al. (1996); Galley and Rothstein (2017) approach will be required.

Finally, let us emphasize that our approach is not limited to dynamical Chern-Simons gravity: dCS is a proof of principle. Any theory with a continuous limit to GR can be treated with the same scheme, and reusing a large fraction of the code. In particular, we will consider EdGB and a class of theories proposed in Endlich et al. (2017). Switching from dCS to another theory will only involve changing source terms that appear on the right hand sides of the differential equations we are solving numerically.

Acknowledgements.

We would like to thank Yanbei Chen, Chad Galley, Luis Lehner, Robert McNees, Frans Pretorius, Thomas Sotiriou, Saul Teukolsky, Helvi Witek, Kent Yagi, and Nico Yunes for many valuable conversations. This work was supported in part by the Sherman Fairchild Foundation, the Brinson Foundation, and NSF grants PHY–1404569 and AST–1333520 at Caltech. MO gratefully acknowledges the support of the Dominic Orr Graduate Fellowship at Caltech. Computations were performed on the Zwicky cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by NSF award PHY–0960291. Some calculations used the computer algebra system Mathematica, in combination with the xAct/xTensor suite xAc ; Martín-García (2008); Brizuela et al. (2009). The figures in this paper were produced with matplotlib Hunter (2007).

Appendix A Scalar field evolution formulation

In this appendix, we discuss the numerical evolution scheme for a (massless) Klein-Gordon field, denoted by the code variable $\Psi$, in greater detail. This is an update of the system described in Holst et al. (2004), which did not include the “$\gamma_{1}\gamma_{2}$” constraint-damping term (see below). The basic equation we are simulating is

[TABLE]

for some prescribed source term $S$ (in this work, the source term is the Pontryagin density ${}^{*}\!RR$).

We first review the 3+1 ADM formalism for the foliation of a spacetime into spatial slices, as used in numerical relativity Baumgarte and Shapiro (2010). We decompose the metric as

[TABLE]

where $g_{ab}$ is the spacetime metric, $n_{a}$ is a timelike unit one-form normal to the spatial slice with $n_{a}n^{a}=-1$, and $\gamma_{ab}$ is the induced spatial metric and projector, with $n^{a}\gamma_{ab}=0$. In ADM variables, the timelike unit normal can be written in terms of a lapse, $\alpha$, and shift $\beta^{i}$, as $n^{a}=(\alpha^{-1},-\alpha^{-1}\beta^{i})$.

We work with the Spectral Einstein Code (SpEC), which uses the generalized harmonic formulation of general relativity, and evolves a symmetric hyperbolic first-order system of metric variables $g_{ab}$, $\Phi_{iab}=\partial_{i}g_{ab}$ and $\Pi_{ab}=-n^{c}\partial_{c}g_{ab}$ Lindblom et al. (2006).

We similarly define a set of first-order variables for the scalar field $\Psi$ as

[TABLE]

From these definitions and the equality of mixed partial derivatives, we can create a system of constraints which vanish in the continuum limit, and which an accurate evolution of the system will satisfy to within some tolerance:

[TABLE]

In Eq. (70) the indices $j,k$ are summed and $[ijk]$ is the completely antisymmetric Levi-Civita symbol, with $[123]=+1$.

The evolution equation (65) thus becomes a set of first-order time evolution equations for $\{\Psi$, $\Phi_{i},\Pi\}$. However, in order to prevent numerical errors in the constraints from making the evolution unstable, we follow what is done in the metric system and add specific multiples of the constraints to the evolution equations. These combinations of constraints are chosen so as to ensure that the system is symmetric hyperbolic and that the constraints are damped out, ensuring a well-posed evolution scheme. The evolution equation for $\Psi$ is thus

[TABLE]

where the first terms come from the definitions of $\Phi_{i}$ and $\Pi$, and the last term is a constraint damping term with coefficient $\gamma_{1}$. The evolution equation for $\Phi_{i}$ is

[TABLE]

where the term with $\gamma_{2}$ is a constraint damping term, and all other terms come from definitions of the first-order variables and equality of mixed partial derivatives. Finally, the evolution equation for $\Pi$ is

[TABLE]

where $K$ is the trace of the extrinsic curvature, $\Gamma^{m}\equiv g^{ab}\Gamma^{m}{}_{ab}$ is a specific contraction of the Christoffel connection coefficients, $S$ is the source term, and the $\gamma_{1}\gamma_{2}$ term is the appropriate constraint-damping term to keep the system symmetric hyperbolic.

This “$\gamma_{1}\gamma_{2}$” term was not included in the previous description Holst et al. (2004), but it is required if both $\gamma_{1}$ and $\gamma_{2}$ are non-zero. The parameters $\gamma_{1}$ and $\gamma_{2}$ play the same role in the damping and characteristic analysis of this Klein-Gordon system as they do in the generalized harmonic system Lindblom et al. (2006). Specifically, in order for the constraint $C_{i}^{(1)}$ to be damped, we must have $\gamma_{2}>0$ (satisfying the constraint $C_{i}^{(1)}$ implies satisfaction of the constraint $C_{i}^{(2)}$). The choice $\gamma_{1}=-1$ makes the system linearly degenerate. In practice we set the values of $\gamma_{1}$ and $\gamma_{2}$ to match those of the generalized harmonic evolution of the metric variables, so that the characteristic speeds of the metric and scalar field systems agree.

The scalar field variables, like the metric variables, are represented spectrally. In order to reduce the amount of numerical noise in the system, we apply the same filters we use for the metric variables to the scalar field system, namely filtering the top $4$ tensor spherical harmonics and using an exponential Chebyshev filter for the radial piece.

In order to assess the accuracy of the simulations, we evaluate the constraints that the generalized harmonic evolution system must satisfy Lindblom et al. (2006), as well as the constraints for the first-order scalar field system given by Eqs. (69) and (70). We combine these constraints, contracting with a Euclidean metric to give a constraint energy as

[TABLE]

Since the code is spectral, we check for exponential convergence of these constraint energies as we increase the number of angular and radial basis functions per subdomain (and hence the resolution). We plot the convergence of the $L^{\infty}$ norm of the constraint energies for the highest spin simulation of this study, which has the greatest level of constraint violation, in Fig. 8. We find that the error decreases exponentially with resolution. The lower spin simulations have similar qualitative behavior.

Appendix B Pontryagin density in 3+1 split

Since numerical relativity computations are formulated in a 3+1 split, we must compute the scalar field’s source term—the Pontryagin density—in terms of 3 dimensional quantities. First, it is straightforward to verify

[TABLE]

where $C_{abcd}$ is the Weyl tensor, and its left dual is ${}^{*}C^{abcd}\equiv\frac{1}{2}\epsilon^{abef}C_{ef}{}^{cd}$. Thus we do not need to consider all of Riemann, but only its trace-free part, Weyl. The Pontryagin density is completely insensitive to the Ricci part of curvature.

In a 4-dimensional numerical relativity simulation, it is especially convenient to decompose Weyl into its electric and magnetic parts, defined as

[TABLE]

The minus sign in (77) follows the conventions of Owen et al. (2011); Nichols et al. (2011) and the implementation in SpEC SpE , though much of the literature has a plus sign. From the symmetries of Weyl, the two tensors $E_{ab}$ and $B_{ab}$ are both symmetric ($E_{ab}=E_{(ab)}$ and $B_{ab}=B_{(ab)}$), purely spatial ($E_{ab}n^{a}=0=B_{ab}n^{a}$), and trace-free ($E^{a}{}_{a}=0=B^{a}{}_{a}$). We may also write an inversion formula for Weyl in terms of $E_{ab}$ and $B_{ab}$ (thanks to Alfonso García-Parrado for bringing this inversion formula to our attention),

[TABLE]

where the operator ${\tiny\young(ac,bd)}$ is a projector that imposes the symmetries of the Riemann tensor ($R_{abcd}=R_{[ab][cd]}=R_{cdab}$). Here we have the induced 3-dimensional volume element,

[TABLE]

For coordinate component calculations, we use the conventions where $\epsilon_{abcd}=+\sqrt{-g}[abcd]$ where $[abcd]$ is the alternating symbol, with $[0123]=+1$ (see e.g. Misner et al. (1973)). We also have $\epsilon^{abcd}=-[abcd]/\sqrt{-g}$, and similar conventions for the 3-dimensional volume element: $\epsilon_{abc}=\sqrt{\gamma}[abc]$ and $\epsilon^{abc}=[abc]/\sqrt{\gamma}$ (this makes use of the identity $\sqrt{-g}=\alpha\sqrt{\gamma}$).

With this above decomposition, it is easy to verify that the Pontryagin density can be expressed simply in terms of the electric and magnetic parts of Weyl,

[TABLE]

Thus all that remains is to compute $E_{ab}$ and $B_{ab}$ from other quantities. Finding these expressions for $E$ and $B$ comes from the standard Gauss-Codazzi-Mainardi (GCM) equations (see Baumgarte and Shapiro (2010) for a didactic explanation). After using the GCM equations, for the electric Weyl tensor we find

[TABLE]

Here ${}^{(3)}\!R_{ab}$ is the spatial 3-Ricci tensor while ${}^{(4)}\!R_{ab}$ is the 4-Ricci tensor, and $K_{ab}$ is the extrinsic curvature of the spacelike hypersurface. The second line of (81) contains 4-Ricci terms which would vanish if the 4-metric was Ricci-flat, for example if it solves the Einstein equations in vacuum. These terms were not included in e.g. Kidder et al. (2005).

Meanwhile, for the magnetic Weyl tensor we find the simple expression

[TABLE]

where $D_{a}$ is the covariant derivative induced on the 3-surface which is compatible with the 3-metric, $D_{a}\gamma_{bc}=0$.