Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences:(I) Inspiral-Merger

Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences:
(I) Inspiral-Merger

Emanuele Berti E. Berti Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA K. Yagi Department of Physics, University of Virginia, Charlottesville, Virginia 22904, USA N. Yunes eXtreme Gravity Institute, Department of Physics, Montana State University, Bozeman, Montana 59717, USA    Kent Yagi E. Berti Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA K. Yagi Department of Physics, University of Virginia, Charlottesville, Virginia 22904, USA N. Yunes eXtreme Gravity Institute, Department of Physics, Montana State University, Bozeman, Montana 59717, USA    Nicolás Yunes E. Berti Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA K. Yagi Department of Physics, University of Virginia, Charlottesville, Virginia 22904, USA N. Yunes eXtreme Gravity Institute, Department of Physics, Montana State University, Bozeman, Montana 59717, USA
Received: date / Accepted: date
Abstract

The observation of the inspiral and merger of compact binaries by the LIGO/Virgo collaboration ushered in a new era in the study of strong-field gravity. We review current and future tests of strong gravity and of the Kerr paradigm with gravitational-wave interferometers, both within a theory-agnostic framework (the parametrized post-Einsteinian formalism) and in the context of specific modified theories of gravity (scalar-tensor, Einstein-dilaton-Gauss-Bonnet, dynamical Chern-Simons, Lorentz-violating, and extra dimensional theories). In this contribution we focus on (i) the information carried by the inspiral radiation, and (ii) recent progress in numerical simulations of compact binary mergers in modified gravity.

Keywords:
Modified Gravity Gravitational Waves Compact Binary Systems

1 Introduction

Why would you modify Einstein’s beautiful theory of gravity, General Relativity (GR)? This question is asked often in the context of experimental relativity. The driving force behind modified gravity is not to discard Einstein’s theory, but rather to address some of the unsolved physics problems that Einstein’s theory has brought to the forefront. The late-time expansion of the Universe and dark energy, the rotation curves of galaxies and dark matter, the matter-antimatter asymmetry in the early Universe, the information loss problem and the existence of singularities, the exponential expansion of the early Universe and inflation are just a few of these problems. One could of course assert that these problems are resolved by (highly fine-tuned) initial conditions or the existence of a new form of some yet-undetected weakly interacting particle. Another avenue is to postulate that they are resolved because Nature is better described by a modified gravity theory that reduces to Einstein’s in some well-studied scenarios (like in the Solar System or in binary pulsars), but that introduces key modifications in other scenarios (like on large scales, in the early Universe or near black holes).

One way to classify these theory-inspired modifications to gravity is through the fundamental pillars of GR that they violate or deform. Einstein’s theory can be thought of as descending from three organizing principles: diffeomorphism invariance, the (strong) equivalence principle, and spacetime as a four-dimensional continuum. Indeed, one can argue that these principles essentially require spacetime to be curved, and that the need for the theory to reduce to Newtonian gravity in the weak-field limit and the Bianchi identities then lead to the Einstein equations. The resulting theory is necessarily parity invariant (on any spacelike hypersurface) and gravity is described by the curvature of a rank-2 tensor field that is free to propagate at the speed of light, i.e. in a quantized picture it would be described by a massless, spin-2 particle. Given these pillars, one can imagine breaking or deforming them to address some of the unsolved problems mentioned above. For example, the observed matter-antimatter asymmetry that develops in the radiation-dominated era requires additional sources of parity violation by the Sakharov conditions Sakharov:1967dj (), which could be addressed by introducing gravitational parity-violating interactions in the Einstein-Hilbert action Alexander:2004xd (); Alexander:2009tp (). Another example is the late-time acceleration of the Universe, which instead of being described by a cosmological constant may perhaps be explained by an additional tensor field that interacts with the gravitational field, as in bi-gravity theories Pilo:2011zz (); Paulos:2012xe ().

Given so many modified gravity alternatives that have been postulated to address these unsolved problems, how do we decide which one, if any, is the best description of Nature? Experiments and observations are the fool-proof method of choice. Many of these modified theories can be straightforwardly ruled out with current Solar System Will:2014kxa (), binary pulsar Stairs:2003eg () and cosmological observations. For example, the Brans-Dicke-Jordan-Fierz flavor of scalar-tensor theories Brans:1961sx (); Fierz:1956zz (); Jordan:1959eg () is incompatible with Shapiro time-delay observations by the Cassini spacecraft Bertotti:2003rm (), unless the coupling constant is tuned to be sufficiently small. These observations, together with the cosmological evolution of scalar fields in the scalar-tensor theories studied by Damour and Esposito-Farèse Damour:1992we (); Damour:1993hw (); Damour:1996ke (); Damour:1992kf (); Damour:1993id (), prohibit the spontaneous scalarization of neutron stars Sampson:2014qqa (); Anderson:2016aoi (), unless one fine tunes the cosmological initial conditions of the evolution of the scalar field.

But even after imposing the requirement that modified theories must pass all current experimental tests, there still remain a large group of them that is only weakly constrained today. One example are effective field theories of gravity that introduce a scalar field sourced by singularities in the spacetime curvature and not by matter Yunes:2009hc (); Alexander:2009tp (); Yunes:2011we (). Such theories will pass Solar System tests because of the weak curvature of spacetime in the Solar System, which thus prevents the scalar field from being strongly excited AliHaimoud:2011bk (); AliHaimoud:2011fw (). These theories will also pass binary pulsar tests because neutron stars do not possess singularities, and thus, the scalar field is not strongly activated Yagi:2011xp (); Yagi:2015oca (). One must then rely on new observations, observations that sample the extreme gravity regime, where the spacetime is strongly curved and highly dynamical. One example are gravitational wave (GW) observations of the coalescence of binary compact objects, like black holes and neutron stars.

GW tests are quite different from other tests of gravity that have been carried out to date. As argued above, GWs are unique probes of the extreme gravity regime, being sensitive both to the propagation and the generation of these waves. Moreover, GWs are weakly interacting, and thus, they travel essentially unimpeded to our detectors on Earth, without being affected by intervening matter. For example, the accretion disk of black holes introduces a very small modification to the GWs emitted when black holes collide, except perhaps for the most massive disks Kocsis:2011dr (); Yunes:2011ws (); Hayasaki:2012qn (); Barausse:2014tra (). Finally, GW tests are localized in spacetime, sampling the gravitational interaction everywhere inside the wave’s lightcone. One could thus imagine that as more GW observations are made, one can begin to build constraint maps that represent the verification of the pillars of GR everywhere in the sky. That is, one could envision a GW sky map (similar to that which depicts anisotropies in the cosmic microwave background) that depicts the sky locations of thousands of GW detections and the constraints placed on a given deformation of GR at that particular location, perhaps even allowing constraint contours in sky location.

The era of GW physics has begun, heralded by the first detections of Advanced LIGO (aLIGO) and Virgo Abbott:2016blz (); Abbott:2016nmj (); Abbott:2017vtc (); Abbott:2017oio (); TheLIGOScientific:2017qsa (); Abbott:2017gyy (). These observations have already confirmed that indeed GWs exist, that they are predominantly quadrupolar in nature, and that they propagate at essentially the speed of light in a Lorentz-invariant fashion. Such tests therefore confirm some of the pillars of GR, but the era of precision experimental relativity with GWs is only beginning. In the next few years, the LIGO/Virgo collaboration will continue to make detections at an ever increasing rate with larger and larger signal-to-noise ratios. In the near future new detectors (such as KAGRA in Japan and LIGO-India in India) will join the network, which will allow for precise probes of the polarization content of GWs. Work has already begun to consider upgrades to third-generation detectors, and in the early 2030s space-based detectors will be launched Audley:2017drz (). At this stage, the stacking of hundreds or thousands of GW observations will yield unprecedented tests of Einstein’s theory of gravity in the mostly unexplored extreme gravity regime.

This contribution is organized as follows. In Section 2 we briefly review models of modified gravity that have been explored in some detail in the context of strong-gravity tests, including scalar-tensor theories, Einstein-dilaton Gauss-Bonnet gravity, dynamical Chern-Simons gravity, Lorentz-violating theories, and theories involving extra dimensions. In Section 3 we present generic and theory-specific tests of modified gravity with compact binary inspirals, reviewing first current constraints, and then future prospects. In Section 4 we discuss the current state of numerical simulations of compact binary mergers in modified gravity. We conclude in Section 5 with a list of important open problems. In Appendix A we derive the scalar charge in decoupled dynamical Gauss-Bonnet gravity for black holes with arbitrary rotation. The final expression was shown in Yunes:2016jcc (), but the details of the derivation (and the expression for the quadrupole scalar charge) are presented here for the first time. In Appendix B we briefly review a “dynamical no-hair theorem” for black holes in scalar-tensor gravity.

2 Modified Theories of Gravity

2.1 Scalar-tensor Theories

One of the simplest extensions of GR is scalar-tensor gravity. In this class of theories, one or more scalar degrees of freedom are included in the gravitational sector through a non-minimal coupling: the Ricci scalar in the Einstein-Hilbert action is multiplied by a function of the scalar field(s). These theories are well motivated: scalar fields with non-minimal couplings to gravity appear (e.g.) in string theory Polchinski:1998rq (), in Kaluza-Klein-like theories Duff:1994tn () and in braneworld scenarios Randall:1999ee (); Randall:1999vf (). These models have also been extensively discussed in a cosmological context Clifton:2011jh (). Due to their simplicity, scalar-tensor theories are a good framework to study the strong-field dynamics of modified gravity. There are many extensive reviews on the subject Damour:1992we (); Chiba:1997ms (); 2003sttg.book…..F (); Faraoni:2004pi (); Sotiriou:2015lxa (); Berti:2015itd (), so we will only sketch the main ingredients of these theories – first in their original form, and then considering generalizations that have attracted some attention in the past few years.

2.1.1 “Bergmann-Wagoner” scalar-tensor theories

The most general action of scalar-tensor theory with one scalar field which is at most quadratic in derivatives of the fields was studied by Bergmann and Wagoner Bergmann:1968ve (); Wagoner:1970vr (). With an appropriate field redefinition, the so-called Jordan-frame action for the theory can be cast in the form:

(1)

where and are arbitrary functions of the scalar field , and is the action of the matter fields . When is constant and , the theory reduces to (Jordan-Fierz-)Brans-Dicke gravity Jordan:1959eg (); Fierz:1956zz (); Brans:1961sx ().

The Bergmann-Wagoner theory (1) can be expressed in a different form through a scalar field redefinition and a conformal transformation of the metric . In particular, fixing , the action (1) transforms into the Einstein-frame action111A straightforward generalization consists of coupling gravity with more than one scalar field. Then the action (1) has the more general form Damour:1992we () where are functions of the scalar fields (). The scalar fields live on a manifold (sometimes called the target space) with metric . This action is invariant not only under space-time diffeomorphisms, but also under target-space diffeomorphisms, i.e. scalar field redefinitions. The geometry of the target space can affect the dynamics and the structure of compact objects Horbatsch:2015bua ().

(2)

where and are the determinant and Ricci scalar of , respectively, and the potential . In the Einstein frame the scalar field is minimally coupled in the gravitational sector, but there is non-minimal coupling in the matter sector of the action: particle masses and fundamental constants depend on the scalar field.

The actions (1) and (2) are different representations of the same theory Flanagan:2004bz (); Sotiriou:2007zu (), and the frame choice is often dictated by computational convenience: for instance, in vacuum the Einstein-frame action (2) formally reduces to GR with a minimally coupled scalar field. However, since the scalar field is minimally coupled to matter in the Jordan frame, test particles follow geodesics of the Jordan-frame metric.

Jordan-frame and Einstein-frame quantities are related by , , where . The theory is fixed by choosing the function – or, equivalently, – and the form of the scalar potential. Neglecting the scalar potential corresponds to neglecting the cosmological term, the mass of the scalar field and any possible scalar self-interaction. In an asymptotically flat spacetime the scalar field tends to a constant at spatial infinity, corresponding to a minimum of the potential. Taylor expanding around yields, to lowest orders, a cosmological constant and a mass term for the scalar field Wagoner:1970vr (); Alsing:2011er ().

Scalar-tensor theory with a vanishing scalar potential is characterized by a single function . The expansion of this function around the asymptotic value can be written in the form

(3)

As mentioned above, the choice  constant (i.e.,  constant) corresponds to Brans-Dicke theory. A more general formulation, proposed by Damour and Esposito-Farèse, is parametrized by and  Damour:1993hw (); Damour:1996ke (). Another simple variant is massive Brans-Dicke theory, in which is constant, but the potential is non-vanishing and has the form , so that the scalar field has a mass .

In the Einstein frame, the field equations are

(4a)
(4b)

where is the Einstein-frame stress-energy tensor of matter fields and is its trace Damour:1992we (). Equation (4b) shows that determines the strength of the coupling of the scalar field to matter Damour:1995kt ().

Astrophysical observations set bounds on the parameter space of scalar-tensor theories. In the case of Brans-Dicke theory, the best observational bound () comes from the Cassini measurement of the Shapiro time delay. In the more general case with , current constraints on have been obtained by observations of binary neutron star and neutron star-white dwarf binary systems Freire:2012mg (); Shao:2017gwu (), as well as from observations of the Shapiro time delay upon accounting for the cosmological evolution of the scalar field Sampson:2014qqa (); Anderson:2016aoi (). Observations of compact binary systems also constrain massive Brans-Dicke theory, leading to exclusion regions in the plane Alsing:2011er ().

2.1.2 Horndeski gravity

The most general scalar-tensor theory with second-order field equations (and one scalar field) is Horndeski gravity Horndeski:1974wa (). The action of Horndeski gravity can be written in terms of Galileon interactions Deffayet:2011gz () as

(5)

where and the ’s () are functions of the scalar field and of its kinetic term , and are derivatives of with respect to the kinetic term . For a particular choice of these functions Kobayashi:2011nu (); Maselli:2015yva () this theory coincides with Gauss-Bonnet gravity, which will be discussed next.

It was recently understood that Horndeski gravity is a subclass of all higher-order scalar-tensor theories that contain a single scalar mode. The crucial ingredient that singles out higher-order theories with a single scalar degree of freedom is the degeneracy of their Lagrangian, and now these theories are commonly referred to as Degenerate Higher-Order Scalar-Tensor (DHOST) theories Langlois:2015skt (); Langlois:2017mdk (). Theories similar to Horndeski gravity have also been constructed for vector (Proca) fields Heisenberg:2014rta (); Heisenberg:2016eld ().

2.2 Einstein-dilaton Gauss-Bonnet Gravity

Einstein-dilaton Gauss-Bonnet (EdGB) gravity is a string-inspired theory that acquires a quadratic-curvature correction to the Einstein-Hilbert action . Motivated from heterotic superstring theory, the action for EdGB gravity is given by Maeda:2009uy ()

(6)

where is the matter action while and are the quadratic-curvature correction and the kinetic term for the scalar field , respectively, given by

(7)

Here is the determinant of the metric, and are coupling constants while

(8)

with , and representing the Ricci scalar, Ricci tensor and Riemann tensor, respectively. This specific combination of curvature quantities ensures derivatives in the field equations to be at most of second order, but the well-posedness of this theory is still an open question. We assume that the scalar field is dimensionless in geometrical units, so has units of length squared, while is dimensionless. For simplicity, we set the scalar field potential to zero. Current constraints on come from the orbital decay rate of low-mass X-ray binaries with a black hole kent-LMXB (), the existence of stellar-mass black holes Pani:2009wy () and the maximum mass of neutron stars Pani:2011xm (). Black hole solutions have been constructed numerically for static configurations Kanti:1995vq (); Torii:1996yi (); Alexeev:1996vs () and rotating configurations Pani:2009wy (); Kleihaus:2011tg (); Kleihaus:2014lba () (see also Kokkotas:2017ymc () for approximate analytic solutions for static black holes).

One can treat this theory as an effective field theory since we neglect terms higher than cubic order in curvature – or terms of – in the action. Namely, one can assume that is much smaller than . This can be done by perturbing the action about and retaining terms up to linear order. Since can be written in a total derivative form upon variation of the action, the leading term where does not couple to the scalar field does not contribute to the field equations. Thus, we define the decoupled dynamical Gauss-Bonnet (DGB) gravity Yagi:2015oca (), where is replaced with

(9)

This theory is shift-invariant, in the sense that the field equations are invariant under a constant shift in the scalar field. In fact, this theory is a special case of Horndeski gravity Kobayashi:2011nu (), the most generic class of scalar-tensor theories with second-order field equations Horndeski:1974wa (); Deffayet:2009mn (); Deffayet:2011gz () as already explained in Sec. 2.1.2. More specifically it falls within the class of shift-symmetric Horndeski theories, and it is free of Ostrogradski ghosts Crisostomi:2017ugk (). The well-posedness of the theory is still an open question: Papallo Papallo:2017ddx () proved that DGB gravity fails to be strongly hyperbolic in any generalized harmonic gauge. Black hole solutions in DGB gravity have been constructed analytically for non-rotating configurations Mignemi:1992nt (); Mignemi:1993ce (); Yunes:2011we (); Sotiriou:2014pfa () and slowly rotating ones Pani:2011gy (); Ayzenberg:2014aka (); Maselli:2015tta ().

Deriving (monopole) scalar charges, also known as “sensitivities” Will:1989sk (), is crucial to understanding the emission of scalar radiation. One can extract such scalar charges by studying the asymptotic behavior of the scalar field at spatial infinity. So far, scalar charges have only been discussed in DGB gravity (but see Yagi:2015oca () for approximate calculations of stellar scalar charges in EdGB gravity). Scalar charges for black holes were extracted for nonrotating configurations Yagi:2011xp (); Sotiriou:2014pfa () and one can easily extract charges for slowly rotating black holes: see e.g. Ayzenberg:2014aka (). In Appendix A, we derive the scalar charge for black holes with arbitrary rotation. The final expression was shown in Yunes:2016jcc (), but the details of the derivation (together with the quadrupole scalar charge expression) are presented here for the first time. On the other hand, scalar charges for ordinary stars vanish Yagi:2011xp (); Yagi:2015oca (). This can be shown by integrating the scalar field equation

(10)

over the entire spacetime. The left-hand side gives the scalar charge, while the right-hand side generates a topological term and boundary terms after applying the generalized Gauss-Bonnet-Chern theorem Alty:1994xj (); Gilkey:2014wca (). Both of these vanish for simply connected, stationary spacetimes (i.e., for stars). This means that when a star collapses to a black hole, it suddenly acquires a scalar charge in a nontrivial manner: this is the opposite of what happens in scalar-tensor theories, where stars lose their scalar charge during collapse to a black hole. Such a nontrivial process was confirmed numerically under an Oppenheimer-Snyder collapse for nonrotating configurations Benkel:2016rlz () (see also Benkel:2016kcq () for related work). One can in fact show that scalar charges for stars vanish more in general for shift-symmetric Horndeski theories Barausse:2015wia (); Barausse:2017gip (); Lehebel:2017fag (), which include DGB gravity as a special class. This argument fails in the absence of shift symmetry Silva:2017uqg (); Doneva:2017duq (); Antoniou:2017acq (); Antoniou:2017hxj ().

2.3 Dynamical Chern-Simons Gravity

Dynamical Chern-Simons (dCS) gravity is an effective field theory model that modifies the Einstein-Hilbert action through a quadratic-curvature correction Alexander:2009tp (). Motivated from the Green-Schwarz anomaly-canceling mechanism in field theory and its higher-dimensional generalization in string theory Alexander:2004xd (), as well as from loop quantum gravity upon scalarization of the Barbero-Immirzi parameter Taveras:2008yf (); Calcagni:2009xz () and from effective field theories of inflation Weinberg:2008hq (), the dCS action is

(11)

where is the matter action, while and are a quadratic-curvature interaction term and the kinetic term for the scalar field respectively, namely

(12)

Here is a coupling constant and

(13)

is the Pontryagin density, with

(14)

the dual Riemann tensor. We assume here that the scalar field is dimensionless in geometrical units, which forces to have units of length squared. For simplicity, we assume we can Taylor expand to leading non-vanishing order. The constant term in the Taylor expansion does not contribute, because the Pontryagin density is a topological invariant, i.e. its integral over the manifold is the Pontryagin number, which is related to the manifold’s winding number. This means that the Pontryagin density can be written as the four-divergence of a four-current, which becomes a boundary term upon integration by parts and does not contribute to the field equations. The resulting theory with is then naturally shift-invariant, and in order to preserve this symmetry one typically sets the scalar field potential to zero.

The dCS field equations can be obtained by varying the action with respect to the metric and the scalar field. The resulting field equations, however, must be understood as effective, given that they arise from an effective action. In particular, this means that must be much smaller than during any process considered; if this were not the case, then higher-order curvature terms would have to be included in the action, which would then modify the field equations upon variation. As an effective field theory, it is trivial to prove that dCS gravity is well-posed Delsate:2014hba (). If one insists in treating the theory as exact, however, the appearance of third-order time derivatives in the field equations probably renders this exact version ill-posed Delsate:2014hba (). Spinning black hole solutions in dCS gravity have been found analytically for slowly rotating black holes Yunes:2009hc (); Konno:2009kg (); Yagi:2012ya () and for black holes in the extremal limit Stein:2014xba (); McNees:2015srl (), both sets of which differ from the Kerr metric family. Non-spinning black holes, and in fact, any static and spherically symmetric matter configuration have zero dCS modification because the Pontryagin density vanishes identically Jackiw:2003pm (); Yunes:2007ss (); Grumiller:2007rv (). Current constraints on come from Solar System and table top experiments, but they are extremely weak due to the weak fields in the Solar System alihaimoud-chen (); Yagi:2012ya ().

Just as in the case of EdGB gravity, the calculation of (monopole) scalar charges (or sensitivities) is also crucial here to understand the emission of scalar radiation. The asymptotic behavior of the scalar field at spatial infinity has revealed that black holes have a rather large scalar charge Yunes:2009hc (); Yagi:2012ya (), while the scalar charge of neutron stars is greatly suppressed Yagi:2011xp (); Yagi:2013mbt (); the derivation of the latter result follows closely the explanation in EdGB gravity presented in the previous section. As in the EdGB case, this means that stars have a tiny scalar field, but this grows upon gravitational collapse, until the charge asymptotes that of isolated black holes.

2.4 Lorentz-violating Gravity

Lorentz symmetry is one of the fundamental building blocks of GR. Some modified theories of gravity break this symmetry in order to provide a power-counting renormalizable completion of GR in the ultra-violet regime, as proposed by Hořava Horava:2009uw () (see Wang:2017brl () for a recent review). Lorentz symmetry breaking in the matter sector has been constrained very stringently from, e.g., particle physics experiments Kostelecky:2003fs (); Kostelecky:2008ts (); Mattingly:2005re (); Jacobson:2005bg (). A model-independent framework called the Standard Model Extension (SME) Colladay:1998fq (); Kostelecky:1998id (); Kostelecky:1999rh () was developed to map various observations to bounds on Lorentz symmetry breaking. This framework is efficient to probe such a breakage in the matter sector Kostelecky:2008ts () or in the sector where matter directly couples with gravity Kostelecky:2010ze (). On the other hand, Lorentz symmetry breaking in the gravity sector has not been well constrained, and different mechanisms exist to prevent the breakage in the matter sector to percolate into the gravity sector (see e.g. Liberati:2013xla (); Pospelov:2010mp ()). Bounds on gravitational Lorentz symmetry within the SME context can be found in Bailey:2006fd (); Bailey:2009me (); Bailey:2013oda ().

One example of Lorentz-violating gravity is Einstein-Æther (EA) theory Jacobson:2000xp (); Jacobson:2008aj (), which breaks the gravitational Lorentz symmetry by introducing a preferred time direction at each point in spacetime via a timelike unit vector . One can consider this theory as a low-energy effective theory of some unknown dynamics at high energy Blas:2009qj (). This theory is the most general vector-tensor model whose action only depends on a unit timelike vector field (Æther field) and its first derivative, and is quadratic in the latter. Such an action is given by , where

(15)

and

(16)

Here , is the bare gravitational constant in the theory, while are coupling constants. Solar System experiments essentially reduce the parameter space from four to two,  Jacobson:2008aj (). These remaining parameters have been constrained from binary pulsar observations Yagi:2013qpa (); Yagi:2013ava () and from the recent coincident GW and electromagnetic observation of merging neutron stars Monitor:2017mdv (). The binary pulsar bounds are derived by calculating the sensitivities of neutron stars, which are obtained by constructing a slowly moving neutron star solutions with respect to the Æther field (black hole sensitivities have not been calculated yet). The gravitational wave bounds come from the speed of the propagating degrees of freedom in EA theory, the two tensor modes, the two vector modes and the one scalar mode Jacobson:2004ts (); Foster:2005dk (); Foster:2006az (); Foster:2007gr (). Tensor, vector and scalar modes propagate at speeds , and respectively, which can be expressed in terms of the coupling constants, as summarized in Table 1. Here, we defined and .

Theory mode propagation speed
Einstein-Æther tensor
vector
scalar
khronometric tensor
scalar
Table 1: Propagation speed of tensor, vector and scalar modes in EA and khronometric theory. Vector modes are absent in the latter.

Another example of Lorentz-violating gravity that is related to EA theory is khronometric gravity Blas:2010hb (). This theory corresponds to the low-energy limit of Hořava gravity. The action is the same as the EA case, but now the Æther field is hypersurface orthogonal. In other words, this theory introduces a preferred time foliation of spacetime (or a global preferred time) and the Æther field is normal to the hypersurface of constant preferred time. This reduces the number of coupling constants from four to three: , and . Imposing the Solar System bounds essentially eliminates . One can further constrain the remaining parameters from binary pulsar Yagi:2013qpa (); Yagi:2013ava (), Big Bang nucleosynthesis Audren:2013dwa () and gravitational-wave Monitor:2017mdv (); Gumrukcuoglu:2017ijh () observations. Regarding the propagation degrees of freedom in khronometric gravity, there are two tensor modes and one scalar mode Blas:2011zd (), whose speeds are given by and , respectively. Their explicit forms are also summarized in Table 1.

2.5 Extra Dimensional Theories

String theory predicts that we live in a higher-dimensional spacetime, and extra dimensions are compactified. One simple example of such a compactification is the Kaluza-Klein compactification. Particle physics experiments place bounds on the size of the extra dimension: cm. Arkani-Hamed, Dimopoulos and Dvali (ADD) ArkaniHamed:1998rs (); ArkaniHamed:1998nn () proposed a braneworld model with a tensionless brane (on which we live) embedded in a flat and compact bulk spacetime. In the ADD model matter is localized on the brane, so that only gravitons can propagate through the bulk. The size of the extra dimensions can be relatively large in this model, since bounds on the gravity sector are much weaker than those on the matter sector. Furthermore, this model can naturally explain the hierarchy problem between the electroweak and Planck scale.

Another braneworld model was proposed by Randall and Sundrum. The first model (RS-I) Randall:1999ee () introduces one positive-tension and one negative-tension brane in a five-dimensional anti-de Sitter (AdS) bulk spacetime. With this model, one can choose the separation of two branes such that the fundamental five-dimensional Planck mass is TeV, while the four-dimensional Planck mass becomes GeV. In their second model (RS-II) Randall:1999vf (), they take the negative-tension brane to infinity, thus the model effectively includes only one brane. The remarkable feature of this model is that Newtonian gravity is reproduced in the low-energy limit Garriga:1999yh () even if the size of the extra dimension is relatively large. Table-top experiments place the bound Adelberger:2006dh ().

Since the bulk of the RS-II model is AdS, one can apply the AdS/CFT conjecture Maldacena:1997re (); Aharony:1999ti (), which states that the four-dimensional super Yang-Mills theory on the AdS boundary can be recovered from gravity in the spacetime. In particular, one can apply this conjecture to brane-localized black holes. In the CFT picture, due to the large number of CFT degrees of freedom Aharony:1999ti (), four-dimensional black holes on the brane evaporate significantly faster than in GR. In the AdS picture, this is seen as classical evaporation of a black hole of mass  emparan-conj (); tanaka-conj (), with evaporation rate given by Emparan:2002jp ()

(17)

Later, static brane-localized black hole solutions were found numerically Figueras:2011gd (); Abdolrahimi:2012qi (), which questions the validity of this conjecture. Nevertheless, we consider placing bounds on assuming that the conjecture is correct. This is because one can easily map this mass loss effect of black holes in the RS-II model to a similar effect due to (e.g.) phantom energy accretion in GR Babichev:2004yx (); Babichev:2005py (); Babichev:2014lda (). In terms of GWs, a change in the black hole mass is similar to a change in the gravitational constant  Yunes:2009bv (), which is predicted in many modified theories of gravity.

3 Inspiral Tests of Modified Gravity

3.1 Generic Tests

The plethora of modified gravity models that have been proposed and the extreme difficulty in constructing sufficiently accurate GW models for data analysis suggests that theories ought not to be treated on a case-by-case basis. Indeed, it took the gravity theory community approximately 50 years to obtain a sufficiently accurate (third post-Newtonian (PN) order) model of the GWs emitted in the inspiral of compact binaries in GR. Instead of carrying out similar calculations on a theory-by-theory basis, it is much more appealing to develop generic tests of Einstein’s theory given the available data. Indeed, a similar approach was successfully pursued when carrying out tests with Solar System observations, which led to the development of the parameterized PN framework of Will and Nordtvedt Nordtvedt:1968qs (); 1971ApJ…163..611W (); 1972ApJ…177..757W (); 1972ApJ…177..775N ().

The first attempt at such a generic test consisted of verifying the PN structure of the waveform phase Arun:2006yw (). The idea was to decompose the Fourier-domain waveform model into a frequency-dependent amplitude and a frequency-dependent phase, and to then rewrite the phase as222The terms and contain contributions that depend on , which the authors treat as constant in Arun:2006yw (). In their follow-up papers Arun:2006hn (); mishra (), they modified Eq. (18) by adding further terms of the form .

(18)

where are PN coefficients, which in GR are known functions of the parameters of the binary (to be more precise , the individual masses and for non-spinning black hole binaries), and is the orbital velocity, with the binary’s total mass. The proposal was then to treat all of these coefficients as independent and find the best-fit values by comparing the above template waveform with the data. One can then draw error regions of each coefficient in the - plane assuming GR is correct to check for consistency, namely to check if there is a region where all of error regions overlap. Later the authors only considered three out of eight coefficients, so that correlations among parameters could be reduced and one could carry out a stronger test by shrinking the error regions Arun:2006hn (); mishra (). This procedure resembles binary pulsar tests in the parameterized post-Keplerian formalism Damour:1991rd (); Stairs:2003eg ().

Although feasible in principle, the above test has a few limitations. First, it has the strong bias of assuming Nature follows the same exact functional structure of the PN approximation in GR, i.e. that the Fourier phase can be expressed as a series in integer powers of velocity, with the leading-order term starting at . Indeed, many examples of modified gravity effects and modified gravity theories exist which do not admit this structure; examples of this include dipole emission (), variability of the fundamental constants (), parity violation in eccentric binaries (), and massive gravitons in eccentric binaries (), to name a few. Second, the framework does not allow for tests of modified gravity theories that lead predominantly to amplitude modifications, without affecting the phase evolution much; examples of this include gravitational birefringence Alexander:2007kv (); Yunes:2010yf (); Yagi:2017zhb (). Third, the framework assumes that polynomials in velocity are a good basis to expand the Fourier phase during the entire inspiral, including right up to plunge and merger. Today, we know that this is not the case, with the series requiring arctangent corrections Husa:2015iqa (); Khan:2015jqa ().

An extension and generalization of this method that resolves all of the above problems is the parameterized post-Einsteinian (ppE) approach Yunes:2009ke (). In this framework, one extends the GR waveform model via

(19)

where and are the Fourier amplitude and Fourier phase in the most accurate GR model developed at that stage in time. The quantities are ppE constants that control the magnitude of deviations from GR, while are real numbers that determine the type of deviation that is being constrained. This parameterization resolves the issues listed above because (i) it allows for deformations from GR that enter at pre-Newtonian order, i.e. in negative powers of velocity relative to the leading-order PN term in GR, (ii) it allows for both amplitude and phase deformations, and (iii) it allows the deformation of the most sophisticated waveform model in GR, such that when the ppE parameters are zero, one recovers exactly the most accurate model. Different variants of this approximation have been and can be used, with for example different scalings in the ppE terms, working entirely in the time domain, or working with multiple amplitude and phase ppE parameters at the same time Cornish:2011ys (); Chatziioannou:2012rf (); Sampson:2013jpa (); Sampson:2013lpa ().

The idea of the ppE framework is then to use the above model and let the data decide (the posteriors for) the magnitude of the ppE parameters for whatever choice of modification one wishes to study. Once these posteriors have been constructed, one can then map them to posteriors on coupling parameters of specific theories. This mapping step is crucial because it is what allows us to draw inferences on different specific aspects of theoretical physics from the observations, as done recently with the first GW observations Yunes:2016jcc (). This step, of course, requires that one first calculate predictions for the GWs emitted by compact binaries in specific theories. These predictions have been obtained in many theories, as we detail in the next subsection. The ppE parameterization also relaxes the a priori assumption that GR is correct, allowing the data itself to select the theory that best supports it, thus minimizing what has been dubbed as theoretical bias Yunes:2009ke (); Vallisneri:2013rc (). Fortunately, all GW observations to date seem consistent with GR, implying that theoretical bias is not a problem right now.

Inferences drawn from constraints on ppE parameters can be classified into two classes, depending on the sector of the theory they constrain: generation and propagation. The generation sector deals with the way the particular theory sources GWs, and any additional degree of freedom, and how these evolve in time and back-react on the evolution of the binary system. The propagation sector deals with how the waves travel away from the binary system to Earth, according to the wave’s dispersion relation. One can, of course, devise more generic parameterizations that account for more exotic phenomena during the generation and propagation of GWs. For example, one could imagine that an additional massive degree of freedom is present in the theory, leading to modifications to the dynamics only once the binary radiates at the Compton wavelength of the new degree of freedom. This would lead to a sharp turn-on and turn-off of GR modifications that are not captured by the parameterization above. However, data analysis investigations have revealed that even if Nature were that cruel, the ppE parameterization presented above would still be able to detect an anomaly in the signal, although in a suboptimal way Sampson:2013lpa (). Similarly, one could imagine that in a yet-to-be-discovered modified theory GWs propagate differently depending on the direction in which they travel Tso:2016mvv (). Such a direction-dependent dispersion relation would not fit perfectly within the ppE parameterization, but a signal of that type can be recovered with the simplest ppE model written above, although sub-optimally.

A recent idea has been put forth to draw generation-type inferences on the nature of the black holes that generate GWs in a binary. The idea is to allow for the quadrupole moment of each black hole in a binary to be a new parameter in the waveform model, and then to check whether the measured parameter agrees with the GR expectation per the no-hair theorems Krishnendu:2017shb (). A generic quadrupole in the waveform will introduce a second PN order correction to the GW phase, which thus maps in a one-to-one fashion to the ppE parameterization described above. This correction, however, will be degenerate with the spins of the black holes, so unless these can be separately extracted, it would be difficult to implement this test. The independent extraction of the spins can be achieved if the signal and the waveform model include spin-orbit and spin-spin precession effects, which induce amplitude modulations that can in principle break this degeneracy Chatziioannou:2014coa (); Chatziioannou:2014bma (); Chatziioannou:2016ezg (); Chatziioannou:2017tdw ().

More in general, the binary dynamics of hypothetical black hole alternatives in binary systems is driven by their so-called tidal Love numbers. Tidal Love numbers encode the deformability of a self-gravitating object immersed in a tidal environment, and they depend both on the object’s internal structure and on the dynamics of the gravitational field. In classical general relativity, the tidal Love numbers of black holes are exactly zero. Recent work computed the tidal Love numbers of various exotic compact objects, including boson stars, gravastars and wormholes, as well as black holes in various theories of modified gravity (Einstein-Maxwell, Brans-Dicke and Chern-Simons gravity) Cardoso:2017cfl (). In general, these calculations showed that the tidal Love numbers of exotic objects depend logarithmically on the location of the putative “surface” replacing the black hole’s horizon. LIGO-like detectors, third-generation Earth-based detectors and LISA can impose interesting constraints on the tidal Love numbers of boson stars, while a LISA-like detector could probe even extremely compact objects Maselli:2017cmm (), as long as systematic errors in general relativity are under control.

Theory
scalar-tensor Jacobson:1999vr (); Horbatsch:2011ye (); Yunes:2016jcc ()
EdGB, DGB Yagi:2011xp ()
dCS Yagi:2012vf ()
EA Hansen:2014ewa ()
khronometric Hansen:2014ewa ()
extra dimension yagi:brane ()
varying  Yunes:2009bv ()
mod. disp. rel. Mirshekari:2011yq ()
Table 2: Mapping of ppE parameters to those in each theory for a black hole binary. In scalar-tensor theories, black holes acquires a scalar charge for a cosmologically evolving scalar field Jacobson:1999vr (); Horbatsch:2011ye (). Such a scalar charge is proportional to . is related to the black hole scalar charge in DGB in Eq. (37) as . The dimensionless coupling constant in quadratic-curvature theories is defined by . Propagation speeds in Lorentz-violating theories are summarized in Table 1. can be calculated from Eq. (17). , where is the Planck constant. The distance is defined in Eq. (22).

3.2 Mapping to Specific Theories

The mapping between the ppE parameters and the coupling constants in various modified theories of gravity is summarized in Table 2. Since matched filtering data analysis is more sensitive to deviations in the waveform phase than the amplitude, we only show the mapping for the phase. The mapping in the table corresponds to GWs from black hole binaries with component masses and dimensionless spins (). We also introduce the total mass , symmetric mass ratio and chirp mass , as well as two (symmetric and antisymmetric) spin parameters and .

Of all the theories considered in the inspiral of compact objects, scalar-tensor theories are by far the most well-studied. Pioneering studies of compact binaries in scalar-tensor theory were carried out by Eardley 1975ApJ…196L..59E (), Will and Zaglauer Will:1977wq (); Will:1989sk () and Damour and Esposito-Farèse Damour:1992we (); Damour:1996ke (); Damour:1998jk (), among others. More recently, the post-Minkowskian technique of direct integration of the relaxed Einstein equations (DIRE) Will:1996zj (); Pati:2000vt (); Pati:2002ux () was used to compute the equations of motion for a system of compact objects at 2.5PN order, as well as the gravitational waveform and energy flux at 2PN (relative) order for binaries on generic orbits Mirshekari:2013vb (); Lang:2013fna (); Lang:2014osa (). Sennett et al. Sennett:2016klh () specialized the results of Mirshekari:2013vb (); Lang:2013fna (); Lang:2014osa () to non-spinning binary systems on quasi-circular orbits in scalar-tensor gravity at 2PN relative order333By imposing the stringent constraints set by current astrophysical observations (cf. Table II of Sennett:2016klh ()), they find that dipolar radiation is subdominant to quadrupolar radiation for most prospective GW sources: in the absence of spontaneous scalarization, the dipole term can dominate only at frequencies Hz in binary neutron star or neutron-star/stellar-mass-black-hole systems, and at frequencies Hz in neutron-star/intermediate-mass-black-hole systems. Therefore, ground- and space-based GW detectors would only observe binary systems whose inspiral is driven by the next-to-leading order flux.. Such waveforms introduce PN corrections to the mapping presented in Table 2. Julié and Deruelle Julie:2017pkb (); Julie:2017ucp () use these higher order PN results to begin to extend the effective-one-body (EOB) formalism of Buonanno and Damour Buonanno:1998gg () to scalar-tensor gravity. Such resummed waveform models cannot be analytically mapped to the ppE waveforms directly.

Theory
massive gravity Will:1997bb (); Rubakov:2008nh (); Hinterbichler:2011tt (); deRham:2014zqa () 0
multifractional spacetime Calcagni:2009kc (); Calcagni:2011kn (); Calcagni:2011sz (); Calcagni:2016zqv () (timelike spacetime) 2–3
(spacelike spacetime)
double Special Relativity AmelinoCamelia:2000ge (); Magueijo:2001cr (); AmelinoCamelia:2002wr (); AmelinoCamelia:2010pd () 3
extra dimension Sefiedgar:2010we () 4
SME Kostelecky:2016kfm () (even )
(odd )
Hořava Horava:2008ih (); Horava:2009uw (); Vacaru:2010rd (); Blas:2011zd () 4
Table 3: Mapping between modified dispersion relation parameters for the graviton in Eq. (20) and the parameters of each theory. The meaning of the parameters is as follows. : the graviton mass; : the characteristic length scale above which spacetime is discrete; : the characteristic observer-independent length scale; : the square of the Planck length in extra dimensional theories; and : parameters controlling the Lorentz-violation operators in SME in the rotation-invariant limit; : a parameter related to the bare gravitational constant; : a parameter related to the deformation in the “detailed balance” conditions in Hořava gravity.

All the mappings in Table 2 (except for the last one) originate from non-GR effects created at the level of generation of GWs, while such waves in general acquire modifications also at the level of their propagation. The dispersion relation of the graviton in non-GR theories can be expressed in terms of two parameters and as Mirshekari:2011yq ()

(20)

We present the mapping between these two parameters and coupling constants in each modified theory of gravity in Table 3. From the above dispersion relation, one finds the group velocity of the graviton

(21)

The ppE parameters and with the modified dispersion relation of the graviton in Eq. (20) are given by the last line in Table 2. Here, the distance parameter is defined by Mirshekari:2011yq ()

(22)

where is the local Hubble parameter, represents the redshift, and and are the energy density of dark matter and dark energy, respectively.

3.3 Current and Future Tests

3.3.1 Current Tests

The LIGO/Virgo Collaboration (LVC) used the observed data to perform various model-independent tests of gravity in the extreme field regime TheLIGOScientific:2016src (). The first test that they performed was to estimate the residual signal-to-noise ratio between the GW150914 data and the GR waveform template. They concluded that GR violations in GW150914 that cannot be absorbed into a redefinition of binary parameters have been constrained to less than 4%. The second test was to check the consistency in the measurement of the final black hole mass and spin using the inspiral and post-inspiral data. The third test was to look for non-tensorial polarization modes. The presence of such polarization modes was inconclusive for the first three events due to the near-alignment of the two aLIGO detectors. Since Virgo also detected signals for the fourth event GW170814 Abbott:2017oio (), the LVC could carry out a meaningful test of GR with GW polarizations for the first time. They concluded that the data favors purely tensor polarizations over purely scalar (vector) polarizations with a Bayes factor larger than 1000 (200), respectively. The fourth test was to constrain deviations in the waveform phase away from GR at each PN order.

Figure 1: [Adapted from Yunes:2016jcc ().] 90%-confidence constraints on the ppE parameter entering at th PN order. The green crosses represent the bounds reported in TheLIGOScientific:2016pea (); TheLIGOScientific:2016src () using a Bayesian analysis of event GW150914. The red (blue) dots and line represent bounds from GW150914 (GW151226) estimated using a Fisher analysis. The dashed black line and the cyan star corresponds to bounds from binary pulsar observations Yunes:2010qb () and Solar System experiments Sampson:2013wia (). Observe that the GW bounds are stronger than binary pulsar bounds on entering at a positive PN order.

The last test mentioned in the previous paragraph can easily be mapped to a test with ppE waveforms. Green crosses in Fig. 1 show the GW150914 bound on at each PN order that is obtained by mapping the results in TheLIGOScientific:2016src () to ppE parameters. These results constrain non-GR corrections arising from GW generation mechanisms, as the non-GR corrections are introduced only in the inspiral part of the waveform. For the merger-ringdown phase, one needs to carry out more merger simulations in non-GR theories to construct valid parameterized waveforms. Observe that bounds on with GW150914 are much stronger than those from the orbital decay rate of the double binary pulsar J0737-3039 Yunes:2010qb () for positive PN corrections. This is because such positive PN corrections become larger in the stronger field regime, precisely where GW observations are most probing. The cyan star indicates the bound on at 1PN order obtained from Solar System experiments. Although such a bound is stronger than the GW bound, these two bounds have a different meaning: the former probes weak, non-dynamical gravity, while the latter probes the strong and dynamical regime.

Reference Yunes:2016jcc () extended this result by deriving GW bounds also for negative PN corrections444The LVC derived bounds on the PN term with GW170814 Abbott:2017oio ().. Instead of carrying out a full Bayesian analysis using the actual data, the authors carried out a simpler calculation using a Fisher analysis. They injected a GR waveform whose parameters are consistent with those reported in TheLIGOScientific:2016pea (), and recovered them with ppE templates. Bounds on for both positive and negative PN corrections with this method for GW150914 are shown by the red curve and circles. Observe first that the bounds on the positive PN corrections agree nicely with those obtained in TheLIGOScientific:2016src (). Such an agreement demonstrates the validity of using a Fisher analysis and a GR waveform for injection. Observe also that bounds for negative PN corrections are much weaker than the binary pulsar bounds. A stronger bound on does not necessarily mean a stronger bound on a specific modified theory of gravity. This is because depends not only on theoretical coupling constants, but also on binary parameters like masses and spins. Moreover, there are theories like DGB gravity where non-GR corrections from neutron star binaries are highly suppressed, and thus one needs to use observations from binaries with at least one black hole to constrain the theories more efficiently. The blue dashed curve and triangles present bounds from GW151226. The bounds are stronger than GW150914, especially for negative PN corrections. This is because the masses of GW151226 are smaller than those of GW150914, which means that the velocity of binary components are smaller at a fixed frequency, making the corrections larger at negative PN orders.

Theory GR Pillar PN Repr. Parameters GW150914 Other Bounds
EdGB, DGB SEP [km]  Amendola:2007ni (), 2 Kanti:1995vq (); kent-LMXB (); Pani:2009wy ()
scalar-tensor [1/sec]  Horbatsch:2011ye ()
dCS SEP, PI [km]  alihaimoud-chen (); Yagi:2012ya ()
Einstein-Æther SEP, LI  Yagi:2013qpa (); Yagi:2013ava ()
khronometric  Yagi:2013qpa (); Yagi:2013ava ()
Extra Dimensions 4D [m] 10–10 Johannsen:2008tm (); johannsen2 (); Adelberger:2006dh (); psaltis-RS (); gnedin ()
Time-Varying SEP [10/yr] 0.1–1 Manchester:2015mda (); 2011Icar..211..401K (); Hofmann (); Copi:2003xd (); Bambi:2005fi ()
Massive graviton [eV]  TheLIGOScientific:2016src ()  talmadge (); sutton (); goldhaber (); Hare:1973px (); Brito:2013wya ()
Multifractional LI [eV] (time)
[eV] (space)  Kiyota:2015dla ()
double Special Rel. LI
 Kiyota:2015dla ()
Extra Dimensions 4D  Kiyota:2015dla ()
Stand. Model Ext. LI  Kostelecky:2010ze (); Kiyota:2015dla ()
[cm]  Kostelecky:2016kfm ()  Kostelecky:2010ze (); Kiyota:2015dla ()
[cm]
[cm]  Kostelecky:2010ze (); Kiyota:2015dla ()
[cm]
Hořava-Lifshitz LI [1/eV]
Einstein-Æther LI  Blas:2016qmn ()  Yagi:2013qpa (); Yagi:2013ava ()
Table 4: Various bounds on example theories that violate certain fundamental pillars in GR. The meaning of each column is as follows. 1st: names of modified theories of gravity; 2nd: GR fundamental pillars that each theory break; 3rd: the leading PN order in the gravitational waveform at which the non-GR effect enters; 4th: representative parameters in each theory; 5th: bounds on each non-GR parameter from GW150914 derived mostly in Yunes:2016jcc (); 6th: other bounds on each theory. Theories in the top (bottom) half of the table modifies the waveform at the level of generation (propagation). One finds similar bounds from GW151226 Yunes:2016jcc (). This table is taken and edited from Yunes:2016jcc ().

Although Fig. 1 is useful as the bounds presented are model independent, it is unclear what fundamental pillars of GR we are testing. Reference Yunes:2016jcc () addressed such a question by mapping the bound on to that on each example non-GR theory that violates certain fundamental aspects of GR, such as the strong equivalence principle (SEP), Lorentz invariance (LI), parity invariance (PI), 4 dimensional spacetime (4D) and a massless graviton . The first half of Table 4 summarizes the results by presenting GW150914 bounds on each theory that are mapped from Fig. 1, together with other bounds from (e.g.) Solar System experiments or binary pulsar observations. We also list the fundamental pillars of GR violated by each theory, and the PN order at which dominant corrections appear. Observe first that we do not show any GW150914 bounds on the first three theories. This is because the bounds are so weak that they violate the small-coupling approximation used to derive corrections to the waveform. On the other hand, one can place meaningful constraints on the other four theories considered that break SEP, LI or the assumption of a four-dimensional spacetime. Although these GW bounds are much weaker than other existing bounds, they are the first bounds obtained in the extreme gravity regime. The potential of GW observations in probing GR is currently limited by the lack of our knowledge of non-GR effects in the merger-ringdown regime.

One can carry out a similar test for corrections generated at the level of GW propagation. For such a case, one can include corrections not only in the inspiral phase, but also in the merger-ringdown part of the waveform. GW150914 bounds on the mass of the graviton, which are three times stronger than Solar System bounds talmadge (), were obtained in this way TheLIGOScientific:2016src (). The second half of Table 4 summarizes bounds on each example theory that modifies the dispersion relation of the graviton. Unlike bounds on GW generation, those on GW propagation are complementary to cosmic ray bounds obtained from the absence of the gravitational Cherenkov radiation Kiyota:2015dla (). In other words, some of the parameter space in each theory has been constrained for the first time using GW150914. Some theories (like EA theory) are difficult to constrain from deviations in the waveform phase due to degeneracies with other parameters, like the time of coalescence. Thus the bound on such a theory in the last line of Table 4 is obtained from bounds on the propagation speed of GWs due to the arrival time difference between the Hanford and Livingston detectors Blas:2016qmn () (see also Cornish:2017jml ()).

Observed GW events were used to constrain other modified theories of gravity, such as non-commutative spacetime Kobakhidze:2016cqh (). These events were also used to constrain graviton oscillations in bigravity Max:2017flc (), a phenomenon similar to neutrino oscillations, first proposed in DeFelice:2013nba ().

Very recently, the LVC detected yet another GW signal, that is consistent with the source being a binary neutron star merger (GW170817) TheLIGOScientific:2017qsa (). Unlike the binary black hole events, where confirmed electromagnetic wave counterparts were absent, GW170817-associated counterpart signals were detected by gamma-rays, X-rays, ultraviolet, optical, infrared and radio waves GBM:2017lvd (). This historic observation marks the dawn of the era of the multi-messenger astronomy. Thanks to the simultaneous detection of GWs and electromagnetic waves, the LVC carried out new tests of GR Monitor:2017mdv () (see also Boran:2017rdn (); Shoemaker:2017nqv (); Wang:2017rpx (); Wei:2017nyl ()). For example, they measured the propagation speed of GWs from the arrival time difference between GWs and gamma-rays. Conservatively assuming that photons were emitted within 10 seconds compared to the graviton’s emission time, they obtained the bound  Monitor:2017mdv (). As pointed out e.g. in Lombriser:2015sxa (); Lombriser:2016yzn (), this rules out many models in modified theories of gravity that aim to explain the current accelerating expansion of the Universe, including the quartic, quintic and covariant Galileons Sakstein:2017xjx (); Creminelli:2017sry (); Ezquiaga:2017ekz (); Baker:2017hug (); Arai:2017hxj (). The observations constrain the parameter in Einstein-Æther theory and in khronometric gravity to be of order  Hansen:2014ewa (); Yunes:2016jcc (); Baker:2017hug (), improving over previous bounds by roughly 13 orders of magnitude. Such a bound on the propagation speed of the graviton can be used to probe gravitational Lorentz violation through the SME Kostelecky:2016kfm (). The LVC placed bounds on gravitational SME parameters which are generally a few to 10 times more stringent than previous bounds Monitor:2017mdv (). GW170817 also probes violations of the equivalence principle. This was done by testing whether gravitons and photons feel the same gravitational potential as they propagate. The difference between the graviton’s and photon’s parameterized PN parameter ( and respectively) was bounded as  Monitor:2017mdv () via Shapiro delay measurements. GW170817 has been used to set stringent constraints on the Vainshtein mechanism Crisostomi:2017lbg (); Langlois:2017dyl (); Dima:2017pwp () and on several modifications of GR, including theories with extra dimensions Visinelli:2017bny (), Hor̂ava gravity Gumrukcuoglu:2017ijh () models Nojiri:2017hai (); Lee:2017dox () and massive gravity Heisenberg:2017qka ().

3.3.2 Future Tests

GW tests will benefit in the future from three types of accomplishments: (i) multiple detections, (ii) low-frequency detections, and (iii) multi-wavelength observations. It is easy to understand how multiple detections can yield improved tests: one can either stack the events to enhance the power in the signal, or simply combine the events by multiplying posteriors together to enhance constraints. This will strengthen inferences on both the generation and propagation of waves, roughly by a factor if is the number of comparable signal-to-noise ratio events. In reality, the enhancement factor will be dominated by the loudest events (see e.g. Berti:2011jz ()). The second accomplishment refers specifically to observations with space-borne detectors, which will be sensitive to waves in the milli-Hz range. These observations are unique because they will allow for very high signal-to-noise ratios (in the hundreds to thousands) and very large distances (roughly Gpc and beyond), and they can probe systems at much larger separations, or much lower orbital frequency. Such observations will strengthen inferences both on the generation and propagation sectors, although propagation bounds will benefit the most, due to the long baseline of the measurements. The third accomplishment refers to observations that could be done first by space-borne detectors at deci-Hz frequencies, when the binary system is widely separated, and then again by ground-based detectors at hecto-Hz frequencies, when the same binary merges. This will allow for precise tests that metaphorically tie the theory at both ends: during the early inspiral and during the merger simultaneously.

Figure 2: [From Chamberlain:2017fjl ().] Projected constraints on modified gravity effects as a function of ppE PN order at which they first enter, for a variety of space-based (left) and ground-based (right) detectors and a variety of systems. Anything above the regions is projected to be ruled out. The shaded regions are bounded by the highest and lowest constraints that can be placed at a given PN order for all instruments studied. For comparison, we also include the constraints that have already been placed by aLIGO with the GW150914 detection Abbott:2016blz (); Abbott:2016nmj () (thin cyan line), as well as constraints that can be placed with binary pulsars Yunes:2010qb () (dashed black line). Observe that the magnitudes of the projected constraints with space-based and ground-based instruments are comparable at positive PN orders, with space-based constraints being better by roughly orders of magnitude at negative PN order.

Given these expected advancements, one may wonder how much more stringent future constraints will become in the future Barausse:2016eii (); Chamberlain:2017fjl (); Samajdar:2017mka (). Fig. 2 shows projected constraints on the parameter as a function of the PN order at which they enter – i.e., a term of PN order is associated with a ppE correction proportional to – for a variety of single GW observations from the inspiral of compact binaries. The shaded regions correspond to variations in the constraints due to using different future instruments, including aLIGO at design sensitivity, A+, Voyager, Cosmic Explorer, the Einstein Telescope, and different incarnations of LISA. Observe that future constraints will be many orders of magnitude more stringent than current constraints (represented here with the aLIGO observation of GW150914 in cyan). These constraints would be enhanced by a factor of roughly given observations, where one should keep in mind that ground-based and space-based detectors are not expected to see the same number of events; although this is strongly dependent on the uncertain event rate, one expects to see roughly sources with ground-based instruments, and roughly sources with space-based instruments (see e.g. Berti:2016lat ()).

Figure 3: [From Chamberlain:2017fjl ().] Projected fractional improvement of constraints on modified gravity effects as a function of ppE PN order at which they first enter for a variety of detectors (relative to aLIGO at design sensitivity). The shaded regions are the same as in Fig. 2. Observe that 3G detectors improve constraints by more than an order of magnitude even for single detectors.

A related question is how the strength of the constraint changes with different third-generation (3G) detector configuration. Fig. 3 shows the fractional improvement of projected constraints as a function of ppE PN order at which modified gravity effects first enter Chamberlain:2017fjl () for a variety of detectors (relative to aLIGO at design sensitivity). Although minor upgrades, like A+ and Voyager, will only lead to modest improvements in constraints, 3G detectors can achieve improvements that are better than an order of magnitude. The fractional improvement dramatically increases at negative PN order in the ET case, simply because of this detector’s greatly improved sensitivity at low frequencies Chamberlain:2017fjl ().

Figure 4: [From Chamberlain:2017fjl ().] Projected constraints on dipole emission (left) and the mass of the graviton (right) for a variety of sources and detectors. Anything above the points is projected to be ruled out. The horizontal dashed lines correspond to current constraints. Observe that the magnitude of the projected constraints with space-based and ground-based instruments are always a great improvement over what aLIGO at design sensitivity will be able to achieve.

Given these more stringent constraints, the natural question to ask is: what new physics will be probed in the future? Multi-wavelength observations with space- and ground-based instruments will allow for constraints on violations of the strong equivalence principle that are 8 orders of magnitude more stringent than all current bounds Barausse:2016eii (). Single observations with future instruments will allow for constraints on the size of a large extra-dimension (in Randall-Sundrum type models) that are 5 orders of magnitude more stringent than current bounds with aLIGO Chamberlain:2017fjl (). Similarly, constraints on the mass of the graviton from propagation effects in the dispersion relation will be about 5 orders of magnitude better than current bounds Chamberlain:2017fjl (). These constraints would begin to approach the natural value of the mass of the graviton in eV that one would expect if such a mass is somehow connected to a solution to the dark-energy problem. Fig. 4 shows how constraints on dipole emission and the mass of the graviton improve for a variety of sources and detectors Chamberlain:2017fjl ().

4 Merger Tests of Modified Gravity

There are very few merger simulations in non-GR theories. This is because simulations are possible only for theories where a gauge has been found in which the theory is well-posed. For example, dCS gravity is likely to be ill-posed when treated as an exact theory Delsate:2014hba (); Crisostomi:2017ugk (), although it becomes well-posed if one treats it as an effective field theory and solves the field equations order by order Delsate:2014hba (), as done in Okounkova:2017yby (). The well-posedness of Lovelock gravity and of a certain subclass of Horndeski theories has been studied in Papallo:2017qvl (), including also DGB gravity Papallo:2017ddx (). The latter theory was found to be not strongly hyperbolic in a specific class of (generalized harmonic) gauges. Therefore the well-posedness of EdGB and DGB gravity in general is still an open question. Recently Cayuso et alCayuso:2017iqc () proposed a new framework to make a truncated theory well-behaved following the Israel-Stuart formalism for relativistic viscous hydrodynamics, and they applied this framework to toy models for linearized non-commutative geometry and dCS gravity. The application of this framework to the full theories and to other non-GR theories is an active research area.

In contrast to the theories discussed above, proving well-posedness is almost trivial for the Bergmann-Wagoner scalar-tensor theories555Since theories are equivalent to scalar-tensor gravity, they are also well-posed LanahanTremblay:2007sg (); Paschalidis:2011ww (). discussed in section 2.1.1 Salgado:2008xh (). We will first review analytical work on the dynamics of binary systems in these theories, and then summarize some numerical simulations of compact binary systems that have been carried out in these theories.

4.1 Compact Binaries in Scalar-Tensor Theories: Analytical Results

As stated earlier, scalar-tensor gravity is the most studied modified theory of gravity. It was studied in depth already in the 1950s Brans:1961sx (); Fierz:1956zz (); Jordan:1959eg (), and most variants of Bergmann-Wagoner theories have been constrained to be extremely close to GR, at least in the weak-field limit. These theories are well-posed, so compact binary systems can be evolved throughout merger. The question is: are there variants of scalar-tensor gravity which are still compatible with weak-field and cosmological observations, but at the same time make predictions in the strong-gravity regime that are sufficiently different from GR to be testable with GW interferometers?

Weak-field observations imply that in Eq. (3), so that deviations from GR are generally small. However there are cases where scalar-tensor gravity may lead to observable differences from GR in the strong-field regime targeted by GW detectors. Three such smoking guns of scalar-tensor gravity are:

(i) The emission of dipolar gravitational radiation from compact binary systems. Dipolar gravitational radiation is “pre-Newtonian,” i.e. it occurs at lower PN order than quadrupole radiation, and it does not exist in GR 1975ApJ…196L..59E (); Will:1989sk (). So far, LIGO/Virgo binary black hole detections outnumber the detections of binary systems containing neutron stars. Unfortunately, the phenomenology of scalar-tensor theory in vacuum spacetimes (such as binary black hole spacetimes) is less interesting than in the presence of matter. This is because, when the matter action can be neglected, the Einstein-frame formulation of the theory is equivalent to GR minimally coupled to a scalar field. Isolated black holes in Bergmann-Wagoner theories666Some classes of Horndeski theory (for example, those that can be shown to be equivalent to Einstein-dilaton-Gauss-Bonnet gravity through integration by parts) are such that these no-hair theorem can be circumvented Sotiriou:2013qea (); Sotiriou:2014pfa (); Maselli:2015yva (); Silva:2017uqg (), so that stationary black hole solutions can be different from GR. satisfy the same no-hair theorem as in GR, and thus the stationary black hole solutions in the two theories coincide777There are some proposal to circumvent these no-hair theorems involving time-dependent scalar fields Herdeiro:2014goa (); Herdeiro:2015gia (); Herdeiro:2015waa (). Recent evidence shows that the resulting solutions are unstable Ganchev:2017uuo (), but the instability is astrophysically irrelevant in some regions of the parameter space Degollado:2018ypf (). Heusler:1995qj (); Sotiriou:2011dz (). Dipolar radiation is proportional to the difference in sensitivity between the two binary members [see Eq. (41) below], so it vanishes identically for binary black hole systems. In fact, as we explain in Appendix B, dynamical (vacuum) black hole binary spacetimes satisfy what we could call a generalized no-hair theorem: the dynamics of a black hole binary system in Bergmann-Wagoner theory with vanishing potential in asymptotically flat spacetimes are the same as in GR up to at least PN order for generic mass ratios Damour:1992we (); Mirshekari:2013vb (), and at any PN order in the extreme mass-ratio limit Yunes:2011aa ().

(ii) Spontaneous scalarization. Sizable strong-field deviations from GR can occur in the presence of non-perturbative neutron star (or black hole) solutions for which the scalar field amplitude is finite even when . This spontaneous scalarization phenomenon Damour:1993hw (); Damour:1996ke () could significantly affect the mass and radius of a neutron star, and therefore the orbital motion of a compact binary system, even far from coalescence. This mechanism relies on the presence of matter, so (once again) it does not apply to black holes; however, it has recently been pointed out that there are some scalar-tensor theories where black hole solutions can scalarize Doneva:2017bvd (); Silva:2017uqg (); Antoniou:2017acq (); Antoniou:2017hxj ().

(iii) Superradiance in rotating black hole spacetimes. Another class of non-perturbative mechanisms involves rotating black holes and the phenomenon known as superradiance (see Brito:2015oca () for a review). The coupling of massive scalar fields to matter in orbit around rotating black holes leads to a surprising effect: because of superradiance, matter could in principle hover into “floating orbits” for which the net gravitational energy loss at infinity is entirely provided by the black hole’s rotational energy Cardoso:2011xi (). This phenomenon could in principle occur in nature, but it is possible that radiation reaction Zimmerman:2015hua () and finite-size effects Fujita:2016yav () will destabilize floating orbits. Tachyonic instabilities similar to spontaneous scalarization can also be produced when “ordinary” rotating black holes in GR are surrounded by matter, as long as the trace of the stress-energy tensor changes sign Cardoso:2013opa (); Cardoso:2013fwa ().

Two complementary approaches are used to model these effects in compact binary systems: analytical calculations (usually PN expansions) and numerical relativity simulations. Below we will review numerical studies of compact binaries in scalar-tensor gravity and other modified theories of gravity.

4.2 Compact Binaries in Scalar-Tensor Theories: Numerical Simulations

Even within GR, obtaining numerically stable and accurate time evolutions of the Einstein equations took almost 50 years. A stable evolution of binary black holes and neutron stars in numerical relativity requires an understanding of many complex issues, such as the well-posedness of the evolution system, the construction of initial data and gauge conditions Baumgarte:2010ndz (). These same questions naturally arise also in modified theories of gravity, and at present they remain unanswered for most of these theories. Scalar-tensor theories are a notable exception, because they can be formulated in close analogy to GR. As discussed in Section 2.1, the action of scalar-tensor theories in the Einstein frame is the same as the Einstein-Hilbert action, except for a minimal coupling with the scalar field in the gravitational sector. A non-minimal coupling with the scalar field only appears in the matter sector. The resulting field equations are similar to the field equations of GR, and the evolution of the scalar field is dictated by a hyperbolic wave equation. Therefore scalar-tensor theories can be evolved using relatively minor generalizations of the numerical codes developed for GR. As shown by Salgado et al. Salgado:2005hx (); Salgado:2008xh (), a strongly hyperbolic formulation can be obtained also in the physical (Jordan) frame. However, the Einstein frame is exceptionally convenient for applications to binary black hole spacetimes, because in vacuum the evolution equations are independent of the coupling function . In this sense, a single numerical evolution represents a whole class of theories characterized by different functional forms of for a given potential . Different choices of the function result in different physical predictions, but all of these predictions can be calculated by post-processing data from a single numerical simulation. This would not be possible in the Jordan frame, where the coupling function appears explicitly in the field equations.

Black hole binaries.

Because of the “dynamical no-hair theorem” summarized in Appendix  B, interesting black hole binary dynamics that are significantly different from GR requires somewhat contrived scenarios. The “dynamical no-hair theorem” of Appendix B relies on the following assumptions: (1) vacuum spacetime, (2) the scalar-tensor action is truncated at second order in a derivative expansion, (3) the potential vanishes, and (4) the metric is asymptotically flat, with an asymptotically constant scalar field.

Deviations from GR in black hole binaries can occur if we violate any of these assumptions. What happens when we violate hypothesis (1) will be discussed in the next section. Violating hypothesis (2) by introducing higher-order derivatives in the action would lead to substantially more complicated equations, whose well-posedness is presently unclear (see e.g. Berti:2013gfa ()).

Healy et al. Healy:2011ef () introduced nontrivial dynamics by placing the black holes inside a scalar field “bubble” which in some cases includes a non-vanishing scalar field potential, thus violating hypothesis (3). As the bubble collapses, the black holes accrete the scalar field and grow in mass. This mass growth affects the binary dynamics and the emitted gravitational radiation, and the inclusion of a potential term introduces longer lived dynamics in the scalar field mode. Their work supports the view that an evolving scalar field is required to introduce interesting dynamics in black hole binaries. However, for the effects to be observable the merging black holes must accrete enough scalar field to appreciably change their masses and modify the binary evolution.

Horbatsch & Burgess Horbatsch:2011ye () suggested violating hypothesis (4) instead: if the scalar field is time-dependent at the boundary, the black holes in a binary could retain scalar hair Jacobson:1999vr () and emit dipole radiation, as long as their masses are not exactly equal. Numerical evolutions implementing this idea were carried out in Berti:2013gfa () by introducing non-asymptotically flat or constant boundary conditions. The main motivation for relaxing this assumption comes from cosmology: indeed, inhomogeneous scalar fields have been considered as an alternative to dark matter Sahni:1999qe (); Hu:2000ke () and as models of supermassive boson stars Macedo:2013qea (). For scalar-field profiles that vary on a length scale much larger than the black hole binary orbital separation, the scalar-field gradient can be considered approximately constant. Ref. Berti:2013gfa () studied the quasi-circular inspiral of a non-spinning black hole binary (with mass ratio ) in a scalar-field gradient perpendicular to the orbital angular momentum vector. The lowest multipoles of the radiation are effectively indistinguishable from their GR counterparts. However, a non-vanishing scalar field gradient can lead to (mostly dipolar) emission of scalar radiation, which is not present in GR at twice the orbital frequency. At first glance this may appear surprising, but a simple calculation reveals that this feature is a consequence of the interaction of the orbital motion with a background field with an azimuthal component : cf. the discussion around Eqs. (36)-(38) of Berti:2013gfa ().

In summary, the simulations of Berti:2013gfa () and Healy:2011ef () demonstrate that non-asymptotically flat boundary conditions provide a mechanism to generate scalar radiation in black hole inspirals in scalar-tensor theories of gravity, at least in principle. Unfortunately, this scalar radiation will be unobservable in the near future for cosmologically realistic values of the scalar-field gradients.

Neutron star binaries.

The dynamics of scalar-tensor theories of gravity in the presence of matter sources can be very different from GR. The theory can violate the strong equivalence principle, so that self-gravitating objects follow trajectories that depend on their internal composition/structure: this is the well-known “Nordtvedt effect” Nordtvedt:1968qr (); Roll:1964rd (); 1975ApJ…196L..59E ().

For Bergmann-Wagoner theories, the dimensionless coupling between the scalar field and matter can be expanded as in Eq. (3), where and are dimensionless constants, and is the asymptotic value of the scalar field. The leading-order term in this expansion, i.e. , is severely constrained by Solar System experiments: , or  Will:2014kxa (). Furthermore , otherwise spontaneous scalarization would affect the dynamics in ways that are severely constrained by binary pulsar data (note that these constraints are mildly dependent on, and slightly degenerate with, the equation of state of nuclear matter Shibata:2013pra (); Silva:2014fca (); Shao:2017gwu ()).

Recently, Barausse et al. Barausse:2012da () and Palenzuela et al. Palenzuela:2013hsa () (using numerical simulations and semi-analytical arguments, respectively) discovered a phenomenon similar to spontaneous scalarization in the late stages of the evolution of binary neutron star binaries, that they called “dynamical scalarization” (see also Shibata:2013pra (); Sennett:2016rwa ()). Even when the individual neutron stars would not spontaneously scalarize in isolation, the scalar field inside each star can grow in amplitude when the binary separation decreases to about , affecting the binary dynamics and speeding up the merger with respect to GR. The resulting gravitational waveforms are significantly different from GR at frequencies , and deviations at even lower frequencies are possible for certain binary systems and theory parameters. Therefore, the effects of dynamical scalarization are in principle detectable (at least in some cases) by LIGO/Virgo like detectors for values of the coupling parameters and that are still allowed by Solar System and binary pulsar tests Shao:2017gwu (): see Sampson:2014qqa () for an extensive discussion.

Unfortunately, the coupling parameters needed in order to obtain any type of scalarization seem to be incompatible with Solar System observations when accounting for the cosmological evolution of the scalar field Damour:1992kf (); Damour:1993id (); Sampson:2014qqa (); Anderson:2016aoi (). This is because for scalarization to occur, typically must be negative (or positive and very large Mendes:2016fby (), which leads to instabilities Palenzuela:2015ima ()), but this also forces the scalar field to evolve cosmologically away from General Relativity, thus leading to maximal deviations by the present epoch in the Solar System. That is, when accounting for the cosmological evolution of the scalar field, the observation of the Shapiro time delay and its consistency with GR essentially require to be positive (and scalarization not to occur) in the simplest models. One way out of this is to allow the scalar field to have a mass or to fine-tune the cosmological initial conditions, so that the scalar field does not evolve cosmologically at all.

Deviations from GR may also be observable in the electromagnetic signal from binaries of magnetized NSs Ponce:2014hha (), and therefore they could be tested by multi-messenger observations similar to GW170817 TheLIGOScientific:2017qsa (). This work was recently generalized to certain models – namely, – which are equivalent to massive scalar-tensor theories. The inspiral/merger dynamics of neutron star binaries in these theories, as observed by GW detectors, can set stringent bounds on attractive finite-range scalar forces Sagunski:2017nzb ().

4.3 Compact binaries in other theories of gravity

As discussed above, our ability to simulate compact binary mergers in modified theories of gravity that go beyond “plain vanilla” scalar-tensor theories is limited by our poor understanding of their well-posedness.

Jai-akson et al. Jai-akson:2017ldo () estimated the outcome of a merger in Einstein-Maxwell-dilaton gravity (whose action arises in the low-energy limit of string theory) using a “geodesic analogy” approach inspired by Buonanno:2007sv (), rather than numerical simulations. Their qualitative conclusions were confirmed by full-blown numerical simulations of isolated and binary black hole systems Hirschmann:2017psw (). For the coupling parameters considered in their work, the dilaton can largely be ignored, and GW signals from binary black hole systems in these theories are hard to distinguish from their GR counterparts.

Another line of inquiry tries to understand (or bypass) well-posedness issues, such as the study of the initial value problem in Lovelock gravity and in a subclass of Horndeski theories Papallo:2017qvl () including EdGB gravity Papallo:2017ddx (), and the Israel-Stewart Israel:1976tn (); Israel:1979wp () type approach Cayuso:2017iqc (), as we already reviewed. Another idea to cure these pathologies involves an iterative strategy, where the zeroth-order solution provided by GR is employed to evaluate higher-order “corrections” Endlich:2017tqa (). Okounkova et al. Okounkova:2017yby () performed binary black hole merger simulations in dCS gravity using a well-posed perturbation scheme for numerically integrating beyond-GR theories that have a continuous limit to GR, and working to linear order in the perturbation parameter. They computed gravitational waveforms in GR and energy fluxes of the dCS pseudo-scalar field, finding good agreement with analytic predictions at early times Yagi:2011xp () (including the absence of pseudo-scalar dipole radiation) and discovering new phenomenology, including in particular a burst of dipole radiation during merger. Perhaps unsurprisingly, they found that LIGO observations could place bounds on the new dCS length scale that are approximately comparable to the size of the black hole horizon, i.e.  km.

More interesting merger dynamics could occur in the presence of black hole scalarization phenomena, which have recently been studied in Einstein-Maxwell-dilaton Julie:2017rpw () and scalar-Gauss-Bonnet Silva:2017uqg () gravity models.

5 Outlook and Discussion

We end this review by discussing some important future directions that need to be pursued. In order to improve our ability to probe extreme gravity with GWs, one needs to prepare parametrized template waveforms that capture non-GR effects not only in the inspiral, but also in the merger-ringdown phase. Ideally one would achieve this goal through black hole merger simulations in various non-GR theories, which are necessary to understand qualitatively and quantitatively how non-GR modifications affect the merger-ringdown phase, and to construct complete inspiral-merger-ringdown parameterized waveforms. These simulations require a preliminary understanding of the well-posedness properties of modified theories of gravity. Some recent attempts to build parametrized frameworks for the post-merger/ringdown waveform (see e.g. Glampedakis:2017dvb (); Tattersall:2017erk ()) will be discussed in the second part of this contribution.

Current work linking GW observations with the “fundamental” parameters of various modified theories of gravity Yunes:2016jcc () can and should be improved in various ways. First, one needs to use the actual data and carry out a Bayesian analysis to derive bounds on modified theories of gravity with GW observations. Second, one needs to study how the bounds on each theory improve by combining all the observed GW events. Bounds on generic non-GR parameters in the waveform with multiple events have been derived by the LIGO/Virgo Collaboration in TheLIGOScientific:2016pea (); Abbott:2017vtc (); Abbott:2017oio (), though it is difficult to map such combined bounds on generic parameters to those on theoretical coupling constants in each theory. This is because the former depends not only on the latter, but also on source parameters (like masses) that are different from one source to another. Thus, to derive combined bounds on each theory, one may need to carry out a theory-based (rather than model-independent) analysis. Third, one needs to repeat and extend these calculations taking into account the binary neutron star merger event GW170817.

Other ways to improve bounds on non-GR theories include deriving subleading PN corrections to the gravitational waveform. Higher-order corrections, while subdominant in certain theories (such as Brans-Dicke theory Yunes:2016jcc ()), may be important in other theories. For example, the leading PN correction in EdGB gravity is suppressed for nearly equal-mass and equal-spin black hole binaries (unless scalarization occurs Silva:2017uqg ()), and thus the subleading 0PN correction may dominate the PN contribution.

Another avenue for future work includes deriving black hole sensitivities in theories like EA and khronometric gravity. These sensitivities can be calculated by constructing a slowly moving black hole solution relative to the vector field, following previous work on neutron star sensitivities Yagi:2013qpa (); Yagi:2013ava (). Reference Yunes:2016jcc () used the 0PN correction to the waveform for these theories as the dominant contribution, but the PN correction due to scalar and vector radiation (that depends on the sensitivities) may dominate over the 0PN correction.

One also needs to construct non-GR waveforms for generic binaries with precessing and eccentric orbits. For example, the dCS correction to the waveform in Table 2 is only valid for spin-aligned systems, and thus it is important to derive the waveform for precessing binaries in this theory, as spins are crucial to constrain parity violation in gravity. However in order to construct parameterized non-GR template waveforms for generic orbits, one needs to first construct such waveforms with sufficient reliability within GR.

Another important avenue to pursue is to reveal the importance of corrections to GW generation when probing corrections to GW propagation. Reference Yunes:2016jcc () showed that generation effects are negligible compared to the propagation effects for a specific model of massive gravity, but this may not be true in other modified gravity theories.

Acknowledgements.
E.B. was supported by NSF Grants No. PHY-1607130 and AST-1716715. N.Y. acknowledges support through the NSF CAREER Grant PHY-1250636 and NASA Grants NNX16AB98G and 80NSSC17M0041.

Appendix A Derivation of the Black Hole Scalar Charge in decoupled dynamical Gauss-Bonnet Gravity

The goal of this appendix is to derive the scalar charge in DGB gravity for a stationary black hole, valid to arbitrary order in spin. We closely follow the calculation in dCS gravity in Yagi:2012vf (). The scalar charge can be read off from the coefficient in the asymptotic behavior of the scalar field at spatial infinity as , where is the black hole mass. Since we work within the small-coupling approximation, we can take the background metric to be Kerr, and the above equation becomes

(23)

where we work in the rescaled radial coordinate and with representing the dimensionless Kerr parameter, and

(24)

with .

In order to solve the above field equation using Green’s functions, we decompose the scalar field and the source term as Teukolsky:1973ap ()

(25)
(26)

where is normalized as

(27)

Inverting Eq. (26), one obtains

(28)

Eq. (23) can be split into radial and angular parts as

(29)
(30)

The solution to the second equation is nothing but the mode of the spherical harmonics .

Let us first derive the scalar monopole charge by concentrating on the mode. The solution to Eq. (29) consists of homogeneous and particular solutions. Let us first study the former. Modulo overall integration constants, homogeneous solutions for the mode of Eq. (29) are given by

(31)

The asymptotic behavior of at spatial infinity and at the Kerr horizon is given by

(33)
(34)

Since EdGB gravity is a shift-symmetric theory, one can set without loss of generality (namely no contribution from