Gravitational-Wave Tests of General Relativity with Ground-Based Detectors and Pulsar-Timing Arrays

Abstract

This review is focused on tests of Einstein’s theory of general relativity with gravitational waves that are detectable by ground-based interferometers and pulsar-timing experiments. Einstein’s theory has been greatly constrained in the quasi-linear, quasi-stationary regime, where gravity is weak and velocities are small. Gravitational waves will allow us to probe a complimentary, yet previously unexplored regime: the non-linear and dynamical strong-field regime. Such a regime is, for example, applicable to compact binaries coalescing, where characteristic velocities can reach fifty percent the speed of light and gravitational fields are large and dynamical. This review begins with the theoretical basis and the predicted gravitational-wave observables of modified gravity theories. The review continues with a brief description of the detectors, including both gravitational-wave interferometers and pulsar-timing arrays, leading to a discussion of the data analysis formalism that is applicable for such tests. The review ends with a discussion of gravitational-wave tests for compact binary systems.

1 Introduction

1.1 The importance of testing

The era of precision gravitational-wave astrophysics is at our doorstep. With it, a plethora of previously unavailable information will flood in, allowing for unprecedented astrophysical measurements and tests of fundamental theories. Nobody would question the importance of more precise astrophysical measurements, but one may wonder whether fundamental tests are truly necessary, considering the many successes of Einstein’s theory of general relativity (GR). Indeed, GR has passed many tests, including solar system ones, binary pulsar ones and cosmological ones (for a recent review, see [438, 359]).

What all of these tests have in common is that they sample the quasi-stationary, quasi-linear weak field regime of GR. That is, they sample the regime of spacetime where the gravitational field is weak relative to the mass-energy of the system, the characteristic velocities of gravitating bodies are small relative to the speed of light, and the gravitational field is stationary or quasi-stationary relative to the characteristic size of the system. A direct consequence of this is that gravitational waves emitted by weakly-gravitating, quasi-stationary sources are necessarily extremely weak. To make this more concrete, let us define the gravitational compactness as a measure of the strength of the gravitational field:

$${\mathcal C} = {{\mathcal M} \over {\mathcal R}},$$

(1)

where $\mathcal{M}$ is the characteristic mass of the system, $\mathcal{R}$ is the characteristic length scale associated with gravitational radiation, and henceforth we set G = c = 1. For binary systems, the orbital separation serves as this characteristic length scale. The strength of gravitational waves and the mutual gravitational interaction between bodies scale linearly with this compactness measure. Let us also define the characteristic velocities of such a system $\mathcal{V}$ as a quantity related to the rate of change of the gravitational field in the center of mass frame. We can then more formally define the weak field as the region of spacetime where the following two conditions are simultaneously satisfied:

$${\bf{weak}} \,\, {\bf{Field}}:\quad {\mathcal C} \ll 1,\quad {\mathcal V} \ll 1.$$

(2)

By similarity, the strong field is defined as the region of spacetime where both conditions in Eq. (2) are not valid simultaneously¹.

Let us provide some examples. For the Earth-Sun system, $\mathcal{M}$ is essentially the mass of the sun, while $\mathcal{R}$ is the orbital separation, which leads to $\mathcal{C}\approx 9.8\times 10^{-9}$ and $\mathcal{V}\approx 9.9\times 10^{-5}$. Even if an object were in a circular orbit at the surface of the sun, its gravitational compactness would be $\mathcal{O}(10^{-6})$ and its characteristic velocity $\mathcal{O}(10^{-3})$. Thus, we conclude that all solar-system experiments are necessarily sampling the weak field regime of gravity. Similarly, for the binary pulsar J0737−3039 [298, 274], $\mathcal{C}\approx 6\times 10^{-6}$ and $\mathcal{V}\approx 2\times 10^{-3}$, where we have set the characteristic length $\mathcal{R}$ to the orbital separation via $\mathcal{R}\approx [MP^{2}/(4\pi^{2})]^{1/3}\approx 10^{6}$ km, where P = 0.1 days is the orbital period and M ≈ 3 M_⊙ is the total mass. Although neutron stars are sources of strong gravity (the ratio of their mass to their radius is of order one tenth), binary pulsars are most sensitive to the quasi-static part of the post-Newtonian effective potential or to the leading-order (Newtonian piece) of the radiation-reaction force. On the other hand, in compact binary coalescence the gravitational compactness and the characteristic speed can reach values much closer to unity. Therefore, although in much of the pulsar-timing literature binary pulsar timing is said to allow for strong-field tests of gravity, gravitational information during compact binary coalescence would be a much-stronger—field test.

Even though current data does not give us access to the full non-linear and dynamical regime of GR, solar-system tests and binary-pulsar observations have served (and will continue to serve) an invaluable role in testing Einstein’s theory. Solar-system tests effectively cured an outbreak of modified gravity theories in the 1970s and 1980s, as summarized for example in [438]. Binary pulsars were crucial as the first indirect detectors of gravitational waves, and later to kill certain theories, like Rosen’s bimetric gravity [365], and heavily constrain others that predict dipolar energy loss, as we see in Sections 2 and 5. Binary pulsars are probes of GR in a certain sector of the strong field: in the dynamical but quasi-linear sector, verifying that compact objects move as described by a perturbative, post-Newtonian analysis to leading order. Binary pulsars can be used to test GR in the “strong field” only in the sense that they probe non-linear stellar-structure effects, but they say very little to nothing about non-linear radiative effects. Similarly, future electromagnetic observations of black-hole—accretion disks may probe GR in another strong-field sector: the non-linear but fully stationary regime, verifying that black holes are described by the Kerr metric. As of this writing, only gravitational waves will allow for tests of GR in the full strong-field regime, where gravity is both heavily non-linear and inherently dynamical.

No experiments exist to date that validate Einstein’s theory of GR in the highly-dynamical, strong-field region. Due to previous successes of GR, one might consider such validation unnecessary. However, as most scientists would agree, the role of science is to predict and verify and not to assume without proof. Moreover, the incompleteness of GR in the quantum regime, together with the somewhat unsatisfactory requirement of the dark sector of cosmology (including dark energy and dark matter), have prompted more than one physicist to consider deviations from GR more seriously. Gravitational waves will soon allow us to verify Einstein’s theory in a regime previously inaccessible to us, and as such, these tests are invaluable.

However, in many areas of physics GR is so ingrained that questioning its validity (even in a regime where Einstein’s theory has not yet been validated) is synonymous with heresy. Dimensional arguments are usually employed to argue that any quantum gravitational correction will necessarily and unavoidably be unobservable with gravitational waves, as the former are expected at a (Planck) scale that is inaccessible to gravitational-wave detectors. This rationalization is dangerous, as it introduces a theoretical bias in the analysis of new observations of the universe, thus reducing the potential for new discoveries. For example, if astrophysicists had followed such a rationalization when studying supernova data, they would not have discovered that the universe is expanding. Dimensional arguments suggest that the cosmological constant is over 100 orders of magnitude larger than the value required to agree with observations. When observing the universe for the first time in a completely new way, it seems more conservative to remain agnostic about what is expected and what is not, thus allowing the data itself to guide our efforts toward theoretically understanding the gravitational interaction.

1.2 Testing general relativity versus testing alternative theories

When testing GR, one considers Einstein’s theory as a null hypothesis and searches for generic deviations. On the other hand, when testing alternative theories one starts from a particular modified gravity model, develops its equations and solutions and then predicts certain observables that then might or might not agree with experiment. Similarly, one may define a bottom-up approach versus a top-down approach. In the former, one starts from some observables in an attempt to discover fundamental symmetries that may lead to a more complete theory, as was done when constructing the standard model of elementary particles. On the other hand, a top-down approach starts from some fundamental theory and then derives its consequence.

Both approaches possess strengths and weaknesses. In the top-down approach one has complete control over the theory under study, being able to write down the full equations of motion, answer questions about well-posedness and stability of solutions, and predict observables. But, as we see in Section 2, carrying out such an approach can be quixotic within any one model. What is worse, the lack of a complete and compelling alternative to GR makes choosing a particular modified theory difficult.

Given this, one might wish to attempt a bottom-up approach, where one considers a set of principles one wishes to test without explicit mention of any particular theory. One usually starts by assuming GR as a null-hypothesis and then considers deformations away from GR. The hope is that experiments will be sensitive to such deformations, thus either constraining the size of the deformations or pointing toward a possible inconsistency. But if experiments do confirm a GR deviation, a bottom-up approach fails at providing a given particular action from which to derive such a deformation. In fact, there can be several actions that lead to similar deformations, all of which can be consistent with the data within its experimental uncertainties.

Nonetheless, both approaches are complementary. The bottom-up approach draws inspiration from particular examples carried out in the top-down approach. Given a plausible measured deviation from GR within a bottom-up approach, one will still need to understand what plausible top-down theories can lead to such deviations. From this standpoint, then, both approaches are intrinsically intertwined and worth pursuing.

1.3 Gravitational-wave tests versus other tests of general relativity

Gravitational-wave tests differ from other tests of GR in many ways. Perhaps one of the most important differences is the spacetime regime gravitational waves sample. Indeed, as already mentioned, gravitational waves have access to the most extreme gravitational environments in nature. Moreover, gravitational waves travel essentially unimpeded from their source to Earth, and thus, they do not suffer from issues associated with obscuration. Gravitational waves also exist in the absence of luminous matter, thus allowing us to observe electromagnetically dark objects, such as black-hole inspirals.

This last point is particularly important as gravitational waves from inspiral—black-hole binaries are one of the cleanest astrophysical systems in nature. In the last stages of inspiral, when such gravitational waves would be detectable by ground-based interferometers, the evolution of a blackhole binary is essentially unaffected by any other matter or electromagnetic fields present in the system. As such, one does not need to deal with uncertainties associated with astrophysical matter. Unlike other tests of GR, such as those attempted with accretion-disk observations, black-hole—binary gravitational-wave tests may well be the cleanest probes of Einstein’s theory.

Of course, what is an advantage here, can also be a huge disadvantage in another context. Gravitational waves from compact binaries are intrinsically transient (they turn on for a certain amount of time and then shut off). This is unlike binary pulsar systems, for which astrophysicists have already collected tens of years of data. Moreover, gravitational wave tests rely on specific detections that cannot be anticipated beforehand. This is in contrast to Earth-based laboratory experiments, where one has complete control over the experimental setup. Finally, the intrinsic weakness of gravitational waves makes detection a very difficult task that requires complex data-analysis algorithms to extract signals from the noise. As such, gravitational-wave tests are limited by the signal-to-noise ratio and affected by systematics associated with the modeling of the waves, issues that are not as important in other loud astrophysical systems.

1.4 Ground-based vs space-based detectors and interferometers vs pulsar timing

This review article focuses only on ground-based detectors, by which we mean both gravitational-wave interferometers, such as the Laser Interferometer Gravitational Observatory (LIGO) [3, 2, 217], Virgo [5, 6] and the Einstein Telescope (ET) [361, 377], as well as pulsar-timing arrays (for a recent review of gravitational-wave tests of GR with space-based detectors, see [183, 446]). Ground-based detectors have the limitation of being contaminated by man-made and nature-made noise, such as ground and air traffic, logging, earthquakes, ocean tides and waves, which are clearly absent in space-based detectors. Ground-based detectors, however, have the clear benefit that they can be continuously upgraded and repaired in case of malfunction, which is obviously not possible with space-based detectors.

As far as tests of GR are concerned, there is a drastic difference in space-based and ground-based detectors: the gravitational-wave frequencies these detectors are sensitive to. For various reasons that we will not go into, space-based interferometers are likely to have million kilometer long arms, and thus, be sensitive in the milli-Hz band. On the other hand, ground-based interferometers are bound to the surface and curvature of the Earth, and thus, they have kilometer-long arms and are sensitive in the deca- and hecta-Hz band. Different types of interferometers are then sensitive to different types of gravitational-wave sources. For example, when considering binary coalescences, ground-based interferometers are sensitive to late inspirals and mergers of neutron stars and stellar-mass black holes, while space-based detectors will be sensitive to supermassive—black-hole binaries with masses around 10⁵ M_⊙.

The impact of a different population of sources in tests of GR depends on the particular modified gravity theory considered. When studying quadratic gravity theories, as we see in Section 2, the Einstein—Hilbert action is modified by introducing higher-order curvature operators, which are naturally suppressed by powers of the inverse of the radius of curvature. Thus, space-based detectors will not be ideal at constraining these theories, as the radius of curvature of supermassive black holes is much larger than that of stellar-mass black holes at merger. Moreover, space-based detectors will not be sensitive to neutron-star-binary coalescences; they are sensitive to supermassive black-hole/neutron-star coalescences, where the radius of curvature of the system is controlled by the supermassive black hole.

On the other hand, space-based detectors are unique in their potential to probe the spacetime geometry of supermassive black holes through gravitational waves emitted during extreme-massratio inspirals. These inspirals consist of a stellar-mass compact object in a generic decaying orbit around a supermassive black hole. Such inspirals produce millions of cycles of gravitational waves in the sensitivity band of space-based detectors (in fact, they can easily out-live the observation time!). Therefore, even small changes to the radiation-reaction force, or to the background geometry, can lead to noticeable effects in the waveform observable and thus strong tests of GR, albeit constrained to the radius of curvature of the supermassive black hole. For recent work on such systems and tests, see [23, 370, 371, 263, 196, 50, 289, 182, 390, 471, 31, 297, 184, 116, 93, 183].

Space-based detectors also have the advantage of range, which is particularly important when considering theories where gravitons do not travel at light speed [316]. Space-based detectors have a horizon distance much larger than ground-based detectors; the former can see black-hole mergers to redshifts of order 10 if there are any at such early times in the universe, while the latter are confined to events within redshift 1. Gravitational waves emitted from distant regions in spacetime need a longer time to propagate from the source to the detectors. Thus, theories that modify the propagation of gravitational waves will be best constrained by space-based type systems. Of course, such theories are also likely to modify the generation of gravitational waves, which ground-based detectors should also be sensitive to.

Another important difference between detectors is in their response to an impinging gravitational wave. Ground-based detectors, as we see in Section 3, cannot separate between the two possible scalar modes (the longitudinal and the breathing modes) of metric theories of gravity, due to an intrinsic degeneracy in the response functions. Space-based detectors in principle also possess this degeneracy, but they may be able to break it through Doppler modulation if the interferometer orbits the Sun. Pulsar-timing arrays, on the other hand, lack this degeneracy altogether, and thus, they can in principle constrain the existence of both modes independently.

Pulsar-timing arrays differ from interferometers in their potential to test GR mostly by the frequency space they are most sensitive to. The latter can observe the late inspiral and merger of compact binaries, while the former is restricted to the very early inspiral. This is why pulsar timing arrays do not need very accurate waveform templates that account for the highly-dynamical and non-linear nature of gravity to detect gravitational waves; leading-order quadrupole waveforms are sufficient [120]. In turn, this implies that pulsar timing arrays cannot constrain theories that only deviate significantly from GR in the late inspiral, while they are exceptionally well-suited for constraining low-frequency deviations.

Therefore, we see a complementarity emerging: different detectors can test GR in different complementary regimes:

• Ground-based detectors are best at constraining higher-curvature type modified theories that deviate from GR the most in the late inspiral and merger phase.

• Space-based detectors are best at constraining modified graviton dispersion relations and the geometry of supermassive compact objects.

• Pulsar-timing arrays are best at independently constraining the existence of both scalar modes and any deviation from GR that dominates at low orbital frequencies.

Through the simultaneous implementation of all these tests, GR can be put on a much firmer footing in all phases of the strong-field regime.

1.5 Notation and conventions

We mainly follow the notation of [318], where Greek indices stand for spacetime coordinates and spatial indices in the middle of the alphabet (i, j, k, …) for spatial indices. Parenthesis and square brackets in index lists stand for symmetrization and antisymmetrization respectively, e.g., A_(μν) = (A _μν + A _νμ )/2 and A_[μν] = (A _μν − A _νμ )/2. Partial derivatives with respect to spacetime and spatial coordinates are denoted ∂ _μ A = A, _μ and ∂ _i A = A, _i respectively. Covariant differentiation is denoted ∇_μA = A_;μ, multiple covariant derivatives ∇^μν… = ∇^μ∇^ν …, and the curved spacetime D’Alembertian □A = ∇ _μ ∇^μ A. The determinant of the metric g _μν is g, R _μνδσ is the Riemann tensor, R _μν is the Ricci tensor, R is the Ricci scalar and g _μν is the Einstein tensor. The Levi-Civita tensor and symbol are ∈^μνδσ and ∈^−μνδσ respectively, with ${\bar \varepsilon ^{0123}} = + 1$ in an orthonormal, positively-oriented frame. We use geometric units (G = c =1) and the Einstein summation convention is implied.

We will mostly be concerned with metric theories, where gravitational radiation is only defined much farther than a gravitational-wave wavelength from the source. In this far or radiation zone, the metric tensor can be decomposed as

$${g_{\mu \nu}} = {\eta _{\mu \nu}} + {h_{\mu \nu}},$$

(3)

with η _μν the Minkowski metric and h _μν the metric perturbation. If the theory considered has additional fields ϕ, these can also be decomposed in the far zone as

$$\phi = {\phi _0} + \psi ,$$

(4)

with ϕ₀ the background value of the field and Ψ a perturbation. With such a decomposition, the field equations for the metric will usually be wave equations for the metric perturbation and for the field perturbation, in a suitable gauge.

2 Alternative Theories of Gravity

In this section, we discuss the many possible alternative theories that have been studied so far in the context of gravitational-wave tests. We begin with a description of the theoretically desirable properties that such theories must have. We then proceed with a review of the theories so far explored as far as gravitational waves are concerned. We will leave out the description of many theories in this chapter, especially those which currently lack a gravitational-wave analysis. We will conclude with a brief description of unexplored theories as possible avenues for future research.

2.1 Desirable theoretical properties

The space of possible theories is infinite, and thus, one is tempted to reduce it by considering a subspace that satisfies a certain number of properties. Although the number and details of such properties depend on the theorist’s taste, there is at least one fundamental property that all scientists would agree on:

• 1. Precision Tests. The theory must produce predictions that pass all solar system, binary pulsar, cosmological and experimental tests that have been carried out so far.

This requirement can be further divided into the following:

• 1.a General Relativity Limit. There must exist some limit, continuous or discontinuous, such as the weak-field one, in which the predictions of the theory are consistent with those of GR within experimental precision.

• 1.b Existence of Known Solutions [426]. The theory must admit solutions that correspond to observed phenomena, including but not limited to (nearly) flat spacetime, (nearly) Newtonian stars, and cosmological solutions.

• 1.c Stability of Solutions [426]. The special solutions described in property (1.b) must be stable to small perturbations on timescales smaller than the age of the universe. For example, perturbations to (nearly) Newtonian stars, such as impact by asteroids, should not render such solutions unstable.

Of course, these properties are not all necessarily independent, as the existence of a weak-field limit usually also implies the existence of known solutions. On the other hand, the mere existence of solutions does not necessarily imply that these are stable.

In addition to these fundamental requirements, one might also wish to require that any new modified gravity theory possesses certain theoretical properties. These properties will vary depending on the theorist, but the two most common ones are listed below:

• 2. Well-motivated from Fundamental Physics. There must be some fundamental theory or principle from which the modified theory (effective or not) derives. This fundamental theory would solve some fundamental problem in physics, such as late-time acceleration or the incompatibility between quantum mechanics and GR.

• 3. Well-posed Initial Value Formulation [426]. A wide class of freely specifiable initial data must exist, such that there is a uniquely determined solution to the modified field equations that depends continuously on this data.

The second property goes without saying at some level, as one expects modified-gravity-theory constructions to be motivated from some (perhaps yet incomplete) quantum-gravitational description of nature. As for the third property, the continuity requirement is necessary because otherwise the theory would lose predictive power, given that initial conditions can only be measured to a finite accuracy. Moreover, small changes in the initial data should not lead to solutions outside the causal future of the data; that is, causality must be preserved. Section 2.2 expands on this well-posedness property further.

One might be concerned that Property (2) automatically implies that any predicted deviation to astrophysical observables will be too small to be detectable. This argument usually goes as follows. Any quantum gravitational correction to the action will “naturally” introduce at least one new scale, and this, by dimensional analysis, must be the Planck scale. Since this scale is usually assumed to be larger than 1 TeV in natural units (or 10⁻³⁵ m in geometric units), gravitational-wave observations will never be able to observe quantum-gravitational modifications (see, e.g., [155] for a similar argument). Although this might be true, in our view such arguments can be extremely dangerous, since they induce a certain theoretical bias in the search for new phenomena. For example, let us consider the supernova observations of the late-time expansion of the universe that led to the discovery of the cosmological constant. The above argument certainly fails for the cosmological constant, which on dimensional arguments is over 100 orders of magnitude too small. If the supernova teams had respected this argument, they would not have searched for a cosmological constant in their data. Today, we try to explain our way out of the failure of such dimensional arguments by claiming that there must be some exquisite cancellation that renders the cosmological constant small; but this, of course, came only after the constant had been measured. One is not trying to argue here that cancellations of this type are common and that quantum gravitational modifications are necessarily expected in gravitational-wave observations. Rather, we are arguing that one should remain agnostic about what is expected and what is not, and allow oneself to be surprised without suppressing the potential for new discoveries that will accompany the new era of gravitational-wave astrophysics.

One last property that we wish to consider for the purposes of this review is:

• 4. Strong Field Inconsistency. The theory must lead to observable deviations from GR in the strong-field regime. Many modified gravity models have been proposed that pose infrared or cosmological modifications to GR, aimed at explaining certain astrophysical or cosmological observables, like the late expansion of the universe. Such modified models usually reduce to GR in the strong-field regime, for example via a Vainshtein-like mechanism [413, 140, 45] in a static spherically-symmetric context. Extending this mechanism to highly-dynamical strong-field scenarios has not been fully worked out yet [137, 138]. Gravitational-wave tests of GR, however, are concerned with modified theories that predict deviations in the strong-field, precisely where cosmological modified models do not. Clearly, Property (4) is not necessary for a theory to be a valid description of nature. This is because a theory might be identical to GR in the weak and strong fields, yet different at the Planck scale, where it would be unified with quantum mechanics. However, Property (4) is a desirable feature if one is to test this theory with gravitational wave observations.

2.2 Well-posedness and effective theories

Property (3) not only requires the existence of an initial value formulation, but also that it be well posed, which is not necessarily guaranteed. For example, the Cauchy-Kowalewski theorem states that a system of n partial differential equations for n unknown functions ϕ _i of the form ϕ_i,tt = F _i (x^μ; ϕ_{j, μ}; ϕ_j,ti; ϕ_j,ik), with F _i analytic functions has an initial value formulation (see, e.g., [425]). However, this theorem does not guarantee continuity or the causal conditions described above. For this, one has to rely on more general energy arguments, for example constructing a suitable energy measure that obeys the dominant energy condition and using it to show well-posedness (see, e.g., [225, 425]). One can show that second-order, hyperbolic partial differential equations, i.e., equations of the form

$${\nabla ^\mu}{\nabla _\mu}\phi + {A^\mu}{\nabla _\mu}\phi + B\phi + C = 0,$$

(5)

where A^μ is an arbitrary vector field and (B, C) are smooth functions, have a well-posed initial value formulation. Moreover, the Leray theorem proves that any quasilinear, diagonal, second-order hyperbolic system also has a well-posed initial value formulation [425].

Proving the well-posedness of an initial-value formulation for systems of higher-than-second-order, partial differential equations is much more difficult. In fact, to our knowledge, no general theorems exist of the type described above that apply to third, fourth or higher-order, partial, nonlinear and coupled differential equations. Usually, one resorts to the Ostrogradski theorem [337] to rule out (or at the very least cast serious doubt on) theories that lead to such higher-order field equations. Ostrogradski’s theorem states that Lagrangians that contain terms with higher-than-first-time derivatives possess a linear instability in the Hamiltonian (see, e.g., [443] for a nice review).² As an example, consider the Lagrangian density

$${\mathcal L} = {m \over 2}{\dot q^2} - {{m{\omega ^2}} \over 2}{q^2} - {{gm} \over {2{\omega ^2}}}{\ddot q^2},$$

(6)

whose equations of motion,

$$\ddot q + {\omega ^2}q = - {g \over {{\omega ^2}}}{\overset{\ldots.}{q}},$$

(7)

obviously contain higher derivatives. The exact solution to this differential equation is

$$q = {A_1}\cos {k_1}t + {B_1}\sin {k_1}t + {A_2}\cos {k_2}t + {B_2}\sin {k_2}t,$$

(8)

where (A _i , B _i ) are constants and $k_{1,2}^2/{\omega ^2} = (1 \mp \sqrt {1 - 4g)/(2g)}$. The on-shell Hamiltonian is then

$$H = {m \over 2}\sqrt {1 - 4g} k_1^2(A_1^2 + B_1^2) - {m \over 2}\sqrt {1 - 4g} k_2^2(A_2^2 + B_2^2),$$

(9)

from which it is clear that mode 1 carries positive energy, while mode 2 carries negative energy and forces the Hamiltonian to be unbounded from below. The latter implies that dynamical degrees of freedom can reach arbitrarily negative energy states. If interactions are present, then an “empty” state would instantaneously decay into a collection of positive and negative energy particles, which cannot describe the universe we live in [443].

However, the Ostrogradski theorem [337] can be evaded if the Lagrangian in Eq. (6) describes an effective theory, i.e., a theory that is a truncation of a more general or complete theory. Let us reconsider the particular example above, assuming now that the coupling constant g is an effective theory parameter and Eq. (6) is only valid to linear order in g. One approach is to search for perturbative solutions of the form q _pert = x₀ + gx₁ + …, which leads to the system of differential equations

$${\ddot x_n} + {\omega ^2}{x_n} = - {1 \over {{\omega ^2}}}{\overset{\ldots.}{x}}_{n - 1},$$

(10)

with x₋₁ =0. Solving this set of n differential equations and resumming, one finds

$${q_{{\rm{pert}}}} = {A_1}\cos {k_1}t + {B_1}\sin {k_1}t.$$

(11)

Notice that q _pert contains only the positive (well-behaved) energy solution of Eq. (8), i.e., perturbation theory acts to retain only the well-behaved, stable solution of the full theory in the g → 0 limit. One can also think of the perturbative theory as the full theory with additional constraints, i.e., the removal of unstable modes, which is why such an analysis is sometimes called perturbative constraints [117, 118, 466].

Another way to approach effective field theories that lead to equations of motion with higherorder derivatives is to apply the method of order reduction. In this method, one substitutes the low-order derivatives of the field equations into the high-order derivative part, thus rendering the resulting new theory usually well posed. One can think of this as a series resummation, where one changes the non-linear behavior of a function by adding uncontrolled, higher-order terms. Let us provide an explicit example by reconsidering the theory in Eq. (6). To lowest order in g, the equation of motion is that of a simple harmonic oscillator,

$$\ddot q + {\omega ^2}q ={\mathcal O}(g),$$

(12)

which is obviously well posed. One can then order-reduce the full equation of motion, Eq. (7), by substituting Eq. (12) into the right-hand side of Eq. (7). Doing so, one obtains the order-reduced equation of motion

$$\ddot q + {\omega ^2}q = g\ddot q + {\mathcal O}({g^2}),$$

(13)

which now clearly has no high-order derivatives and is well posed, provided g ≪ 1. The solution to this order-reduced differential equation is q _pert once more, but with k₁ linearized in g ≪ 1. Therefore, the solutions obtained with a perturbative decomposition and with the order-reduced equation of motion are the same to linear order in g. Of course, since an effective field theory is only defined to a certain order in its perturbative parameter, both treatments are equally valid, with the unstable mode effectively removed in both cases.

However, such a perturbative analysis can say nothing about the well-posedness of the full theory from which the effective theory derives, or of the effective theory if treated as an exact one (i.e., not as a perturbative expansion). In fact, a well-posed full theory may have both stable and unstable solutions. The arguments presented above only discuss the stability of solutions in an effective theory, and thus, they are self-consistent only within their perturbative scheme. A full theory may have non-perturbative instabilities, but these can only be studied once one has a full (non-truncated in g) theory, from which Eq. (6) derives as a truncated expansion. Lacking a full quantum theory of nature, quantum gravitational models are usually studied in a truncated low-energy expansion, where the leading-order piece is GR and higher-order pieces are multiplied by a small coupling constant. One can perturbatively explore the well-behaved sector of the truncated theory about solutions to the leading-order theory. However, such an analysis is incapable of answering questions about well-posedness or non-linear stability of the full theory.

2.3 Explored theories

In this subsection we briefly describe the theories that have so far been studied in some depth as far as gravitational waves are concerned. In particular, we focus only on those theories that have been sufficiently studied so that predictions of the expected gravitational waveforms (the observables of gravitational-wave detectors) have been obtained for at least a typical source, such as the quasi-circular inspiral of a compact binary.

2.3.1 Scalar-tensor theories

Scalar-tensor theories in the Einstein frame [82, 129, 166, 165, 181, 197] are defined by the action (where we will restore Newton’s gravitational constant G in this section)

$$S_{{\rm{ST}}}^{({\rm{E}})} = {1 \over {16\pi G}}\int {{d^4}x\sqrt {- g} [R - 2{g^{\mu \nu}}({\partial _\mu}\varphi)({\partial _\nu}\varphi) - V(\varphi)] + {S_{{\rm{mat}}}}[{\psi _{{\rm{mat}}}},{A^2}(\varphi){g_{\mu \nu}}],}$$

(14)

where φ is a scalar field, A(φ) is a coupling function, V(φ) is a potential function, Ψ_mat represents matter degrees of freedom and G is Newton’s constant in the Einstein frame. For more details on this theory, we refer the interested reader to the reviews [438, 435]. Of course, one can consider more complicated scalar-tensor theories, for example by including multiple scalar fields, but we will ignore such generalizations here.

The Einstein frame is not the frame where the metric governs clocks and rods, and thus, it is convenient to recast the theory in the Jordan frame through the conformal transformation ${\tilde g_{\mu \nu}} = {A^2}(\varphi){\tilde g_{\mu \nu}}$:

$$S_{{\rm{ST}}}^{({\rm{J}})} = {1 \over {16\pi G}}\int {{d^4}x\sqrt {- \tilde g} \left[ {\phi \tilde R - {{\omega (\phi)} \over \phi}{{\tilde g}^{\mu \nu}}({\partial _\mu}\phi)({\partial _\nu}\phi) - {\phi ^2}V} \right] + {S_{{\rm{mat}}}}[{\psi _{{\rm{mat}}}},{{\tilde g}_{\mu \nu}}],}$$

(15)

where ${\tilde g_{\mu \nu}}$ is the physical metric, the new scalar field ϕ is defined via ϕ = A⁻², the coupling field is ω(ϕ) = (α⁻² − 3)/2 and α = A,φ/A. When cast in the Jordan frame, it is clear that scalar-tensor theories are metric theories (see [438] for a definition), since the matter sector depends only on matter degrees of freedom and the physical metric (without a direct coupling of the scalar field). When the coupling ω(ϕ) = ω_bd is constant, then Eq. (15) reduces to the massless version of Jordan-Fierz-Brans-Dicke theory [82].

The modified field equations in the Einstein frame are

$$\begin{array}{*{20}c} {\square\varphi = {1 \over 4}{{dV} \over {d\varphi}} - 4\pi G{{\delta {S_{{\rm{mat}}}}} \over {\delta \varphi}},}\\ {{G_{\mu \nu}} = 8\pi G\left({T_{\mu \nu}^{{\rm{mat}}} + T_{\mu \nu}^{(\varphi)}} \right),}\\ \end{array}$$

(16)

where

$$T_{\mu \nu}^{(\varphi)} = {1 \over {4\pi}}\left[ {{\varphi _{,\mu}}{\varphi _{,\nu}} - {1 \over 2}{g_{\mu \nu}}{\varphi _{,\delta}}{\varphi ^{,\delta}} - {1 \over 4}{g_{\mu \nu}}V(\varphi)} \right]$$

(17)

is a stress-energy tensor for the scalar field. The matter stress-energy tensor is not constructed from the Einstein-frame metric alone, but by the combination A(φ)²g_μν. In the Jordan frame and neglecting the potential, the modified field equations are [435]

$$\begin{array}{*{20}c} {\tilde \square\phi = {1 \over {3 + 2\omega (\phi)}}\left({8\pi {T^{{\rm{mat}}}} - {{d\omega} \over {d\phi}}{{\tilde g}^{\mu \nu}}{\phi _{,\mu}}{\phi _{,\nu}}} \right),\quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ {{{\tilde G}_{\mu \nu}} = {{8\pi G} \over \phi}T_{\mu \nu}^{{\rm{mat}}} + {\omega \over {{\phi ^2}}}\left({{\phi _{,\mu}}{\phi _{,\nu}} - {1 \over 2}{{\tilde g}_{\mu \nu}}{{\tilde g}^{\sigma \rho}}{\phi _{,\sigma}}{\phi _{,\rho}}} \right) + {1 \over \phi}({\phi _{,\mu \nu}} - {{\tilde g}_{\mu \nu}}\tilde \square\phi),}\\ \end{array}$$

(18)

where T^mat is the trace of the matter stress-energy tensor $T_{\mu \nu}^{{\rm{mat}}}$ constructed from the physical metric ${\tilde g_{\mu \nu}}$. The form of the modified field equations in Jordan frame suggest that in the weak-field limit one may consider scalar-tensor theories as modifying Newton’s gravitational constant via G → G(ϕ) = G/ϕ.

Using the decompositions of Eqs. (3)–(4), the field equations of massless Jordan-Fierz-Brans-Dicke theory can be linearized in the Jordan frame to find (see, e.g., [441])

$${\square_\eta}{\theta ^{\mu \nu}} = - 16\pi {\tau ^{\mu \nu}},\quad {\square_\eta}\psi = - 16\pi S,$$

(19)

where □ _η is the D’Alembertian operator of flat spacetime, we have defined a new metric perturbation

$${\theta ^{\mu \nu}} = {h^{\mu \nu}} - {1 \over 2}{\eta ^{\mu \nu}}h - {\psi \over {{\phi _0}}}{\eta ^{\mu \nu}},$$

(20)

i.e., the metric perturbation in the Einstein frame, with h the trace of the metric perturbation and

$${\tau ^{\mu \nu}} = \phi _0^{- 1}T_{{\rm{mat}}}^{\mu \nu} + {t^{\mu \nu}},$$

(21)

$$S = - {1 \over {6 + 4{\omega _{BD}}}}\left({{T^{{\rm{mat}}}} - 3\phi {{\partial {T^{{\rm{mat}}}}} \over {\partial \phi}}} \right)\left({1 - {\theta \over 2} - {\psi \over {{\phi _0}}}} \right) - {1 \over {16\pi}}\left({{\psi _{,\mu \nu}}{\theta ^{\mu \nu}} + {1 \over {{\phi _0}}}{\phi _{,\mu}}{\psi ^{,\mu}}} \right),$$

(22)

with cubic remainders in either the metric perturbation or the scalar perturbation. The quantity ∂T^mat/∂ϕ arises in an effective point-particle theory, where the matter action is a functional of both the Jordan-frame metric and the scalar field. The quantity t^μν is a function of quadratic or higher order in θ^μν or Ψ. These equations can now be solved given a particular physical system, as done for quasi-circular binaries in [441, 374, 336]. Given the above evolution equations, Jordan-Fierz-Brans-Dicke theory possesses a scalar (spin-0) mode, in addition to the two transverse-traceless (spin-2) modes of GR, i.e., Jordan-Fierz-Brans-Dicke theory is of Type N₃ in the E(2) classification [161, 438].

Let us now discuss whether scalar-tensor theories satisfy the properties discussed in Section 2.1. Massless Jordan-Fierz-Brans-Dicke theory agrees with all known experimental tests provided ω_BD > 4 × 10⁴, a bound imposed by the tracking of the Cassini spacecraft through observations of the Shapiro time delay [73]. Massive Jordan-Fierz-Brans-Dicke theory has been recently constrained to ω_BD > 4 × 10⁴ and m_s < 2.5 × 10⁻²⁰ eV, with m_s the mass of the scalar field [348, 20]. Of course, these bounds are not independent, as when m_s → 0 one recovers the standard massless constraint, while when m_s → ∞, ω_BD cannot be bounded as the scalar becomes non-dynamical. Observations of the Nordtvedt effect with Lunar Laser Ranging observations, as well as observations of the orbital period derivative of white-dwarf/neutron-star binaries, yield similar constraints [131, 132, 20, 177]. Neglecting any homogeneous, cosmological solutions to the scalar-field evolution equation, it is clear that in the limit ω → ∞ one recovers GR, i.e., scalar-tensor theories have a continuous limit to Einstein’s theory, but see [164] for caveats for certain spacetimes. Moreover, [375, 278, 425] have verified that scalar-tensor theories with minimal or non-minimal coupling in the Jordan frame can be cast in a strongly-hyperbolic form, and thus, they possess a well-posed initial-value formulation. Therefore, scalar-tensor theories possess both Properties (1) and (3).

Scalar-tensor theories also possess Property (2), since they can be derived from the low-energy limit of certain string theories. The integration of string quantum fluctuations leads to a higher-dimensional string theoretical action that reduces locally to a field theory similar to a scalar-tensor one [189, 176], the mapping being ϕ = e^−2Ψ, with Ψ one of the string moduli fields [133, 134]. Moreover, scalar-tensor theories can be mapped to f(R) theories, where one replaces the Ricci scalar by some functional of R. In particular, one can show that f(R) theories are equivalent to Brans-Dicke theory with ω_BD = 0, via the mapping ϕ = df(R)/dR and V(ϕ) = R df (R)/dR − f(R) [104, 396]. For a recent review on this topic, see [135].

Black holes and stars continue to exist in scalar-tensor theories. Stellar configurations are modified from their GR profile [441, 131, 214, 215, 410, 132, 394, 139, 393, 235], while black holes are not, provided one neglects homogeneous, cosmological solutions to the scalar field evolution equation. Indeed, Hawking [224, 159, 222, 98, 244, 363] has proven that Brans-Dicke black holes that are stationary and the endpoint of gravitational collapse are identical to those of GR. This proof has recently been extended to a general class of scalar-tensor models [398]. That is, stationary black holes radiate any excess “hair”, i.e., additional degrees of freedom, after gravitational collapse, a result sometimes referred to as the no-hair theorem for black holes in scalar-tensor theories. This result has recently been extended even further to allow for quasi-stationary scenarios in generic scalar-tensor theories through the study of extreme-mass-ratio inspirals [465] (small black hole in orbit around a much larger one), post-Newtonian comparable-mass inspirals [315] and numerical simulations of comparable-mass black-hole mergers [230, 67].

Damour and Esposito-Farése [129, 130] proposed a different type of scalar-tensor theory, one that can be defined by the action in Eq. (15) but with the conformal factor $A(\varphi) = {e^{\alpha \varphi + \beta {\varphi ^2}/2}}$ or the coupling function ω(ϕ) = −3/2 − 2πG/(β log ϕ), where α and β are constants. When β = 0 one recovers standard Brans-Dicke theory. When β ≲ −4, non-perturbative effects that develop if the gravitational energy is large enough can force neutron stars to spontaneously acquire a non-trivial scalar field profile, to spontaneously scalarize. Through this process, a neutron-star binary that initially had no scalar hair in its early inspiral would acquire it before merger, when the binding energy exceeded some threshold [51]. Binary pulsar observations have constrained this theory in the (α,β) space; very roughly speaking β > −4 and α < 10⁻² [131, 132, 177]

As for Property (4), scalar tensor theories are not built with the aim of introducing strong-field corrections to GR.³ Instead, they naturally lead to modifications of Einstein’s theory in the weak-field, modifications that dominate in scenarios with sufficiently weak gravitational interactions. Although this might seem strange, it is natural if one considers, for example, one of the key modifications introduced by scalar-tensor theories: the emission of dipolar gravitational radiation. Such dipolar emission dominates over the general relativistic quadrupolar emission for systems in the weak to intermediate field regime, such as in binary pulsars or in the very early inspiral of compact binaries. Therefore, one would expect scalar-tensor theories to be best constrained by experiments or observations of weakly-gravitating systems, as it has recently been explicitly shown in [465].

2.3.2 Massive graviton theories and Lorentz violation

Massive graviton theories are those in which the gravitational interaction is propagated by a massive gauge boson, i.e., a graviton with mass m _g ≠ 0 or Compton wavelength λ _g = h/(m _g c) < ∞. Einstein’s theory predicts massless gravitons and thus gravitational propagation at light speed, but if this were not the case, then a certain delay would develop between electromagnetic and gravitational signals emitted simultaneously at the source. Fierz and Pauli [169] were the first to write down an action for a free massive graviton, and ever since then, much work has gone into the construction of such models. For a detailed review, see, e.g., [232].

Gravitational theories with massive gravitons are somewhat well-motivated from a fundamental physics perspective, and thus, one can say they possess Property (2). Indeed, in loop quantum cosmology [42, 77], the cosmological extension to loop quantum gravity, the graviton dispersion relation acquires holonomy corrections during loop quantization that endow the graviton with a mass [78] m _g = Δ^−1/2γ⁻¹(ρ/ρ _c ), with γ the Barbero-Immirzi parameter, Δ the area operator, and ρ and ρ _c the total and critical energy density respectively. In string-theory-inspired effective theories, such as Dvali’s compact, extra-dimensional theory [157], such massive modes also arise

Massive graviton modes also occur in many other modified gravity models. In Rosen’s bimetric theory [365], for example, photons and gravitons follow null geodesics of different metrics [438, 435]. In Visser’s massive graviton theory [424], the graviton is given a mass at the level of the action through an effective perturbative description of gravity, at the cost of introducing a non-dynamical background metric, i.e., a prior geometry. A recent re-incarnation of this model goes by the name of bigravity, where again two metric tensors are introduced [349, 346, 219, 220]. In Bekenstein’s Tensor-Vector-Scalar (TeVeS) theory [54], the existence of a scalar and a vector field lead to subluminal gravitational-wave propagation.

Massive graviton theories have a theoretical issue, the van Dam-Veltman-Zakharov (vDVZ) discontinuity [418, 475], which is associated with Property 1.a, i.e., a GR limit. The problem is that certain predictions of massive graviton theories do not reduce to those of GR in the mg → 0 limit. This can be understood qualitatively by studying how the 5 spin states of the graviton behave in this limit. Two of them become the two GR helicity states of the massless graviton. Another two become helicity states of a massless vector that decouples from the tensor perturbations in the m _g → 0 limit. However, the last state, the scalar mode, retains a finite coupling to the trace of the stress-energy tensor in this limit. Therefore, massive graviton theories in the m _g → 0 limit do not reduce to GR, since the scalar mode does not decouple.

However, the vDVZ discontinuity can be evaded, for example, by carefully including non-linearities. Vainshtein [413, 269, 140, 45] showed that around any spherically-symmetric source of mass M, there exists a certain radius $r < {r_V} \equiv {({r_S}\lambda _g^4)^{1/5}}$, with r _S the Schwarzschild radius, where linear theory cannot be trusted. Since r _V → ∞ as m _g → 0, this implies that there is no radius at which the linear approximation (and thus vDVZ discontinuity) can be trusted. Of course, to determine then whether massive graviton theories have a continuous limit to GR, one must include non-linear corrections to the action (see also an argument by [34]), which are more difficult to uniquely predict from fundamental theory. Recently, there has been much activity in the development of new, non-linear massive gravity theories [60, 136, 211, 61, 137, 138].

Lacking a particular action for massive graviton theories that modifies the strong-field regime and is free of non-linear and radiatively-induced ghosts, it is difficult to ascertain many of its properties, but this does not prevent us from considering certain phenomenological effects. If the graviton is truly massive, whatever the action may be, two main modifications to Einstein’s theory will be introduced:

• (i) Modification to Newton’s laws;

• (ii) Modification to gravitational wave propagation.

Modifications of class (i) correspond to the replacement of the Newtonian potential by a Yukawa type potential (in the non-radiative, near-zone of any body of mass M): V = (M/r) → (M/r) exp (− r/λ _g ), where r is the distance to the body [437]. Tests of such a Yukawa interaction have been proposed through observations of bound clusters, tidal interactions between galaxies [200] and weak gravitational lensing [106], but such tests are model dependent.

Modifications of class (ii) are in the form of a non-zero graviton mass that induces a modified gravitational-wave dispersion relation. Such a modification to the dispersion relation was originally parameterized via [437]

$${{v_g^2} \over {{c^2}}} = 1 - {{m_g^2{c^4}} \over {{E^2}}},$$

(23)

where υ _g and m _g are the speed and mass of the graviton, while E is its energy, usually associated to its frequency via the quantum mechanical relation E = hf. This modified dispersion relation is inspired by special relativity, a more general version of which, inspired by quantum gravitational theories, is [316]

$${{v_g^2} \over {{c^2}}} = 1 - {\lambda ^\alpha},$$

(24)

where α is now a parameter that depends on the theory and λ represents deviations from light-speed propagation. For example, in Rosen’s bimetric theory [365], the graviton does not travel at the speed of light, but at some other speed partially determined by the prior geometry. In metric theories of gravity, $\lambda = Am_g^2{c^4}/{E^2}$, where A is some amplitude that depends on the metric theory (see discussion in [316]). Either modification to the dispersion relation has the net effect of slowing gravitons down, such that there is a difference in the time of arrival of photons and gravitons. Moreover, such an energy-dependent dispersion relation would also affect the accumulated gravitational-wave phase observed by gravitational-wave detectors, as we discuss in Section 5. Given these modifications to the dispersion relation, one would expect the generation of gravitational waves to also be greatly affected in such theories, but again, lacking a particular healthy action to consider, this topic remains today mostly unexplored.

From the structure of the above phenomenological modifications, it is clear that GR can be recovered in the m _g → 0 limit, avoiding the vDVZ issue altogether by construction. Such phenomenological modifications have been constrained by several types of experiments and observations. Using the modification to Newton’s third law and precise observations of the motion of the inner planets of the solar system together with Kepler’s third law, [437] found a bound of λ _g > 2.8 × 10¹² km. Such a constraint is purely static, as it does not sample the radiative sector of the theory. Dynamical constraints, however, do exist: through observations of the decay of the orbital period of binary pulsars, [174] found a bound of λ _g > 1.6 × 10¹⁰ km;⁴ by investigating the stability of Schwarzschild and Kerr black holes, [88] placed the constraint λ _g > 2.4 × 10¹³ km in Fierz-Pauli theory [169]. New constraints that use gravitational waves have been proposed, including measuring a difference in time of arrival of electromagnetic and gravitational waves [126, 266], as well as direct observation of gravitational waves emitted by binary pulsars (see Section 5).

Although massive gravity theories unavoidably lead to a modification to the graviton dispersion relation, the converse is not necessarily true. A modification of the dispersion relation is usually accompanied by a modification to either the Lorentz group or its action in real or momentum space. Such Lorentz-violating effects are commonly found in quantum gravitational theories, including loop quantum gravity [78] and string theory [107, 403], as well as other effective models [58, 59]. In doubly-special relativity [26, 300, 27, 28], the graviton dispersion relation is modified at high energies by modifying the law of transformation of inertial observers. Modified graviton dispersion relations have also been shown to arise in generic extra-dimensional models [381], in Hořava-Lifshitz theory [233, 234, 412, 76] and in theories with non-commutative geometries [186, 187, 188]. None of these theories necessarily requires a massive graviton, but rather the modification to the dispersion relation is introduced due to Lorentz-violating effects.

One might be concerned that the mass of the graviton and subsequent modifications to the graviton dispersion relation should be suppressed by the Planck scale. However, Collins, et al. [111, 110] have suggested that Lorentz violations in perturbative quantum field theories could be dramatically enhanced when one regularizes and renormalizes them. This is because terms that vanish upon renormalization due to Lorentz invariance do not vanish in Lorentz-violating theories, thus leading to an enhancement [185]. Whether such an enhancement is truly present cannot currently be ascertained.

2.3.3 Modified quadratic gravity

Modified quadratic gravity is a family of models first discussed in the context of black holes and gravitational waves in [473, 447]. The 4-dimensional action is given by

$$\begin{array}{*{20}c} {S \equiv \int {{d^4}x\sqrt {- g} \left\{{\kappa R + {\alpha _1}{f_1}(\vartheta){R^2} + {\alpha _2}{f_2}(\vartheta){R_{\mu \nu}}{R^{\mu \nu}} + {\alpha _3}{f_3}(\vartheta){R_{\mu \nu \delta \sigma}}{R^{\mu \nu \delta \sigma}}} \right.}}\\ {\left. {+ {\alpha _4}{f_4}(\vartheta){R_{\mu \nu \delta \sigma}}^{\ast}{R^{\mu \nu \delta \sigma}} - {\beta \over 2}[{\nabla _\mu}\vartheta {\nabla ^\mu}\vartheta + 2V(\vartheta)] + {{\mathcal L}_{{\rm{mat}}}}} \right\}.\quad \quad \quad \quad \quad \quad \quad \quad}\\ \end{array}$$

(25)

The quantity *R^μ _νδσ = (1/2) _∈δσ ^αβ R^μ _ναβ is the dual to the Riemann tensor. The quantity $\mathcal{L}_{\mathrm{mat}}$ is the external matter Lagrangian, while f _i (·) are functionals of the field ϑ, with (α _i ,β) coupling constants and κ = (16πG)−¹. Clearly, the two terms second to last in Eq. (25) represent a canonical kinetic energy term and a potential. At this stage, one might be tempted to set β = 1 or the α _i = 1 via a rescaling of the scalar field functional, but we shall not do so here.

The action in Eq. (25) is well-motivated from fundamental theories, as it contains all possible quadratic, algebraic curvature scalars with running (i.e., non-constant) couplings. The only restriction here is that all quadratic terms are assumed to couple to the same field, which need not be the case. For example, in string theory some terms might couple to the dilaton (a scalar field), while others couple to the axion (a pseudo scalar field). Nevertheless, one can recover well-known and motivated modified gravity theories in simple cases. For example, dynamical Chern-Simons modified gravity [17] is recovered when α₄ = −α _CS /4 and all other α _i = 0. Einstein-Dilaton-Gauss-Bonnet gravity [343] is obtained when α₄ = 0 and (α₁, α₂, α₃) = (1, −4, 1)α_EDGB.⁵ Both theories unavoidably arise as low-energy expansions of heterotic string theory [203, 204, 12, 89]. As such, modified quadratic gravity theories should be treated as a class of effective field theories. Moreover, dynamical Chern-Simons gravity also arises in loop quantum gravity [43, 366] when the Barbero-Immirzi parameter is promoted to a field in the presence of fermions [41, 16, 406, 311, 192].

One should make a clean and clear distinction between the theory defined by the action of Eq. (25) and that of f(R) theories. The latter are defined as functionals of the Ricci scalar only, while Eq. (25) contains terms proportional to the Ricci tensor and Riemann tensor squared. One could think of the subclass of f(R) theories with f(R) = R² as the limit of modified quadratic gravity with only α₁ = 0 and f₁(ϑ) = 1. In that very special case, one can map quadratic gravity theories and f(R) gravity to a scalar-tensor theory. Another important distinction is that f(R) theories are usually treated as exact, while the action presented above is to be interpreted as an effective theory [89] truncated to quadratic order in the curvature in a low-energy expansion of a more fundamental theory. This implies that there are cubic, quartic, etc. terms in the Riemann tensor that are not included in Eq. (25) and that presumably depend on higher powers of α _i . Thus, when studying such an effective theory one should also order-reduce the field equations and treat all quantities that depend on perturbatively, the small-coupling approximation. One can show that such an order reduction removes any additional polarization modes in propagating metric perturbations [390, 400] that naturally arise in f(R) theories. In analogy to the treatment of the Ostrogradski instability in Section 2.1, one would also expect that order-reduction would lead to a theory with a well-posed initial-value formulation.

This family of theories is usually simplified by making the assumption that coupling functions f _i (·) admit a Taylor expansion: ${f_i}(\vartheta) = {f_i}(0) + f_i^{{\prime}}(0)\vartheta + \mathcal{O}({\vartheta ^2})$ for small ϑ, where f _i (0) and $f_i^{{\prime}}(0)$ are constants and ϑ is assumed to vanish at asymptotic spatial infinity. Reabsorbing f _i (0) into the coupling constants $\alpha _i^{(0)} \equiv {\alpha _i}{f_i}(0)$ into the constants $f_i^{{\prime}}(0)$, Eq. (25) becomes S = S_GR + S₀ + S₁ with

$${S_{{\rm{GR}}}} \equiv \int {{d^4}x\sqrt {- g} \{\kappa R + {{\mathcal L}_{{\rm{mat}}}}\} ,}$$

(26a)

$${S_{\rm{0}}} \equiv \int {{d^4}x\sqrt {- g} \left\{{\alpha _1^{(0)}{R^2} + \alpha _2^{(0)}{R_{\mu \nu}}{R^{\mu \nu}} + \alpha _3^{(0)}{R_{\mu \nu \delta \sigma}}{R^{\mu \nu \delta \sigma}}} \right\},}$$

(26b)

$$\begin{array}{*{20}c} {{S_{\rm{1}}} \equiv \int {{d^4}x\sqrt {- g} \left\{{\alpha _1^{(1)}\vartheta {R^2} + \alpha _2^{(1)}\vartheta {R_{\mu \nu}}{R^{\mu \nu}} + \alpha _3^{(1)}\vartheta {R_{\mu \nu \delta \sigma}}{R^{\mu \nu \delta \sigma}}} \right.}}\\ {\left. {+ \alpha _4^{(1)}\vartheta {R_{\mu \nu \delta \sigma}}^{\ast}{R^{\mu \nu \delta \sigma}} - {\beta \over 2}[{\nabla _\mu}\vartheta {\nabla ^\mu}\vartheta + 2V(\vartheta)]} \right\}.\quad \quad \quad}\\ \end{array}$$

(26c)

Here, S_GR is the Einstein-Hilbert plus matter action, while S₀ and S₁ are corrections. The former is decoupled from υ, where the omitted term proportional to $\alpha _4^{(0)}$ does not affect the classical field equations since it is topological, i.e., it can be rewritten as the total 4-divergence of some 4-current. Similarly, if the $\alpha _i^{(0)}$ were chosen to reconstruct the Gauss-Bonnet invariant, $(\alpha _1^{(0)},\alpha _2^{(0)},\alpha _3^{(0)}) = (1, - 4,1){\alpha _{{\rm{GB}}}}$, then this combination would also be topological and not affect the classical field equations. On the other hand, S₁ is a modification to GR with a direct (non-minimal) coupling to ϑ such that as the field goes to zero, the modified theory reduces to GR.

Another restriction one usually makes to simplify modified gravity theories is to neglect the $\alpha _i^{(0)}$ terms and only consider the S₁ modification, the restricted modified quadratic gravity. The $\alpha _i^{(0)}$ terms represent corrections that are non-dynamical. The term proportional to $\alpha _1^{(0)}$ resembles a certain class of f(R) theories. As such, it can be mapped to a scalar-tensor theory with a complicated potential, which has been heavily constrained by torsion-balance Eöt-Wash experiments to $\alpha _1^{(0)} < 2 \times {10^{- 8}}{{\rm{m}}^2}$ [237, 259, 62]. Moreover, these theories have a fixed coupling constant that does not run with energy or scale. In restricted modified gravity, the scalar field is effectively forcing the running of the coupling.

Then, let us concentrate on restricted modified quadratic gravity and drop the superscript in $\alpha _i^{(1)}$. The modified field equations are

$$\begin{array}{*{20}c} {{G_{\mu \nu}} + {{{\alpha _1}\vartheta} \over \kappa}{\mathcal H}_{\mu \nu}^{(0)} + {{{\alpha _2}\vartheta} \over \kappa}{\mathcal I}_{\mu \nu}^{(0)} + {{{\alpha _3}\vartheta} \over \kappa}{\mathcal J}_{\mu \nu}^{(0)} + {{{\alpha _1}} \over \kappa}{\mathcal H}_{\mu \nu}^{(1)} + {{{\alpha _2}} \over \kappa}{\mathcal I}_{\mu \nu}^{(1)} + {{{\alpha _3}} \over \kappa}{\mathcal J}_{\mu \nu}^{(1)} + {{{\alpha _4}} \over \kappa}{\mathcal K}_{\mu \nu}^{(1)}}\\ {= {1 \over {2\kappa}}\left({T_{\mu \nu}^{{\rm{mat}}} + T_{\mu \nu}^{(\vartheta)}} \right),\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\\end{array}$$

(27)

where we have defined

$${\mathcal H}_{\mu \nu}^{(0)} \equiv 2R{R_{\mu \nu}} - {1 \over 2}{g_{\mu \nu}}{R^2} - 2{\nabla _{\mu \nu}}R + 2{g_{\mu \nu}}\square R,$$

(28a)

$${\mathcal I}_{\mu \nu}^{(0)} \equiv \square{R_{\mu \nu}} + 2{R_{\mu \delta \nu \sigma}}{R^{\delta \sigma}} - {1 \over 2}{g_{\mu \nu}}{R^{\delta \sigma}}{R_{\delta \sigma}} + {1 \over 2}{g_{\mu \nu}}\square R - {\nabla _{\mu \nu}}R,$$

(28b)

$${\mathcal J}_{\mu \nu}^{(0)} \equiv 8{R^{\delta \sigma}}{R_{\mu \delta \nu \sigma}} - 2{g_{\mu \nu}}{R^{\delta \sigma}}{R_{\delta \sigma}} + 4\square {R_{\mu \nu}} - 2R\;{R_{\mu \nu}} + {1 \over 2}{g_{\mu \nu}}{R^2} - 2{\nabla _{\mu \nu}}R,$$

(28c)

$${{\mathcal H}_{\mu \nu}^{(1)} \equiv - 4({\nabla _{(\mu}}\vartheta){\nabla _{\nu)}}R - 2R{\nabla _{\mu \nu}}\vartheta + {g_{\mu \nu}}[2R\square\vartheta + 4({\nabla ^\delta}\vartheta){\nabla _\delta}R],}$$

(28d)

$$\begin{array}{*{20}c} {{\mathcal J}_{\mu \nu}^{(1)} \equiv - ({\nabla _{(\mu}}\vartheta){\nabla _{\nu)}}R - 2{\nabla ^\delta}\vartheta {\nabla _{(\mu}}{R_{\nu)\delta}} + 2{\nabla ^\delta}\vartheta {\nabla _\delta}{R_{\mu \nu}} + {R_{\mu \nu}}\square\vartheta \,\quad} \\ {- 2{R_{\delta (\mu}}{\nabla ^\delta}{\nabla _{\nu)}}\vartheta + {g_{\mu \nu}}\left({{\nabla ^\delta}\vartheta {\nabla _\delta}R + {R^{\delta \sigma}}{\nabla _{\delta \sigma}}\vartheta} \right),\quad \,\,\,} \\ \end{array}$$

(28e)

$${\mathcal J}_{\mu \nu}^{(1)} \equiv - 8({\nabla ^\delta}\vartheta)\left({{\nabla _{(\mu}}{R_{\nu)\delta}} - {\nabla _\delta}{R_{\mu \nu}}} \right) + 4{R_{\mu \delta \nu \sigma}}{\nabla ^{\delta \sigma}}\vartheta ,$$

(28f)

$${\mathcal K}_{\mu \nu}^{(1)} \equiv - 4({\nabla ^\delta}\vartheta){\epsilon _{\delta \sigma \chi (\mu}}{\nabla ^\chi}R_{\nu)}^\sigma + 4{({\nabla _{\delta \sigma}}\vartheta)^{\ast}}{R_{(\mu}}{\delta _{\nu)}}\sigma .$$

(28g)

The ϑ stress-energy tensor is

$$T_{\mu \nu}^{(\vartheta)} = \beta \left[ {({\nabla _\mu}\vartheta)({\nabla _\nu}\vartheta) - {1 \over 2}{g_{\mu \nu}}\left({{\nabla _\delta}\vartheta {\nabla ^\delta}\vartheta - 2V(\vartheta)} \right)} \right].$$

(29)

The field equations for the scalar field are

$$\beta \square\vartheta - \beta {{dV} \over {d\vartheta}} = - {\alpha _1}{R^2} - {\alpha _2}{R_{\mu \nu}}{R^{\mu \nu}} - {\alpha _3}{R_{\mu \nu \delta \sigma}}{R^{\mu \nu \delta \sigma}} - {\alpha _4}{R_{\mu \nu \delta \sigma}}^{\ast}{R^{\mu \nu \delta \sigma}}.$$

(30)

Notice that unlike traditional scalar-tensor theories, the scalar field is here sourced by the geometry and not by the matter distribution. This directly implies that black holes in such theories are likely to be hairy.

From the structure of the above equations, it should be clear that the dynamics of ϑ guarantee that the modified field equations are covariantly conserved exactly. That is, one can easily verify that the covariant divergence of Eq. (27) identically vanishes upon imposition of Eq. (30). Such a result had to be so, as the action is diffeomorphism invariant. If one neglected the kinetic and potential energies of ϑ in the action, as was originally done in [245], the theory would possess preferred-frame effects and would not be covariantly conserved. Moreover, such a theory requires an additional constraint, i.e., the right-hand side of (30) would have to vanish, which is an unphysical consequence of treating ϑ as prior structure [470, 207].

One last simplification that is usually made when studying modified quadratic gravity theories is to ignore the potential V(ϑ), i.e., set V(ϑ) = 0. This potential can in principle be non-zero, for example if one wishes to endow ϑ with a mass or if one wishes to introduce a cosine driving term, like that for axions in field and string theory. However, reasons exist to restrict the functional form of such a potential. First, a mass for ϑ will modify the evolution of any gravitational degree of freedom only if this mass is comparable to the inverse length scale of the problem under consideration (such as a binary system). This could be possible if there is an incredibly large number of fields with different masses in the theory, such as perhaps in the string axiverse picture [40, 268, 303]. However, in that picture the moduli fields are endowed with a mass due to shift-symmetry breaking by non-perturbative effects; such masses are not expected to be comparable to the inverse length scale of binary systems. Second, no mass term may appear in a theory with a shift symmetry, i.e., invariance under the transformation ϑ → ϑ + const. Such symmetries are common in four-dimensional, low-energy, effective string theories [79, 204, 203, 92, 89], such as dynamical Chern—Simons and Einstein-Dilaton-Gauss-Bonnet theory. Similar considerations apply to other more complicated potentials, such as a cosine term.

Given these field equations, one can linearize them about Minkowski space to find evolution equations for the perturbation in the small-coupling approximation. Doing so, one finds [447]

$$\begin{array}{*{20}c} {{\square_\eta}\vartheta = - {{{\alpha _1}} \over \beta}{{\left({{1 \over {2\kappa}}} \right)}^2}T_{{\rm{mat}}}^2 - {{{\alpha _2}} \over \beta}{{\left({{1 \over {2\kappa}}} \right)}^2}T_{{\rm{mat}}}^{\mu \nu}T_{\mu \nu}^{{\rm{mat}}}\quad \quad} \\ {- {{2{\alpha _3}} \over \beta}({h_{\alpha \beta ,\mu \nu}}{h^{\alpha [\beta ,\mu ]\nu}} + {h_{\alpha \beta ,\mu \nu}}{h^{\mu [\nu ,\alpha ]\beta}})} \\ {- {{2{\alpha _4}} \over \beta}{\epsilon ^{- \alpha \beta \mu \nu}}{h_{\alpha \delta ,\gamma \beta}}{h_\nu}^{[\gamma ,\delta ]}\mu ,\quad \quad \quad \quad} \\ \end{array}$$

(31)

where we have order-reduced the theory where possible and used the harmonic gauge condition (which is preserved in this class of theories [390, 400]). The corresponding equation for the metric perturbation is rather lengthy and can be found in Eqs. (17)–(24) in [447]. Since these theories are to be considered effective, working always to leading order in α _i , one can show that they are perturbatively of type N₂ in the E(2) classification [161], i.e., in the far zone, the only propagating modes that survive are the two transverse-traceless (spin-2) metric perturbations [390]. However, in the strong-field region it is possible that additional modes are excited, although they decay rapidly as they propagate to future null infinity.

Lastly, let us discuss what is known about whether modified quadratic gravities satisfy the requirements discussed in Section 2.1. As it should be clear from the action itself, this modified gravity theory satisfies the fundamental requirement, i.e., passing all precision tests, provided the couplings α _i are sufficiently small. This is because such theories have a continuous limit to GR as α _i → 0.⁶ Dynamical Chern—Simons gravity is constrained only weakly at the moment, $\xi _4^{1/4} < {10^8}$ km, where ${\xi _4} \equiv \alpha _4^2/(\beta \kappa)$, only through observations of Lense—Thirring precession in the solar system [19]. The Einstein-Dilaton-Gauss—Bonnet gravity coupling constant ${\xi _3} \equiv \alpha _3^2/(\beta \kappa)$, on the other hand, has been constrained by several experiments: solar system observations of the Shapiro time delay with the Cassini spacecraft placed the bound $\xi _3^{1/4} < 1.3 \times {10^7}$ km [73, 29]; the requirement that neutron stars still exist in this theory placed the constraint $\xi _3^{1/4} \underset{\sim}{<}26$ km [342], with the details depending somewhat on the central density of the neutron star; observations of the rate of change of the orbital period in the low-mass X-ray binary A0620−00 [358, 255] has led to the constraint $\xi _3^{1/4} < 1.9$ km [445].

However, not all sub-properties of the fundamental requirement are known to be satisfied. One can show that certain members of modified quadratic gravity possess known solutions and these are stable, at least in the small-coupling approximation. For example, in dynamical Chern—Simons gravity, spherically-symmetric vacuum solutions are given by the Schwarzschild metric with constant ϑ to all orders in α _i [245, 470]. Such a solution is stable to small perturbations [319, 190], as also are non-spinning black holes and branes in anti de Sitter space [144]. On the other hand, spinning solutions continue to be elusive, with approximate solutions in the slow-rotation/small-coupling limit known both for black holes [466, 272, 345, 455] and stars [19, 342]; nothing is currently known about the stability of these spinning solutions. In Einstein-Dilaton-Gauss—Bonnet theory even spherically-symmetric solutions are modified [473, 345] and these are stable to axial perturbations [343].

The study of modified quadratic gravity theories as effective theories is valid provided one is sufficiently far from its cut-off scale, i.e., the scale beyond which higher-order curvature terms cannot be neglected anymore. One can estimate the magnitude of this scale by studying the size of loop corrections to the quadratic curvature terms in the action due to n-point interactions [455]. Simple counting requires that the number of scalar and graviton propagators, P _s and P _g , satisfy the following relation in terms of the number of vertices V:

$${P_s} = {V \over 2},\quad {P_g} = (n - 1){V \over 2}.$$

(32)

Thus, loop corrections are suppressed by factors of $\alpha _i^VM_{{\rm{pl}}}^{(2 - n)V}{\Lambda ^{nV}}$, with M_p1 the Planck mass and Λ the energy scale introduced by dimensional arguments. The cut-off scale above which the theory cannot be treated as an effective one can be approximated as the value of Λ at which the suppression factor becomes equal to unity:

$${\Lambda _c} \equiv M_{{\rm{pl}}}^{1 - 2/n}\alpha _i^{1/n},$$

(33)

This cut-off scale automatically places a constraint on the magnitude of α _i above which higher-curvature corrections must be included. Setting the largest value of Λ _c to be equal to $\mathcal{O}(10\mu\mathrm{m})$, thus saturating bounds from table-top experiments [259], and solving for α _i , we find

$$\alpha _i^{1/2} < {\mathcal O}({10^8}\;{\rm{km}}).$$

(34)

Current solar system bounds on α _i already require the coupling constant to be smaller than 10⁸ km, thus justifying the treatment of these theories as effective models.

As for the other requirements discussed in Section 2.1, it is clear that modified quadratic gravity is well-motivated from fundamental theory, but it is not clear at all whether it has a well-posed initial-value formulation. From an effective point of view, a perturbative treatment in α _i naturally leads to stable solutions and a well-posed initial-value problem, but this is probably not the case when it is treated as an exact theory. In fact, if one were to treat such a theory as exact (to all orders in α _i ), then the evolution system would likely not be hyperbolic, as higher-than-second time derivatives now drive the evolution. Although no proof exists, it is likely that such an exact theory is not well-posed as an initial-value problem. Notice, however, that this says nothing about the fundamental theories that modified quadratic gravity derives from. This is because even if the truncated theory were ill posed, higher-order corrections that are neglected in the truncated version could restore well-posedness.

As for the last requirement (that the theory modifies the strong field), modified quadratic theories are ideal in this respect. This is because they introduce corrections to the action that depend on higher powers of the curvature. In the strong-field, such higher powers could potentially become non-negligible relative to the Einstein—Hilbert action. Moreover, since the curvature scales inversely with the mass of the objects under consideration, one expects the largest deviations in systems with small total mass, such as stellar-mass black-hole mergers. On the other hand, deviations from GR should be small for small compact objects spiraling into a supermassive black hole, since here the spacetime curvature is dominated by the large object, and thus it is small, as discussed in [390].

2.3.4 Variable G theories and large extra dimensions

Variable G theories are defined as those where Newton’s gravitational constant is promoted to a spacetime function. Such a modification breaks the principle of equivalence (see [438]) because the laws of physics now become local position dependent. In turn, this implies that experimental results now depend on the spacetime position of the laboratory frame at the time of the experiment.

Many known alternative theories that violate the principle of equivalence, and in particular, the strong equivalence principle, predict a varying gravitational constant. A classic example is scalar-tensor theory [435], which, as explained in Section 2.3.1, modifies the gravitational sector of the action by multiplying the Ricci scalar by a scalar field (in the Jordan frame). In such theories, one can effectively think of the scalar as promoting the coupling between gravity and matter to a field-dependent quantity G → G(ϕ), thus violating local position invariance when ϕ varies. Another example are bimetric theories, such as that of Lightman—Lee [293], where the gravitational constant becomes time-dependent even in the absence of matter, due to possibly time-dependent cosmological evolution of the prior geometry. A final example are higher-dimensional, brane-world scenarios, where enhanced Hawking radiation inexorably leads to a time-varying effective 4D gravitational constant [141], whose rate of change depends on the curvature radius of extra dimensions [255].

One can also construct f(R)-type actions that introduce variability to Newton’s constant. For example, consider the f(R) model [180]

$$S = \int {{d^4}x\sqrt {- g} \;\kappa R} \left[ {1 + {\alpha _0}\ln \left({{R \over {{R_0}}}} \right)} \right] + {S_{{\rm{mat}}}},$$

(35)

where κ = (16πG)⁻¹, α₀ is a coupling constant and R₀ is a curvature scale. This action is motivated by certain renormalization group flow arguments [180]. The field equations are

$${G_{\mu \nu}} = {1 \over {2\bar \kappa}}T_{\mu \nu}^{{\rm{mat}}} - {{{\alpha _0}} \over {\bar \kappa}}{R_{\mu \nu}} - 2{\kappa \over {\bar \kappa}}{{{\alpha _0}} \over {{R^2}}}{\nabla _{(\mu}}R{\nabla _{\nu)}}R - {1 \over 2}{{{\alpha _0}\kappa} \over {\bar \kappa}}{g_{\mu \nu}}\square R,$$

(36)

where we have defined the new constant

$$\bar \kappa : = \kappa \left[ {1 + {{{\alpha _0}} \over \kappa}\ln \left({{R \over {{R_0}}}} \right)} \right].$$

(37)

Clearly, the new coupling constant $\bar \kappa$ depends on the curvature scale involved in the problem, and thus, on the geometry, forcing G to run to zero in the ultraviolet limit.

An important point to address is whether variable G theories can lead to modifications to a vacuum spacetime, such as a black-hole—binary inspiral. In Einstein’s theory, G appears as the coupling constant between geometry, encoded by the Einstein tensor g _μν , and matter, encoded by the stress energy tensor $T_{\mu \nu}^{{\rm{mat}}}$. When considering vacuum spacetimes, $T_{\mu \nu}^{{\rm{mat}}} = 0$ and one might naively conclude that a variable G would not introduce any modification to such spacetimes. In fact, this is the case in scalar-tensor theories (without homogeneous, cosmological solutions to the scalar field equation), where the no-hair theorem establishes that black-hole solutions are not modified. On the other hand, scalar-tensor theories with a cosmological, homogeneous scalar field solution can violate the no-hair theorem, endowing black holes with time-dependent hair, which in turn would introduce variability into G even in vacuum spacetimes [246, 236, 67].

In general, Newton’s constant plays a much more fundamental role than merely a coupling constant: it defines the relationship between energy and length. For example, for the vacuum Schwarzschild solution, G establishes the relationship between the radius R of the black hole and the rest-mass energy E of the spacetime via R = 2 GE/c⁴. Similarly, in a black-hole—binary spacetime, each black hole introduces an energy scale into the problem that is quantified by a specification of Newton’s constant. Therefore, one can treat variable G modifications as induced by some effective theory that modifies the mapping between the curvature scale and the energy scale of the problem, as is done for example in theories with extra dimensions.

An explicit example of this idea is realized in braneworld models. Superstring theory suggests that physics should be described by 4 large dimensions, plus another 6 that are compactified and very small [354, 355]. The size of these extra dimensions is greatly constrained by particle theory experiments. However, braneworld models, where a certain higher-dimensional membrane is embedded in a higher-dimensional bulk spacetime, can evade this constraint as only gravitons can interact with the bulk. The ADD model [32, 33] is a particular example of such a braneworld, where the bulk is flat and compact and the brane is tensionless with ordinary fields localized on it. Since gravitational-wave experiments have not yet constrained deviations from Einstein’s theory in the strong field, the size of these extra dimensions is constrained to micrometer scales only by table-top experiments [259, 7].

What is relevant to gravitational-wave experiments is that in many of these braneworld models black holes may not remain static [163, 405]. The argument goes roughly as follows: a five-dimensional black hole is dual to a four-dimensional one with conformal fields on it by the ADS/CFT conjecture [301, 9], but since the latter must evolve via Hawking radiation, the black hole must be losing mass. The Hawking mass loss rate is here enhanced by the large number of degrees of freedom in the conformal field theory, leading to an effective modification to Newton’s laws and to the emission of gravitational radiation. Effectively, one can think of the black-hole mass loss as due to the black hole being stretched away from the brane into the bulk due to a universal acceleration, that essentially reduces the size of the brane-localized black hole. For black-hole binaries, one can then draw an analogy between this induced time dependence in the black-hole mass and a variable G theory, where Newton’s constant decays due to the presence of black holes. Of course, this is only analogy, since large extra dimensions would not predict a time-evolving mass in neutron-star binaries.

Recently, however, Figueras et al. [170, 172, 171] numerically found stable solutions that do not require a radiation component. If such solutions were the ones realized in nature as a result of gravitational collapse on the brane, then the black hole mass would be time independent, up to quantum correction due to Hawking evaporation, a negligible effect for realistic astrophysical systems. Unfortunately, we currently lack numerical simulations of the dynamics of gravitational collapse in such scenarios.

Many experiments have been carried out to measure possible deviations from a constant G value, and they can broadly be classified into two groups: (a) those that search for the present or nearly present rate of variation (at redshifts close to zero); (b) those that search for secular variations over long time periods (at very large redshifts). Examples of experiments or observations of the first class include planetary radar ranging [350], surface temperature observations of low-redshift millisecond pulsars [249, 362], lunar ranging observations [442] and pulsar timing observations [260, 143], the latter two being the most stringent. Examples of experiments of the second class include the evolution of the sun [208] and Big-Bang Nucleosynthesis (BBN) calculations [119, 47], again with the latter being more stringent. For either class, the strongest constraints are about Ġ/G ≲ 10⁻¹³ yr⁻¹, varying somewhat from experiment to experiment.

Lacking a particularly compelling action to describe variable G theories, one is usually left with a phenomenological model of how such a modification to Einstein’s theory would impact gravitational waves. Given that the part of the waveform that detectors are most sensitive to is the gravitational wave phase, one can model the effect of variable G theories by studying how the rate of change of its frequency would be modified. Assuming a Taylor expansion for Newton’s constant one can derive the modification to the evolution equation for the gravitational wave frequency, given whichever physical scenario one is considering. Solving such an evolution equation then leads to a modification in the accumulated gravitational-wave phase observed at detectors on Earth. In Section 5 we will provide an explicit example of this for a compact binary system.

Let us discuss whether such theories satisfy the criteria defined in Section 2.1. The fundamental property can be satisfied if the rate of change of Newton’s constant is small enough, as variable G theories usually have a continuous limit to GR (as all derivatives of G go to zero). Whether variable theories are well-motivated from fundamental physics (Property 2) depends somewhat on the particular effective model or action that one considers. But in general, Property 2 is usually satisfied, considering that such variability naturally arises in theories with extra dimensions, and the latter are also natural in all string theories. However, variable theories usually fail at introducing modifications in the strong-field region. Usually, such variability is parameterized as a Taylor expansion about some initial point with constant coefficients. That is, the variability of is not usually constructed so as to become stronger closer to merger.

The well-posed property and the sub-properties of the fundamental property depend somewhat on the particular effective theory used to describe varying G modifications. In the f(R) case, one can impose restrictions on the functional form f(·) such that no ghosts (f′ > 0) or instabilities (f″ > 0) arise [180]. This, of course, does not guarantee that this (or any other such) theory is well posed. A much more detailed analysis would be required to prove well-posedness of the class of theories that lead to a variable Newton’s constant, but such is currently lacking.

2.3.5 Non-commutative geometry

Non-commutative geometry is a gravitational theory that generalizes the continuum Riemannian manifold of Einstein’s theory with the product of it with a tiny, discrete, finite non-commutative space, composed of only two points. Although the non-commutative space has zero spacetime dimension, as the product manifold remains four dimensional, its internal dimensions are 6 to account for Weyl and chiral fermions. This space is discrete to avoid the infinite tower of massive particles that would otherwise be generated, as in string theory. Through this construction, one can recover the standard model of elementary particles, while accounting for all (elementary particle) experimental data to date. Of course, the simple non-commutative space described above is expected to be replaced by a more complex model at Planckian energies. Thus, one is expected to treat such non-commutative geometry models as effective theories. Essentially nothing is currently known about the full non-commutative theory of which the theories described in this section are an effective low-energy limit.

Before proceeding with an action-principle description of non-commutative geometry theories, we must distinguish between the spectral geometry approach championed by Connes [114], and Moyal-type non-commutative geometries [389, 206, 322]. In the former, the manifold is promoted to a non-commutative object through the product of a Riemann manifold with a non-commutative space. In the latter, instead, a non-trivial set of commutation relations is imposed between operators corresponding to position. These two theories are in principle unrelated. In this review, we will concentrate only on the former, as it is the only type of non-commutative GR extension that has been studied in the context of gravitational-wave theory.

The effective action for spectral non-commutative geometry theories (henceforth, non-commutative geometries for short) is

$$S = \int {{d^4}x\sqrt {- g} (\kappa R + {\alpha _0}{C_{\mu \nu \delta \sigma}}{C^{\mu \nu \delta \sigma}} + {\tau _0}{R^{\ast}}{R^{\ast}} - {\xi _0}R\vert H\vert ^{2}) + {S_{{\rm{mat}}}},}$$

(38)

where H is related to the Higgs field, C _μνδσ is the Weyl tensor, (α₀, τ₀, ξ₀) are couplings constants and we have defined the quantity

$${R^{\ast}}{R^{\ast}}: = {1 \over 4}{\epsilon ^{\mu \nu \rho \sigma}}{\epsilon _{\alpha \beta \gamma \delta}}{R_{\mu \nu}}^{\alpha \beta}{R_{\rho \sigma}}^{\gamma \delta}.$$

(39)

Notice that this term integrates to the Euler characteristic, and since τ₀ is a constant, it is topological and does not affect the classical field equations. The last term of Eq. (38) is usually ignored as H is assumed to be relevant only in the early universe. Finally, the second term can be rewritten in terms of the Riemann and Ricci tensors as

$${C_{\mu \nu \delta \sigma}}{C^{\mu \nu \delta \sigma}} = {1 \over 3}{R^2} - 2{R_{\mu \nu}}{R^{\mu \nu}} + {R_{\mu \nu \delta \sigma}}{R^{\mu \nu \delta \sigma}}.$$

(40)

Notice that this corresponds to the modified quadratic gravity action of Eq. (26) with all $\alpha _i^{(1)} = 0$ and $\alpha _1^{(0)},\alpha _2^{(0)},\alpha _3^{(0)} = (1/3, - 2,1)$, which is not the Gauss—Bonnet invariant. Notice also that this model is not usually studied in modified quadratic gravity theory, as one usually concentrates on the terms that have an explicit scalar field coupling.

The field equations of this theory can be read directly from Eq. (27), but we repeat them here for completeness:

$${G_{\mu \nu}} - {{2{\alpha _0}} \over \kappa}[2{\nabla ^{\kappa \lambda}} + {R^{\lambda \kappa}}]\;{C_{\mu \lambda \nu \kappa}} = {1 \over {2\kappa}}T_{\mu \nu}^{{\rm{mat}}}.$$

(41)

One could in principle rewrite this in terms of the Riemann and Ricci tensors, but the expressions become quite complicated, as calculated explicitly in Eqs. (2) and (3) of [473]. Due to the absence of a dynamical degree of freedom coupling to the modifications to the Einstein—Hilbert action, this theory is not covariantly conserved in vacuum. By this we mean that the covariant divergence of Eq. (41) does not vanish in vacuum, thus violating the weak-equivalence principle and leading to additional equations that might over-constrain the system. In the presence of matter, the equations of motion will not be given by the vanishing of the covariant divergence of the matter stress-energy alone, but now there will be additional geometric terms.

Given these field equations, one can linearize them about a flat background to find the evolution equations for the metric perturbation [326, 325]

$$(1 - {\beta ^{- 2}}{\square_\eta})\;{\square_\eta}{h_{\mu \nu}} = - 16\pi T_{\mu \nu}^{{\rm{mat}}},$$

(42)

where the term proportional to β² = (−32πα₀)⁻¹ acts like a mass term. Here, one has imposed the transverse-traceless gauge (a refinement of Lorenz gauge), which can be shown to exist [326, 325]. Clearly, even though the full non-linear equations are not covariantly conserved, its linearized version is, as one can easily show that the divergence of the left-hand side of Eq. (42) vanishes. Because of these features, if one works perturbatively in β⁻¹, then such a theory will only possess the two usual transverse-traceless (spin-2) polarization modes, i.e., it is perturbatively of type N₂ in the E(2) classification [161].

Let us now discuss whether such a theory satisfies the properties discussed in Section 2.1. Non-commutative geometry theories clearly possess the fundamental property, as one can always take α₀ → 0 (or equivalently β⁻² → 0) to recover GR. Therefore, there must exist a sufficiently small α₀ such that all precision tests carried out to date are satisfied. As for the existence and stability of known solutions, [326, 325] have shown that Minkowski spacetime is stable only for α₀ < 0, as otherwise a tachyonic term appears in the evolution of the metric perturbation, as can be seen from Eq. (42). This then automatically implies that β must be real.

Current constraints on Weyl terms of this form come mostly from solar system experiments. Ni [328] recently studied an action of the form of Eq. (38) minimally coupled to matter in light of solar system experiments. He calculated the relativistic Shapiro time-delay and light deflection about a massive body in this theory and found that observations of the Cassini satellite place constraints on |α₀|^1/2 < 5.7 km [328]. This is currently the strongest bound we are aware of on α₀.

Many solutions of GR are preserved in non-commutative geometries. Regarding black holes, all solutions that are Ricci flat (vacuum solutions of the Einstein equations) are also solutions to Eq. (41). This is because by the second Bianchi identity, one can show that

$${\nabla ^{\kappa \lambda}}{R_{\mu \lambda \nu \kappa}} = {\nabla ^\kappa}_\nu {R_{\mu \kappa}} - \square{R_{\mu \nu}},$$

(43)

and the right-hand side vanishes in vacuum, forcing the entire left-hand side of Eq. (41) to vanish. However, this is not so for neutron stars where the equations of motion are likely to be modified, unless they are static [324]. Moreover, as of now there has been no stability analysis of black-hole or stellar solutions and no study of whether the theory is well posed as an initial-value problem, even as an effective theory. Thus, except for the fundamental property, it is not clear that non-commutative geometries satisfy any of the other criteria listed in Section 2.1.

2.3.6 Gravitational parity violation

Parity, the symmetry transformation that flips the sign of the spatial triad, has been found to be broken in the standard model of elementary interactions. Only the combination of charge conjugation, parity transformation and time inversion (CPT) still remains a true symmetry of the standard model. Experimentally, it is curious that the weak interaction exhibits maximal parity violation, while other fundamental forces seem to not exhibit any. Theoretically, parity violation unavoidably arises in the standard model [55, 8, 21], as there exist one-loop chiral anomalies that give rise to parity-violating terms coupled to lepton number [428]. In certain sectors of string theory, such as in heterotic and Type I superstring theories, parity violation terms are also generated through the Green—Schwarz gauge anomaly-canceling mechanism [204, 355, 12]. Finally, in loop quantum gravity [41], the scalarization of the Barbero—Immirzi parameter coupled to fermions leads to an effective action that contains parity-violating terms [406, 90, 311, 192]. Even without a particular theoretical model, one can show that effective field theories of inflation generically contain non-vanishing, second-order, parity-violating curvature corrections to the Einstein—Hilbert action [429]. Alternatively, phenomenological parity-violating extensions of GR have been proposed through a scalarization of the fundamental constants of nature [115].

One is then naturally led to ask whether the gravitational interaction is parity invariant in the strong field. A violation of parity invariance would occur if the Einstein—Hilbert action were modified through a term that involved a Levi-Civita tensor and parity invariant tensors or scalars. Let us try to construct such terms with only single powers of the Riemann tensor and a single scalar field ϑ:

$$\begin{array}{*{20}c} {({\rm{ia}})\,{R_{\alpha \beta \gamma \delta}}\,{\epsilon ^{\alpha \beta \gamma \delta}},\,\,\quad \quad \,\,} & {({\rm{ib}})\,{R_{\alpha \beta \gamma \mu}}\,{\epsilon ^{\alpha \beta \gamma \nu}}{\nabla ^\mu}_\nu \vartheta ,\,\quad} \\ {({\rm{ic}})\,{R_{\alpha \beta \gamma \mu}}\,{\epsilon ^{\alpha \beta \delta \nu}}\,{\nabla ^{\mu \gamma}}_{\nu \delta}\vartheta ,} & {({\rm{id}})\,{R_{\alpha \zeta \gamma \mu}}\,{\epsilon ^{\alpha \beta \delta \nu}}\nabla _{\,\,\,\,\,\,\beta \nu{\delta ^\zeta}\vartheta}^{\mu \gamma}.} \\ \end{array}$$

Option (ia) and (ib) vanish by the Bianchi identities. Options (ic) and (id) include the commutator of covariant derivatives, which can be rewritten in terms of a Riemann tensor, and thus it leads to terms that are at least quadratic in the Riemann tensor. Therefore, no scalar can be constructed that includes contractions with the Levi-Civita tensor from a single Riemann curvature tensor and a single field. One can try to construct a scalar from the Ricci tensor

$$({\rm{iia}})\;{R_{\alpha \beta}}\;{\epsilon^{\alpha \beta \gamma \delta}}\;{\nabla _{\gamma \delta}}\vartheta ,\quad ({\rm{iib}})\;{R_{\alpha \beta}}\;{\epsilon^{\alpha \mu \gamma \delta}}\;{\nabla _{\gamma \delta \mu}}^\beta \vartheta ,$$

(44)

but again (iia) vanishes by the symmetries of the Ricci tensor, while (iib) involves the commutator of covariant derivatives, which introduces another power of the curvature tensor. Obviously, the only term one can write with the Ricci scalar would lead to a double commutator of covariant derivatives, leading to extra factors of the curvature tensor.

One is then forced to consider either theories with two mutually-independent fields or theories with quadratic curvature tensors. Of the latter, the only combination that can be constructed and that does not vanish by the Bianchi identities is the Pontryagin density, i.e., R* R, and therefore, the action [245, 17]

$$S = \int {{d^4}x\sqrt {- g} \left({\kappa R + {\alpha \over 4}\vartheta \;{R^{\ast}}R} \right),}$$

(45)

is the most general, quadratic action with a single scalar field that violates parity invariance, where we have rescaled the α prefactor to follow historical conventions. This action defines non-dynamical Chern—Simons modified gravity, initially proposed by Jackiw and Pi [245, 17]. Notice that this is the same as the term proportional to α₄ in the quadratic gravity action of Eq. (26), except that here ϑ is prior geometry, i.e., it does not possess self-consistent dynamics or an evolution equation. Such a term violates parity invariance because the Pontryagin density is a pseudo-scalar, while ϑ is assumed to be a scalar.

The field equations for this theory are⁷

$${G_{\mu \nu}} + {\alpha \over {4\kappa}}{\mathcal K}_{\mu \nu}^{(1)} = {1 \over {2\kappa}}T_{\mu \nu}^{{\rm{mat}}},$$

(46)

which is simply Eq. (27) with (α₁, α₂, α₃) set to zero and no stress-energy for ϑ. Clearly, these field equations are not covariantly conserved in vacuum, i.e., taking the covariant divergence one finds the constraint

$$\alpha {R^{\ast}}R = 0.$$

(47)

This constraint restricts the space of allowed solutions, for example disallowing the Kerr metric [207]. Therefore, it might seem that the evolution equations for the metric are now overconstrained, given that the field equations provide 10 differential conditions for the 10 independent components of the metric tensor, while the constraint adds one additional, independent differential condition. Moreover, unless the Pontryagin constraint, Eq. (47), is satisfied, matter fields will not evolve according to ${\nabla ^\mu}T_{\mu \nu}^{{\rm{mat}}} = 0$, thus violating the equivalence principle.

From the field equations, we can derive an evolution equation for the metric perturbation when linearizing about a flat background, namely

$${\square_\eta}{h_{\mu \nu}} + {\alpha \over \kappa}({\vartheta _{,\gamma}}{\epsilon_{(\mu}}^{\gamma \delta \chi}{\square_\eta}{h_{\nu)\delta ,\chi}} - {\vartheta _{,\gamma}}^\zeta {\epsilon_{(\mu}}^{\gamma \delta \chi}{h_{\vert \delta \zeta \vert ,\nu)\chi}} + {\vartheta _{,\gamma}}^\zeta {\epsilon_{(\mu}}^{\gamma \delta \chi}{h_{\nu)\delta ,\chi \zeta}}) = - {2 \over \kappa}T_{\mu \nu}^{{\rm{mat}}}$$

(48)

in a transverse-traceless gauge, which can be shown to exist in this theory [11, 460]. The constraint of Eq. (47) is identically satisfied to second order in the metric perturbation. However, without further information about ϑ one cannot proceed any further, except for a few general observations. As is clear from Eq. (48), the evolution equation for the metric perturbation can contain third time derivatives, which generically will lead to instabilities. In fact, as shown in [13] the general solution to these equations will contain exponentially growing and decaying modes. However, the theory defined by Eq. (45) is an effective theory, and thus, there can be higher-order operators not included in this action that may stabilize the solution. Regardless, when studying this theory order-reduction is necessary if one is to consider it an effective model.

Let us now discuss the properties of such an effective theory. Because of the structure of the modification to the field equations, one can always choose a sufficiently small value for α such that all solar system tests are satisfied. In fact, one can see from the equations in this section that in the limit α → 0, one recovers GR. Non-dynamical Chern—Simons gravity leads to modifications to the non-radiative (near-zone) metric in the gravitomagnetic sector, leading to corrections to Lense—Thirring precession [14, 15]. This fact has been used to constrain the theory through observations of the orbital motion of the LAGEOS satellites [388] to $(\alpha/\kappa){\dot \vartheta ^{- 1}} < 2 \times {10^4}$ km, or equivalently $(\kappa/\alpha){\dot \vartheta ^{- 1}} \underset{\sim}{>}{10^{- 14}}{\rm{eV}}$. However, much better constraints can be placed through observations of the binary pulsar [472, 18]: $({\alpha _4}/\kappa)\dot \vartheta < 0.8$ km.

Some of the sub-properties of the fundamental requirement are satisfied in non-dynamical Chern—Simons gravity. On the one hand, all spherically-symmetric metrics that are solutions to the Einstein equations are also solutions in this theory for a “canonical” scalar field (θ ∞ t) [207]. On the other hand, axisymmetric solutions to the Einstein equations are generically not solutions in this theory. Moreover, although spherically-symmetric solutions are preserved, perturbations of such spacetimes that are solutions to the Einstein equations are not generically solutions to the modified theory [470]. What is perhaps worse, the evolution of perturbations to non-spinning black holes have been found to be generically overconstrained [470]. This is a consequence of the lack of scalar field dynamics in the modified theory, which, via Eq. (47), tends to overconstrain it. Such a conclusion also suggests that this theory does not posses a well-posed initial-value problem.

One can argue that non-dynamical Chern—Simons gravity is well-motivated from fundamental theories [17], except that in the latter, the scalar field is always dynamical, instead of having to be prescribed a priori. Thus, perhaps the strongest motivation for such a model is as a phenomenological proxy to test whether the gravitational interaction remains parity invariant in the strong field, a test that is uniquely suited to this modified model.

2.4 Currently unexplored theories in the gravitational-wave sector

The list of theories we have described here is by no means exhaustive. In fact, there are many fascinating theories that we have chosen to leave out because they have not yet been analyzed in the gravitational wave context in detail. Examples of these include the following:

• Einstein-Aether Theory [247] and Hořava—Lifshitz Theory [234];

• Standard Model Extension [109];

• Eddington-inspired Born—Infeld gravity [48];

• New Massive Gravity [60, 136] and Bi-Gravity Theories [349, 346, 219, 220].

We will update this review with a description of these theories, once a detailed gravitational-wave study for compact binaries or supernovae sources is carried out and the predictions for the gravitational waveform observables are made for any physical system plausibly detectable by current or near future gravitational-wave experiments.

3 Detectors

3.1 Gravitational-wave interferometers

Kilometer-scale gravitational-wave interferometers have been in operation for over a decade. These types of detectors use laser interferometry to monitor the locations of test masses at the ends of the arms with exquisite precision. Gravitational waves change the relative length of the optical cavities in the interferometer (or equivalently, the proper travel time of photons) resulting in a strain

$$h = {{\Delta L} \over L},$$

where ΔL is the path length difference between the two arms of the interferometer.

Fractional changes in the difference in path lengths along the two arms can be monitored to better than 1 part in 10²⁰. It is not hard to understand how such precision can be achieved. For a simple Michelson interferometer, a difference in path length of order the size of a fringe can easily be detected. For the typically-used, infrared lasers of wavelength λ ∼ 1 μm, and interferometer arms of length L = 4 km, the minimum detectable strain is

$$h \sim {\lambda \over L} \sim 3 \times {10^{- 10}}.$$

This is still far off the 10⁻²⁰ mark. In principle, however, changes in the length of the cavities corresponding to fractions of a single fringe can also be measured provided we have a sensitive photodiode at the dark port of the interferometer, and enough photons to perform the measurement. This way we can track changes in the amount of light incident on the photodiode as the lengths of the arms change and we move over a fringe. The rate at which photons arrive at the photodiode is a Poisson process and the fluctuations in the number of photons is ∼ N¹/², where N is the number of photons. Therefore we can track changes in the path length difference of order

$$\Delta L\sim{\lambda \over {{N^{1/2}}}}.$$

The number of photons depends on the laser power P, and the amount of time available to perform the measurement. For a gravitational wave of frequency f, we can collect photons for a time t ∼ 1/f, so the number of photons is

$$N\sim{P \over {f{h_p}\nu}},$$

where h _p is Planck’s constant and ν = c/λ is the laser frequency. For a typical laser power P ∼ 1 W, a gravitational-wave frequency f = 100 Hz, and λ ∼ 1 μm the number of photons is

$$N\sim{10^{16}},$$

so that the strain we are sensitive to becomes

$$h\sim{10^{- 18}}.$$

The sensitivity can be further improved by increasing the effective length of the arms. In the LIGO instruments, for example, each of the two arms forms a resonant Fabry-Pérot cavity. For gravitational-wave frequencies smaller than the inverse of the light storage time, the light in the cavities makes many back and forth trips in the arms, while the wave is traversing the instrument. For gravitational waves of frequencies around 100 Hz and below, the light makes about a thousand back and forth trips while the gravitational wave is traversing the interferometer, which results in a three-orders-of-magnitude improvement in sensitivity,

$$h\sim{10^{- 21}}.$$

For frequencies larger than 100 Hz the number of round trips the light makes in the Fabry-Pérot cavities while the gravitational wave is traversing the instrument is reduced and the sensitivity is degraded.

The proper light travel time of photons in interferometers is controlled by the metric perturbation, which can be expressed as a sum over polarization modes

$${h_{ij}}(t,\vec x) = \sum\limits_A {h_{ij}^A} (t,\vec x),$$

(49)

where A labels the six possible polarization modes in metric theories of gravity. The metric perturbation for each mode can be written in terms of a plane wave expansion,

$$h_{ij}^A(t,\vec x) = \int\nolimits_{- \infty}^\infty {df} \int\nolimits_{{S^2}} {d\hat \Omega} \;{e^{i2\pi f(t - \hat \Omega \cdot \vec x)}}{\tilde h^A}(f,\hat \Omega)\epsilon _{ij}^A(\hat \Omega){.}$$

(50)

Here f is the frequency of the gravitational waves, $\vec k = 2\pi f\hat \Omega$ is the wave vector, Ω is a unit vector that points in the direction of propagation of the gravitational waves, Ω is the Ath polarization tensor, with i,j = x,y,z spatial indices. The metric perturbation due to mode A from the direction $\hat \Omega$ can be written by integrating over all frequencies,

$$h_{ij}^A(t - \hat \Omega \cdot \vec x) = \int\nolimits_{- \infty}^\infty {df} \;{e^{i2\pi f(t - \hat \Omega \cdot\vec x)}}{\tilde h^A}(f,\hat \Omega)\epsilon_{ij}^A(\hat \Omega){.}$$

(51)

By integrating Eq. (50) over all frequencies we have an expression for the metric perturbation from a particular direction $\hat \Omega$, i.e., only a function of $t - \hat \Omega \cdot \vec x$. The full metric perturbation due to a gravitational wave from a direction $\hat \Omega$ can be written as a sum over all polarization modes

$${h_{ij}}(t - \hat \Omega \cdot \vec x) = \sum\limits_A {{h^A}} (t - \hat \Omega \cdot \vec x)\epsilon_{ij}^A(\hat \Omega).$$

(52)

The response of an interferometer to gravitational waves is generally referred to as the antenna pattern response, and depends on the geometry of the detector and the direction and polarization of the gravitational wave. To derive the antenna pattern response of an interferometer for all six polarization modes we follow the discussion in [329] closely. For a gravitational wave propagating in the z direction, the polarization tensors are as follows

$$\begin{array}{*{20}c} {\epsilon _{ij}^ + = \left({\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & {- 1} & 0 \\ 0 & 0 & 0 \\ \end{array}} \right),\epsilon _{ij}^ \times = \left({\begin{array}{*{20}c} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array}} \right),} \\ {\epsilon _{ij}^x = \left({\begin{array}{*{20}c} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ \end{array}} \right),\epsilon _{ij}^y = \left({\begin{array}{*{20}c} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{array}} \right),} \\ {\epsilon _{ij}^b = \left({\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{array}} \right),\epsilon _{ij}^\ell = \left({\begin{array}{*{20}c} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \\ \end{array}} \right),} \\ \end{array}$$

(53)

where the superscripts +, ×, x, y, b, and ℓ correspond to the plus, cross, vector-x, vector-y, breathing, and longitudinal modes.

Suppose that the coordinate system for the detector is $\hat x = (1,0,0),\,\hat y = (0,1,0),\,\hat z = (0,0,1)$, as in Figure 1. Relative to the detector, the gravitational-wave coordinate system is rotated by angles $(\theta ,\phi),\,\hat x{\prime} = (\cos \theta \cos \phi ,\cos \theta \sin \phi , - \sin \theta),\,\hat y{\prime} = (\sin \phi ,\cos \phi ,0)$, and $\hat z{\prime}(\sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \theta)$. We still have the freedom to perform a rotation about the gravitational-wave propagation direction, which introduces the polarization angle ψ,

$$\begin{array}{*{20}c} {\hat m = \hat x{\prime}\cos \psi + \hat y{\prime}\sin \psi ,}\\ {\hat n = - \hat x{\prime}\sin \psi + \hat y{\prime}\cos \psi ,}\\ {\hat \Omega = \hat z{\prime}.\quad \quad \quad \quad \quad \quad \;\;}\\ \end{array}$$

(54)

The coordinate systems $(\hat x,\hat y,\hat z)$ and $(\hat m,\hat n,\hat \Omega)$ are also shown in Figure 1. To generalize the polarization tensors in Eq. (53) to a wave coming from a direction $\hat \Omega$, we use the unit vectors $\hat m,\hat n$, and $\hat \Omega$ as follows

$$\begin{array}{*{20}c} {{\epsilon ^ +} = \hat m \otimes \hat m - \hat n \otimes \hat n,}\\ {{\epsilon ^ \times} = \hat m \otimes \hat n + \hat n \otimes \hat m,}\\ {{\epsilon ^x} = \hat m \otimes \hat \Omega + \hat \Omega \otimes \hat m,}\\ {{\epsilon ^y} = \hat n \otimes \hat \Omega + \hat \Omega \otimes \hat n,}\\ {{\epsilon ^b} = \hat m \otimes \hat m + \hat n \otimes \hat n,}\\ {{\epsilon ^\ell} = \hat \Omega \otimes \hat \Omega .\quad \quad \quad}\\ \end{array}$$

(55)

For LIGO and VIRGO the arms are perpendicular so that the antenna pattern response can be written as the difference of projection of the polarization tensor onto each of the interferometer arms,

$${F^A}(\hat \Omega ,\psi) = {1 \over 2}({\hat x^i}{\hat x^j} - {\hat y^i}{\hat y^j})\epsilon _{ij}^A(\hat \Omega ,\psi){.}$$

(56)

This means that the strain measured by an interferometer due to a gravitational wave from direction $\hat \Omega$ and polarization angle ψ takes the form

$$h(t) = \sum\limits_A {{h_A}(t - \hat \Omega \cdot x){F^A}} (\hat \Omega ,\psi){.}$$

(57)

Explicitly, the antenna pattern functions are,

$$\begin{array}{*{20}c} {{F^ +}(\theta ,\phi ,\psi) = {1 \over 2}(1 + {{\cos}^2}\theta)\cos 2\phi \cos 2\psi - \cos \theta \sin 2\phi \sin 2\psi ,}\\ {{F^ \times}(\theta ,\phi ,\psi) = - {1 \over 2}(1 + {{\cos}^2}\theta)\cos 2\phi \sin 2\psi - \cos \theta \sin 2\phi \cos 2\psi ,}\\ {{F^x}(\theta ,\phi ,\psi) = \sin \theta (\cos \theta \cos 2\phi \cos \psi - \sin 2\phi \sin \psi),\quad \quad \quad \quad \;\;}\\ {{F^y}(\theta ,\phi ,\psi) = - \sin \theta (\cos \theta \cos 2\phi \sin \psi + \sin 2\phi \cos \psi),\quad \quad \quad \quad}\\ {{F^b}(\theta ,\phi) = - {1 \over 2}{{\sin}^2}\theta \cos 2\phi ,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ {{F^\ell}(\theta ,\phi) = {1 \over 2}{{\sin}^2}\theta \cos 2\phi .\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ \end{array}$$

(58)

The dependence on the polarization angles ψ reveals that the + and × polarizations are spin-2 tensor modes, the x and y polarizations are spin-1 vector modes, and the b and ℓ polarizations are spin-0 scalar modes. Note that for interferometers the antenna pattern responses of the scalar modes are degenerate. Figure 2 shows the antenna patterns for the various polarizations given in Eq. (58) with ψ = 0. The color indicates the strength of the response with red being the strongest and blue being the weakest.

Figure 1

Detector coordinate system and gravitational-wave coordinate system.

Figure 2

Antenna pattern response functions of an interferometer (see Eqs. (58)) for ψ = 0. Panels (a) and (b) show the plus (∣F₊∣) and cross (∣F_×∣) modes, panels (c) and (d) the vector x and vector y modes (∣F _x ∣ and ∣F _y ∣), and panel (e) shows the scalar modes (up to a sign, it is the same for both breathing and longitudinal). Color indicates the strength of the response with red being the strongest and blue being the weakest. The black lines near the center give the orientation of the interferometer arms.

3.2 Pulsar timing arrays

Neutron stars can emit powerful beams of radio waves from their magnetic poles. If the rotational and magnetic axes are not aligned, the beams sweep through space like the beacon on a lighthouse. If the line of sight is aligned with the magnetic axis at any point during the neutron star’s rotation the star is observed as a source of periodic radio-wave bursts. Such a neutron star is referred to as a pulsar. Due to their large moment of inertia pulsars are very stable rotators, and their radio pulses arrive on Earth with extraordinary regularity. Pulsar timing experiments exploit this regularity: gravitational waves are expected to cause fluctuations in the time of arrival of radio pulses from pulsars.

The effect of a gravitational wave on the pulses propagating from a pulsar to Earth was first computed in the late 1970s by Sazhin and Detweiler [378, 145]. Gravitational waves induce a redshift in the pulse train

$$z(t,\hat \Omega) = {1 \over 2}{{{{\hat p}^i}{{\hat p}^j}} \over {1 + \hat \Omega \cdot \hat p}}\Delta {h_{ij}},$$

(59)

where $\hat p$ is a unit vector that points in the direction of the pulsar, $\hat \Omega$ is the unit vector of gravitational wave propagation, and

$$\Delta {h_{ij}} \equiv {h_{ij}}({t_{\rm{e}}},\hat \Omega) - {h_{ij}}({t_{\rm{p}}},\hat \Omega),$$

(60)

is the difference in the metric perturbation at the pulsar when the pulse was emitted and the metric perturbation on Earth when the pulse was received. The inner product in Eq. (59) is computed with the Euclidean metric.

In pulsar timing experiments it is not the redshift, but rather the timing residual that is measured. The times of arrival of pulses are measured and the timing residual is produced by subtracting off a model that includes the rotational frequency of the pulsar, the spin-down (frequency derivative), binary parameters if the pulsar is in a binary, sky location and proper motion, etc. The timing residual induced by a gravitational wave, R(t), is just the integral of the redshift

$$R(t) \equiv \int\nolimits_0^t {dt{\prime}z(t{\prime})}.$$

(61)

Times-of-arrival (TOAs) are measured a few times a year over the course of several years allowing for gravitational waves in the nano-Hertz band to be probed. Currently, the best timed pulsars have residual root-mean-squares (RMS) of a few 10 s of ns over a few years.

The equations above ((59)ff) can be used to estimate the strain sensitivity of pulsar timing experiments. For gravitational waves of frequency f the expected induced residual is

$$R \sim {h \over f},$$

so that for pulsars with RMS residuals R ∼ 100 ns, and gravitational waves of frequency f ∼ 10⁻⁸ Hz, gravitational waves with strains

$$h\sim Rf\sim {10^{- 5}}$$

would produce a measurable effect.

To find the antenna pattern response of the pulsar-Earth system, we are free to place the pulsar on the z-axis. The response to gravitational waves of different polarizations can then be written as

$${F^A}(\hat \Omega ,\psi) = {1 \over 2}{{{{\hat z}^i}{{\hat z}^j}} \over {1 + \cos \theta}}\epsilon _{ij}^A(\hat \Omega ,\psi),$$

(62)

which allows us to express the Fourier transform of (59) as

$$\tilde z(f,\hat \Omega) = \left({1 - {e^{- 2\pi ifL(1 + \hat \Omega \cdot \hat p)}}} \right)\sum\limits_A {{{\tilde h}_A}(f,\hat \Omega){F^A}(\hat \Omega),}$$

(63)

where the sum is over all possible gravitational-wave polarizations: A = +, ×,x,y,b,l, and L is the distance to the pulsar.

Explicitly,

$${F^ +}(\theta ,\psi) = {\sin ^2}{\theta \over 2}\cos 2\psi ,$$

(64)

$${F^ \times}(\theta ,\psi) = - {\sin ^2}{\theta \over 2}\sin 2\psi ,$$

(65)

$${F^x}(\theta ,\psi) = - {1 \over 2}{{\sin 2\theta} \over {1 + \cos \theta}}\cos \psi ,$$

(66)

$${F^y}(\theta ,\psi) = {1 \over 2}{{\sin 2\theta} \over {1 + \cos \theta}}\sin \psi ,$$

(67)

$${F^b}(\theta) = {\sin ^2}{\theta \over 2},$$

(68)

$${F^\ell}(\theta) = {1 \over 2}{{{{\cos}^2}\theta} \over {1 + \cos \theta}}.$$

(69)

Just like for the interferometer case, the dependence on the polarization angle ψ, reveals that the + and × polarizations are spin-2 tensor modes, the x and y polarizations are spin-1 vector modes, and the b and ℓ polarizations are spin-0 scalar modes. Unlike interferometers, the antenna pattern responses of the pulsar-Earth system do not depend on the azimuthal angle of the gravitational wave, and the scalar modes are not degenerate.

In the literature, it is common to write the antenna pattern response by fixing the gravitational-wave direction and changing the location of the pulsar. In this case the antenna pattern responses are [284, 22, 99]

$${\tilde F^ +}({\theta _p},{\phi _p}) = {\sin ^2}{{{\theta _p}} \over 2}\cos 2{\phi _p},$$

(70)

$${\tilde F^ \times}({\theta _p},{\phi _p}) = {\sin ^2}{{{\theta _p}} \over 2}\sin 2{\phi _p},$$

(71)

$${\tilde F^x}({\theta _p},{\phi _p}) = {1 \over 2}{{\sin 2{\theta _p}} \over {1 + \cos {\theta _p}}}\cos {\phi _p},$$

(72)

$${\tilde F^y}({\theta _p},{\phi _p}) = {1 \over 2}{{\sin 2{\theta _p}} \over {1 + \cos {\theta _p}}}\sin {\phi _p},$$

(73)

$${\tilde F^b}({\theta _p}) = {\sin ^2}{{{\theta _p}} \over 2},$$

(74)

$${\tilde F^\ell}({\theta _p}) = {1 \over 2}{{{{\cos}^2}{\theta _p}} \over {1 + \cos {\theta _p}}},$$

(75)

where θ _p and ϕ _p are the polar and azimuthal angles, respectively, of the vector pointing to the pulsar. Up to signs, these expressions are the same as Eq. (69) taking θ → θ _p and ψ → ϕ _p . This is because fixing the gravitational-wave propagation direction while allowing the pulsar location to change is analogous to fixing the pulsar position while allowing the direction of gravitational-wave propagation to change — there is degeneracy in the gravitational-wave polarization angle and the pulsar’s azimuthal angle ϕ _p . For example, changing the polarization angle of a gravitational wave traveling in the z-direction is the same as performing a rotation about the z-axis that changes the pulsar’s azimuthal angle. Antenna patterns for the pulsar-Earth system using Eqs. (75) are shown in Figure 3. The color indicates the strength of the response, red being the largest and blue the smallest.

Figure 3

Antenna patterns for the pulsar-Earth system. The plus mode is shown in (a), breathing modes in (b), the vector-x mode in (c), and longitudinal modes in (d), as computed from Eq. (75). The cross mode and the vector-y mode are rotated versions of the plus mode and the vector-x mode, respectively, so we did not include them here. The gravitational wave propagates in the positive z-direction with the Earth at the origin, and the antenna pattern depends on the pulsar’s location. The color indicates the strength of the response, red being the largest and blue the smallest.

4 Testing Techniques

4.1 Coalescence analysis

Gravitational waves emitted during the inspiral, merger and ringdown of compact binaries are the most studied in the context of data analysis and parameter estimation. In this section, we will review some of the main data analysis techniques employed in the context of parameter estimation and tests of GR. We begin with a discussion of matched filtering and Fisher theory (for a detailed review, see [173, 103, 125, 174, 248]). We then continue with a discussion of Bayesian parameter estimation and hypothesis testing (for a detailed review, see [387, 205, 123, 294]).

4.1.1 Matched filtering and Fisher’s analysis

When the detector noise n(t) is Gaussian and stationary, and when the signal s(t) is known very well, the optimal detection strategy is matched filtering. For any given realization, such noise can be characterized by its power spectral density S _n (f), defined via

$$\langle \tilde n(f)\tilde n^\ast (f{\prime})\rangle = {1 \over 2}{S_n}(f)\delta (f - f{\prime}),$$

(76)

where the tilde stands for the Fourier transform, the asterisk for complex conjugation and the brackets for the expectation value.

The detectability of a signal is determined by its signal-to-noise ratio or SNR, which is defined via

$${\rho ^2} = {{(s\vert h)} \over {\sqrt {(h\vert n)(n\vert h)}}},$$

(77)

where h is a template with parameters λⁱ and we have defined the inner product

$$(A\vert B) \equiv 4\Re \int\nolimits_0^\infty {{{\tilde A\ast\tilde B} \over {{S_n}}}df.}$$

(78)

If the templates do not exactly match the signal, then the SNR is reduced by a factor of $\mathcal{M}$, called the match:

$$\bar {\mathcal M} \equiv {{(s\vert h)} \over {\sqrt {(s\vert s)(h\vert h)}}},$$

(79)

where $1 - \bar{\mathcal{M}} = \mathcal{M}\mathcal{M}$ is the mismatch.

For the noise assumptions made here, the probability of measuring s(t) in the detector output, given a template h, is given by

$$p \propto {e^{- (s - h\vert s - h)/2}},$$

(80)

and thus the waveform h that best fits the signal is that with best-fit parameters such that the argument of the exponential is minimized. For large SNR, the best-fit parameters will have a multivariate Gaussian distribution centered on the true values of the signal ${\hat \lambda ^i}$, and thus, the waveform parameters that best fit the signal minimize the argument of the exponential. The parameter errors δλⁱ will be distributed according to

$$p(\delta {\lambda ^i}) \propto {e^{- {1 \over 2}{\Gamma _{ij}}\delta {\lambda ^i}\delta {\lambda ^j}}},$$

(81)

where Γ _ij is the Fisher matrix

$${\Gamma _{ij}} \equiv \left({{{\partial h} \over {\partial {\lambda ^i}}}\left\vert {{{\partial h} \over {\partial {\lambda ^j}}}} \right.} \right).$$

(82)

The root-mean-squared (1σ) error on a given parameter λ^ī is then

$$\sqrt {\langle {{(\delta {\lambda ^{\bar i}})}^2}\rangle} = \sqrt {{\Sigma ^{\overline {ii}}}} ,$$

(83)

where Σ^ij ≡ (Γ _ij )⁻¹ is the variance-covariance matrix and summation is not implied in Eq. (83) (λ^ī denotes a particular element of the vector λⁱ). This root-mean-squared error is sometimes referred to as the statistical error in the measurement of λ^ī. One can use Eq. (83) to estimate how well modified gravity parameters can be measured. Put another way, if a gravitational wave were detected and found consistent with GR, Eq. (83) would provide an estimate of how close to zero these modified gravity parameters would have to be.

The Fisher method to estimate projected constraints on modified gravity theory parameters is as follows. First, one constructs a waveform model in the particular modified gravity theory one wishes to constrain. Usually, this waveform will be similar to the GR one, but it will contain an additional parameter, κ, such that the template parameters are now λⁱ plus κ. Let us assume that as κ → 0, the modified gravity waveform reduces to the GR expectation. Then, the accuracy to which κ can be measured, or the accuracy to which we can say κ is zero, is approximately (Σ^κκ)^1/2, where the Fisher matrix associated with this variance-covariance matrix must be computed with the non-GR model evaluated at the GR limit (κ → 0). Such a method for estimating how well modified gravity theories can be constrained was pioneered by Will in [436, 353], and since then, it has been widely employed as a first-cut estimate of the accuracy to which different theories can be constrained.

The Fisher method described above can dangerously lead to incorrect results if abused [414, 415]. One must understand that this method is suitable only if the noise is stationary and Gaussian and if the SNR is sufficiently large. How large an SNR is required for Fisher methods to work depends somewhat on the signals considered, but usually for applications concerning tests of GR, one would be safe with ρ ≳ 30 or so. In real data analysis, the first two conditions are almost never satisfied. Moreover, the first detections that will be made will probably be of low SNR, i.e., ρ ∼ 8, for which again the Fisher method fails. In such cases, more sophisticated parameter estimation methods need to be employed.

4.1.2 Bayesian theory and model testing

Bayesian theory is ideal for parameter estimation and model selection. Let us then assume that a detection has been made and that the gravitational wave signal in the data can be described by some model $\mathcal{M}$, parameterized by the vector λⁱ. Using Bayes’ theorem, the posterior distribution function (PDF) or the probability density function for the model parameters, given data d and model $\mathcal{M}$, is

$$p({\lambda ^i}\vert d,{\mathcal M}) = {{p(d\vert {\lambda ^i},{\mathcal M})p({\lambda ^i}\vert {\mathcal M})} \over {p(d\vert {\mathcal M})}}.$$

(84)

Obviously, the global maximum of the PDF in the parameter manifold gives the best fit parameters for that model. The prior probability density $p({\lambda ^i}\vert \mathcal{M})$ represents our prior beliefs of the parameter range in model $\mathcal{M}$. The marginalized likelihood or evidence, is the normalization constant

$$p(d\vert {\mathcal M}) = \int {d{\lambda ^i}p} (d\vert {\lambda ^i},{\mathcal M})\;p({\lambda ^i}\vert {\mathcal M}),$$

(85)

which clearly guarantees that the integral of Eq. (84) integrates to unity. The quantity $p(d\vert {\lambda ^i},\mathcal{M})$ is the likelihood function, which is simply given by Eq. (80), with a given normalization. In that equation we used slightly different notation, with s being the data d and h the template associated with model $\mathcal{M}$ and parameterized by λⁱ. The marginalized PDF, which represents the probability density function for a given parameter λ^ī (recall that λ^ī is a particular element of λⁱ), after marginalizing over all other parameters, is given by

$$p({\lambda ^{\bar i}}\vert d,{\mathcal M}) = \int\nolimits_{i \neq \bar i} {d{\lambda ^i}p({\lambda ^i}\vert {\mathcal M})p(d\vert {\lambda ^i},{\mathcal M})} ,$$

(86)

where the integration is not to be carried out over ī.

Let us now switch gears to model selection. In hypothesis testing, one wishes to determine whether the data is more consistent with hypothesis A (e.g., that a GR waveform correctly models the signal) or with hypothesis B (e.g., that a non-GR waveform correctly models the signal). Using Bayes’ theorem, the PDF for model A given the data is

$$p(A\vert d) = {{p(d\vert A)p(A)} \over {p(d)}}.$$

(87)

As before, p(A) is the prior probability of hypothesis A, namely the strength of our prior belief that hypothesis A is correct. The normalization constant p(d) is given by

$$p(d) = \int {d{\mathcal M}\;p(d\vert {\mathcal M})\;p({\mathcal M}),}$$

(88)

where the integral is to be taken over all models. Thus, it is clear that this normalization constant does not depend on the model. Similar relations hold for hypothesis B by replacing A → B in Eq. (87).

When hypothesis A and B refer to fundamental theories of nature we can take different viewpoints regarding the priors. If we argue that we know nothing about whether hypothesis A or B better describes nature, then we would assign equal priors to both hypotheses. If, on the other hand, we believe GR is the correct theory of nature, based on all previous experiments performed in the solar system and with binary pulsars, then we would assign p(A) > p(B). This assigning of priors necessarily biases the inferences derived from the calculated posteriors, which is sometimes heavily debated when comparing Bayesian theory to a frequentist approach. However, this “biasing” is really unavoidable and merely a reflection of our state of knowledge of nature (for a more detailed discussion on such issues, please refer to [294]).

The integral over all models in Eq. (88) can never be calculated in practice, simply because we do not know all models. Thus, one is forced to investigate relative probabilities between models, such that the normalization constant p(d) cancels out. The odds ratio is defined by

$${{\mathcal O}_{A,B}} = {{p(A\vert d)} \over {p(B\vert d)}} = {{p(A)} \over {p(B)}}{{\mathcal B}_{A,B}},$$

(89)

where ${\mathcal{B}_{A,B}} \equiv p(d\vert A)/p(d\vert B)$ is the Bayes factor and the prefactor p(A)/p(B) is the prior odds. The odds-ratio is a convenient quantity to calculate because the evidence p(d), which is difficult to compute, actually cancels out. Recently, Vallisneri [416] has investigated the possibility of calculating the odds-ratio using only frequentist tools and without having to compute full evidences. The odds-ratio should be interpreted as the betting-odds of model A over model B. For example, an odds-ratio of unity means that both models are equally supported by the data, while an odds-ratio of 10² means that there is a 100 to 1 odds that model A better describes the data than model B.

The main difficulty in Bayesian inference (both in parameter estimation and model selection) is sampling the PDF sufficiently accurately. Several methods have been developed for this purpose, but currently the two main workhorses in gravitational-wave data analysis are Markov Chain Monte Carlo and Nested Sampling. In the former, one samples the likelihood through the Metropolis-Hastings algorithm [314, 221, 122, 367]. This is computationally expensive in high-dimensional cases, and thus, there are several techniques to improve the efficiency of the method, e.g., parallel tempering [402]. Once the PDF has been sampled, one can then calculate the evidence integral, for example via thermodynamic integration [420, 167, 419]. In Nested Sampling, the evidence is calculated directly by laying out a fixed number of points in the prior volume, which are then allowed to move and coalesce toward regions of high posterior probability. With the evidence in hand, one can then infer the PDF. As in the previous case, Nested Sampling can be computationally expensive in high-dimensional cases.

Del Pozzo et al. [142] were the first to carry out a Bayesian implementation of model selection in the context of tests of GR. Their analysis focused on tests of a particular massive graviton theory, using the gravitational wave signal from quasi-circular inspiral of non-spinning black holes. Cornish et al. [124, 376] extended this analysis by considering model-independent deviations from GR, using the parameterized post-Einsteinian (ppE) approach (Section 5.3.4) [467]. Recently, this was continued by Li et al. [290, 291], who carried out a similar analysis on a large statistical sample of Advanced LIGO (aLIGO) detections using simulated data and a restricted ppE model. All of these studies suggest that Bayesian tests of GR are possible, given sufficiently-high SNR events. Of course, whether deviations from GR are observable will depend on the strong-field character and strength of the deviation, as well as the availability of sufficiently-accurate GR waveforms.

4.1.3 Systematics in model selection

The model selection techniques described above are affected by other systematics present in data analysis. In general, we can classify these into the following [417]:

• Mismodeling Systematic, caused by inaccurate models of the gravitational-wave template.

• Instrumental Systematic, caused by inaccurate models of the gravitational-wave response.

• Astrophysical Systematic, caused by inaccurate models of the astrophysical environment.

Mismodeling systematics are introduced due to the lack of an exact solution to the Einstein equations from which to extract an exact template, given a particular astrophysical scenario. Inspiral templates, for example, are approximated through post-Newtonian theory and become increasingly less accurate as the binary components approach each other. Cutler and Vallisneri [127] were the first to carry out a formal and practical investigation of such a systematic in the context of parameter estimation from a frequentist approach.

Mismodeling systematics will prevent us from testing GR effectively with signals that we do not understand sufficiently well. For example, when considering signals from black hole coalescences, if the the total mass of the binary is sufficiently high, the binary will merge in band. The higher the total mass, the fewer the inspiral cycles that will be in band, until eventually only the merger is in band. Since the merger phase is the least understood phase, it stands to reason that our ability to test GR will deteriorate as the total mass increases. Of course, we do understand the ringdown phase very well, and tests of the no-hair theorem would be allowed during this phase, provided a sufficiently-high SNR [65]. On the other hand, for neutron star binaries or very-low-mass black-hole binaries, the merger phase is expected to be essentially out of band for aLIGO (above 1 kHz), and thus, the noise spectrum itself may shield us from our ignorance.

Instrumental systematics are introduced by our ignorance of the transfer function, which connects the detector output to the incoming gravitational waves. Through sophisticated calibration studies with real data, one can approximate the transfer function very well [4, 1]. However, this function is not time-independent, because the noise in the instrument is not stationary or Gaussian. Thus, un-modeled drifts in the transfer function can introduce systematics in parameter estimation that are as large as 10% in the amplitude and the phase [4].

Instrumental systematics can affect tests of GR, if these are performed with a single instrument. However, one expects multiple detectors to be online in the future and for gravitational-wave detections to be made in several of them simultaneously. Instrumental systematics should be present in all such detections, but since the noise will be mostly uncorrelated between different instruments, one should be able to ameliorate its effects through cross-correlating outputs from several instruments.

Astrophysical systematics are induced by our lack of a priori knowledge of the gravitational wave source. As explained above, matched filtering requires knowledge of a waveform template with which to filter the data. Usually, we assume the sources are in a perfect vacuum and isolated. For example, when considering inspiral signals, we ignore any third bodies, electric or magnetic fields, neutron star hydrodynamics, the expansion of the universe, etc. Fortunately, however, most of these effects are expected to be small: the probability of finding third bodies sufficiently close to a binary system is very small [463]; for low redshift events, the expansion of the universe induces an acceleration of the center of mass, which is also very small [468]; electromagnetic fields and neutron-star hydrodynamic effects may affect the inspiral of black holes and neutron stars, but not until the very last stages, when most signals will be out of band anyways. For example, tidal deformation effects enter a neutron-star-binary inspiral waveform at 5 post-Newtonian order, which therefore affects the signal outside of the most sensitive part of the aLIGO sensitivity bucket.

Perhaps the most dangerous source of astrophysical systematics is due to the assumptions made about the astrophysical systems we expect to observe. For example, when considering neutron-star-binary inspirals, one usually assumes the orbit will have circularized by the time it enters the sensitivity band. Moreover, one assumes that any residual spin angular momentum that the neutron stars may possess is very small and aligned or counter-aligned with the orbital angular momentum. These assumptions certainly simplify the construction of waveform templates, but if they happen to be wrong, they would introduce mismodeling systematics that could also affect parameter estimation and tests of GR.

4.2 Burst analyses

In alternative theories of gravity, gravitational-wave sources such as core collapse supernovae may result in the production of gravitational waves in more than just the plus and cross-polarizations [384, 380, 216, 334, 333, 369]. Indeed, the near-spherical geometry of the collapse can be a source of scalar breathing-mode gravitational waves. However, the precise form of the waveform is unknown because it is sensitive to the initial conditions.

When searching for un-modeled bursts in alternative theories of gravity, a general approach involves optimized linear combinations of data streams from all available detectors to form maximum likelihood estimates of the waveforms in the various polarizations, and the use of null streams. In the context of ground-based detectors and GR, these ideas were first explored by Gürsel and Tinto [212] and later by Chatterji et al. [101] with the aim of separating false-alarm events from real detections. The main idea was to construct a linear combination of data streams received by a network of detectors, so that the combination contained only noise. Of course, in GR one need only include h₊. and h_× polarizations, and thus a network of three detectors suffices. This concept can be extended to develop null tests of GR, as originally proposed by Chatziioannou et al. [102] and recently implemented by Hayama et al. [228].

Let us consider a network of D ≥ 6 detectors with uncorrelated noise and a detection by all D detectors. For a source that emits gravitational waves in the direction $\hat \Omega$, a single data point (either in the time-domain, or a time-frequency pixel) from an array of D detectors (either pulsars or interferometers) can be written as

$$d = Fh + n.$$

(90)

Here

$$d \equiv \left[ {\begin{array}{*{20}c} {{d_1}} \\ {{d_2}} \\ \vdots \\ {{d_D}} \\ \end{array}} \right],\quad h \equiv \left[ {\begin{array}{*{20}c} {{h_ +}} \\ {{h_ \times}} \\ {{h_x}} \\ {{h_y}} \\ {{h_b}} \\ {{h_\ell}} \\ \end{array}} \right],\quad n \equiv \left[ {\begin{array}{*{20}c} {{n_1}} \\ {{n_2}} \\ \vdots \\ {{n_D}} \\ \end{array}} \right],$$

(91)

where n is a vector with the noise. The antenna pattern functions are given by the matrix,

$$[{F^ +}{F^ \times}{F^x}{F^y}{F^b}{F^\ell}] \equiv \left[ {\begin{array}{*{20}c} {F_1^ +} & {F_1^ \times} & {F_1^x} & {F_1^y} & {F_1^b} & {F_1^\ell}\\ {F_2^ +} & {F_2^ \times} & {F_2^x} & {F_2^y} & {F_2^b} & {F_2^\ell}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {F_D^ +} & {F_D^ \times} & {F_D^x} & {F_D^y} & {F_D^b} & {F_D^\ell}\\ \end{array}} \right].$$

(92)

For simplicity we have suppressed the sky-location dependence of the antenna pattern functions. These can either be the interferometric antenna pattern functions in Eqs. (58), or the pulsar response functions in Eqs. (69). For interferometers, since the breathing and longitudinal antenna pattern response functions are degenerate, and even though F is a 6 × D matrix, there are only five linearly-independent vectors [81, 80, 102, 228].

If we do not know the form of the signal present in our data, we can obtain maximum likelihood estimators for it. For simplicity, let us assume the data are Gaussian and of unit variance (the latter can be achieved by whitening the data). Just as we did in Eq. (80), we can write the probability of obtaining datum d, in the presence of a gravitational wave h as

$$P(d\vert h) = {1 \over {{{(2\pi)}^{D/2}}}}\exp \left[ {- {1 \over 2}\vert d - Fh\vert ^{2}} \right].$$

(93)

The logarithm of the likelihood ratio, i.e., the logarithm of the ratio of the likelihood when a signal is present to that when a signal is absent, can then be written as

$$L \equiv \ln {{P(d\vert h)} \over {P(d\vert 0)}} = {1 \over 2}\left[ {\vert d\vert ^{2} - \vert d - Fh\vert ^{2}} \right].$$

(94)

if we treat the waveform values for each datum as free parameters, we can maximize the likelihood ratio

$$0 = {\left. {{{\partial L} \over {\partial h}}} \right\vert _{h = {h_{{\rm{MAX}}}}}},$$

(95)

and obtain maximum likelihood estimators for the gravitational wave,

$${h_{{\rm{MAX}}}} = {({F^T}F)^{- 1}}{F^T}d.$$

(96)

We can further substitute this solution into the likelihood, to obtain the value of the likelihood at the maximum,

$${E_{{\rm{SL}}}} \equiv 2L(h_{{\rm{MAX}}}) = {d^T}{P^{{\rm{GW}}}}d,$$

(97)

where

$${P^{{\rm{GW}}}} \equiv F{({F^T}F)^{- 1}}{F^T}.$$

(98)

The maximized likelihood can be thought of as the power in the signal, and can be used as a detection statistic. P^GW is a projection operator that projects the data into the subspace spanned by F. An orthogonal projector can also be constructed,

$${P^{{\rm{null}}}} \equiv (I - {P^{{\rm{GW}}}}),$$

(99)

which projects the data onto a sub-space orthogonal to F. Thus one can construct a certain linear combination of data streams that has no component of a certain polarization by projecting them to a direction orthogonal to the direction defined by the beam pattern functions of this polarization mode

$${d^{{\rm{null}}}} = {P^{{\rm{null}}}}d.$$

(100)

This is called a null stream and, in the context of GR, it was introduced as a means of separating false-alarm events due, say, to instrumental glitches from real detections [212, 101].

With just three independent detectors, we can choose to eliminate the two tensor modes (the plus and cross-polarizations) and construct a GR null stream: a linear combination of data streams that contains no signal consistent within GR, but could contain a signal in another gravitational theory, as illustrated in Figure 4. With such a GR null stream, one can carry out null tests of GR and study whether such a stream contains any statistically-significant deviations from noise. Notice that this approach does not require a template; if one were parametrically constructed, such as in [102], more powerful null tests could be applied. In the future, we expect several gravitational wave detectors to be online: the two aLIGO ones in the United States, Advanced VIRGO (adVirgo) in Italy, LIGO-India in India, and KAGRA in Japan. Given a gravitational-wave observation that is detected by all five detectors, one can then construct three GR null streams, each with power in a signal direction.

Figure 4

Schematic diagram of the projection of the data stream d orthogonal to the GR subspace spanned by F⁺ and F^×, along with a perpendicular subspace, for 3 detectors to build the GR null stream.

For pulsar timing experiments where one is dealing with data streams of about a few tens of pulsars, waveform reconstruction for all polarization states, as well as numerous null streams, can be constructed.

4.3 Stochastic background searches

Much work has been done on the response of ground-based interferometers to non-tensorial polarization modes, stochastic background detection prospects, and data analysis techniques [299, 323, 191, 329, 121]. In the context of pulsar timing, the first work to deal with the detection of such backgrounds in the context of alternative theories of gravity is due to Lee et al. [284], who used a coherence statistic approach to the detection of non-Einsteinian polarizations. They studied the number of pulsars required to detect the various extra polarization modes, and found that pulsar timing arrays are especially sensitive to the longitudinal mode. Alves and Tinto [22] also found enhanced sensitivity to longitudinal and vector modes. Here we follow the work in [329, 99] that deals with the LIGO and pulsar timing cases using the optimal statistic, a cross-correlation that maximizes the SNR.

In the context of the optimal statistic, the derivations of the effect of extra polarization states for ground-based instruments and pulsar timing are very similar. We begin with the metric perturbation written in terms of a plane wave expansion, as in Eq. (50). If we assume that the background is unpolarized, isotropic, and stationary, we have that

$$\langle \tilde h_A^{\ast}(f,\hat \Omega){\tilde h_{A{\prime}}}(f{\prime},\hat \Omega {\prime})\rangle = {\delta ^2}(\hat \Omega ,\hat \Omega {\prime}){\delta _{AA{\prime}}}\delta (f - f{\prime}){H_A}(f),$$

(101)

where H _A (f) is the gravitational-wave power spectrum for polarization A. H _A (f) is related to the energy density in gravitational waves per logarithmic frequency interval for that polarization through

$${\Omega _A}(f) \equiv {1 \over {{\rho _{{\rm{crit}}}}}}{{d{\rho _A}} \over {d\ln f}},$$

(102)

where ${\rho _{{\rm{crit}}}} = 3H_0^2/8\pi$ is the closure density of the universe, and

$${\rho _A} = {1 \over {32\pi}}\langle {\dot h_{Aij}}(t,\vec x)\dot h_A^{ij}(t,\vec x)\rangle$$

(103)

is the energy density in gravitational waves for polarization A. It follows from the plane wave expansion in Eq. (51), along with Eqs. (101) and (102) in Eq. (103), that

$${H_A}(f) = {{3H_0^2} \over {16{\pi ^3}}}\vert f\vert ^{{- 3}}{\Omega _A}(\vert f\vert),$$

(104)

and therefore

$$\langle \tilde h_A^{\ast}(f,\hat \Omega){\tilde h_{A{\prime}}}(f{\prime},\hat \Omega {\prime})\rangle = {{3H_0^2} \over {16{\pi ^3}}}{\delta ^2}(\hat \Omega ,\hat \Omega {\prime}){\delta _{AA{\prime}}}\delta (f - f{\prime})\vert f\vert ^{{- 3}}{\Omega _A}(\vert f\vert).$$

(105)

For both ground-based interferometers and pulsar-timing experiments, an isotropic stochastic background of gravitational waves appears in the data as correlated noise between measurements from different instruments. The data set from the a^th instrument is of the form

$${d_a}(t) = {s_a}(t) + {n_a}(t),$$

(106)

where s _a (t) corresponds to the gravitational-wave signal and n _a (t) to noise. The noise is assumed in this case to be stationary and Gaussian, and uncorrelated between different detectors,

$$\langle {n_a}(t)\rangle = 0,$$

(107)

$$\langle {n_a}(t){n_b}(t)\rangle = 0,$$

(108)

for a ≠ b.

Since the gravitational-wave signal is correlated, we can use cross-correlations to search for it. The cross-correlation statistic is defined as

$${S_{ab}} = \int\nolimits_{- T/2}^{T/2} {dt} \int\nolimits_{- T/2}^{T/2} {dt} {\prime}{d_a}(t){d_b}(t{\prime}){Q_{ab}}(t - t{\prime}),$$

(109)

where Q _ab (t − t′) is a filter function to be determined. Henceforth, no summation is implied on the detector indices (a, b, …). At this stage it is not clear why Q _ab (t − t′) depends on the pair of data sets being correlated. We will show how this comes about later. The optimal filter is determined by maximizing the expected SNR

$${\rm{SNR =}}{{{\mu _{ab}}} \over {{\sigma _{ab}}}}.$$

(110)

Here μ _ab is the mean 〈S _ab ,〉 and σ _ab is the square root of the variance $\sigma _{ab}^2 = \left\langle {S_{ab}^2} \right\rangle - {\left\langle {S_{ab}^2} \right\rangle ^2}$.

The expressions for the mean and variance of the cross-correlation statistic, μ _ab and ${\mu _{ab}}$ respectively, take the same form for both pulsar timing and ground-based instruments. In the frequency domain, Eq. (109) becomes

$${S_{ab}} = \int\nolimits_{- \infty}^\infty {df} \int\nolimits_{- \infty}^\infty {d{f\prime}} {\delta_T}(f - f{\prime})\tilde d_a^{\ast}(f){\tilde d_b}(f{\prime}){\tilde Q_{ab}}(f{\prime}),$$

(111)

by the convolution theorem, and the mean μ is then

$${\mu _{ab}} \equiv \langle {S_{ab}}\rangle = \int\nolimits_{- \infty}^\infty {df} \int\nolimits_{- \infty}^\infty {d{f\prime}} {\delta _T}(f - f{\prime})\langle \tilde s_a^{\ast}(f){\tilde s_b}(f{\prime})\rangle {\tilde Q_{ab}}(f{\prime}),$$

(112)

where δ _T is the finite time approximation to the delta function, δ _T (f) = sinπft/(πf). With this in hand, the mean of the cross-correlation statistic is

$${\mu _{ab}}{{3H_0^2} \over {16{\pi ^3}}}T\sum\limits_A {\int\nolimits_{- \infty}^\infty {df}} \vert f\vert ^{{- 3}}{\tilde Q_{ab}}(f){\Omega _A}(f)\Gamma _{ab}^A(f),$$

(113)

and the variance in the weak signal limit is

$$\begin{array}{*{20}c} {\sigma _{ab}^2 \equiv \langle S_{ab}^2\rangle - {{\langle {S_{ab}}\rangle}^2} \approx \langle S_{ab}^2\rangle \quad \quad \quad \quad \quad \quad \;}\\ {\approx {T \over 4}\int\nolimits_{- \infty}^\infty {df} {P_a}(\vert f\vert){P_b}(\vert f\vert){{\left\vert {{{\tilde Q}_{ab}}(f)} \right\vert}^2},}\\ \end{array}$$

(114)

where the one-sided power spectra of the noise are defined by

$$\langle \tilde n_a^{\ast}(f){\tilde n_a}(f{\prime})\rangle = {1 \over 2}\delta (f - f{\prime}){P_a}(\vert f\vert),$$

(115)

in analogy to Eq. (76), where P _a plays here the role of S _n (f).

The mean and variance can be rewritten more compactly if we define a positive-definite inner product using the noise power spectra of the two data streams

$${(A,B)_{ab}} \equiv \int\nolimits_{- \infty}^\infty {df} A^\ast(f)B(f){P_a}(\vert f\vert){P_b}(\vert f\vert),$$

(116)

again in analogy to the inner product in Eq. (78), when considering inspirals. Using this definition

$${\mu _{ab}} = {{3H_0^2} \over {16{\pi ^3}}}T{\left({{{\tilde Q}_{ab}},{{{\Sigma _A}{\Omega _A}(\vert f\vert)\Gamma _{ab}^A(\vert f\vert)} \over {\vert f\vert ^{3}{P_a}(\vert f\vert){P_b}(\vert f\vert)}}} \right)_{ab}},$$

(117)

$$\sigma _{ab}^2 \approx {T \over 4}{\left({\tilde Q,\tilde Q} \right)_{ab}},$$

(118)

where we recall that the capital Latin indices (A, B, …) stand for the polarization content. From the definition of the SNR and the Schwartz’s inequality, it follows that the optimal filter is given by

$${\tilde Q_{ab}}(f) = N{{{\Sigma _A}{\Omega _A}(\vert f\vert)\Gamma _{ab}^A(\vert f\vert)} \over {\vert f\vert ^{3}{P_a}(\vert f\vert){P_b}(\vert f\vert)}},$$

(119)

where N is an arbitrary normalization constant, normally chosen so that the mean of the statistic gives the amplitude of the stochastic background.

The differences in the optimal filter between interferometers and pulsars arise only from differences in the overlap reduction functions, $\Gamma _{ab}^A(f)$. For ground-based instruments, the signal data s _a are the strains given by Eq. (57). The overlap reduction functions are then given by

$$\Gamma _{ab}^A(f) = \int\nolimits_{{S^2}} {d\hat \Omega F_a^A(\hat \Omega)F_b^A(\hat \Omega)} {e^{2\pi if\hat \Omega \cdot ({{\vec x}_a} - {{\vec x}_b})}},$$

(120)

where ${\vec x_a}$ and ${\vec x_b}$ are the locations of the two interferometers. The integrals in this case all have solutions in terms of spherical Bessel functions [329], which we do not summarize here for brevity.

For pulsar timing arrays, the signal data s _a are the redshifts z _a , given by Eq. (63). The overlap reduction functions are then given by

$$\Gamma _{ab}^A(f) = {3 \over {4\pi}}\int\nolimits_{{S^2}} {d\hat \Omega \left({{e^{i2\pi f{L_a}(1 + \hat \Omega \cdot {{\hat p}_a})}} - 1} \right)\left({{e^{- i2\pi f{L_b}(1 + \hat \Omega \cdot {{\hat p}_b})}} - 1} \right)F_a^A(\hat \Omega)F_b^A(\hat \Omega)} ,$$

(121)

where L _a and L _b are the distances to the two pulsars. For all transverse modes pulsar timing experiments are in a regime where the exponential factors in Eq. (121) can be neglected [30, 99], and the overlap reduction functions effectively become frequency independent. For the + and × mode the overlap reduction function becomes

$$\Gamma _{ab}^ + = 3\left\{{{1 \over 3} + {{1 - \cos {\xi _{ab}}} \over 2}\left[ {\ln \left({{{1 - \cos {\xi _{ab}}} \over 2}} \right) - {1 \over 6}} \right]} \right\},$$

(122)

where ${\xi _{ab}} = {\rm{co}}{{\rm{s}}^{- 1}}({\hat p_a} \cdot {\hat p_b})$ is the angle between the two pulsars. This quantity is proportional to the Hellings and Downs curve [231]. For the breathing mode, the overlap reduction function takes the closed form expression [284]:

$$\Gamma _{ab}^b = {1 \over 4}(3 + \cos {\xi _{ab}}){.}$$

(123)

For the vector and longitudinal modes the overlap reduction functions remain frequency dependent and there are no general analytic solutions.

The result for the combination of cross-correlation pairs to form an optimal network statistic is also the same in both ground-based interferometer and pulsar timing cases: a sum of the cross-correlations of all detector pairs weighted by their variances. The detector network optimal statistic is,

$${S_{{\rm{opt}}}} = {{{\Sigma _{ab}}\sigma _{ab}^{- 2}{S_{ab}}} \over {{\Sigma _{ab}}\sigma _{ab}^{- 2}}},$$

(124)

where Σ _ab is a sum over all detector pairs.

In order to perform a search for a given polarization mode one first needs to compute the overlap reduction functions (using either Eq. (120) or (121)) for that mode. With that in hand and a form for the stochastic background spectrum Ω _A (f), one can construct optimal filters for all pairs in the detector network using Eq. (119), and perform the cross-correlations using either Eq. (109) (or equivalently Eq. (111)). Finally, we can calculate the overall network statistic Eq. (124), by first finding the variances using Eq. (114).

It is important to point out that the procedure outlined above is straightforward for ground-based interferometers. However, pulsar timing data are irregularly sampled, and have a pulsar-timing model subtracted out. This needs to be accounted for, and generally, a time-domain approach is more appropriate for these data sets. The procedure is similar to what we have outlined above, but power spectra and gravitational-wave spectra in the frequency domain need to be replaced by auto-covariance and cross-covariance matrices in the time domain that account for the model fitting (for an example of how to do this see [162]).

Interestingly, Nishizawa et al. [329] show that with three spatially-separated detectors the tensor, vector, and scalar contributions to the energy density in gravitational waves can be measured independently. Lee et al. [284] and Alves and Tinto [22] showed that pulsar timing experiments are especially sensitive to the longitudinal mode, and to a lesser extent the vector modes. Chamberlin and Siemens [99] showed that the sensitivity of the cross-correlation to the longitudinal mode using nearby pulsar pairs can be enhanced significantly compared to that of the transverse modes. For example, for the NANOGrav pulsar timing array, two pulsar pairs separated by 3° result in an enhancement of 4 orders of magnitude in sensitivity to the longitudinal mode relative to the transverse modes. The main contribution to this effect is due to gravitational waves that are coming from roughly the same direction as the pulses from the pulsars. In this case, the induced redshift for any gravitational-wave polarization mode is proportional to fL, the product of the gravitational-wave frequency and the distance to the pulsar, which can be large. When the gravitational waves and the pulse direction are exactly parallel, the redshift for the transverse and vector modes vanishes, but it is proportional to f L for the scalar-longitudinal mode.

Lee et al. [285] studied the detectability of massive gravitons in pulsar timing arrays through stochastic background searches. They introduced a modification to Eq. (59) to account for graviton dispersion, and found the modified overlap reduction functions (i.e., modifications to the Hellings-Downs curves Eq. (122)) for various values of the graviton mass. They conclude that a large number of stable pulsars (≥ 60) are required to distinguish between the massive and massless cases, and that future pulsar timing experiments could be sensitive to graviton masses of about 10−²² eV (∼ 10¹³ km). This method is competitive with some of the compact binary tests described later in Section 5.3.1 (see Table 2). In addition, since the method of Lee et al. [285] only depends on the form of the overlap reduction functions, it is advantageous in that it does not involve matched filtering (and therefore prior knowledge of the waveforms), and generally makes few assumptions about the gravitational-wave source properties.

5 Compact Binary Tests

In this section, we discuss gravitational wave tests of GR with signals emitted by compact binary systems. We begin by explaining the difference between direct and generic tests. We then proceed to describe the many direct or top-down tests and generic or bottom-up tests that have been proposed once gravitational waves are detected, including tests of the no-hair theorems. We concentrate here only on binaries composed of compact objects, such as neutron stars, black holes or other compact exotica. We will not discuss tests one could carry out with electromagnetic information from binary (or double) pulsars, as these are already described in [438]. We will also not review tests of GR with accretion disk observations, for which we refer the interested reader to [359].

5.1 Direct and generic tests

Gravitational-wave tests of Einstein’s theory can be classed into two distinct subgroups: direct tests and generic tests. Direct tests employ a top-down approach, where one starts from a particular modified gravity theory with a known action, derives the modified field equations and solves them for a particular gravitational wave-emitting system. On the other hand, generic tests adopt a bottom-up approach, where one takes a particular feature of GR and asks what type of signature its absence would leave on the gravitational-wave observable; one then asks whether the data presents a statistically-significant anomaly pointing to that particular signature.

Direct tests have by far been the traditional approach to testing GR with gravitational waves. The prototypical examples here are tests of Jordan-Fierz-Brans-Dicke theory. As described in Section 2, one can solve the modified field equations for a binary system in the post-Newtonian approximation to find a prediction for the gravitational-wave observable, as we will see in more detail later in this section. Other examples of direct tests include those concerning modified quadratic gravity models and non-commutative geometry theories.

The main advantage of such direct tests is also its main disadvantage: one has to pick a particular modified gravity theory. Because of this, one has a well-defined set of field equations that one can solve, but at the same time, one can only make predictions about that modified gravity model. Unfortunately, we currently lack a particular modified gravity theory that is particularly compelling; many modified gravity theories exist, but none possess all the criteria described in Section 2, except perhaps for the subclass of scalar-tensor theories with spontaneous scalarization. Lacking a clear alternative to GR, it is not obvious which theory one should pick. Given that the full development (from the action to the gravitational wave observable) of any particular theory can be incredibly difficult, time and computationally consuming, carrying out direct tests of all possible modified gravity models once gravitational waves are detected is clearly unfeasible.

Given this, one is led to generic tests of GR, where one asks how the absence of specific features contained in GR could impact the gravitational wave observable. For example, one can ask how such an observable would be modified if the graviton had a mass, if the gravitational interaction were Lorentz or parity violating, or if there existed large extra dimensions. From these general considerations, one can then construct a “meta”-observable, i.e., one that does not belong to a particular theory, but that interpolates over all known possibilities in a well-defined way. This model has come to be known as the parameterized post-Einsteinian framework, in analogy to the parameterized post-Newtonian scheme used to test GR in the solar system [438]. Given such a construction, one can then ask whether the data points to a statistically-significant deviation from GR.

The main advantage of generic tests is precisely that one does not have to specify a particular model, but instead one lets the data select whether it contains any statistically-significant deviations from our canonical beliefs. Such an approach is, of course, not new to physics, having most recently been successfully employed by the WMAP team [57]. The intrinsic disadvantage of this method is that, if a deviation is found, there is no one-to-one mapping between it and a particular action, but instead one has to point to a class of possible models. Of course, such a disadvantage is not that limiting, since it would provide strong hints as to what type of symmetries or properties of GR would have to be violated in an ultra-violet completion of Einstein’s theory.

5.2 Direct tests

5.2.1 Scalar-tensor theories

Let us first concentrate on Jordan-Fierz-Brans-Dicke theory, where black holes and neutron stars have been shown to exist. In this theory, the gravitational mass depends on the value of the scalar field, as Newton’s constant is effectively promoted to a function, thus leading to violations of the weak-equivalence principle [160, 434, 441]. The usual prescription for the modeling of binary systems in this theory is due to Eardley [160].⁸ He showed that such a scalar-field effect can be captured by replacing the constant inertial mass by a function of the scalar field in the distributional stress-energy tensor and then Taylor expanding about the cosmological constant value of the scalar field at spatial infinity, i.e.,

$${m_a} \rightarrow {m_a}(\phi) = {m_a}({\phi _0})\left\{{1 + {s_a}{\psi \over {{\phi _0}}} - {1 \over 2}(s_a{\prime} - s_a^2 + {s_a}){{\left({{\psi \over {{\phi _0}}}} \right)}^2} + {\mathcal O}\left[ {{{\left({{\psi \over {{\phi _0}}}} \right)}^3}} \right]} \right\},$$

(125)

where the subscript a stands for a different sources, while ψ ≡ ϕ − ϕ₀ ≪ 1 and the sensitivities s _a and s′ _a are defined by

$${s_a} \equiv - {\left[ {{{\partial (\ln {m_a})} \over {\partial (\ln G)}}} \right]_0},\quad \quad {s_a}\prime \equiv - {\left[ {{{{\partial ^2}(\ln {m_a})} \over {\partial {{(\ln G)}^2}}}} \right]_0},$$

(126)

where we remind the reader that G = 1/ϕ, the derivatives are to be taken with the baryon number held fixed and evaluated at ϕ = ϕ₀. These sensitivities encode how the gravitational mass changes due to a non-constant scalar field; one can think of them as measuring the gravitational binding energy per unit mass. The internal gravitational field of each body leads to a non-trivial variation of the scalar field, which then leads to modifications to the gravitational binding energies of the bodies. In carrying out this expansion, one assumes that the scalar field takes on a constant value at spatial infinity ϕ → ϕ₀, disallowing any homogeneous, cosmological solution to the scalar field evolution equation [Eq. (19)].

With this at hand, one can solve the massless Jordan-Fierz-Brans-Dicke modified field equations [Eq. (19)] for the non-dynamical, near-zone field of N compact objects to obtain [441]

$${\psi \over {{\phi _0}}} = {1 \over {2 + {\omega _{{\rm{BD}}}}}}\sum\limits_a {(1 - 2{s_a})} {{{m_a}} \over {{r_a}}} + \ldots ,$$

(127)

$${g_{00}} = - 1 + \sum\limits_a {\left({1 - {{{s_a}} \over {2 + {\omega _{{\rm{BD}}}}}}} \right)} {{2{m_a}} \over {{r_a}}} + \ldots ,$$

(128)

$${g_{0i}} = - 2(1 + \gamma)\sum\limits_a {{{{m_a}} \over {{r_a}}}v_a^i} + \ldots ,$$

(129)

$${g_{ij}} = {\delta _{ij}}\left[ {1 + 2\gamma \sum\limits_a {\left({1 + {{{s_a}} \over {1 + {\omega _{{\rm{BD}}}}}}} \right){{{m_a}} \over {{r_a}}}} + \ldots} \right],$$

(130)

where a runs from 1 to N, we have defined the spatial field point distance ${r_a} \equiv \vert {x^i} - x_a^i\vert$, the parameterized post-Newtonian quantity γ = (1 + ω_BD)(2 + ω_BD)−¹ and we have chosen units in which G = c = 1. This solution is obtained in a post-Newtonian expansion [75], where the ellipses represent higher-order terms in v _a /c and m _a /r _a . From such an analysis, one can also show that compact objects follow geodesics of such a spacetime, to leading order in the post-Newtonian approximation [160], except that Newton’s constant in the coupling between matter and gravity is replaced by $G \to {{\mathcal G}_{12}} = 1 - ({s_1} + {s_2} - 2{s_1}{s_2}){(2 + {\omega _{{\rm{BD}}}})^{- 1}}$, in geometric units.

As is clear from the above analysis, black-hole and neutron-star solutions in this theory generically depend on the quantities ω_BD and s _a . The former determines the strength of the correction, with the theory reducing to GR in the ω_BD → ∞ limit [164]. The latter depends on the compact object that is being studied. For neutron stars, this quantity can be computed as follows. First, neglecting scalar corrections to neutron-star structure and using the Tolman-Oppenheimer-Volkoff equation, one notes that the mass m ∝ N ∝ G^−3/2, for a fixed equation of state and central density, with N the total baryon number. Thus, using Eq. (126), one has that

$${s_a} \equiv {3 \over 2}\left[ {1 - {{\partial (\ln {m_a})} \over {\partial {{(\ln N)}_G}}}} \right],$$

(131)

where the derivative is to be taken holding G fixed. In this way, given an equation of state and central density, one can compute the gravitational mass as a function of baryon number, and from this, obtain the neutron star sensitivities. Eardley [160], Will and Zaglauer [441], and Zaglauer [474] have shown that these sensitivities are always in the range s _a ∈ (0.19, 0.3) for a soft equation of state and s _a ∈ (0.1,0.14) for a stiff one, in both cases monotonically increasing with mass in m _a ∈ (1.1,1.5) M_⊙. Recently, Gralla [202] has found a more general method to compute sensitivities is generic modified gravity theories.

What is the sensitivity of black holes in generic scalar-tensor theories? Will and Zaglauer [474] have argued that the no-hair theorems require s _a = 1/2 for all black holes, no matter what their mass or spin is. As already explained in Section 2, stationary black holes that are the byproduct of gravitational collapse (i.e., with matter that satisfies the energy conditions) in a general class of scalar-tensor theories are identical to their GR counterparts [224, 408, 159, 398].⁹ This is because the scalar field satisfies a free wave equation in vacuum, which forces the scalar field to be constant in the exterior of a stationary, asymptotically-flat spacetime, provided one neglects a homogeneous, cosmological solution. If the scalar field is to be constant, then by Eq. (127), s _a = 1/2 for a single black-hole spacetime.

Such an argument formally applies only to stationary scenarios, so one might wonder whether a similar argument holds for binary systems that are in a quasi-stationary arrangement. Will and Zaglauer [474] and Mirshekari and Will [315] extended this discussion to quasi-stationary spacetimes describing black-hole binaries to higher post-Newtonian order. They argued that the only possible deviations from ψ = 0 are due to tidal deformations of the horizon due to the companion, which are known to arise at very high order in post-Newtonian theory, $\psi = \mathcal{O}[{({m_a}/{r_a})^5}]$. Recently, Yunes et al. [465] extended this argument further by showing that to all orders in post-Newtonian theory, but in the extreme mass-ratio limit, black holes cannot have scalar hair in generic scalar-tensor theories. Finally, Healy et al. [230] have carried out a full numerical simulation of the non-linear field equations, confirming this argument in the full non-linear regime.

The activation of dynamics in the scalar field for a vacuum spacetime requires either a non-constant distribution of initial scalar field (violating the constant cosmological scalar field condition at spatial infinity) or a pure geometrical source to the scalar field evolution equation. The latter would lead to the quadratic modified gravity theories discussed in Section 2.3.3. As for the former, Horbatsch and Burgess [235] have argued that if, for example, one lets ψ = μt, which clearly satisfies □ψ = 0 in a Minkowski background,¹⁰ then a Schwarzschild black hole will acquire modifications that are proportional to μ. Alternatively, scalar hair could also be induced by spatial gradients in the scalar field [67], possibly anchored in matter at galactic scales. Such cosmological hair, however, is likely to be suppressed by a long time scale; in the example above μ must have units of inverse time, and if it is to be associated with the expansion of the universe, then it would be natural to assume $\mu = \mathcal{O}(H)$, where H is the Hubble parameter. Therefore, although such cosmological hair might have an effect on black holes in the early universe, it should not affect black hole observations at moderate to low redshifts.

Scalar field dynamics can be activated in non-vacuum spacetimes, even if initially the stars are not scalarized provided one considers a more general scalar-tensor theory, like the one introduced by Damour and Esposito-Farèse [129, 130]. As discussed in Section 2.3.1, when the conformal factor takes on a particular functional form, non-linear effects induced when the gravitational energy exceeds a certain threshold can spontaneously scalarize merging neutron stars, as demonstrated recently by Barausse, et al [51]. Therefore, neutron stars in binaries are likely to have hair in generic scalar-tensor theories, even if they start their inspiral unscalarized.

What do gravitational waves look like in Jordan-Fierz-Brans-Dicke theory? As described in Section 2.3.1, both the scalar field perturbation ψ and the new metric perturbation θ^μν satisfy a sourced wave equation [Eq. (19)], whose leading-order solution for a two-body inspiral is [436]

$${\theta ^{ij}} = 2(1 + \gamma){\mu \over R}\left({v_{12}^{ij} - {{\mathcal G}_{12}}m{{{x^{ij}}} \over {{r^3}}}} \right),$$

(132)

$${\psi \over {{\phi _0}}} = (1 - \gamma){\mu \over R}\left[ {\Gamma {{({n_i}v_{12}^i)}^2} - {{\mathcal G}_{12}}\Gamma {m \over {{r^3}}}{{({n_i}{x^i})}^2} - {m \over r}({{\mathcal G}_{12}}\Gamma + 2\Lambda) - 2S{n_i}v_{12}^i} \right],$$

(133)

where R is the distance to the detector, nⁱ is a unit vector pointing toward the detector, r is the magnitude of relative position vector ${x^i} \equiv x_1^i - x_2^i$, with $x_a^i$ the trajectory of body a, μ = m₁m₂/m is the reduced mass and m = m₁ + m₂ is the total mass, $v_{12}^i \equiv v_1^i - v_2^i$ is the relative velocity vector and we have defined the shorthands

$$\Gamma \equiv 1 - 2{{{m_1}{s_2} + {m_2}{s_1}} \over m},\quad S \equiv {s_2} - {s_1},$$

(134)

$$\Lambda \equiv {\mathcal{G}_{12}}(1 - {s_1} - {s_2}) - {(2 + {\omega _{{\text{BD}}}})^{ - 1}}[(1 - 2{s_1}){s'_2}_{} + (1 - 2{s_2}){s'_1}_{}].$$

(135)

We have also introduced multi-index notation here, such that A^ij… = AⁱA^j …. Such a solution is derived using the Lorenz gauge condition θ^μν, ν = 0 and in a post-Newtonian expansion, where we have left out subleading terms of relative order $v_{12}^2$ or m/r.

Given the new metric perturbation θ^ij, one can reconstruct the gravitational wave h^ij metric perturbation, and from this, the response function, associated with the quasi-circular inspiral of compact binaries. After using Kepler’s third law to simplify expressions $[\omega = {({\mathcal{G}_{12}}m/{r^3})^{1/2}}$, where ω is the orbital angular frequency and m is the total mass and r is the orbital separation], one finds for a ground-based L-shaped detector [102]:

$$\begin{array}{*{20}c} {h(t) = - {{{{\mathcal M}_c}} \over R}{{(2\pi {{\mathcal M}_c}F)}^{2/3}}{e^{- 2i\Phi}}\left\{{[{F_ +}(1 + {{\cos}^2}\iota) + 2i{F_ \times}\cos \iota ]\left[ {1 - {{1 - \gamma} \over 2}\left({1 + {4 \over 3}{S^2}} \right)} \right]} \right.}\\ {\left. {- {{1 - \gamma} \over 2}\Gamma {F_{\rm{b}}}{{\sin}^2}\iota} \right\} - {\eta ^{1/5}}{{{{\mathcal M}_c}} \over R}{{(2\pi {{\mathcal M}_c}F)}^{1/3}}{e^{- i\Phi}}S(1 - \gamma){F_{\rm{b}}}\sin \iota \quad \quad \quad \quad}\\ {- {{{{\mathcal M}_c}} \over R}{{(2\pi {{\mathcal M}_c}F)}^{2/3}}{{1 - \gamma} \over 2}{F_{\rm{b}}}(\Gamma + 2\Lambda),\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ \end{array}$$

(136)

where η ≡ μ/m is the symmetric mass ratio, ${\mathcal{M}_c} \equiv {\eta ^{3/5}}m$ is the chirp mass, ı is the inclination angle, and where we have used the beam-pattern functions in Eq. (58). In Eq. (136) and henceforth, we linearize all expressions in 1 − γ ≪ 1. Jordan-Fierz-Brans-Dicke theory predicts the generic excitation of three polarizations: the usual plus and cross polarizations, and a breathing, scalar mode. We see that the latter contributes to the response at two, one and zero times the orbital frequency. One should note that all of these corrections arise during the generation of gravitational waves, and not due to a propagation effect. In fact, gravitational waves travel at the speed of light (and the graviton remains massless) in standard Jordan-Fierz-Brans-Dicke theory.

The quantities Φ and F are the orbital phase and frequency respectively, which are to be found by solving the differential equation

$${{dF} \over {dt}} = (1 - \gamma){S^2}{{{\eta ^{2/5}}} \over \pi}{\mathcal M}_c^{- 2}{(2\pi {{\mathcal M}_c}F)^3} + {{48} \over {5\pi}}{\mathcal M}_c^{- 2}{(2\pi {{\mathcal M}_c}F)^{11/3}}\left[ {1 - {{1 - \gamma} \over 2}\left({1 - {{{\Gamma ^2}} \over 6} + {4 \over 3}{S^2}} \right)} \right] \ldots ,$$

(137)

where the ellipses stand for higher-order terms in the post-Newtonian approximation. In this expression, and henceforth, we have kept only the leading-order dipole term and all known post-Newtonian, GR terms. If one wished to include higher post-Newtonian-order Brans-Dicke terms, one would have to include monopole contributions as well as post-Newtonian corrections to the dipole term. The first term in Eq. (137) corresponds to dipole radiation, which is activated by the scalar mode. That is, the scalar field carries energy away from the system modifying the energy balance law to [436, 379, 440]

$${\dot E_{{\rm{BD}}}} = - {2 \over 3}{\mathcal G}_{12}^2{\eta ^2}{{{m^4}} \over {{r^4}}}(1 - \gamma){S^2} - {{32} \over 5}{\mathcal G}_{12}^2{\eta ^2}{\left({{m \over r}} \right)^5}\left[ {1 - {{1 - \gamma} \over 2}\left({1 - {{{\Gamma ^2}} \over 6}} \right)} \right] + \ldots ,$$

(138)

where the ellipses stand again for higher-order terms in the post-Newtonian approximation. Solving the frequency evolution equation perturbatively in 1/ω_BD ≪ 1, one finds

$${{256} \over 5}{{{t_c} - t} \over {{{\mathcal M}_c}}} = {u^{- 8}}\left[ {1 - {1 \over {12}}(1 - \gamma){S^2}{\eta ^{2/5}}{u^{- 2}} + \ldots} \right],$$

(139)

$$\Phi = - {1 \over {64\pi}}{\left({{{256} \over 5}{{{t_c} - t} \over {{{\mathcal M}_c}}}} \right)^{5/8}}\left[ {1 - {5 \over {224}}(1 - \gamma){S^2}{\eta ^{2/5}}{{\left({{{256} \over 5}{{{t_c} - t} \over {{{\mathcal M}_c}}}} \right)}^{1/4}} + \ldots} \right],$$

(140)

where we have defined $u \equiv {(2\pi {\mathcal{M}_c}F)^{1/3}}$. In deriving these equations, we have neglected the last term in Eq. (137), as this is a constant that can be reabsorbed into the chirp mass. Notice that since the two definitions of chirp mass differ only by a term of $\mathcal{O}(\omega _{{\rm{BD}}}^{- 1})$, the first term of Eq. (137) is not modified.

One of the main ingredients that goes into parameter estimation is the Fourier transform of the response function. This can be estimated in the stationary-phase approximation, for a simple, non-spinning, quasi-circular inspiral. In this approximation, one assumes the phase is changing much more rapidly than the amplitude [56, 125, 153, 457]. One finds [102]

$$\tilde h(f) = {{\mathcal A}_{{\rm{BD}}}}{(\pi {{\mathcal M}_c}f)^{- 7/6}}\left[ {1 - {5 \over {96}}{{{S^2}} \over {{\omega _{{\rm{BD}}}}}}{\eta ^{2/5}}{{(\pi {{\mathcal M}_c}f)}^{- 2/3}}} \right]{e^{- i\psi _{{\rm{BD}}}^{(2)}}} + {\gamma _{{\rm{BD}}}}{(\pi {{\mathcal M}_c}f)^{- 3/2}}{e^{- i\Psi _{{\rm{BD}}}^{(1)}}}$$

(141)

where we have defined the amplitudes

$${{\mathcal A}_{{\rm{BD}}}} \equiv {\left({{{5\pi} \over {96}}} \right)^{1/2}}{{{\mathcal M}_c^2} \over R}{\left[ {F_ + ^2{{(1 + {{\cos}^2}\iota)}^2} + 4F_ \times ^2{{\cos}^2}\iota - {F_ +}{F_{\rm{b}}}(1 - {{\cos}^4}\iota){\Gamma \over {{\omega _{{\rm{BD}}}}}}} \right]^{1/2}},$$

(142)

$${\gamma _{{\rm{BD}}}} \equiv - {\left({{{5\pi} \over {48}}} \right)^{1/2}}{{{\mathcal M}_c^2} \over R}{\eta ^{1/5}}{S \over {{\omega _{{\rm{BD}}}}}}{F_{\rm{b}}}\sin \iota ,$$

(143)

and the Fourier phase

$$\begin{array}{*{20}c} {\Psi _{{\rm{BD}}}^{(\ell)} = - 2\pi f{t_c} + \ell \Phi _c^{(\ell)} + {\pi \over 4} - {{3\ell} \over {256}}{{\left({{{2\pi {{\mathcal M}_c}f} \over \ell}} \right)}^{- 5/3}}\sum\limits_{n = 0}^7 {{{\left({{{2\pi {{\mathcal M}_c}f} \over \ell}} \right)}^{n/3}}(c_n^{{\rm{PN}}} + l_n^{{\rm{PN}}}\ln f)}}\\ {+ {{5\ell} \over {7168}}{{{S^2}} \over {{\omega _{{\rm{BD}}}}}}{\eta ^{2/5}}{{\left({{{2\pi {{\mathcal M}_c}f} \over \ell}} \right)}^{- 7/3}},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ \end{array}$$

(144)

where the Brans-Dicke correction is kept only to leading order in $\omega _{{\rm{BD}}}^{- 1}$ and ν, while $(c_n^{{\rm{PN}}},l_n^{{\rm{PN}}})$ are post-Newtonian GR coefficients (see, e.g., [265]). In writing the Fourier response in this way, we had to redefine the phase of coalescence via

$$\Phi _c^{(\ell)} = {\Phi _c} - {\delta _{\ell ,2}}\left\{{\arctan \left[ {{{2\cos \iota \;{F_ \times}} \over {(1 + {{\cos}^2}\iota){F_ +}}}} \right] + {\Gamma \over {{\omega _{{\rm{BD}}}}}}{{\cos \iota (1 - {{\cos}^2}\iota){F_ \times}{F_{\rm{b}}}} \over {{{(1 + {{\cos}^2}\iota)}^2}F_ + ^2 + 4{{\cos}^2}\iota F_ \times ^2}}} \right\},$$

(145)

where δ_ℓ,m is the Kronecker delta and Φ _c is the GR phase of coalescence (defined as an integration constant when the frequency diverges). Of course, in this calculation we have neglected amplitude corrections that arise purely in GR, if one were to carry out the post-Newtonian approximation to higher order.

Many studies have been carried out to determine the level at which such corrections to the waveform could be measured or constrained once a gravitational wave has been detected. The first such study was carried out by Will [436], who determined that given a LIGO detection at SNR ρ = 10 of a (1.4, 3) M_⊙ black-hole/neutron-star non-spinning, quasi-circular inspiral, one could constrain ω_BD > 10³. Scharre and Will [379] carried out a similar analysis but for a LISA detection with ρ = 10 of a (1.4, 10³) M_⊙ intermediate-mass black-hole/neutron-star, non-spinning, quasi-circular inspiral, and found that one could constrain ω_BD > 2.1 × 10⁴. Such an analysis was then repeated by Will and Yunes [440] but as a function of the classic LISA instrument. They found that the bound is independent of the LISA arm length, but inversely proportional to the LISA position noise error, if the position error noise dominates over laser shot noise. All such studies considered an angle-averaged signal that neglected the spin of either body, assumptions that were relaxed by Berti et al. [63, 64]. They carried out Monte-Carlo simulations over all signal sky positions that included spin-orbit precession to find that the projected bound with LISA deteriorates to ω_BD > 0.7 × 10⁴ for the same system and SNR. This was confirmed and extended by Yagi et al. [450], who in addition to spin-orbit precession allowed for non-circular (eccentric) inspirals. In fact, when eccentricity is included, the bound deteriorates even further to ω_BD > 0.5 × 10⁴. The same authors also found that similar gravitational-wave observations with the next-generation detector DECIGO could constrain ω_BD > 1.6 × 10⁶. Similarly, for a non-spinning neutron-star/black-hole binary, the future ground-based detector, the Einstein Telescope (ET) [361], could place constraints about 5 times stronger than the Cassini bound, as shown in [38].

All such projected constraints are to be compared with the current solar system bound of ω_BD > 4 × 10⁴ placed through the tracking of the Cassini spacecraft [73]. Table 1 presents all such bounds for ease of comparison,¹¹ normalized to an SNR of 10. As should be clear, it is unlikely that LIGO observations will be able to constrain ω_BD better than current solar system bounds. In fact, even LISA would probably not be able to do better than the Cassini bound. Table 1 also shows that the inclusion of more complexity in the waveform seems to dilute the level at which ω_BD can be constrained. This is because the inclusion of eccentricity and spin forces one to introduce more parameters in the waveform, without these modifications truly adding enough waveform complexity to break the induced degeneracies. One would then expect that the inclusion of amplitude modulation due to precession and higher harmonics should break such degeneracies, at least partially, as was found for massive black-hole binary [279, 280]. However, even then it seems reasonable to expect that only third-generation detectors will be able to constrain ω_BD beyond solar-system levels.

Table 1 Comparison of proposed tests of scalar-tensor theories.

The main reason that solar-system constraints of Jordan-Fierz-Brans-Dicke theory cannot be beaten with gravitational-wave observations is that the former are particularly well-suited to constrain weak-field deviations of GR. One might have thought that scalar-tensor theories constitute strong-field tests of Einstein’s theory, but this is not quite true, as argued in Section 2.3.1. One can see this clearly by noting that scalar-tensor theory predicts dipolar radiation, which dominates at low velocities over the GR prediction (precisely the opposite behavior that one would expect from a strong-field modification to Einstein’s theory).

However, one should note that all the above analysis considered only the inspiral phase of coalescence, usually truncating their study at the innermost stable-circular orbit. The merger and ringdown phases, where most of the gravitational wave power resides, have so far been mostly neglected. One might expect that an increase in power will be accompanied by an increase in SNR, thus allowing us to constrain ω_BD further, as this scales with 1/SNR [262]. Moreover, during merger and ringdown, dynamical strong-field gravity effects in scalar-tensor theories could affect neutron star parameters and their oscillations [395], as well as possibly induce spontaneous scalarization [51]. All of these non-linear effects could easily lead to a strengthening of projected bounds. However, to date no detailed analysis has attempted to determine how well one could constrain scalar-tensor theories using full information about the entire coalescence of a compact binary.

The subclass of scalar-tensor models described by Jordan-Fierz-Brans-Dicke theory is not the only type of model that can be constrained with gravitational-wave observations. In the extreme-mass-ratio limit, for binaries consisting of a stellar-mass compact object spiraling into a supermassive black hole, Yunes et al. [465] have recently shown that generic scalar-tensor theories reduce to either massless or massive Jordan-Fierz-Brans-Dicke theory. Of course, in this case the sensitivities need to be calculated from the equations of structure within the full scalar-tensor theory. The inclusion of a scalar field mass leads to an interesting possibility: floating orbits [94]. Such orbits arise when the small compact object experiences superradiance, leading to resonances in the scalar flux that can momentarily counteract the gravitational-wave flux, leading to a temporarily-stalled orbit that greatly modifies the orbital-phase evolution. These authors showed that if an extreme mass-ratio inspiral is detected with a template consistent with GR, this alone allows us to rule out a large region of (m _s ,ω _BD ) phase space, where m _s is the mass of the scalar (see Figure 1 in [465]). This is because if such an inspiral had gone through a resonance, a GR template would be grossly different from the signal. Such bounds are dramatically stronger than the current most stringent bound ω_BD > 4 × 10⁴ and m _s < 2.5 × 10⁻²⁰ eV obtained from Cassini measurements of the Shapiro time-delay in the solar system [20]. Even if resonances are not hit, Berti et al. [71] have estimated that second-generation ground-based interferometers could constrain the combination m _s /(ω_BD)¹/² ≲ 10⁻¹⁵ eV with the observation of gravitational waves from neutron-star/binary inspirals at an SNR of 10. These bounds can also be stronger than current constraints, especially for large scalar mass.

Lastly one should mention possible gravitational-wave constraints on other types of scalar tensor theories. Let us first consider Brans-Dicke type scalar-tensor theories, where the coupling constant is allowed to vary. Will [436] has argued that the constraints described in Table 1 go through, with the change

$${{2{{\mathcal G}_{1,2}}} \over {2 + {\omega _{{\rm{BD}}}}}} \to {{2{{\mathcal G}_{1,2}}} \over {2 + {\omega _{{\rm{BD}}}}}}{\left[ {1 + {{2{{\omega \prime}_{{\rm{BD}}}}} \over {{{(3 + 2{\omega _{{\rm{BD}}}})}^2}}}} \right]^2},$$

(146)

where ω′_BD ≡ dω_BD/dϕ,/dϕ. In the ω_BD ≫ 1 limit, this implies the replacement ${\omega _{{\text{BD}}}} \to {\omega _{{\text{BD}}}}[1 + {\omega '_{BD}}/2\omega _{{\text{BD}}}^2){]^{ - 2}}$. Of course, this assumes that there is neither a potential nor a geometric source driving the evolution of the scalar field, and is not applicable for theories where spontaneous scalarization is present [129].

Another interesting scalar-tensor theory to consider is that studied by Damour and Esposito-Farèse [129, 130]. As explained in Section 2.3.1, this theory is defined by the action of Eq. (14) with the conformal factor $A(\psi) = {e^{\beta {\psi ^2}}}$. In standard Brans-Dicke theory, only mixed binaries composed of a black hole and a neutron star lead to large deviations from GR due to dipolar emission. This is because dipole emission is proportional to the difference in sensitivities of the binary components. For neutron-star binaries with similar masses, this difference is close to zero, while for black holes it is identically zero (see Eqs. (134) and (144)). However, in the theory considered by Damour and Esposito-Farèse, when the gravitational energy is large enough, as in the very late inspiral, non-linear effects can lead to drastic modifications from the GR expectation, such as spontaneous scalarization [51]. Unfortunately, most of this happens at rather high frequency, and thus, it is not clear whether such effects are observable with current ground-based detectors.

5.2.2 Modified quadratic gravity

Black holes exist in the classes of modified quadratic gravity that have so far been considered. In non-dynamical theories (when β = 0 and the scalar-fields are constant, refer to Eq. (25)), Stein and Yunes [473] have shown that all metrics that are Ricci tensor flat are also solutions to the modified field equations (see also [360]). This is not so for dynamical theories, since then the ϑ field is sourced by curvature, leading to corrections to the field equations proportional to the Riemann tensor and its dual.

In dynamical Chern-Simons gravity, stationary and spherically-symmetric spacetimes are still described by GR solutions, but stationary and axisymmetric spacetimes are not. Instead, they are represented by [466, 272]

$$ds_{{\rm{CS}}}^2 = ds_{{\rm{Kerr}}}^2 + {5 \over 4}{{\alpha _{{\rm{CS}}}^2} \over {\beta \kappa}}{a \over {{r^4}}}\left({1 + {{12} \over 7}{M \over r} + {{27} \over {10}}{{{M^2}} \over {{r^2}}}} \right){\sin ^2}\theta \;d\theta \;dt + {\mathcal O}({a^2}/{M^2}),$$

(147)

with the scalar field

$${\vartheta _{{\rm{CS}}}} = {5 \over 8}{{{\alpha _{{\rm{CS}}}}} \over \beta}{a \over M}{{\cos \theta} \over {{r^2}}}\left({1 + {{2M} \over r} + {{18{M^2}} \over {5{r^2}}}} \right) + {\mathcal O}({a^3}/{M^3}),$$

(148)

where $ds_{{\rm{Kerr}}}^2$ is the line element of the Kerr metric and we recall that α_CS = −4α₄ in the notation of Section 2.3.3. These expressions are obtained in Boyer-Lindquist coordinates and in the small-rotation/small-coupling limit to $\mathcal{O}(a/M)$ in [466, 272] and to $\mathcal{O}({a^2}/{M^2})$ in [455]. The linear-in-spin corrections modify the frame-dragging effect and they are of 3.5 post-Newtonian order. The quadratic-in-spin corrections modify the quadrupole moment, which induces 2 post-Newtonian-order corrections to the binding energy. However, the stability of these black holes has not yet been demonstrated.

In Einstein-Dilaton-Gauss-Bonnet gravity, stationary and spherically-symmetric spacetimes are described, in the small-coupling approximation, by the line element [473]

$$ds_{{\rm{EDGB}}}^2 = - {f_{{\rm{Schw}}}}(1 + h)d{t^2} + f_{{\rm{Schw}}}^{- 1}(1 + k)d{r^2} + {r^2}d{\Omega ^2},$$

(149)

in Schwarzschild coordinates, where dΩ² is the line element on the two-sphere, f_Schw = 1 − 2M/r is the Schwarzschild factor and we have defined

$$h = {{\alpha _3^2} \over {\beta \kappa {M^4}}}{1 \over {3{f_{{\rm{Schw}}}}}}{{{M^3}} \over {{r^3}}}\left({1 + 26{M \over r} + {{66} \over 5}{{{M^2}} \over {{r^2}}} + {{96} \over 5}{{{M^3}} \over {{r^3}}} - 80{{{M^4}} \over {{r^4}}}} \right),$$

(150)

$$k = - {{\alpha _3^2} \over {\beta \kappa {M^4}}}{1 \over {{f_{{\rm{Schw}}}}}}{{{M^2}} \over {{r^2}}}\left[ {1 + {M \over r} + {{52} \over 3}{{{M^2}} \over {{r^2}}} + 2{{{M^3}} \over {{r^3}}} + {{16} \over 5}{{{M^4}} \over {{r^4}}} - {{368} \over 3}{{{M^5}} \over {{r^5}}}} \right],$$

(151)

while the corresponding scalar field is

$${\vartheta _{{\rm{EDGB}}}} = {{{\alpha _3}} \over \beta}{2 \over {Mr}}\left({1 + {M \over r} + {4 \over 3}{{{M^2}} \over {{r^2}}}} \right).$$

(152)

This solution is not restricted just to Einstein-Dilaton-Gauss-Bonnet gravity, but it is also the most general, stationary and spherically-symmetric solution in quadratic gravity. This is because all terms proportional to α₁,₂ are proportional to the Ricci tensor, which vanishes in vacuum GR, while the α₄ term does not contribute in spherical symmetry (see [473] for more details). Linear slow-rotation corrections to this solution have been found in [345]. Although the stability of these black holes has not yet been demonstrated, other dilatonic black hole solutions obtained numerically (equivalent to those in Einstein-Dilaton-Gauss-Bonnet theory in the limit of small fields) [257] have been found to be stable under axial perturbations [258, 409, 343].

Neutron stars also exist in quadratic modified gravity. In dynamical Chern-Simons gravity, the mass-radius relation remains unmodified to first order in the slow-rotation expansion, but the moment of inertia changes to this order [469, 19], while the quadrupole moment and the mass measured at spatial infinity change to quadratic order in spin [448]. This is because the mass-radius relation, to first order in slow-rotation, depends on the spherically-symmetric part of the metric, which is unmodified in dynamical Chern-Simons gravity. In Einstein-Dilaton-Gauss-Bonnet gravity, the mass-radius relation is modified [342]. As in GR, these functions must be solved for numerically and they depend on the equation of state.

Gravitational waves are also modified in quadratic modified gravity. In dynamical Chern-Simons gravity, Garfinkle et al. [190] have shown that the propagation of such waves on a Minkowski background remains unaltered, and thus, all modifications arise during the generation stage. In Einstein-Dilaton-Gauss-Bonnet theory, no such analysis of the propagation of gravitational waves has yet been carried out. Yagi et al. [447] studied the generation mechanism in both theories during the quasi-circular inspiral of comparable-mass, spinning black holes in the post-Newtonian and small-coupling approximations. They found that a standard post-Newtonian analysis fails for such theories because the assumption that black holes can be described by a distributional stress-energy tensor without any further structure fails. They also found that since black holes acquire scalar hair in these theories, and this scalar field is anchored to the curvature profiles, as black holes move, the scalar fields must follow the singularities, leading to dipole scalar-field emission.

During a quasi-circular inspiral of spinning black holes in dynamical Chern-Simons gravity, the total gravitational wave energy flux carried out to spatial infinity (equal to minus the rate of change of a binary’s binding energy by the balance law) is modified from the GR expectation to leading order by [447]

$${{\delta \dot E_{{\rm{spin}}}^{{\rm{CS}}}} \over {{{\dot E}_{{\rm{GR}}}}}} = {\zeta _4}{\eta ^{- 2}}\left\{{{{25} \over {1536}}[\bar \Delta + 2\langle (\bar \Delta \cdot {{\hat v}_{12}})\rangle ] + {{75} \over {16}}{{{a_1}{a_2}} \over {{m^2}}}\left\langle {\hat S_1^i\hat S_2^j\left({2\hat v_{12}^i\hat v_{12}^j - 2\hat n_{12}^i\hat n_{12}^j} \right)} \right\rangle} \right\}v_{12}^4,$$

(153)

due to scalar field radiation and corrections to the metric perturbation that are of magnetic-type, quadrupole form. In this equation, ${\dot E_{{\rm{GR}}}} = (32/5){\eta ^2}v_{12}^{10}$ is the leading-order GR prediction for the total energy flux, ${\zeta _4} = \alpha _4^2/(\beta \kappa {m^4})$ is the dimensionless Chern-Simons coupling parameter, ν₁₂ is the magnitude of the relative velocity with unit vector $\hat v_{12}^i,{\bar \Delta ^i} = ({m_2}/{m_1})({a_1}/m)\hat S_1^i - ({m_1}/{m_2})({a_2}/m)\hat S_2^i$, where a _A is the Kerr spin parameter of the Ath black hole and $\hat S_A^i$ is the unit vector in the direction of the spin angular momentum, the unit vector $\hat n_{12}^i$ points from body one to two, and the angle brackets stand for an average over several gravitational wave wavelengths. If the black holes are not spinning, then the correction to the scalar energy flux is greatly suppressed [447]

$${{\delta \dot E_{{\rm{no - spin}}}^{{\rm{CS}}}} \over {{{\dot E}_{{\rm{GR}}}}}} = {2 \over 3}\delta _m^2{\zeta _4}v_{12}^{14},$$

(154)

where we have defined the reduced mass difference δ _m ≡ (m₁ − m₂)/m. Notice that this is a 7 post-Newtonian-order correction, instead of a 2 post-Newtonian correction as in Eq. (153). In the non-spinning limit, the dynamical Chern-Simons correction to the metric tensor induces a 6 post-Newtonian-order correction to the gravitational energy flux [447], which is consistent with the numerical results of [344].

On the other hand, in Einstein-Dilaton-Gauss-Bonnet gravity, the corrections to the energy flux are [447]

$${{\delta \dot E_{{\rm{no - spin}}}^{{\rm{EDGB}}}} \over {{{\dot E}_{{\rm{GR}}}}}} = {5 \over {96}}{\eta ^{- 4}}\delta _m^2{\zeta _3}v_{12}^{- 2},$$

(155)

which is a −1 post-Newtonian correction. This is because the scalar field ϑ_EDGB behaves like a monopole (see Eq. (152)), and when such a scalar monopole is dragged by the black hole, it emits electric-type, dipole scalar radiation. Any hairy black hole with monopole hair will thus emit dipolar radiation, leading to −1 post-Newtonian corrections in the energy flux carried to spatial infinity.

Such modifications to the energy flux modify the rate of change of the binary’s binding energy through the balance law, Ė = −Ė_b, which in turn modify the rate of change of the gravitational wave frequency and phase, Ḟ = −Ė (dE _b /dF)−¹. For dynamical Chern-Simons gravity (when the spins are aligned with the orbital angular momentum) and for Einstein-Dilaton-Gauss-Bonnet theory (in the non-spinning case), the Fourier transform of the gravitational-wave response function in the stationary phase approximation becomes [447, 454]

$${\tilde h_{{\rm{dCS,EDGB}}}} = {\tilde h_{{\rm{GR}}}}{e^{i{\beta _{{\rm{dCS,EDGB}}}}{u^{{b_{{\rm{dCS,EDGB}}}}}},}}$$

(156)

where ${\tilde h_{{\rm{GR}}}}$ is the Fourier transform of the response in GR, $u \equiv {(\pi {\mathcal{M}_c}f)^{1/3}}$ with f the gravitational wave frequency and [447, 454]

$${\beta _{{\rm{dCS}}}} = {{1549225} \over {11812864}}{{{\zeta _4}} \over {{\eta ^{14/5}}}}\left[ {\left({1 - {{47953} \over {61969}}\eta} \right)\chi _s^2 + \left({1 - {{199923} \over {61969}}\eta} \right)\chi _a^2 - 2{\delta _m}{\chi _s}{\chi _a}} \right],\quad {b_{{\rm{dCS}}}} = - 1,$$

(157)

$${\beta _{{\rm{EDGB}}}} = - {5 \over {7168}}{\zeta _3}{\eta ^{- 18/5}}\delta _m^2,\quad {b_{{\rm{EDGB}}}} = - 7,$$

(158)

where we have defined the symmetric and antisymmetric spin combinations χ _{s, a} ≡ (a₁/m₁ ± a₂/m₂)/2. We have here neglected any possible amplitude correction, but we have included both deformations to the binding energy and Kepler’s third law, in addition to changes in the energy flux, when computing the phase correction. However, in Einstein-Dilaton-Gauss-Bonnet theory the binding energy is modified at higher post-Newtonian order, and thus, corrections to the energy flux control the modifications to the gravitational-wave response function.

From the above analysis, it should be clear that the corrections to the gravitational-wave observable in quadratic modified gravity are always proportional to the quantity ${\zeta _{3,4}} \equiv {\xi _{3,4}}/{m^4} = \alpha _{3,4}^2/(\beta \kappa {m^4})$. Thus, any measurement that is consistent with GR will allow a constraint of the form ζ_3,4 < Nδ, where N is a number of order unity, and δ is the accuracy of the measurement. Solving for the coupling constants of the theory, such a measurement would lead to $\xi _{3,4}^{1/4} < {(N\delta)^{1/4}}m$ [390]. Therefore, constraints on quadratic modified gravity will weaken for systems with larger characteristic mass. This can be understood by noticing that the corrections to the action scale with positive powers of the Riemann tensor, while this scales inversely with the mass of the object, i.e., the smaller a compact object is, the larger its curvature. Such an analysis then automatically predicts that LIGO will be able to place stronger constraints than LISA-like missions on such theories, because LIGO operates in the 100 Hz frequency band, allowing for the detection of stellar-mass inspirals, while LISA-like missions operate in the mHz band, and are limited to supermassive black-holes inspirals.

How well can these modifications be measured with gravitational-wave observations? Yagi et al. [447] predicted, based on the results of Cornish et al. [124], that a sky-averaged LIGO gravitational-wave observation with SNR of 10 of the quasi-circular inspiral of non-spinning black holes with masses (6,12)M_⊙ would allow a constraint of $\xi _3^{1/4} \underset{\sim}{<}20$ km, where we recall that ${\xi _3} = \alpha _3^2/(\beta \kappa)$. A similar sky-averaged, eLISA observation of a quasi-circular, spin-aligned black-hole inspiral with masses (10⁶,3 × 10⁶M_⊙) would constrain $\xi _3^{1/4} < {10^7}$ km [447]. The loss in constraining power comes from the fact that the constraint on ξ₃ will scale with the total mass of the binary, which is six orders of magnitude larger for space-borne sources. These constraints are not stronger than current bounds from the existence of compact objects [342] (ξ₃ < 26 km) and from the change in the orbital period of the low-mass x-ray binary A0620–00 (ξ₃ < 1.9 km) [444], but they are independent of the nature of the object and sample the theory in a different energy scale. In dynamical Chern-Simons gravity, one expects similar projected gravitational-wave constraints on ξ₄, namely $\xi _4^{1/4} < \mathcal{O}(M)$, where M is the total mass of the binary system in kilometers. Therefore, for binaries detectable with ground-based interferometers, one expects constraints of order $\xi _4^{1/4} < 10$ km. In this case, such a constraint would be roughly six orders of magnitude stronger than current LAGEOS bounds [19]. Dynamical Chern-Simons gravity cannot be constrained with binary pulsar observations, since the theory’s corrections to the post-Keplerian observables are too high post-Newtonian order, given the current observational uncertainties [448]. However, the gravitational wave constraint is more difficult to achieve in the dynamical Chern-Simons case, because the correction to the gravitational wave phase is degenerate with spin. However, Yagi et al. [454] argued that precession should break this degeneracy, and if a signal with sufficiently high SNR is observed, such bounds would be possible. One must be careful, of course, to check that the small-coupling approximation is still satisfied when saturating such a constraint [454].

5.2.3 Non-commutative geometry

Black holes exist in non-commutative geometry theories, as discussed in Section 2.3.5. What is more, the usual Schwarzschild and Kerr solutions of GR persist in these theories. This is not because such solutions have vanishing Weyl tensor, but because the quantity ∇^αβC _μανβ happens to vanish for such metrics. Similarly, one would expect that the two-body, post-Newtonian metric that describes a black-hole-binary system should also satisfy the non-commutative geometry field equations, although this has not been proven explicitly. Similarly, although neutron-star spacetimes have not yet been considered in non-commutative geometries, it is likely that if such spacetimes are stationary and satisfy the Einstein equations, they will also satisfy the modified field equations. Much more work on this is still needed to establish all of these concepts on a firmer basis.

Gravitational waves exist in non-commutative gravity. Their generation for a compact binary system in a circular orbit was analyzed by Nelson et al., in [326, 325]. They began by showing that a transverse-traceless gauge exists in this theory, although the transverse-traceless operator is slightly different from that in GR. They then proceeded to solve the modified field equations for the metric perturbation [Eq. (42)] via a Green’s function approach:

$${h^{ik}} = 2\beta \int {{{dt{\prime}} \over {\sqrt {{{(t - t{\prime})}^2} - \vert r\vert ^{2}}}}{{\ddot I}^{ik}}(t{\prime}){{\mathcal J}_1}(\beta \sqrt {{{(t - t{\prime})}^2} - \vert r\vert ^{2}})} ,$$

(159)

where recall that β² = (− 32πα₀)⁻¹ acts like a mass term, the integral is taken over the entire past light cone, ∂₁(·) is the Bessel function of the first kind, |r| is the distance from the source to the observer and the quadrupole moment is defined as usual:

$${I^{ik}} = \int {{d^3}} x\;T_{{\rm{mat}}}^{00}{x^{ik}},$$

(160)

where T⁰⁰ is the time-time component of the matter stress-energy tensor. Of course, this is only the first term in an infinite multipole expansion.

Although the integral in Eq. (159) has not yet been solved in the post-Newtonian approximation, Nelson et al. [326, 325] did solve for its time derivative to find

$${\dot h^{xx}} = - {\dot h^{yy}} = 32\beta \mu r_{12}^2{\Omega ^4}\left[ {\sin (2\phi){f_c}\left({\beta \vert r\vert ,{\Omega \over \beta}} \right) + \cos (2\phi){f_s}\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right)} \right],$$

(161a)

$${\dot h^{xy}} = - 32\beta \mu r_{12}^2{\Omega ^4}\left[ {\sin \left({2\phi - {\pi \over 2}} \right){f_c}\left({\beta \vert r\vert ,{\Omega \over \beta}} \right) + \cos \left({2\phi - {\pi \over 2}} \right){f_s}\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right)} \right],$$

(161b)

where Ω = 2πF is the orbital angular frequency and we have defined

$${f_s}(x,z) = \int\nolimits_0^\infty {{{ds} \over {\sqrt {{s^2} + {x^2}}}}{{\mathcal J}_1}(s)\sin \;(z\sqrt {{s^2} + {x^2}}),}$$

(162)

$${f_c}(x,z) = \int\nolimits_0^\infty {{{ds} \over {\sqrt {{s^2} + {x^2}}}}{{\mathcal J}_1}(s)\sin \;(z\sqrt {{s^2} + {x^2}}).}$$

(163)

and one has assumed that the binary is in the x-y plane and the observer is on the z-axis. However, if one expands these expressions about β = ∞, one recovers the GR solution to leading order, plus corrections that decay faster than 1/r. This then automatically implies that such modifications to the generation mechanism will be difficult to observe for sources at astronomical distances.

Given such a solution, one can compute the flux of energy carried by gravitational waves to spatial infinity. Stein and Yunes [400] have shown that in quadratic gravity theories, this flux is still given by

$$\dot E = {\kappa \over 2}\int {d\Omega {r^2}\left\langle {{{\dot \bar h}_{\mu \nu}}{{\dot \bar h}^{\mu \nu}}} \right\rangle} ,$$

(164)

where ${\bar h_{\mu v}}$ is the trace-reversed metric perturbation, the integral is taken over a 2-sphere at spatial infinity, and we recall that the angle brackets stand for an average over several wavelengths. Given the solution in Eq. (161), one finds that the energy flux is

$$\dot E = {9 \over {20}}{\mu ^2}r_{12}^2{\Omega ^4}{\beta ^2}\left[ {\vert r\vert ^{2}f_c^2\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right) + \vert r\vert ^{2}f_s^2\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right)} \right].$$

(165)

The asymptotic expansion of the term in between square brackets about β = ∞ is

$$\vert r\vert ^{2}\left[ {f_c^2\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right) + f_s^2\left({\beta \vert r\vert ,{{2\Omega} \over \beta}} \right)} \right] \sim \vert r\vert ^{2}\left\{{{1 \over {{\beta ^2}\vert r\vert ^{2}}}\left[ {1 + {\mathcal O}\left({{1 \over {\vert r\vert ^{2}}}} \right)} \right]} \right\},$$

(166)

which then leads to an energy flux identical to that in GR, as any subdominant term goes to zero when the 2-sphere of integration is taken to spatial infinity. In that case, there are no modifications to the rate of change of the orbital frequency. Of course, if one were not to expand about β = ∞, then the energy flux would lead to certain resonances at β = 2Ω, but the energy flux is only well-defined at future null infinity.

The above analysis was used by Nelson et al. [326, 325] to compute the rate of change of the orbital period of binary pulsars, in the hopes of using this to constrain β. Using data from the binary pulsar, they stipulated an order-of-magnitude constraint of β ≥ 10⁻¹³ m⁻¹. However, such an analysis could be revisited to relax a few assumptions used in [326, 325]. First, binary pulsar constraints on modified gravity theories require the use of at least three observables. These observables can be, for example, the rate of change of the period Ṗ, the line of nodes $\dot \Omega$ and the perihelion shift ẇ. Any one observable depends on the parameters (m₁, m₂) in GR or (m₁, m₂, β) in non-commutative geometries, where m_1,2 are the component masses. Therefore, each observable corresponds to a surface of co-dimension one, i.e., a two-dimensional surface or sheet in the three-dimensional space (m₁, m₂, β). If the binary pulsar observations are consistent with Einstein’s theory, then all sheets will intersect at some point, within a certain uncertainty volume given by the observational error. The simultaneous fitting of all these observables is what allows one to place a bound on β. The analysis of [326, 325] assumed that all binary pulsar observables were known, except for β, but degeneracies between (m₁, m₂, β) could potentially dilute constraints on these quantities. Moreover, this analysis should be generalized to eccentric and inclined binaries, since binary pulsars are known to not be on exactly circular orbits.

But perhaps the most important modification that ought to be made has to do with the calculation of the energy flux itself. The expression for Ė in Eq. (164) in terms of derivatives of the metric perturbation derives from the effective gravitational-wave stress-energy tensor, obtained by perturbatively expanding the action or the field equations and averaging over several wavelengths (the Isaacson procedure [241, 242]). In modified gravity theories, the definition of the effective stress-energy tensor in terms of the metric perturbation is usually modified, as found for example in [400]. In the case of non-commutative geometries, Stein and Yunes [400] showed that Eq. (164) still holds, provided one considers fluxes at spatial infinity. However, the analysis of [326, 325] evaluated this energy flux at a fixed distance, instead of taking the r → ∞ limit.

The balance law relates the rate of change of a binary’s binding energy with the gravitational wave flux emitted by the binary, but for it to hold, one must require the following: (i) that the binary be isolated and possess a well-defined binding energy; (ii) the total stress-energy of the spacetime satisfies a local covariant conservation law. If (ii) holds, one can use this conservation law to relate the rate of change of the volume integral of the energy density, i.e., the energy flux, to the volume integral of the current density, which can be rewritten as an integral over the boundary of the volume through Stokes’ theorem. Since in principle one can choose any integration volume, any physically-meaningful result should be independent of the surface of that volume. This is indeed the case in GR, provided one takes the integration 2-sphere to spatial infinity. Presumably, if one included all the relevant terms in Ė, without taking the limit to i⁰, one would still find a result that is independent of the surface of this two-sphere. However, this has not yet been verified. Therefore, the analysis of [326, 325] should be taken as an interesting first step toward understanding possible changes in the gravitational-wave metric perturbation in non-commutative geometries.

Not much beyond this has been done regarding non-commutative geometries and gravitational waves. In particular, one lacks a study of what the final response function would be if the gravitational-wave propagation were modified, which of course depends on the time-evolution of all propagating gravitational-wave degrees of freedom, and whether there are only the two usual dynamical degrees of freedom in the metric perturbation.

5.3 Generic tests

5.3.1 Massive graviton theories and Lorentz violation

Several massive graviton theories have been proposed to later be discarded due to ghosts, nonlinear or radiative instabilities. Thus, little work has gone into studying whether black holes and neutron stars in these theories persist and are stable, and how the generation of gravitational waves is modified. Such questions will depend on the specific massive gravity model considered, and of course, if a Vainshtein mechanism is employed, then there will not be any modifications.

However, a few generic properties of such theories can still be stated. One of them is that the non-dynamical (near-zone) gravitational field will be corrected, leading to Yukawa-like modifications to the gravitational potential [437]

$${V_{{\rm{MG}}}}(r) = {M \over r}{e^{- r/{\lambda _g}}},\quad {\rm{or}}\quad {V_{{\rm{MG}}}}(r) = {M \over r}\left({1 + {\gamma _{{\rm{MG}}}}{e^{- r/{\lambda _g}}}} \right),$$

(167)

where r is the distance from the source to a field point. For example, the latter parameterization arises in gravitational theories with compactified extra dimensions [261]. Such corrections lead to a fifth force, which then in turn allows us to place constraints on m _g through solar system observations [404]. Nobody has yet considered how such modifications to the near-zone metric could affect the binding energy of compact binaries and their associated gravitational waves.

Another generic consequence of a graviton mass is the appearance of additional propagating degrees of freedom in the gravitational wave metric perturbation. In particular, one expects scalar, longitudinal modes to be excited (see, e.g., [148]). This is, for example, the case if the action is of Pauli-Fierz type [169, 148]. Such longitudinal modes arise due to the non-vanishing of the Ψ₂ and Ψ₃ Newman-Penrose scalars, and can be associated with the presence of spin-0 particles, if the theory is of Type N in the E(2) classification [438]. The specific form of the scalar mode will depend on the structure of the modified field equations, and thus, it is not possible to generically predict its associated contribution to the response function.

A robust prediction of massive graviton theories relates to how the propagation of gravitational waves is affected. If the graviton has a mass, its velocity of propagation will differ from the speed of light, as given for example in Eq. (23). Will [437] showed that such a modification in the dispersion relation leads to a correction in the relation between the difference in time of emission Δt _e and arrival Δt _a of two gravitons:

$$\Delta {t_a} = (1 + z)\left[ {\Delta {t_e} + {D \over {2\lambda _g^2}}\left({{1 \over {f_e^2}} + {1 \over {f_e^{\prime2}}}} \right)} \right],$$

(168)

where z is the redshift, λ _g is the graviton’s Compton wavelength, f _e and $f_e^{\prime}$ are the emission frequencies of the two gravitons and D is the distance measure

$$D = {{1 + z} \over {{H_0}}}\int\nolimits_0^z {{{dz{\prime}} \over {{{(1 + {z{\prime}})}^2}{{[{\Omega _M}{{(1 + {z{\prime}})}^3} + {\Omega _\Lambda}]}^{1/2}}}},}$$

(169)

where H₀ is the present value of the Hubble parameter, Ω _M is the matter energy density and Ω_Λ is the vacuum energy density (for a zero spatial-curvature universe).

Even if the gravitational wave at the source is unmodified, the graviton time delay will leave an imprint on the Fourier transform of the response function by the time it reaches the detector [437]. This is because the Fourier phase is proportional to

$$\Psi \propto 2\pi \int\nolimits_{{f_c}}^f {[t(f) - {t_c}]df{\prime},}$$

(170)

where t is now not a constant but a function of frequency as given by Eq. (168). Carrying out the integration, one finds that the Fourier transform of the response function becomes

$${\tilde h_{{\rm{MG}}}} = {\tilde h_{{\rm{GR}}}}{e^{i{\beta _{{\rm{MG}}}}{u^{{b_{{\rm{MG}}}}}}}},$$

(171)

where ${\tilde h_{{\rm{GR}}}}$ is the Fourier transform of the response function in GR, we recall that $u = {(\pi {\mathcal{M}_c}f)^{1/3}}$ and we have defined

$${\beta _{{\rm{MG}}}} = - {{{\pi ^2}D{{\mathcal M}_c}} \over {\lambda _g^2(1 + z)}},\quad {b_{{\rm{MG}}}} = - 3.$$

(172)

Such a correction is of 1 post-Newtonian order relative to the leading-order, Newtonian term in the Fourier phase. Notice also that there are no modifications to the amplitude at all.

Numerous studies have considered possible bounds on λ _g . The most stringent solar system constraint is λ _g > 2.8 × 10¹² km and it comes from observations of Kepler’s third law (mainly Mars’ orbit), which if the graviton had a mass would be modified by the Yukawa factor in Eq. (167). Observations of the rate of decay of the period in binary pulsars [174, 53] can also be used to place the more stringent constraint λ > 1.5 × 10¹⁴ km. Similarly, studies of the stability of Kerr black holes in Pauli-Fierz theory [169] have yielded constraints of λ _g > 2.4 × 10¹³ km [88]. Gravitational-wave observations of binary systems could also be used to constrain the mass of the graviton once gravitational waves are detected. One possible test is to compare the times of arrival of coincident gravitational wave and electromagnetic signals, for example in white-dwarf binary systems. Larson and Hiscock [281] and Cutler et al. [126] estimated that one could constrain λ _g > 3 × 10¹³ km with classic LISA. Will [437] was the first to consider constraints on λ _g from gravitational-wave observations only. He considered sky-averaged, quasi-circular inspirals and found that LIGO observations of 10 M_⊙ equal-mass black holes would lead to a constraint of λ _g > 6 × 10¹² km with a Fisher analysis. Such constraints are improved to λ _g > 6.9 × 10¹⁶ km with classic LISA observations of 10⁷ M_⊙, equal-mass black holes. This increase comes about because the massive graviton correction accumulates with distance traveled (see Eq. (171)). Since classic LISA would have been able to observe sources at Gpc scales with high SNR, its constraints on λ _g would have been similarly stronger than what one would achieve with LIGO observations. Will’s study was later generalized by Will and Yunes [440], who considered how the detector characteristics affected the possible bounds on λ _g . They found that this bound scales with the square-root of the LISA arm length and inversely with the square root of the LISA acceleration noise. The initial study of Will was then expanded by Berti et al. [63], Yagi and Tanaka [450], Arun and Will [39], Stavridis and Will [399] and Berti et al. [70] to allow for non-sky-averaged responses, spin-orbit and spin-spin coupling, higher harmonics in the gravitational wave amplitude, eccentricity and multiple detections. Although the bound deteriorates on average for sources that are not optimally oriented relative to the detector, the bound improves when one includes spin couplings, higher harmonics, eccentricity, and multiple detections as the additional information and power encoded in the waveform increases, helping to break parameter degeneracies. However, all of these studies neglected the merger and ringdown phases of the coalescence, an assumption that was relaxed by Keppel and Ajith [262], leading to the strongest projected bounds λ _g > 4 × 10¹⁷ km. Moreover, all studies until then had computed bounds using a Fisher analysis prescription, an assumption relaxed by del Pozzo et al. [142], who found that a Bayesian analysis with priors consistent with solar system experiments leads to bounds stronger than Fisher ones by roughly a factor of two. All of these results are summarized in Table 2, normalizing everything to an SNR of 10. In summary, projected constraints on λ _g are generically stronger than current solar system or binary pulsar constraints by several orders of magnitude, given a LISA observation of massive black-hole mergers. Even an aLIGO observation would do better than current solar system constraints by a factor between a few [142] to an order of magnitude [262], depending on the source.

Table 2 Comparison of proposed tests of massive graviton theories. Ang. Ave. stands for an angular average over all sky locations.

Before proceeding, we should note that the correction to the propagation of gravitational waves due to a non-zero graviton mass are not exclusive to binary systems. In fact, any gravitational wave that propagates a significant distance from the source will suffer from the time delays described in this section. Binary inspirals are particularly useful as probes of this effect because one knows the functional form of the waveform, and thus, one can employ matched filtering to obtain a strong constraint. But, in principle, one could use gravitational-wave bursts from supernovae or other sources.

We have so far concentrated on massive graviton theories, but, as discussed in Section 2.3.2, there is a strong connection between such theories and Lorentz violation. Modifications to the dispersion relation are usually a result of a modification of the Lorentz group or its action in real or momentum space. For this reason, it is interesting to consider generic Lorentz-violating-inspired, modified dispersion relations of the form of Eq. (24), or more precisely [316]

$${{v_g^2} \over {{c^2}}} = 1 - A{E^{{\alpha _{{\rm{LV}}}} - 2}},$$

(173)

where α_LV controls the structure of the modification and A its amplitude. When α_LV = 0 and $A = m_g^2{c^2}$ one recovers the standard modified dispersion relation of Eq. (23). Eq. (173) introduces a generalized time delay between subsequent gravitons of the form [316]

$$\Delta {t_a} = (1 + z)\left[ {\Delta {t_e} + {{{D_{{\alpha _{{\rm{LV}}}}}}} \over {2\lambda _a^{2 - {\alpha _{{\rm{LV}}}}}}}\left({{1 \over {f_e^{2 - {\alpha _{{\rm{LV}}}}}}} - {1 \over {f_e^{\prime2 - {\alpha _{{\rm{LV}}}}}}}} \right)} \right],$$

(174)

where we have defined ${\lambda _A} \equiv {h_p}{A^{1/({\alpha _{LV}} - 2)}}$, with h _p Planck’s constant, and the generalized distance measure [316]

$${D_{{\alpha _{{\rm{LV}}}}}} = {{{{(1 + z)}^{1 - {\alpha _{{\rm{LV}}}}}}} \over {{H_0}}}\int\nolimits_0^z {{{{{(1 + z{\prime})}^{{\alpha _{{\rm{LV}}}} - 2}}} \over {{{[{\Omega _M}{{(1 + z{\prime})}^3} + {\Omega _\Lambda}]}^{1/2}}}}dz{\prime}.}$$

(175)

Such a modification then leads to the following correction to the Fourier transform of the response function [316]

$${\tilde h_{{\rm{LV}}}} = {\tilde h_{{\rm{GR}}}}{e^{i{\beta _{{\rm{LV}}}}{u^{{b_{{\rm{LV}}}}}}}},$$

(176)

where ${\tilde h_{{\rm{GR}}}}$ is the Fourier transform of the response function in GR and we have defined [316]

$$\beta _{{\rm{LV}}}^{{\alpha _{{\rm{LV}}}} \neq 1} = - {{{\pi ^{2 - {\alpha _{{\rm{LV}}}}}}} \over {1 - {\alpha _{{\rm{LV}}}}}}{{{D_{{\alpha _{{\rm{LV}}}}}}} \over {\lambda _A^{2 - {\alpha _{{\rm{LV}}}}}}}{{{\mathcal M}_c^{1 - {\alpha _{{\rm{LV}}}}}} \over {{{(1 + z)}^{1 - {\alpha _{{\rm{LV}}}}}}}},\quad b_{{\rm{LV}}}^{{\alpha _{{\rm{LV}}}} \neq 1} = 3({\alpha _{{\rm{LV}}}} - 1).$$

(177)

The case α_LV = 1 is special leading to the Fourier phase correction [316]

$$\delta {\Psi _{{\alpha _{{\rm{LV}}}} = 1}} = {{3\pi {D_1}} \over {{\lambda _A}}}\ln u.$$

(178)

The reason for this is that when α_LV = 1 the Fourier phase is proportional to the integral of 1/f, which then leads to a natural logarithm.

Different α_LV limits deserve further discussion here. Of course, when α_LV = 0, one recovers the standard massive graviton result with the mapping $\lambda _g^{- 2} \rightarrow \lambda _g^{- 2} + \lambda _A^{- 2}$. When α_LV = 2, the dispersion relation is identical to that in Eq. (23), but with a redefinition of the speed of light, and should thus be unobservable. Indeed, in this limit the correction to the Fourier phase in Eq. (176) becomes linear in frequency, and this is 100% degenerate with the time of coalescence parameter in the standard GR Fourier phase. Finally, relative to the standard GR terms that arise in the post-Newtonian expansion of the Fourier phase, the new corrections are of (1 + 3α_LV/2) post-Newtonian order. Then, if LIGO gravitational-wave observations were incapable of discerning between a 4 post-Newtonian and a 5 post-Newtonian waveform, then such observations would not be able to see the modified dispersion effect if α_LV > 2. Mirshekari et al. [316] confirmed this expectation with a Fisher analysis of non-spinning, comparable-mass quasi-circular inspirals. They found that for α_LV = 3, one can place very weak bounds on λ _A , namely A < 10⁻⁷ eV⁻¹ with a LIGO observation of a (1.4, 1.4) M_⊙ neutron star inspiral, A < 0.2 eV⁻¹ with an enhanced-LISA or NGO observation of a (10⁵, 10⁵) M_⊙ black-hole inspiral, assuming a SNR of 10 and 100 respectively. A word of caution is due here, though, as these analyses neglect any Lorentz-violating correction to the generation of gravitational waves, including the excitation of additional polarization modes. One would expect that the inclusion of such effects would only strengthen the bounds one could place on Lorentz-violating theories, but this must be done on a theory by theory basis.

5.3.2 Variable G theories and large extra dimensions

The lack of a particular Lagrangian associated with variable G theories, excluding scalar-tensor theories, or extra dimensions, makes it difficult to ascertain whether black-hole or neutron-star binaries exist in such theories. Whether this is so will depend on the particular variable G model considered. In spite of this, if such binaries do exist, the gravitational waves emitted by such systems will carry some generic modifications relative to the GR expectation.

Most current tests of the variability of Newton’s gravitational constant rely on electromagnetic observations of massive bodies, such as neutron stars. As discussed in Section 2.3.4, scalar-tensor theories can be interpreted as variable-G theories, where the variability of G is really a variation in the coupling between gravity and matter. However, Newton’s constant serves the more fundamental role of defining the relationship between geometry or length and energy, and such a relationship is not altered in most scalar-tensor theories, unless the scalar fields are allowed to vary on a cosmological scale (background, homogeneous scalar solution).

For this reason, one might wish to consider a possible temporal variation of Newton’s constant in pure vacuum spacetimes, such as in black-hole-binary inspirals. Such temporal variation would encode Ġ/G at the time and location of the merger event. Thus, once a sufficiently large number of gravitational wave events has been observed and found consistent with GR, one could reconstruct a constraint map that bounds Ġ/G along our past light cone (as a function of redshift and sky position). Since our past-light cone with gravitational waves would have extended to roughly redshift 10 with classic LISA (limited by the existence of merger events at such high redshifts), such a constraint map would have been much more complete than what one can achieve with current tests at redshift almost zero. Big Bang nucleosynthesis constraints also allow us to bound a linear drift in Ġ/G from z ≫ 10³ to zero, but these become degenerate with limits on the number of relativistic species. Moreover, these bounds exploit the huge lever-arm provided by integrating over cosmic time, but they are insensitive to local, oscillatory variations of G with periods much less than the cosmic observation time. Thus, gravitational-wave constraint maps would test one of the pillars of GR: local position invariance. This principle (encoded in the equivalence principle) states that the laws of physics (and thus the fundamental constants of nature) are the same everywhere in the universe.

Let us then promote G to a function of time of the form [468]

$$G(t,x,y,z) \approx {G_{\rm{c}}} + {\dot G_{\rm{c}}}({t_c} - t),$$

(179)

where G _c = G(t _c , x _c , y _c , z _c ) and Ġ _c = (∂G/∂t)(t _c , x _c , y _c , z _c ) are constants, and the sub-index c means that these quantities are evaluated at coalescence. Clearly, this is a Taylor expansion to first order in time and position about the coalescence event $({t_c},x_c^i)$, which is valid provided the spatial variation of G is much smaller than its temporal variation, i.e., |∇ⁱG| ≪ Ġ, and the characteristic period of the temporal variation is longer than the observation window (at most, T_obs ≤ 3 years for classic LISA), so that Ġ_cT_obs ≪ G_c. Similar parameterization of G(t) have been used to study deviations from Newton’s second law in the solar system [149, 430, 427, 411]. Thus, one can think of this modification as the consequence of some effective theory that could represent the predictions of several different alternative theories.

The promotion of Newton’s constant to a function of time changes the rate of change of the orbital frequency, which then directly impacts the gravitational-wave phase evolution. To leading order, Yunes et al. [468] find

$$\dot F = {\dot F_{{\rm{GR}}}} + {{195} \over {256\pi}}{\mathcal M}_c^{- 2}{x^3}{\eta ^{3/5}}({\dot G_c}{{\mathcal M}_c}),$$

(180)

where Ḟ_GR is the rate of change of the orbital frequency in GR, due to the emission of gravitational waves and x = (2πMF)¹/³. Such a modification to the orbital frequency evolution leads to the following modification [468] to the Fourier transform of the response function in the stationary-phase approximation [56, 125, 153, 457]

$$\tilde h = {\tilde h_{{\rm{GR}}}}(1 + {\alpha _{\dot G}}{u^{{a_{\dot G}}}}){e^{i{\beta _{\dot G}}{u^{{b_{\dot G}}}}}},$$

(181)

where we recall again that $u = {(\pi {\mathcal{M}_c}f)^{1/3}}$ and have defined the constant parameters [468]

$${\alpha _{\dot G}} = - {5 \over {512}}{{{{\dot G}_c}} \over {{G_c}}}({G_c}{{\mathcal M}_z}),\quad {\beta _{\dot G}} = - {{25} \over {65536}}{{{{\dot G}_c}} \over {{G_c}}}({G_c}{{\mathcal M}_z}),\quad a = - 8,\quad b = - 13,$$

(182)

to leading order in the post-Newtonian approximation. We note that this corresponds to a correction of −4 post-Newtonian order in the phase, relative to the leading-order term, and that the corrections are independent of the symmetric mass ratio, scaling only with the redshifted chirp mass ${\mathcal{M}_z}$. Due to this, one expects the strongest effects to be seen in low-frequency gravitational waves, such as those one could detect with LISA or DECIGO/BBO.

Given such corrections to the gravitational-wave response function, one can investigate the level to which a gravitational-wave observation consistent with GR would allow us to constrain Ġ _c . Yunes et al. [468] carried out such a study and found that for comparable-mass black-hole inspirals of total redshifted mass m _z = 10⁶ M_⊙ with LISA, one could constrain Ġ _c /G _c ≲ 10⁻⁹ yr⁻¹ or better to redshift 10 (assuming an SNR of 10³). Similar constraints are possible with observations of extreme mass-ratio inspirals. The constraint is strengthened when one considers intermediate-mass black-hole inspirals, where one would be able to achieve a bound of Ġ _c /G _c ≲ 10⁻¹¹ yr⁻¹. Although this is not as stringent as the strongest constraints from other observations (see Section 2.3.4), we recall that gravitational-wave constraints would measure local variations at the source, as opposed to local variations at zero redshift or integrated variations from the very early universe.

The effect of promoting Newton’s constant to a function of time is degenerate with several different effects. One such effect is a temporal variability of the black hole masses, i.e., if ṁ ≠ 0. Such time-variation could be induced by gravitational leakage into the bulk in certain brane-world scenarios [255], as explained in Section 2.3.4. For a black hole of mass M, the rate of black hole evaporation is given by

$${{dM} \over {dt}} = - 2.8 \times {10^{- 7}}{\left({{{1{M_ \odot}} \over M}} \right)^2}{\left({{\ell \over {10\mu {\rm{m}}}}} \right)^2}{M_ \odot}\;{\rm{y}}{{\rm{r}}^{- 1}},$$

(183)

where ℓ is the size of the large extra dimension. As expected, such a modification to a black-hole-binary inspiral will lead to a correction to the Fourier transform of the response function that is identical in structure to that of Eq. (181), but the parameters (β _Ġ , b _Ġ ) → (β_ED, b_ED) change to [449]

$${\beta _{{\rm{ED}}}} = - 8.378 \times {10^{- 8}}{\left({{\ell \over {{{\mathcal M}_c}}}} \right)^4}\left({1 - {{26} \over 3}\eta + 34{\eta ^2}} \right),\quad {b_{{\rm{ED}}}} = - 13.$$

(184)

A similar expression is found for a neutron-star/black-hole inspiral, except that the η-dependent factor in between parenthesis is corrected.

Given a gravitational-wave detection consistent with GR, one could then, in principle, place an upper bound on ℓ. Yagi et al. [449] carried out a Fisher analysis and found that a 1-year LISA detection would constrain ℓ ≤ 10³ μm with a (10, 10⁵) M_Q binary inspiral at an SNR of 100. This constraint is roughly two orders of magnitude weaker than current table-top experiment constraints [7]. Moreover, the constraint weakens somewhat for more generic inspirals, due to degeneracies between ℓ and eccentricity and spin. However, a similar observation with the third generation detector DECIGO/BBO should be able to beat current constraints by roughly one order of magnitude. Such a constraint could be strengthened by roughly one order of magnitude further, if one included the statistical enhancement in parameter estimation due to detection of order 10⁵ sources by DECIGO/BBO.

Another way to place a constraint on ℓ is to consider the effect of mass loss in the orbital dynamics [308]. When a system loses mass, the evolution of its semi-major axis a will acquire a correction of the form $\dot a = - (\dot M/M)a$, due to conservation of specific orbital angular momentum. There is then a critical semi-major axis a _c at which this correction balances the semi-major decay rate due to gravitational wave emission. McWilliams [308] argues that systems with a < a _c are then gravitational-wave dominated and will thus inspiral, while systems with a > a _c will be mass-loss dominated and will thus outspiral. If a gravitational wave arising from an inspiraling binary is detected at a given semi-major axis, then ℓ is automatically constrained to about $\mathcal{O}(20\mu {\rm{m}})$. Yagi et al. [449] extended this analysis to find that such a constraint is weaker than what one could achieve via matched filtering with a waveform in the form of Eq. (181), using the DECIGO detector.

The Ġ correction to the gravitational-wave phase evolution is also degenerate with cosmological acceleration. That is, if a gravitational wave is generated at high-redshift, its phase will be affected by the acceleration of the universe. To zeroth-order, the correction is a simple redshift of all physical scales. However, if one allows the redshift to be a function of time

$$z \sim {z_c} + {\dot z_c}(t - {t_c}) \sim {z_c} + {H_0}\left[ {{{(1 + {z_c})}^2} - {{(1 + {z_c})}^{5/2}}\Omega _M^{1/2}} \right](t - {t_c}),$$

(185)

then the observed waveform at the detector becomes structurally identical to Eq. (181) but with the parameters

$${\beta _{\dot z}} = {{25} \over {32768}}{\dot z_c}{{\mathcal M}_z},\quad {b_{\dot z}} = - 13.$$

(186)

However, using the measured values of the cosmological parameters from the WMAP analysis [271, 156], one finds that this effect is roughly 10−³ times smaller than that of a possible Ġ correction at the level of the possible bounds quoted above [468]. Of course, if one could in the future constrain Ġ better by 3 orders of magnitude, possible degeneracies with ż would become an issue.

A final possible degeneracy arises with the effect of a third body [463], accretion disk migration [267, 462] and the interaction of a binary with a circumbinary accretion disk [229]. All of these effects introduce corrections to the gravitational-wave phase of negative PN order, just like the effect of a variable gravitational constant. However, degeneracies of this type are only expected to affect a small subset of black-hole-binary observations, namely those with a third body sufficiently close to the binary, or a sufficiently massive accretion disk.

5.3.3 Parity violation

As discussed in Section 2.3.6 the simplest action to model parity violation in the gravitational interaction is given in Eq. (45). Black holes and neutron stars exist in this theory, albeit non-rotating. A generic feature of this theory is that parity violation imprints onto the propagation of gravitational waves, an effect that has been dubbed amplitude birefringence. Such birefringence is not to be confused with optical or electromagnetic birefringence, in which the gauge boson interacts with a medium and is doubly-refracted into two separate rays. In amplitude birefringence, right-(left)-circularly polarized gravitational waves are enhanced or suppressed (suppressed or enhanced) relative to the GR expectation as they propagate [245, 295, 11, 460, 17, 464].

One can understand amplitude birefringence in gravitational wave propagation due to a possible non-commutativity of the parity operator and the Hamiltonian. The Hamiltonian is the generator of time evolution, and thus, one can write [464]

$$\left({\begin{array}{*{20}c} {{h_{+ ,k}}(t)}\\ {{h_{\times ,k}}(t)}\\ \end{array}} \right) = {e^{- ift}}\left({\begin{array}{*{20}c} {{u_c}\;iv}\\ {- iv\;{u_c}}\\ \end{array}} \right)\;\left({\begin{array}{*{20}c} {{h_{+ ,k}}(0)}\\ {{h_{\times ,k}}(0)}\\ \end{array}} \right),$$

(187)

where f is the gravitational-wave angular frequency, t is time, and h_+,×,k are the gravitational wave Fourier components with wavenumber k. The quantity u _c models possible background curvature effects, with u _c = 1 for propagation on a Minkowski metric, and υ proportional to redshift for propagation on a Friedman-Robertson-Walker metric [277]. The quantity υ models possible parity-violating effects, with υ = 0 in GR. One can rewrite the above equation in terms of right and left-circular polarizations, ${h_{{\rm{R,L}}}} = ({h_ +} \pm i{h_ \times})/\sqrt 2$ to find

$$\left({\begin{array}{*{20}c} {{h_{{\rm{R}},k}}(t)}\\ {{h_{{\rm{L}},k}}(t)}\\ \end{array}} \right) = {e^{- ift}}\left({\begin{array}{*{20}c} {{u_c} + v} & 0\\ {\quad \;\;0} & {{u_c} - v}\\ \end{array}} \right)\;\left({\begin{array}{*{20}c} {{h_{{\rm{R}},k}}(0)}\\ {{h_{{\rm{L}},k}}(0)}\\ \end{array}} \right).$$

(188)

Amplitude birefringence has the effect of modifying the eigenvalues of the diagonal propagator matrix for right and left-polarized waves, with right modes amplified or suppressed and left modes suppressed or enhanced relative to GR, depending on the sign of υ. In addition to these parity-violating propagation effects, parity violation should also leave an imprint in the generation of gravitational waves. However, such effects need to be analyzed on a theory by theory basis. Moreover, the propagation-distance-independent nature of generation effects should make them easily distinguishable from the propagation effects we consider here.

The degree of parity violation, υ, can be expressed entirely in terms of the waveform observables via [464]

$$v = {1 \over 2}\left({{{{h_{\rm{R}}}} \over {h_{\rm{R}}^{{\rm{GR}}}}} - {{{h_{\rm{L}}}} \over {h_{\rm{L}}^{{\rm{GR}}}}}} \right) = {i \over 2}(\delta {\phi _{\rm{L}}} - \delta {\phi _{\rm{R}}}),$$

(189)

where $h_{{\rm{R,L}}}^{{\rm{GR}}}$ is the GR expectation for a right or left-polarized gravitational wave. In the last equality we have also introduced the notation δϕ≡ ϕ − ϕ^GR, where ϕ^GR is the GR gravitational-wave phase and

$${h_{{\rm{R,L}}}} = {h_{0,{\rm{R,L}}}}{e^{- i[\phi (\eta) - {\kappa _i}{\chi ^i}]}},$$

(190)

where h_0,R,L is a constant factor, κ is the conformal wave number and (η, χⁱ) are conformal coordinates for propagation in a Friedmann-Robertson-Walker universe. The precise form of υ will depend on the particular theory under consideration. For example, in non-dynamical Chern-Simons gravity with a field ϑ = ϑ(t), and in an expansion about z ≪ 1, one finds [464]

$$v = {\alpha \over \kappa}\pi fz\left({{{\dot \vartheta}_0} - {{{{\ddot \vartheta}_0}} \over {{H_0}}}} \right) = {\alpha \over \kappa}\pi fD\left({{H_0}{{\dot \vartheta}_0} - {{\ddot \vartheta}_0}} \right),$$

(191)

where ϑ₀ is the Chern-Simons scalar field at the detector, with α the Chern-Simons coupling constant [see, e.g., Eq. (45)], z is redshift, D is the comoving distance and H₀ is the value of the Hubble parameter today and f is the observed gravitational-wave frequency. When considering propagation on a Minkowski background, one obtains the above equation in the limit as $\dot a \rightarrow 0$, so the second term dominates, where a is the scale factor. To leading-order in a curvature expansion, the parity-violating coefficient υ will always be linear in frequency, as shown in Eq. (191). For more general parity violation and flat-spacetime propagation, υ will be proportional to (fD)f^aα, where α is a coupling constant of the theory (or a certain derivative of a coupling field) with units of [Length]^a (in the previous case, a = 0, so the correction was simply proportional to fDα, where $\alpha \propto \ddot \vartheta$).

How does such parity violation affect the waveform? By using Eq. (188) one can easily show that the Fourier transform of the response function becomes [11, 460, 464]

$${\tilde h_{{\rm{PV}}}} = ({F_ +} - i\;v\;{F_ \times}){\tilde h_ +} + ({F_ \times} + i\;v\;{F_ +}){\tilde h_ \times}.$$

(192)

Of course, one can rewrite this in terms of a real amplitude correction and a real phase correction. Expanding in υ ≪ 1 to leading order, we find [464]

$${\tilde h_{{\rm{PV}}}} = {\tilde h^{{\rm{GR}}}}(1 + v\;\delta {Q_{{\rm{PV}}}}){e^{i{v^2}\delta {\psi _{{\rm{PV}}}}}},$$

(193)

where ${\tilde h_{{\rm{GR}}}}$ is the Fourier transform of the response function in GR and we have defined

$${Q_{{\rm{GR}}}} = \sqrt {F_ + ^2{{(1 + {{\cos}^2}\iota)}^2} + 4{{\cos}^2}\iota F_ \times ^2} ,$$

(194)

$$\delta {Q_{{\rm{PV}}}} = {{2(1 + {{\cos}^2}\iota)\cos \iota \;(F_ + ^2 + F_ \times ^2)} \over {Q_{{\rm{GR}}}^2}},$$

(195)

$$\delta \psi _{\rm{PV}} = {2\cos \iota (1 + \cos^{2}\iota)(1 - \cos^{2}\iota)^{2}(F_{+}^{2} + F_{\times}^{2})F_{+}F_ {\times} \over Q_{\rm{GR}}^{4}}.$$

(196)

We see then that amplitude birefringence modifies both the amplitude and the phase of the response function. Using the non-dynamical Chern-Simons expression for υ in Eq. (191), we can rewrite Eq. (193) as [464]

$${\tilde h_{{\rm{PV}}}} = {\tilde h^{{\rm{GR}}}}(1 + {\alpha _{{\rm{PV}}}}{u^{{a_{{\rm{PV}}}}}}){e^{i{\beta _{{\rm{PV}}}}{u^{{b_{{\rm{PV}}}}}}}},$$

(197)

where we have defined the coefficients

$${\alpha _{{\rm{PV}}}} = \left({{D \over {{{\mathcal M}_c}}}} \right)\left[ {{{2(1 + {{\cos}^2}\iota)\cos \iota (F_ + ^2 + F_ \times ^2)} \over {Q_{{\rm{GR}}}^2}}} \right]{\alpha \over \kappa}\left({{H_0}{{\dot \vartheta}_0} - {{\ddot \vartheta}_0}} \right),\;\;\;{a_{{\rm{PV}}}} = 3,$$

(198)

$${\beta _{{\rm{PV}}}} = {\left({{D \over {{{\mathcal M}_c}}}} \right)^2}\left[ {{{2\cos \iota (1 + {{\cos}^2}\iota){{(1 - {{\cos}^2}\iota)}^2}(F_ + ^2 + F_ \times ^2){F_ +}{F_ \times}} \over {Q_{{\rm{GR}}}^2}}} \right]{\alpha \over \kappa}{\left({{H_0}{{\dot \vartheta}_0} - {{\ddot \vartheta}_0}} \right)^2},\;\;\;{b_{{\rm{PV}}}} = 6,$$

(199)

where we recall that $u = {(\pi {\mathcal{M}_c}f)^{1/3}}$. The phase correction corresponds to a term of 5.5 post-Newtonian order relative to the Newtonian contribution, and it scales quadratically with the Chern-Simons coupling field ϑ, which is why it was left out in [464]. The amplitude correction, on the other hand, is of 1.5 post-Newtonian order relative to the Newtonian contribution. Since both of these appear as positive-order, post-Newtonian corrections, there is a possibility of degeneracy between them and standard waveform template parameters.

Given such a modification to the response function, one can ask whether such parity violation is observable with current detectors. Alexander et al. [11, 460] argued that a gravitational wave observation with LISA would be able to constrain an integrated measure of υ, because LISA can observe massive-black-hole mergers to cosmological distances, while amplitude birefringence accumulates with distance traveled. For such an analysis, one cannot Taylor expand ϑ about its present value, and instead, one finds that

$${{1 + v} \over {1 - v}} = {e^{2\pi f\zeta (z)}},$$

(200)

where we have defined

$$\zeta (z) = {{\alpha {H_0}} \over \kappa}\int\nolimits_0^z {dz} {(1 + z)^{5/2}}\left[ {{7 \over 2}{{d\vartheta} \over {dz}} + (1 + z){{{d^2}\vartheta} \over {d{z^2}}}} \right].$$

(201)

We can solve the above equation to find

$$v = {{{e^{2\pi f\zeta (z)}} - 1} \over {1 + {e^{2\pi f\zeta (z)}}}} \sim \pi f\zeta (z),$$

(202)

where in the second equality we have linearized about υ ≪ 1 and fζ ≪ 1. Alexander et al. [11, 460] realized that this induces a time-dependent change in the inclination angle (i.e., the apparent orientation of the binary’s orbital angular momentum with respect to the observer’s line-of-sight), since the latter can be defined by the ratio h_R/h_L. They then carried out a simplified Fisher analysis and found that a LISA observation of the inspiral of two massive black holes with component masses 10⁶ M_⊙(1 + z)⁻¹ at redshift z = 15 would allow us to constrain the integrated dimensionless measure ζ < 10⁻¹⁹ to 1σ. One might worry that such an effect would be degenerate with other standard GR processes that induce similar time-dependencies, such as spin-orbit coupling. However, this time-dependence is very different from that of the parity-violating effect, and thus, Alexander et al. [11, 460] argued that these effects would be weakly correlated.

Another test of parity violation was proposed by Yunes et al. [464], who considered the coincident detection of a gravitational wave and a gamma-ray burst with the SWIFT [193] and GLAST/Fermi [97] gamma-ray satellites, and the ground-based LIGO [2] and Virgo [6] gravitational wave detectors. If the progenitor of the gamma-ray burst is a neutron-star/neutron-star merger, the gamma-ray jet is expected to be highly collimated. Therefore, an electromagnetic observation of such an event implies that the binary’s orbital angular momentum at merger must be pointing along the line of sight to Earth, leading to a strongly-circularly-polarized gravitational-wave signal and to maximal parity violation. If the gamma-ray burst observation were to provide an accurate sky location, one would be able to obtain an accurate distance measurement from the gravitational wave signal alone. Moreover, since GLAST/Fermi observations of gamma-ray bursts occur at low redshift, one would also possess a purely electromagnetic measurement of the distance to the source. Amplitude birefringence would manifest itself as a discrepancy between these two distance measurements. Therefore, if no discrepancy is found, the error ellipse on the distance measurement would allow us to place an upper limit on any possible gravitational parity violation. Because of the nature of such a test, one is constraining generic parity violation over distances of hundreds of Mpc, along the light cone on which the gravitational waves propagate.

The coincident gamma-ray burst/gravitational-wave test compares favorably to the pure LISA test, with the sensitivity to parity violation being about 2–3 orders of magnitude better in the former case. This is because, although the fractional error in the gravitational-wave distance measurement is much smaller for LISA than for LIGO, since it is inversely proportional to the SNR, the parity violating effect also depends on the gravitational-wave frequency, which is much larger for neutron-star inspirals than massive black-hole coalescences. Mathematically, the simplest models of gravitational parity violation will lead to a signature in the response function that is proportional to the gravitational-wave wavelength¹² λ_GW ∝ Df. Although the coincident test requires small distances and low SNRs (by roughly 1–2 orders of magnitude), the frequency is also larger by a factor of 5–6 orders of magnitude for the LIGO-Virgo network.

The coincident gamma-ray burst/gravitational-wave test also compares favorably to current solar system constraints. Using the motion of the LAGEOS satellites, Smith et al. [388] have placed the 1σ bound ${\dot \vartheta _0} < 2000$ km assuming ${\ddot \vartheta _0} = 0$. A similar assumption leads to a 2σ bound of ${\dot \vartheta _0} < 200$ km with a coincident gamma-ray burst/gravitational-wave observation. Moreover, the latter test also allows us to constrain the second time-derivative of the scalar field. Finally, a LISA observation would constrain the integrated history of ϑ along the past light cone on which the gravitational wave propagated. However, these tests are not as stringent as the recently proposed test by Dyda et al. [158], ${\dot \vartheta _0} < {10^{- 7}}$ km, assuming the effective theory cut-off scale is less than 10 eV and obtained by demanding that the energy density in photons created by vacuum decay over the lifetime of the universe not violate observational bounds.

The coincident test is somewhat idealistic in that there are certain astrophysical uncertainties that could hamper the degree to which we could constrain parity violation. One of the most important uncertainties relates to our knowledge of the inclination angle, as gamma-ray burst jets are not necessarily perfectly aligned with the line of sight. If the inclination angle is not known a priori, it will become degenerate with the distance in the waveform template, decreasing the accuracy to which the luminosity could be extracted from a pure gravitational wave observation by at least a factor of two. Even after taking such uncertainties into account, Yunes et al. [464] found that ${\dot \vartheta _0}$ could be constrained much better with gravitational waves than with current solar system observations.

5.3.4 Parameterized post-Einsteinian framework

One of the biggest disadvantages of a top-down or direct approach toward testing GR is that one must pick a particular theory from the beginning of the analysis. However, given the large number of possible modifications to Einstein’s theory and the lack of a particularly compelling alternative, it is entirely possible that none of these will represent the correct gravitational theory in the strong field. Thus, if one carries out a top-down approach, one will be forced to make the assumption that we, as theorists, know which modifications of gravity are possible and which are not [467]. The parameterized post-Einsteinian (ppE) approach is a framework developed specifically to alleviate such a bias by allowing the data to select the correct theory of nature through the systematic study of statistically significant anomalies.

For detection purposes, one usually expects to use match filters that are consistent with GR. But if GR happened to be wrong in the strong field, it is possible that a GR template would still extract the signal, but with the wrong parameters. That is, the best fit parameters obtained from a matched filtering analysis with GR templates will be biased by the assumption that GR is sufficiently accurate to model the entire coalescence. This fundamental bias could lead to a highly distorted image of the gravitational-wave universe. In fact, recent work by Vallisneri and Yunes [417] indicates that such fundamental bias could indeed be present in observations of neutron star inspirals, if GR is not quite the right theory in the strong-field.

One of the primary motivations for the development of the ppE scheme was to alleviate fundamental bias, and one of its most dangerous incarnations: stealth-bias [124]. If GR is not the right theory of nature, yet all our future detections are of low SNR, we may estimate the wrong parameters from a matched-filtering analysis, yet without being able to identify that there is a non-GR anomaly in the data. Thus, stealth bias is nothing but fundamental bias hidden by our limited SNR observations. Vallisneri and Yunes [417] have found that such stealth-bias is indeed possible in a certain sector of parameter space, inducing errors in parameter estimation that could be larger than statistical ones, without us being able to identify the presence of a non-GR anomaly.

5.3.4.1 Historical development

The ppE scheme was designed in close analogy with the parameterized post-Newtonian (ppN) framework, developed in the 1970s to test GR with solar system observations (see, e.g., [438] for a review). In the solar system, all direct observables depend on a single quantity, the metric, which can be obtained by a small-velocity/weak-field post-Newtonian expansion of the field equations of whatever theory one is considering. Thus, Will and Nordtvedt [331, 432, 439, 332, 433] proposed the generalization of the solar system metric into a meta-metric that could effectively interpolate between the predictions of many different alternative theories. This meta-metric depends on the product of certain Green function potentials and ppN parameters. For example, the spatial-spatial components of the meta-metric take the form

$${g_{ij}} = {\delta _{ij}}(1 + 2\gamma U + \ldots),$$

(203)

where δ _ij is the Kronecker delta, U is the Newtonian potential and γ is one of the ppN parameters, which acquires different values in different theories: γ = 1 in GR, $\gamma = (1 + {\omega _{{\rm{BD}}}}){(2 + {\omega _{{\rm{BD}}}})^{- 1}}\sim 1 - \omega _{{\rm{BD}}}^{- 1}$ in Jordan-Fierz-Brans-Dicke theory, etc. Therefore, any solar system observable could then be written in terms of system parameters, such as the masses of the planets, and the ppN parameters. An observation consistent with GR allows for a bound on these parameters, thus simultaneously constraining a large class of modified gravity theories.

The idea behind the ppE framework was to develop a formalism that allowed for similar generic tests but with gravitational waves instead of solar system observations. The first such attempt was by Arun et al. [37, 317], who considered the quasi-circular inspiral of compact objects. They suggested the waveform template family

$${\tilde h_{{\rm{PNT}}}} = {\tilde h^{{\rm{GR}}}}{e^{i{\beta _{{\rm{PNT}}u}}^{{b_{{\rm{PN}}}}}}}.$$

(204)

This waveform depends on the standard system parameters that are always present in GR waveforms, plus one theory parameter β_PNT that is to be constrained. The quantity b_PN is a number chosen by the data analyst and is restricted to be equal to one of the post-Newtonian predictions for the phase frequency exponents, i.e., b_PN = (−5, −3, −2, −1, …).

The template family in Eq. (204) allows for post-Newtonian tests of GR, i.e., consistency checks of the signal with the post-Newtonian expansion. For example, let us imagine that a gravitational wave has been detected with sufficient SNR that the chirp mass and mass ratio have been measured from the Newtonian and 1 post-Newtonian terms in the waveform phase. One can then ask whether the 1.5 post-Newtonian term in the phase is consistent with these values of chirp mass and mass ratio. Put another way, each term in the phase can be thought of as a curve in $({M_c},\eta)$ space. If GR is correct, all these curves should intersect inside some uncertainty box, just like when one tests GR with binary pulsar data. From that standpoint, these tests can be thought of as null-tests of GR and one can ask: given an event, is the data consistent with the hypothesis β_rppE = 0 for the restricted set of frequency exponents b_PN?

A Fisher and a Bayesian data analysis study of how well β_PNT could be constrained given a certain b_PN was carried out in [317, 240, 290]. Mishra et al. [317] considered the quasi-circular inspiral of non-spinning compact objects and showed that aLIGO observations would allow one to constrain β_PNT to 6% up to the 1.5 post-Newtonian order correction (b_PN = −2). Third-generation detectors, such as ET, should allow for better constraints on all post-Newtonian coefficients to roughly 2%. Clearly, the higher the value of b_PN, the worse the bound on β_PN because the power contained in higher frequency exponent terms decreases, i.e., the number of useful additional cycles induced by the ${\beta _{PNT}}{u^{b{\rm{PN}}}}$ term decreases as b_PN increases. Huwyler et al. [240] repeated this analysis but for LISA observations of the quasi-circular inspiral of black hole binaries with spin precession. They found that the inclusion of precessing spins forces one to introduce more parameters into the waveform, which dilutes information and weakens constraints on β_PNT by as much as a factor of 5. Li et al. [290] carried out a Bayesian analysis of the odds-ratio between GR and restricted ppE templates given a non-spinning, quasi-circular compact binary inspiral observation with aLIGO and adVirgo. They calculated the odds ratio for each value of b_PN listed above and then combined all of this into a single probability measure that allows one to quantify how likely the data is to be consistent with GR.

5.3.4.2 The simplest ppE model

One of the main disadvantages of the post-Newtonian template family in Eq. (204) is that it is not rooted on a theoretical understanding of modified gravity theories. To alleviate this problem, Yunes and Pretorius [467] re-considered the quasi-circular inspiral of compact objects. They proposed a more general ppE template family through generic deformations of the ℓ = 2 harmonic of the response function in Fourier space:

$$\tilde h_{{\rm{ppE,insp,1}}}^{(\ell = 2)} = {\tilde h^{{\rm{GR}}}}(1 + {\alpha _{{\rm{ppE}}}}{u^{{a_{{\rm{ppE}}}}}}){e^{i{\beta _{{\rm{ppE}}}}{u^{{b_{{\rm{ppE}}}}}}}},$$

(205)

where now (α_ppE, a_ppE, β_ppE, b_ppE) are all free parameters to be fitted by the data, in addition to the usual system parameters. This waveform family reproduces all predictions from known modified gravity theories: when (α_ppE, β_ppE) = (0,0), the waveform reduces exactly to GR, while for other parameters one reproduces the modified gravity predictions of Table 3.

Table 3 Parameters that define the deformation of the response function in a variety of modified gravity theories. The notation · means that a value for this parameter is irrelevant, as its amplitude is zero.

In Table 3, recall that S is the difference in the square of the sensitivities and ω_BD is the Brans-Dicke coupling parameter (see Section 5.2.1; we have here neglected the scalar mode), ζ₃ is the coupling parameter in Einstein-Dilaton-Gauss-Bonnet theory (see Section 5.2.2), where we have here included both the dissipative and the conservative corrections, D is a certain distance measure and λ _g is the Compton wavelength of the graviton (see Section 5.3.1), λ_LV is a distance scale at which Lorentz-violation becomes important and γ_LV is the graviton momentum exponent in the deformation of the dispersion relation (see Section 5.3.1), Ġ _c is the value of the time derivative of Newton’s constant at coalescence and dM/dt is the mass loss due to enhanced Hawking radiation in extra-dimensional scenarios (see Section 5.3.2), β_dCS is given in Eq. (157) and (α_PV,β_PV) are given in Eqs. (198) and (199) of Section 5.3.3.

Although there are only a few modified gravity theories where the leading-order post-Newtonian correction to the Fourier transform of the response function can be parameterized by post-Newtonian waveforms of Eq. (204), all such predictions can be modeled with the ppE templates of Eq. (205). In fact, only massive graviton theories, certain classes of Lorentz-violating theories and dynamical Chern-Simons gravity lead to waveform corrections that can be parameterized via Eq. (204). For example, the lack of amplitude corrections in Eq. (204) does not allow for tests of gravitational parity violation or non-dynamical Chern-Simons gravity.

However, this does not imply that Eq. (205) can parameterize all possible deformations of GR. First, Eq. (205) can be understood as a single-parameter deformation away from Einstein’s theory. If the correct theory of nature happens to be a deformation of GR with several parameters (e.g., several coupling constants, mass terms, potentials, etc.), then Eq. (205) will only be able to parameterize the one that leads to the most useful cycles. This was recently verified by Sampson et al. [376]. Second, Eq. (205) assumes that the modification can be represented as a power series in velocity, with possibly non-integer values. Such an assumption does not allow for possible logarithmic terms, which are known to arise due to non-linear memory interactions at sufficiently-high post-Newtonian order. It also does not allow for interactions that are screened, e.g., in theories with massive degrees of freedom. Nonetheless, the parameterization in Eq. (205) will still be able to signal that the detection is not a pure Einstein event, at the cost of biasing their true value.

The inspiral ppE model of Eq. (205) is motivated not only from examples of modified gravity predictions, but from generic modifications to the physical quantities that drive the inspiral: the binding energy or Hamiltonian and the radiation-reaction force or the fluxes of the constants of the motion. Yunes and Pretorius [467] and Chatziioannou et al. [102] considered generic modifications of the form

$$E = {\mu \over 2}{m \over r}\left[ {1 + {A_{{\rm{ppE}}}}{{\left({{m \over r}} \right)}^{{p_{{\rm{ppE}}}}}}} \right],\quad \quad \dot E = {\dot E_{{\rm{GR}}}}\left[ {1 + {B_{{\rm{ppE}}}}{{\left({{m \over r}} \right)}^{{q_{{\rm{ppE}}}}}}} \right],$$

(206)

where (p, q) ∈ ℤ, since otherwise one would lose analyticity in the limit of zero velocities for circular inspirals, and where (A, B) are parameters that depend on the modified gravity theory and, in principle, could depend on dimensionless quantities like the symmetric mass ratio. Such modifications lead to the following corrections to the SPA Fourier transform of the ℓ = 2 time-domain response function for a quasi-circular binary inspiral template (to leading order in the deformations and in post-Newtonian theory)

$$\begin{array}{*{20}c} {\tilde h = A{{(\pi {{\mathcal M}_c}f)}^{- 7/6}}{e^{- i{\psi _{{\rm{GR}}}}}}\left[ {1 - {{{B_{{\rm{ppE}}}}} \over 2}{\eta ^{- 2{q_{{\rm{ppE}}}}/5}}} \right.{{(\pi {{\mathcal M}_c}f)}^{2{q_{{\rm{ppE}}}}}}\quad \quad}\\ {\left. {+ {{{A_{{\rm{ppE}}}}} \over 6}(6 + 4{p_{{\rm{ppE}}}} - 5p_{{\rm{ppE}}}^2){\eta ^{- 2{p_{{\rm{ppE}}}}/5}}{{(\pi {{\mathcal M}_c}f)}^{2{p_{{\rm{ppE}}}}}}} \right]{e^{- i\delta {\Psi _{{\rm{ppE}}}}}},}\\ \end{array}$$

(207)

$$\begin{array}{*{20}c} {\delta {\Psi _{{\rm{ppE}}}} = {5 \over {32}}A{{5p_{{\rm{ppE}}}^2 - 2{p_{{\rm{ppE}}}} - 6} \over {(4 - {p_{{\rm{ppE}}}})(5 - {p_{{\rm{ppE}}}})}}{\eta ^{- 2{p_{{\rm{ppE}}}}/5}}{{(\pi {{\mathcal M}_c}f)}^{2{p_{{\rm{ppE}}}} - 5}}\quad \quad}\\ {+ {{15} \over {32}}{{{B_{{\rm{ppE}}}}} \over {(4 - {q_{{\rm{ppE}}}})(5 - 2{q_{{\rm{ppE}}}})}}{\eta ^{- 2{q_{{\rm{ppE}}}}/5}}{{(\pi {{\mathcal M}_c}f)}^{2{q_{{\rm{ppE}}}} - 5}}.}\\ \end{array}$$

(208)

Of course, usually one of these two modifications dominates over the other, depending on whether q > p or p < q. In Jordan-Fierz-Brans-Dicke theory, for example, the radiation-reaction correction dominates as q < p. If, in addition to these modifications in the generation of gravitational waves, one also allows for modifications in the propagation, one is then led to the following template family [102]

$$\tilde h_{{\rm{ppE}},{\rm{insp,2}}}^{(\ell = 2)} = {\mathcal A}\;{(\pi {{\mathcal M}_c}f)^{- 7/6}}{e^{- i{\Psi _{{\rm{GR}}}}}}\left[ {1 + c{\beta _{{\rm{ppE}}}}{{(\pi {{\mathcal M}_c}f)}^{{b_{{\rm{ppE}}}}/3 + 5/3}}} \right]{e^{2i{\beta _{{\rm{ppE}}}}{u^{{b_{{\rm{ppE}}}}}}}}{e^{i{k_{{\rm{ppE}}}}{u^{{k_{{\rm{ppE}}}}}}}}.$$

(209)

Here (b_ppE,β_ppE) and (k_ppE, κ_ppE) are ppE parameters induced by modifications to the generation and propagation of gravitational waves respectively, where still (b_ppE,k_ppE) ∈ ℤ, while c is fully determined by the former set via

$${c_{{\rm{cons}}}} = - {{16} \over {15}}{{(3 - b)(42b + 61 + 5{b^2})} \over {5{b^2} + 46b + 81}},$$

(210)

if the modifications to the binding energy dominate,

$${c_{{\rm{diss}}}} = - {{16} \over {15}}(3 - b)b,$$

(211)

if the modifications to the energy flux dominate, or

$${c_{{\rm{both}}}} = - {{32} \over {15}}{{b(3 - b)(44b + 71 + 5{b^2})} \over {5{b^2} + 46b + 81}},$$

(212)

if both corrections enter at the same post-Newtonian order. Noticing again that if only a single term in the phase correction dominates in the post-Newtonian approximation (or both will enter at the same post-Newtonian order), one can map Eq. (207) to Eq. (205) by a suitable redefinition of constants.

5.3.4.3 More complex ppE models

Of course, one can introduce more ppE parameters to increase the complexity of the waveform family, and thus, Eq. (205) should be thought of as a minimal choice. In fact, one expects any modified theory of gravity to introduce not just a single parametric modification to the amplitude and the phase of the signal, but two new functional degrees of freedom:

$${\alpha _{{\rm{ppE}}}}{u^{{a_{{\rm{ppE}}}}}} \rightarrow \delta {A_{{\rm{ppE}}}}({\lambda ^a},{\theta ^a};u),\quad {\beta _{{\rm{ppE}}}}{u^{{b_{{\rm{ppE}}}}}} \rightarrow \delta {\Psi _{{\rm{ppE}}}}({\lambda ^a},{\theta ^a};u),$$

(213)

where these functions will depend on the frequency u, as well as on system parameters λ^a and theory parameters θ^a. In a post-Newtonian expansion, one expects these functions to reduce to leading-order on the left-hand sides of Eqs. (213), but also to acquire post-Newtonian corrections of the form

$$\delta {A_{{\rm{ppE}}}}({\lambda ^a},{\theta ^a};u) = {\alpha _{{\rm{ppE}}}}({\lambda ^a},{\theta ^a}){u^{{a_{{\rm{ppE}}}}}}\sum\limits_n {{\alpha _{n,{\rm{ppE}}}}} ({\lambda ^a},{\theta ^a}){u^n},$$

(214)

$$\delta {\Psi _{{\rm{ppE}}}}({\lambda ^a},{\theta ^a};u) = {\beta _{{\rm{ppE}}}}({\lambda ^a},{\theta ^a}){u^{{b_{{\rm{ppE}}}}}}\sum\limits_n {{\beta _{n,{\rm{ppE}}}}} ({\lambda ^a},{\theta ^a}){u^n},$$

(215)

where here the structure of the series is assumed to be of the form uⁿ with u > 0. Such a model, also suggested by Yunes and Pretorius [467], would introduce too many new parameters that would dilute the information content of the waveform model. Recently, Sampson et al. [376] demonstrated that the simplest ppE model of Eq. (205) suffices to signal a deviation from GR, even if the injection contains three terms in the phase.

In fact, this is precisely one of the most important differences between the ppE and ppN frameworks. In ppN, it does not matter how many ppN parameters are introduced, because the observations are of very high SNR, and thus, templates are not needed to extract the signal from the noise. On the other hand, in gravitational wave astrophysics, templates are essential to make detections and do parameter estimation. Spurious parameters in these templates that are not needed to match the signal will deteriorate the accuracy to which all parameters can be measured because of an Occam penalty. Thus, in gravitational wave astrophysics and data analysis one wishes to minimize the number of theory parameters when testing GR [124, 376]. One must then find a balance between the number of additional theory parameters to introduce and the amount of bias contained in the templates.

At this junction, one must emphasize that frequency exponents in the amplitude and phase correction were above assumed to be integers, i.e., (a_ppE, b_ppE,n) ∈ ℤ. This must be the case if these corrections arise due to modifications that can be represented as integer powers of the momenta or velocity. We are not aware of any theory that predicts corrections proportional to fractional powers of the velocity for circular inspirals. Moreover, one can show that theories that introduce non-integer powers of the velocity into the equations of motion will lead to issues with analyticity at zero velocity and a breakdown of uniqueness of solutions [102]. In spite of this, modified theories can introduce logarithmic terms, that for example enter at high post-Newtonian order in GR due to non-linear propagation effects (see, e.g., [75] and references therein). Moreover, certain modified gravity theories introduce screened modifications that become “active” only above a certain frequency. Such effects would be modeled through a Heaviside function, for example needed when dealing with massive Brans-Dicke gravity [147, 94, 20, 465]. However, even these non-polynomial injections would be detectable with the simplest ppE model. In essence, one finds similar results as if one were trying to fit a 3-parameter injection with the simplest 1-parameter ppE model [376].

Of course, one can also generalize the inspiral ppE waveform families to more general orbits, for example through the inclusion of spins aligned or counter-aligned with the orbital angular momentum. More general inspirals would still lead to waveform families of the form of Eq. (205) or (209), but where the parameters (α_ppE, β_ppE) would now depend on the mass ratio, mass difference, and the spin parameters of the black holes. With a single detection, one cannot break the degeneracy in the ppE parameters and separately fit for its system parameter dependencies. However, given multiple detections one should be able to break such a degeneracy, at least to a certain degree [124]. Such breaking of degeneracies begins to become possible when the number of detections exceeds the number of additional parameters required to capture the physical parameter dependencies of (α_ppE, β_ppE).

PpE waveforms can be extended to account for the merger and ringdown phases of coalescence. Yunes and Pretorius have suggested the following template family to account for this as well [467]

$$\tilde h_{{\rm{ppE,full}}}^{(\ell = 2)} = \left\{{\begin{array}{*{20}c} {{{\tilde h}_{{\rm{ppE}}}}\quad \quad \quad \quad \quad \quad \quad \quad} & {f < {f_{{\rm{IM}}}},\quad \quad \quad \quad}\\ {\gamma {u^c}{e^{i(\delta + \epsilon u)}}\quad \quad \quad \quad \quad \quad} & {{f_{{\rm{IM}}}} < f < {f_{{\rm{MRD}}}},\quad}\\ {\zeta {\tau \over {1 + 4{\pi ^2}{\tau ^2}\kappa {{(f - {f_{{\rm{RD}}}})}^d}}}} & {f > {f_{{\rm{MRD}}}},\quad \quad \quad \quad}\\ \end{array}} \right.$$

(216)

where the subscripts IM and MRD stand for inspiral merger and merger ringdown, respectively. The merger phase (f_IM < f < f_MRD) is modeled here as an interpolating region between the inspiral and ringdown, where the merger parameters (γ, δ) are set by continuity and differentiability, and the ppE merger parameters (c, ϵ) should be fit for. In the ringdown phase (f > f_MRD), the response function is modeled as a single-mode generalized Lorentzian, with real and imaginary dominant frequencies f_RD and τ, ringdown parameter ζ also set by continuity and differentiability, and the ppE ringdown parameters (κ, d) are to be fit for. The transition frequencies (f_IM, f_MRD) can either be treated as ppE parameters or set via some physical criteria, such as at light-ring frequency and the fundamental ringdown frequency, respectively.

Recently, there has been effort to generalize the ppE templates to allow for the excitation of non-GR gravitational-wave polarizations. Modifications to only the two GR polarizations map to corrections to terms in the time-domain Fourier transform that are proportional to the ℓ = 2 harmonic of the orbital phase. However, Arun suggested that if additional polarizations are present, other terms proportional to the ℓ = 0 and ℓ = 1 harmonic will also arise [36]. Chatziioannou, Yunes and Cornish [102] have found that the presence of such harmonics can be captured through the more complete single-detector template family

$$\begin{array}{*{20}c} {\tilde h_{{\rm{ppE}},{\rm{insp}}}^{{\rm{all}}\;\ell}(f) = {\mathcal A}\;{{(\pi {{\mathcal M}_c}f)}^{- 7/6}}{e^{- i{\Psi ^{(2)}}{\rm{GR}}}}\left[ {1 + c\;{\beta _{{\rm{ppE}}}}{{(\pi {{\mathcal M}_c}f)}^{{b_{{\rm{ppE}}}}/3 + 5/3}}} \right]{e^{2i{\beta _{{\rm{ppE}}}}u_2^{{b_{{\rm{ppE}}}}}}}{e^{2i{k_{{\rm{ppE}}}}u_2^{{\kappa _{{\rm{ppE}}}}}}}}\\ {+ {\gamma _{{\rm{ppE}}}}\;u_1^{- 9/2}{e^{- i{\Psi ^{{{(1)}_{{\rm{GR}}}}}}}}{e^{i{\beta _{{\rm{ppE}}}}u_1^{{b_{{\rm{ppE}}}}}}}{e^{2i{k_{{\rm{ppE}}}}u_1^{{\kappa _{{\rm{ppE}}}}}}},\quad \quad \quad \quad \quad \quad \quad \quad \quad}\\ \end{array}$$

(217)

$$\Psi _{{\rm{GR}}}^{(\ell)} = - 2\pi f{t_c} + \ell \Phi _c^{(\ell)} + {\pi \over 4} - {{3\ell} \over {256u_\ell ^5}}\sum\limits_{n = 0}^7 {u_\ell ^{n/3}(c_n^{{\rm{PN}}} + l_n^{{\rm{PN}}}\ln {u_\ell}),}$$

(218)

where we have defined ${u_\ell} = {(2\pi {M_c}f/\ell)^{1/3}}$.

The ppE theory parameters are now $\vec \theta = ({b_{{\rm{ppE}}}},{\beta _{{\rm{ppE}}}},{k_{{\rm{ppE}}}},{\kappa _{{\rm{ppE}}}},{\gamma _{{\rm{ppE}}}},\Phi _c^{(1)})$. Of course, one may ignore (k_ppE, κ_ppE) altogether, if one wishes to ignore propagation effects. Such a parameterization recovers the predictions of Jordan-Fierz-Brans-Dicke theory for a single-detector response function [102], as well as Arun’s analysis for generic dipole radiation [36].

One might worry that the corrections introduced by the ℓ = 1 harmonic, i.e., terms proportional to γ_ppE in Eq. (217), will be degenerate with post-Newtonian corrections to the amplitude of the ℓ = 2 mode (not displayed in Eq. (217)). However, this is clearly not the case, as the latter scale as ${(\pi {\mathcal{M}_c}f)^{- 7/6 + n/3}}$ with n an integer greater than 0, while the ℓ = 1 mode is proportional to ${(\pi {M_c}f)^{- 3/2}}$, which would correspond to a (−0.5) post-Newtonian order correction, i.e., n = − 1. On the other hand, the ppE amplitude corrections to the ℓ = 2 mode, i.e., terms proportional to β_ppE in the amplitude of Eq. (217), can be degenerate with such post-Newtonian corrections when b_ppE is an integer greater than −4.

5.3.4.4 Applications of the ppE formalism

The two models in Eq. (205) and (209) answer different questions. The latter contains a stronger prior (that ppE frequency exponents be integers), and thus, it is ideal for fitting a particular set of theoretical models. On the other hand, Eq. (205) with continuous ppE frequency exponents allows one to search for generic deviations that are statistically significant, without imposing such theoretical priors. That is, if a deviation from GR is present, then Eq. (205) is more likely to be able to fit it, than Eq. (209). If one prioritizes the introduction of the least number of new parameters, Eq. (205) with (a_ppE, b_ppE) ∈ ℝ can still recover deviations from GR, even if the latter cannot be represented as a correction proportional to an integer power of velocity.

Given these ppE waveforms, how should they be used in a data analysis pipeline? The main idea behind the ppE framework is to match filter or perform Bayesian statistics with ppE enhanced template banks to allow the data to select the best-fit values of θ^a. As discussed in [467, 124] and then later in [290], one might wish to first run detection searches with GR template banks, and then, once a signal has been found, do a Bayesian model selection analysis with ppE templates. The first such Bayesian analysis was carried out by Cornish et al. [124], who concluded that an aLIGO detection at SNR of 20 for a quasi-circular, non-spinning black-hole inspiral would allow us to constrain α_ppE and β_ppE much better than existent constraints for sufficiently strong-field corrections, e.g., b_ppE > −5. This is because for lower values of the frequency exponents, the corrections to the waveform are weak-field and better constrained with binary pulsar observations [461]. The large statistical study of Li et al. [290] uses a reduced set of ppE waveforms and investigates our ability to detect deviations of GR when considering a catalogue of aLIGO/adVirgo detections. Of course, the disadvantage of such a pipeline is that it requires a first detection, and if the gravitational interaction is too different from GR’s prediction, it is possible that a search with GR templates might miss the signal all together; we deem this possibility to be less likely.

A built-in problem with the ppE and the ppN formalisms is that if a non-zero ppE or ppN parameter is detected, then one cannot necessarily map it back to a particular modified gravity action. On the contrary, as suggested in Table 3, there can be more than one theory that predicts structurally-similar corrections to the Fourier transform of the response function. For example, both Jordan-Fierz-Brans-Dicke theory and the dissipative sector of Einstein-Dilaton-Gauss-Bonnet theory predict the same type of leading-order correction to the waveform phase. However, if a given ppE parameter is measured to be non-zero, this could provide very useful information as to the type of correction that should be investigated further at the level of the action. The information that could be extracted is presented in Table 4, which is derived from knowledge of the type of corrections that lead to Table 3.

Table 4 Interpretation of non-zero ppE parameters.

Moreover, if a follow-up search is done with the ppE model in Eq. (209), one could infer whether the correction is one due to modifications to the generation or the propagation of gravitational waves. In this way, a non-zero ppE detection could inform theories of what type of GR modification is preferred by nature.

5.3.4.5 Degeneracies

However, much care must be taken to avoid confusing a ppE theory modification with some other systematic, such as an astrophysical, a mismodeling or an instrumental effect. Instrumental effects can be easily remedied by requiring that several instruments, with presumably unrelated instrumental systematics, independently derive a posterior probability for (α_ppE, β_ppE) that peaks away from zero. Astrophysical uncertainties can also be alleviated by requiring that different events lead to the same posteriors for ppE parameters (after breaking degeneracies with system parameters). However, astrophysically there are a limited number of scenarios that could lead to corrections in the waveforms that are large enough to interfere with these tests. For comparable-mass-ratio inspirals, this is usually not a problem as the inertia of each binary component is too large for any astrophysical environment to affect the orbital trajectory [229]. Magnetohydrodynamic effects could affect the merger of neutron-star binaries, but this usually occurs outside of the sensitivity band of ground-based interferometers. However, in extreme-mass-ratio inspirals the small compact object can be easily nudged away by astrophysical effects, such as the presence of an accretion disk [462, 267] or a third supermassive black hole [463]. However, these astrophysical effects present the interesting feature that they correct the waveform in a form similar to Eq. (205) but with b_ppE < −5. This is because the larger the orbital separation, the stronger the perturbations of the astrophysical environment, either because the compact object gets closer to the third body or because it leaves the inner edge of the accretion disk and the disk density increases with separation. Such effects, however, are not likely to be present in all sources observed, as few extreme-mass-ratio inspirals are expected to be embedded in an accretion disk or sufficiently close to a third body (≲ 0.1 pc) for the latter to have an effect on the waveform.

Perhaps the most dangerous systematic is mismodeling, which is due to the use of approximation schemes when constructing waveform templates. For example, in the inspiral one uses the post-Newtonian approximation series, expanding and truncating the waveform at a given power of orbital velocity. Moreover, neutron stars are usually modeled as test-particles (with a Dirac distributional density profile), when in reality they have a finite radius, which will depend on its equation of state. Such finite-size effects enter at 5 post-Newtonian order (the effacement principle [227, 128]), but with a post-Newtonian coefficient that can be rather large [320, 72, 175]. Ignorance of the post-Newtonian series beyond 3 post-Newtonian order can lead to systematics in the determination of physical parameters and possibly also to confusion when carrying out ppE-like tests. Much more work is needed to determine the systems and SNRs for which such systematics are truly a problem.

5.3.5 Searching for non-tensorial gravitational-wave polarizations

Another way to search for generic deviations from GR is to ask whether any gravitational-wave signal detected contains more than the two traditional polarizations expected in GR. A general approach to answer this question is through null streams, as discussed in Section 4.3. This concept was first studied by Gürsel and Tinto [212] and later by Chatterji et al. [101] with the aim to separate false-alarm events from real detections. Chatziioannou et al. [102] proposed the extension of the idea of null streams to develop null tests of GR, which was proposed using stochastic gravitational wave backgrounds in [329, 330] and recently implemented in [228] to reconstruct the independent polarization modes in time-series data of a ground-based detector network.

Given a gravitational-wave detection, one can ask whether the data is consistent with two polarizations by constructing a null stream through the combination of data streams from 3 or more detectors. As explained in Section 4.3, such a null stream should be consistent with noise in GR, while it would present a systematic deviation from noise if the gravitational wave metric perturbation possessed more than two polarizations. Notice that such a test would not require a template; if one were parametrically constructed, such as in [102], more powerful null tests could be applied to such a null steam. In the future, we expect several gravitational wave detectors to be online: the two aLIGO ones in the United States, adVIRGO in Italy, LIGO-India in India, and KAGRA in Japan. Given a gravitational-wave observation that is detected by all five detectors, one can then construct three enhanced GR null streams, each with power in a signal null direction.

5.3.6 I-Love-Q tests

Neutron stars in the slow-rotation limit can be characterized by their mass and radius (to zeroth-order in spin), by their moment of inertia (to first-order in spin), and by their quadrupole moment and Love numbers (to second-order in spin). One may expect these quantities to be quite sensitive to the neutron star’s internal structure, which can be parameterized by its equation of state, i.e., the relation between its internal pressure and its internal energy density. Since the equation of state cannot be well-constrained at super-nuclear densities in the laboratory, one is left with a variety of possibilities that predict different neutron-star mass-radius relations.

Recently, however, Yagi and Yunes [453, 452] have demonstrated that there are relations between the moment of inertia (I), the Love numbers (λ), and the quadrupole moment (Q), the I-Love-Q relations that are essentially insensitive to the equation of state. Figure 5 shows two of these relations (the normalized I-Love and Q-Love relations — see caption) for a variety of equations of state, including APR [10], SLy [150, 385], Lattimer-Swesty with nuclear incompressibility of 220 MeV (LS220) [283, 335], Shen [382, 383, 335], the latter two with temperature of 0.01 MeV and an electron fraction of 30%, and polytropic equations of state with indices of n = 0.6, 0.8 and 1.¹³ The bottom panels show the difference between the numerical results and the analytical, fitting curve. Observe that all equations of state lead to the same I-Love and Q-Love relations, with discrepancies smaller than 1% for realistic neutron-star masses. These results have recently been verified in [304] through the post-Newtonian-Affine approach [168, 305], which proves the I-Love-Q relations hold not only during the inspiral, but also close to plunge and merger.

Figure 5

Top: Fitting curves (solid curve) and numerical results (points) of the universal I-Love (left) and Q-Love (right) relations for various equations of state, normalized as $\bar I = I/M_{{\rm{NS}}}^3$, ${\bar \lambda ^{({\rm{tid}})}} = {\lambda ^{({\rm{tid}})}}/M_{{\rm{NS}}}^5$ and $\bar Q = - {Q^{({\rm{rot}})}}/[M_{{\rm{NS}}}^3{(S/M_{{\rm{NS}}}^2)^2}]$, M_NS is the neutron-star mass, λ^(tid) is the tidal Love number, Q^(rot) is the rotation-induced quadrupole moment, and S is the magnitude of the neutron-star spin angular momentum. The neutron-star central density is the parameter varied along each curve, or equivalently the neutron-star compactness. The top axis shows the neutron star mass for the APR equation of state, with the vertical dashed line showing M_NS = 1 M_⊙. Bottom: Relative fractional errors between the fitting curve and the numerical results. Observe that these relations are essentially independent of the equation of state, with loss of universality at the 1% level. Image reproduced by permission from [452], copyright by APS.

Given the independent measurement of any two members of the I-Love-Q trio, one could carry out a (null) model-independent and equation-of-state-independent test of GR [453, 452]. For example, assume that electromagnetic observations of the binary pulsar J0737-3039 have measured the moment of inertia to 10% accuracy [282, 273, 274]. The slow-rotation approximation is perfectly valid for this binary pulsar, due to its relatively long spin period. Assume further that a gravitational-wave observation of a neutron-star-binary inspiral, with individual masses similar to that of the primary in J0737-3039, manages to measure the neutron star tidal Love number to 60% accuracy [453, 452]. These observations then lead to an error box in the I-Love plane, which must contain the curve in the left-panel of Figure 5.

A similar test could be carried out by using data from only binary pulsar observations or only gravitational wave detections. In the case of the latter, one would have to simultaneously measure or constrain the value of the quadrupole moment and the Love number, since the moment of inertia is not measurable with gravitational wave observations. In the case of the former, one would have to extract the moment of inertia and the quadrupole moment, the latter of which will be difficult to measure. Therefore, the combination of electromagnetic and gravitational wave observations would be the ideal way to carry out such tests.

Such a test of GR, of course, is powerful only as long as modified gravity theories predict I-Love-Q relations that are not degenerated with the general relativistic ones. Yagi and Yunes [453, 452] investigated such a relation in dynamical Chern-Simons gravity to find that such degeneracy is only present in the limit ζ_CS → 0. That is, for any finite value of ζ_CS, the dynamical Chern-Simons I-Love-Q relation differs from that of GR, with the distance to the GR expectation increasing for larger ζ_CS. Yagi and Yunes [453, 452] predicted that a test similar to the one described above could constrain dynamical Chern-Simons gravity to roughly $\xi _{{\rm{CS}}}^{1/4} < 10{M_{NS}}\sim 15$ km, where recall that ${\xi _{{\rm{CS}}}} = \alpha _{{\rm{CS}}}^2/(\beta \kappa)$.

The test described above, of course, only holds provided the I-Love-Q relations are valid, which in turn depends on the assumptions made in deriving them. In particular, Yagi and Yunes [453, 452] assumed that the neutron stars are uniformly and slowly rotating, as well as only slightly tidally deformed by their rotational velocity or companion. These assumptions would not be valid for newly-born neutron stars, which are probably differentially rotating and doing so quickly. However, the gravitational waves emitted by neutron-star inspirals are expected to have binary components that are old and not rapidly spinning by the time they enter the detector sensitivity band [74]. Some short-period, millisecond pulsars may spin at a non-negligible rate, for which the normalized moment of inertia, quadrupole moment and Love number would not be independent of the rotational angular velocity. However, if then the above tests should still be possible, since binary pulsar observations would also automatically determine the rotational angular velocity, for which a unique I-Love-Q relation should exist in GR.

5.4 Tests of the no-hair theorems

Another important class of generic tests of GR are those that concern the no-hair theorems. Since much work has been done on this area, we have decided to separate this topic from the main generic tests section (5.3). In what follows, we describe what these theorems are and the possible tests one could carry out with gravitational-wave observations emitted by black-hole-binary systems.

5.4.1 The no-hair theorems

The no-hair theorems state that the only stationary, vacuum solution to the Einstein equations that is non-singular outside the event horizon is completely characterized by three quantities: its mass M, its spin S and its charge Q. This conclusion is arrived at by combining several different theorems. First, Hawking [223, 222] proved that a stationary black hole must have an event horizon with a spherical topology and that it must be either static or axially symmetric. Israel [243, 244] then proved that the exterior gravitational field of such static black holes is uniquely determined by M and Q and it must be given by the Schwarzschild or the Reissner-Nordström metrics. Carter [98] constructed a similar proof for uncharged, stationary, axially-symmetric black holes, where this time black holes fall into disjoint families, not deformable into each other and with an exterior gravitational field uniquely determined by M and S. Robinson [363] and Mazur [306] later proved that such black holes must be described by either the Kerr or the Kerr-Newman metric. See also [318, 352] for more details.

The no-hair theorems apply under a restrictive set of conditions. First, the theorems only apply in stationary situations. Black-hole horizons can be tidally deformed in dynamical situations, and if so, Hawking’s theorems [223, 222] about spherical horizon topologies do not apply. This then implies that all other theorems described above also do not apply, and thus, dynamical black holes will generically have hair. Second, the theorems only apply in vacuum. Consider, for example, an axially-symmetric black hole in the presence of a non-symmetrical matter distribution outside the event horizon. One might naively think that this would tidally distort the event horizon, leading to a rotating, stationary black hole that is not axisymmetric. However, Hawking and Hartle [226] showed that in such a case the matter distribution torques the black hole forcing it to spin down, thus leading to a non-stationary scenario. If the black hole is non-stationary, then again the no-hair theorems do not apply by the arguments described at the beginning of this paragraph, and thus non-isolated black holes can have hair. Third, the theorems only apply within GR, i.e., through the use of the Einstein equations. Therefore, it is plausible that black holes in modified gravity theories or in GR with singularities outside any event horizons (naked singularities) will have hair.

The no-hair theorems imply that the exterior gravitational field of isolated, stationary, uncharged and vacuum black holes (in GR and provided the spacetime is regular outside all event horizons) can be written as an infinite sum of mass and current multipole moments, where only two of them are independent: the mass monopole moment M and the current dipole moment S. One can extend these relations to include charge, but astrophysical black holes are expected to be essentially neutral due to charge accretion. If the no-hair theorems hold, all other multipole moments can be determined from [195, 194, 213]

$${M_\ell} + {\rm{i}}{S_\ell} = M{({\rm{i}}a)^\ell},$$

(219)

where M _ℓ and S _ℓ are the ℓth mass and current multipole moments. Even if the black-hole progenitor was not stationary or axisymmetric, the no-hair theorems guarantee that any excess multipole moments will be shed-off during gravitational collapse [356, 357]. Eventually, after the black hole has settled down and reached an equilibrium configuration, it will be described purely in terms of M₀ = M and S₁ = S = Ma², where a is the Kerr spin parameter.

An astrophysical observation of a hairy black hole would not imply that the no-hair theorems are wrong, but rather that one of the assumptions made in deriving these theorems is not appropriate to describe nature. As described above, the three main assumptions are stationarity, vacuum and that GR and the regularity condition hold. Astrophysical black holes will generically be hairy due to a violation of the first two assumptions, since they will neither be perfectly stationary, nor exist in a perfect vacuum. Astrophysical black holes will always suffer small perturbations by other stars, electromagnetic fields, other forms of matter, like dust, plasma or dark matter, etc, which will induce non-zero deviations from Eq. (219) and thus evade the no-hair theorems. However, in all cases of interest such perturbations are expected to be too small to be observable, which is why one argues that even astrophysical black holes should obey the no-hair theorems if GR holds. Put another way, an observation of the violation of the no-hair theorems would be more likely to indicate a failure of GR in the strong-field, than an unreasonably large amount of astrophysical hair.

Tests of the no-hair theorems come in two flavors: through electromagnetic observations [250, 251, 253, 254] and through gravitational wave observations [370, 371, 112, 196, 44, 50, 289, 390, 471, 422, 421, 184, 423, 364]. The former rely on radiation emitted by accelerating particles in an accretion disk around black holes. However, such tests are not clean as they require the modeling of complicated astrophysics, with matter and electromagnetic fields. Gravitational wave tests are clean in that respect, but unlike electromagnetic tests, they cannot be carried out yet due to lack of data. Other electromagnetic tests of the no-hair theorems exist, for example through the observation of close stellar orbits around Sgr A* [312, 313, 373] and pulsar-black-hole binaries [431], but these cannot yet probe the near-horizon, strong-field regime, since electromagnetic observations cannot yet resolve horizon scales. See [359] for reviews on this topic.

5.4.2 Extreme mass-ratio tests of the no-hair theorem

Gravitational wave tests of the no-hair theorems require the detection of either extreme mass-ratio inspirals or the ringdown of comparable-mass black-hole mergers with future space-borne gravitational-wave detectors [25, 24]. Extreme mass-ratio inspirals consist of a stellar-mass compact object spiraling into a supermassive black hole in a generic orbit within astronomical units from the event horizon of the supermassive object [23]. These events outlive the observation time of future detectors, emitting millions of gravitational wave cycles, with the stellar-mass compact object essentially acting as a tracer of the supermassive black hole spacetime [397]. Ringdown gravitational waves are always emitted after black holes merge and the remnant settles down into its final configuration. During the ringdown, the highly-distorted remnant radiates all excess degrees of freedom and this radiation carries a signature of whether the no-hair theorems hold in its quasi-normal mode spectrum (see, e.g., [68] for a recent review).

Both electromagnetic and gravitational wave tests need a metric with which to model accretion disks, quasi-periodic oscillations, or extreme mass-ratio inspirals. One can classify these metrics as direct or generic, paralleling the discussion in Section 5.2. Direct metrics are exact solutions to a specific set of field equations, with which one can derive observables. Examples of such metrics are the Manko-Novikov metric [302] and the slowly-spinning black-hole metric in dynamical Chern-Simons gravity [466]. When computing observables with these metrics, one usually assumes that all radiative and dynamical process (e.g., the radiation-reaction force) are as predicted in GR. Generic metrics are those that parametrically modify the Kerr spacetime, such that for certain parameter choices one recovers identically the Kerr metric, while for others, one has a deformation of Kerr. Generic metrics can be further classified into two subclasses, Ricci-flat versus non-Ricci-flat, depending on whether they satisfy R _μν = 0.

Let us first consider direct metric tests of the no-hair theorem. The most studied direct metric is the Manko-Novikov one, which although an exact, stationary and axisymmetric solution to the vacuum Einstein equations, does not represent a black hole, as the event horizon is broken along the equator by a ring singularity [302]. Just like the Kerr metric, the Manko-Novikov metric possesses an ergoregion, but unlike the former, it also possesses regions of closed time-like curves that overlap the ergoregion. Nonetheless, an appealing property of this metric is that it deviates continuously from the Kerr metric through certain parameters that characterize the higher multiple moments of the solution.

The first geodesic study of Manko-Novikov spacetimes was carried out by Gair et al. [182]. They found that there are two ring-like regions of bound orbits: an outer one where orbits look regular and integrable, as there exist four isolating integrals of the motion; and an inner one where orbits are chaotic and thus ergodic. Gair et al. [182] suggested that orbits that transition from the integrable to the chaotic region would leave a clear observable signature in the frequency spectrum of the emitted gravitational waves. However, they also noted that chaotic regions exist only very close to the central body and are probably not astrophysically accessible. The study of Gair et al. [182] was recently confirmed and followed up by Contopoulos et al. [116]. They studied a wide range of geodesics and found that, in addition to an inner chaotic region and an outer regular region, there are also certain Birkhoff islands of stability. When an extreme mass-ratio inspiral traverses such a region, the ratio of resonant fundamental frequencies would remain constant in time, instead of increasing monotonically. Such a feature would impact the gravitational waves emitted by such a system, and it would signal that the orbit equations are non-integrable and the central object is not a Kerr black hole.

The study of chaotic motion in geodesics of non-Kerr spacetimes is by no means new. Chaos has also been found in geodesics of Zipoy-Voorhees-Weyl and Curzon spacetimes with multiple singularities [391, 392] and in general for Zipoy-Voorhees spacetimes in [296], of perturbed Schwarzschild spacetimes [287], of Schwarzschild spacetimes with a dipolar halo [286, 288, 209] of Erez-Rosen spacetimes [210], and of deformed generalizations of the Tomimatsy-Sato spacetime [154]. One might worry that such chaotic orbits will depend on the particular spacetime considered, but recently Apostolatos et al. [31] and Lukes-Gerakopoulos et al. [297] have argued that the Birkhoff islands of stability are a general feature. Although the Kolmogorov, Arnold, and Moser theorem [270, 35, 321] states that phase orbit tori of an integrable system are only deformed if the Hamiltonian is perturbed, the Poincare-Birkhoff theorem [292] states that resonant tori of integrable systems actually disintegrate, leaving behind a chain of Birkhoff islands. These islands are only characterized by the ratio of winding frequencies that equals a rational number, and thus, they constitute a distinct and generic feature of non-integrable systems [31, 297]. Given an extreme mass-ratio gravitational-wave detection, one can monitor the ratio of fundamental frequencies and search for plateaus in their evolution, which would signal non-integrability. Of course, whether detectors can resolve such plateaus depends on the initial conditions of the orbits and the physical system under consideration (these determine the thickness of the islands), as well as the mass ratio (this determines the radiation-reaction timescale) and the distance and mass of the central black hole (this determines the SNR).

Another example of a direct metric test of the no-hair theorem is through the use of the slowly-rotating dynamical Chern-Simons black hole metric [466]. Unlike the Manko-Novikov metric, the dynamical Chern-Simons one does represent a black hole, i.e., it possesses an event horizon, but it evades the no-hair theorems because it is not a solution to the Einstein equations. Sopuerta and Yunes [390] carried out the first extreme mass-ratio inspiral analysis when the background supermassive black hole object is taken to be such a Chern-Simons black hole. They used a semi-relativistic model [368] to evolve extreme mass-ratio inspirals and found that the leading-order modification comes from a modification to the geodesic trajectories, induced by the non-Kerr modifications of the background. Because the latter correspond to a strong-field modification to GR, modifications in the trajectories are most prominent for zoom-whirl orbits, as the small compact object zooms around the supermassive black hole in a region of unstable orbits, close to the event horizon. These modifications were then found to propagate into the gravitational waves emitted, leading to a dephasing that could be observed or ruled out with future gravitational-wave observations to roughly the horizon scale of the supermassive black hole, as has been recently confirmed by Canizares et al. [93]. However, these studies may be underestimates, given that they treat the black hole background in dynamical Chern-Simons gravity only to first-order in spin.

A final example of a direct metric test of the no-hair theorems is to consider black holes that are not in vacuum. Barausse et al. [52] studied extreme-mass-ratio inspirals in a Kerr-black-hole background that is perturbed by a self-gravitating, homogeneous torus that is compact, massive and close to the Kerr black hole. They found that the presence of this torus impacts the gravitational waves emitted during such inspirals, but only weakly, making it difficult to distinguish the presence of matter. Yunes et al. [462] and Kocsis et al. [267] carried out a similar study, where this time they considered a small compact object inspiraling completely within a geometrically thin, radiation-pressure dominated accretion disk. They found that disk-induced migration can modify the radiation-reaction force sufficiently so as to leave observable signatures in the waveform, provided the accretion disk is sufficiently dense in the radiation-dominated regime and a gap opens up. However, these tests of the no-hair theorem will be rather difficult as most extreme-mass-ratio inspirals are not expected to be in an accretion disk.

Let us now consider generic metric tests of the no-hair theorem. Generic Ricci-flat deformed metrics will lead to Laplace-type equations for the deformation functions in the far-field since they must satisfy R _μν = 0 to linear order in the perturbations. The solution to such an equation can be expanded in a sum of mass and current multipole moments, when expressed in asymptotically Cartesian and mass-centered coordinates [407]. These multipoles can be expressed via [112, 422, 421]

$${M_\ell} + {\rm{i}}{S_\ell} = M{({\rm{i}}a)^\ell} + \delta {M_\ell} + {\rm{i}}\delta {M_\ell},$$

(220)

where δM _ℓ and δS _ℓ are mass and current multipole deformations. Ryan [370, 371] showed that the measurement of three or more multipole moments would allow for a test of the no-hair theorem. Generic non-Ricci flat metrics, on the other hand, will not necessarily lead to Laplace-type equations for the deformation functions in the far field, and thus, the far-field solution and Eq. (220) will depend on a sum of ℓ and m multipole moments.

The first attempt to construct a generic, Ricci-flat metric was by Collins and Hughes [112]: the bumpy black-hole metric. In this approach, the metric is assumed to be of the form

$${g_{\mu \nu}} = g_{\mu \nu}^{({\rm{Kerr}})} + \epsilon \delta {g_{\mu \nu}},$$

(221)

where ϵ < 1 is a bookkeeping parameter that enforces that δg _μν is a perturbation of the Kerr background. This metric is then required to satisfy the Einstein equations linearized in ϵ, which then leads to differential equations for the metric deformation. Collins and Hughes [112] assumed a non-spinning, stationary spacetime, and thus δg _μν only possessed two degrees of freedom, both of which were functions of radius only: ψ₁(r), which must be a harmonic function and which changes the Newtonian part of the gravitational field at spatial infinity; and γ₁(r) which is completely determined through the linearized Einstein equations once ψ₁ is specified. One then has the freedom to choose how to prescribe ψ₁ and Collins and Hughes investigate [112] two choices that correspond physically to point-like and ring-like naked singularities, thus violating cosmic censorship [347]. Vigeland and Hughes [422] and Vigeland [421] then extend this analysis to stationary, axisymmetric spacetimes via the Newman-Janis method [327, 151], showing how such metric deformations modify Eq. (220), and computing how these bumps imprint themselves onto the orbital frequencies and thus the gravitational waves emitted during an extreme-mass-ratio inspiral.

That the bumps represent unphysical matter should not be a surprise, since by the no-hair theorems, if the bumps are to satisfy the vacuum Einstein equations they must either break stationarity or violate the regularity condition. Naked singularities are an example of the latter. A Lorentz-violating massive field coupled to the Einstein tensor is another example [155]. Gravitational wave tests with bumpy black holes must then be understood as null tests: one assumes the default hypothesis that GR is correct and then sets out to test whether the data rejects or fails to reject this hypothesis (a null hypothesis can never be proven). Unfortunately, however, bumpy black hole metrics cannot parameterize spacetimes in modified gravity theories that lead to corrections in the field equations that are not proportional to the Ricci tensor, such as for example in dynamical Chern-Simons or in Einstein-Dilaton-Gauss-Bonnet modified gravity.

Other bumpy black hole metrics have also been recently proposed. Glampedakis and Babak [196] proposed a different type of stationary and axisymmetric bumpy black hole through the Hartle-Thorne metric [218], with modifications to the quadrupole moment. They then constructed a “kludge” extreme mass-ratio inspiral waveform and estimated how well the quadrupole deformation could be measured [44]. However, this metric is valid only when the supermassive black hole is slowly-rotating, as it derives from the Hartle-Thorne ansatz. Recently, Johansen and Psaltis [252] proposed yet another metric to represent bumpy stationary and spherically-symmetric spacetimes. This metric introduces one new degree of freedom, which is a function of radius only and assumed to be a series in M/r. Johansen and Psaltis then rotated this metric via the Newman-Janis method [327, 151] to obtain a new bumpy metric for axially-symmetric spacetimes. However, such a metric possesses a naked ring singularity on the equator, and naked singularities on the poles. As before, none of these bumpy metrics can be mapped to known modified gravity black hole solutions, in the Glampedakis and Babak case [196] because the Einstein equations are assumed to hold to leading order in the spin, while in the Johansen and Psaltis case [252] because a single degree of freedom is not sufficient to model the three degrees of freedom contained in stationary and axisymmetric spacetimes [401, 423].

The only generic non-Ricci-flat bumpy black-hole metric so far is that of Vigeland, Yunes and Stein [423]. They allowed generic deformations in the metric tensor, only requiring that the new metric perturbatively retained the Killing symmetries of the Kerr spacetime: the existence of two Killing vectors associated with stationarity and axisymmetry, as well as the perturbative existence of a Killing tensor (and thus a Carter-like constant), at least to leading order in the metric deformation. Such requirements imply that the geodesic equations in this new background are fully integrable, at least perturbatively in the metric deformation, which then allows one to solve for the orbital motion of extreme-mass-ratio inspirals by adapting previously existing tools. Brink [83, 84, 85, 86, 87] studied the existence of such a second-order Killing tensor in generic, vacuum, stationary and axisymmetric spacetimes in Einstein’s theory and found that these are difficult to construct exactly. By relaxing this exact requirement, Vigeland, Yunes and Stein [423] found that the existence of a perturbative Killing tensor poses simple differential conditions on the metric perturbation that can be analytically solved. Moreover, they also showed how this new bumpy metric can reproduce all known modified gravity black hole solutions in the appropriate limits, provided these have an at least approximate Killing tensor; thus, these metrics are still vacuum solutions even though R ≠ 1, since they satisfy a set of modified field equations. Although unclear at this junction, it seems that the imposition that the spacetime retains the Kerr Killing symmetries leads to a bumpy metric that is well-behaved everywhere outside the event horizon (no singularities, no closed-time-like curves, no loss of Lorentz signature). Recently, Gair and Yunes [184] studied how the geodesic equations are modified for a test-particle in a generic orbit in such a spacetime and showed that the bumps are indeed encoded in the orbital motion, and thus, in the gravitational waves emitted during an extreme-mass-ratio inspiral.

One might be concerned that such no-hair tests of GR cannot constrain modified gravity theories, because Kerr black holes can also be solutions in the latter [360]. This is indeed true provided the modified field equations depend only on the Ricci tensor or scalar. In Einstein-Dilaton-Gauss-Bonnet or dynamical Chern-Simons gravity, the modified field equations depend on the Riemann tensor, and thus, Ricci-flat metric need not solve these modified set [473]. Moreover, just because the metric background is identically Kerr does not imply that inspiral gravitational waves will be identical to those predicted in GR. All studies carried out to date, be it direct metric tests or generic metric tests, assume that the only quantity that is modified is the metric tensor, or equivalently, the Hamiltonian or binding energy. Inspiral motion, of course, does not depend just on this quantity, but also on the radiation-reaction force that pushes the small object from geodesic to geodesic. Moreover, the gravitational waves generated during such an inspiral depend on the field equations of the theory considered. Therefore, all metric tests discussed above should be considered as partial tests. In general, strong-field modified gravity theories will modify the Hamiltonian, the radiation-reaction force and the wave generation.

5.4.3 Ringdown tests of the no-hair theorem

Let us now consider tests of the no-hair theorems with gravitational waves emitted by comparable-mass binaries during the ringdown phase. Gravitational waves emitted during ringdown can be described by a superposition of exponentially-damped sinusoids [69]:

$${h_ +}(t) + i\;{h_ \times}(t) = {M \over r}\sum\limits_{\ell mn} {\left\{{{{\mathcal A}_{\ell mn}}{e^{i({\omega _{\ell mn}}t + {\phi _{\ell mn}})}}{e^{- t/{\tau _{\ell mn}}}}{S_{\ell mn}} + {\mathcal A}_{\ell mn}^\prime{e^{i(- {\omega _{\ell mn}}t + \phi _{\ell mn}^\prime)}}{e^{- t/{\tau _{\ell mn}}}}S_{\ell mn}^{\ast}} \right\},}$$

(222)

where r is the distance from the source to the detector, the asterisk stands for complex conjugation, the real mode amplitudes ${\mathcal{A}_{\ell ,m,n}}$ and $\mathcal{A}_{\ell ,m,n}^{\prime}$ and the real phases ϕ _nℓm and $\phi _{\ell ,m,n}^{\prime}$ depend on the initial conditions, S _ℓmn are spheroidal functions evaluated at the complex quasinormal ringdown frequencies ω _nℓm = 2πf _nℓm + i/τ _nℓm , and the real physical frequency f _nℓm and the real damping times τ _nℓm are both functions of the mass M and the Kerr spin parameter a only, provided the no-hair theorems hold. These frequencies and damping times can be computed numerically or semi-analytically, given a particular black-hole metric (see [68] for a recent review). The Fourier transform of a given (ℓ, m, n) mode is [69]

$$\tilde h_ + ^{(\ell ,m,n)}(\omega) = {M \over r}{\mathcal A}_{\ell m n}^ + \left[ {{e^{i\phi _{\ell mn}^ +}}{S_{\ell mn}}{b_ +}(\omega) + {e^{- i\phi _{\ell mn}^ +}}S_{\ell mn}^{\ast}{b_ -}(\omega)} \right],$$

(223)

$$\tilde h_ \times ^{(\ell ,m,n)}(\omega) = {M \over r}{\mathcal A}_{\ell mn}^\times \left[ {{e^{i\phi _{\ell mn}^\times}}{S_{\ell mn}}{b_ +}(\omega) + {e^{- i\phi _{\ell mn}^\times}}S_{\ell mn}^{\ast}{b_ -}(\omega)} \right],$$

(224)

where we have defined $\mathcal{A}_{\ell mn}^{+ , \times}{e^{i\phi _{\ell mn}^{+ , \times}}} \equiv {\mathcal{A}_{\ell mn}}{e^{i{\phi _{\ell mn}}}} \pm \mathcal{A}{\prime}{e^{- i\phi _{\ell mn}^{\prime}}}$ as well as the Lorentzian functions

$${b_ \pm}(\omega) = {{{\tau _{\ell mn}}} \over {1 + \tau _{\ell mn}^2{{(\omega \pm {\omega _{\ell mn}})}^2}}}.$$

(225)

Ringdown gravitational waves will all be of the form of Eq. (222) provided that the characteristic nature of the differential equation that controls the evolution of ringdown modes is not modified, i.e., provided that one only modifies the potential in the Teukolsky equation or other subdominant terms, which in turn depend on the modified field equations.

Tests of the no-hair theorems through the observation of black-hole ringdown date back to Detweiler [146], and it was recently worked out in detail by Dreyer et al. [152]. Let us first imagine that a single complex mode is detected ${\omega _{{\ell _1}{m_1}{n_1}}}$ and one measures separately its real and imaginary parts. Of course, from such a measurement, one cannot extract the measured harmonic triplet (ℓ₁,m₁,n₁), but instead one only measures the complex frequency ${\omega _{{\ell _1}{m_1}{n_1}}}$. This information is not sufficient to extract the mass and spin angular momentum of the black hole because different quintuplets (M, a, ℓ, m, n) can lead to the same complex frequency ${\omega _{{\ell _1}{m_1}{n_1}}}$. The best way to think of this is graphically: a given observation of $\omega _{{\ell _1}{m_1}{n_1}}^{(1)}$ traces a line in the complex ${\Omega _{{\ell _1}{m_1}{n_1}}} = M\omega _{{\ell _1}{m_1}{n_1}}^{(1)}$ plane; a given (ℓ, m, n) triplet defines a complex frequency ω _ℓmn that also traces a curve in the complex Ω _ℓmn plane; each intersection of the measured line ${\Omega _{{\ell _1}{m_1}{n_1}}}$ with Ω _ℓmn defines a possible doublet (M, a); since different (ℓ, m, n) triplets lead to different ω _ℓmn curves and thus different intersections, one ends up with a set of doublets S₁, out of which only one represents the correct black-hole parameters. We thus conclude that a single mode observation of ringdown gravitational waves is not sufficient to test the no-hair theorem [152, 69]

Let us then imagine that one has detected two complex modes, ${\omega _{{\ell _1}{m_1}{n_1}}}$ and ${\omega _{{\ell _2}{m_2}{n_2}}}$. Each detection leads to a separate line ${\Omega _{{\ell _1}{m_1}{n_1}}}$ and ${\Omega _{{\ell _2}{m_2}{n_2}}}$ the complex plane. As before, each (n, ℓ, m) triplet leads to separate curves Ω _ℓmn which will intersect with both ${\Omega _{{\ell _1}{m_1}{n_1}}}$ and ${\Omega _{{\ell _2}{m_2}{n_2}}}$ in the complex plane. Each intersection between Ω _ℓmn and ${\Omega _{{\ell m n}}}$ leads to a set of doublets S₁, while each intersection between Ωℓmn and ${\Omega _{{\ell _1}{m_1}{n_1}}}$ leads to another set of doublets S₂. However, if the no-hair theorems hold sets S₁ and S₂ must have at least one element in common. Therefore, a two-mode detection allows for tests of the no-hair theorem [152, 69]. However, when dealing with a quasi-circular black-hole-binary inspiral within GR one knows that the dominant mode is ℓ = m = 2. In such a case, the observation of this complex mode by itself allows one to extract the mass and spin angular momentum of the black hole. Then, the detection of the real frequency in an additional mode can be used to test the no-hair theorem [69, 65].

Although the logic behind these tests is clear, one must study them carefully to determine whether all systematic and statistical errors are sufficiently under control so that they are feasible. Berti et al. [69, 65] investigated such tests carefully through a frequentist approach. First, they found that a matched-filtering type analysis with two-mode ringdown templates would increase the volume of the template manifold by roughly three orders of magnitude. A better strategy then is perhaps to carry out a Bayesian analysis, like that of Gossan et al. [256, 201]; through such a study one can determine whether a given detection is consistent with a two-mode or a one-mode hypothesis. Berti et al. [69, 65] also calculated that a SNR of $\mathcal{O}({10^2})$ would be sufficient to detect the presence of two modes in the ringdown signal and to resolve their frequencies, so that no-hair tests would be possible. Strong signals are necessary because one must be able to distinguish at least two modes in the signal. Unfortunately, however, whether the ringdown leads to such strong SNRs and whether the sub-dominant ringdown modes are of a sufficiently large amplitude depends on a plethora of conditions: the location of the source in the sky, the mass of the final black hole, which depends on the rest mass fraction that is converted into ringdown gravitational waves (the ringdown efficiency), the mass ratio of the progenitor, the magnitude and direction of the spin angular momentum of the final remnant and probably also of the progenitor and the initial conditions that lead to ringdown. Thus, although such tests are possible, one would have to be quite fortunate to detect a signal with the right properties so that a two-mode extraction and a test of the no-hair theorems is feasible.

5.4.4 The hairy search for exotica

Another way to test GR is to modify the matter sector of the theory through the introduction of matter corrections to the Einstein-Hilbert action that violate the assumptions made in the no-hair theorems. More precisely, one can study whether gravitational waves emitted by binaries composed of strange stars, like quark stars, or horizonless objects, such as boson stars or gravastars, are different from waves emitted by more traditional neutron-star or black-hole binaries. In what follows, we will describe such hairy tests of the existence of compact exotica.

Boson stars are a classic example of a compact object that is essentially indistinguishable from a black hole in the weak field, but which differs drastically from one in the strong field due to its lack of an event horizon. A boson star is a coherent scalar-field configuration supported against gravitational collapse by its self-interaction. One can construct several Lagrangian densities that would allow for the existence of such an object, including mini-boson stars [178, 179], axially-symmetric solitons [372], and nonsolitonic stars supported by a non-canonical scalar potential energy [113]. Boson stars are well-motivated from fundamental theory, since they are the gravitationally-coupled limit of q-balls [108, 276], a coherent scalar condensate that can be described classically as a non-topological soliton and that arises unavoidably in viable supersymmetric extensions of the standard model [275]. In all studies carried out to date, boson stars have been studied within GR, but they are also allowed in scalar-tensor theories [46].

At this junction, one should point out that the choice of a boson star is by no means special; the key point here is to select a straw-man to determine whether gravitational waves emitted during the coalescence of compact binaries are sensitive to the presence of an event horizon or the evasion of the no-hair theorems induced by a non-vacuum spacetime. Of course, depending on the specific model chosen, it is possible that the exotic object will be unstable to evolution or even to its own rotation. For example, in the case of an extreme mass-ratio inspiral, one could imagine that as the small compact object enters the boson star’s surface, it will accrete the scalar field, forcing the boson star to collapse into a black hole. Alternatively, one can imagine that as two supermassive boson stars merge, the remnant might collapse into a black hole, emitting the usual GR quasinormal modes. What is worse, even when such objects are in isolation, they are unstable under small perturbations if their angular momentum is large, possibly leading to gravitational collapse into a black hole or possibly a scalar explosion [95, 96]. Since most astrophysical black hole candidates are believed to have high spins, such instabilities somewhat limit the interest of horizonless objects. Even then, however, the existence of slowly spinning or non spinning horizonless compact objects cannot be currently ruled out by observation.

Boson stars evade the no-hair theorems within GR because they are not vacuum spacetimes, and thus, their metric and quasinormal mode spectrum cannot be described by just their mass and spin angular momentum; one also requires other quantities intrinsic to the scalar-field energy momentum tensor, scalar hair. Therefore, as before, two types of gravitational wave tests for scalar hair have been proposed: extreme-mass-ratio inspiral tests and ringdown tests. As for the former, several studies have been carried out that considered a supermassive boson star background. Kesden et al. [263] showed that stable circular orbits exist both outside and inside of the surface of the boson star, provided the small compact object interacts with the background only gravitationally. This is because the effective potential for geodesic motion in such a boson-star background lacks the Schwarzschild-like singular behavior at small radius, instead turning over and allowing for a new minimum. Gravitational waves emitted in such a system would then stably continue beyond what one would expect if the background had been a supermassive black hole; in the latter case the small compact object would simply disappear into the horizon. Kesden et al. [263] found that orbits inside the boson star exhibit strong precession, exciting high frequency harmonics in the waveform, and thus allowing one to easily distinguish between such boson stars from black-hole backgrounds.

Just as the inspiral phase is modified by the presence of a boson star, the merger phase is also greatly altered, but this must be treated fully numerically. A few studies have found that the merger of boson stars leads to a spinning bar configuration that either fragments or collapses into a Kerr black hole [339, 338]. Of course, the gravitational waves emitted during such a merger will be drastically different from those produced when black holes merge. Unfortunately, the complexity of such simulations makes predictions difficult for any one given example, and the generalization to other more complicated scenarios, such as theories with modified field equations, is currently not feasible.

Recently, Pani et al. [340, 341] revisited this problem, but instead of considering a supermassive boson star, they considered a gravastar. This object consists of a Schwarzschild exterior and a de Sitter interior, separated by an infinitely thin shell with finite tension [307, 100]. Pani et al. [341] calculated the gravitational waves emitted by a stellar-mass compact object in a quasi-circular orbit around such a gravastar background. In addition to considering a different background, Pani et al. used a radiative-adiabatic waveform generation model to describe the gravitational waves [351, 238, 239, 458, 456, 459], instead of the kludge scheme used by Kesden et al. [49, 44, 456]. Pani el al. [341] concluded that the waves emitted during such inspirals are sufficiently different that they could be used to discern between a Kerr black hole and a gravastar.

On the ringdown side of no-hair tests, several studies have been carried out. Berti and Cardoso [66] calculated the quasi-normal mode spectrum of boson stars. Chirenti and Rezzolla [105] studied the non-radial, axial perturbations of gravastars, and Pani et al. [340] the non-radial, axial and polar oscillations of gravastars. Medved et al. [309, 310] considered the quasinormal ringdown spectrum of skyrmion black holes [386]. In all cases, it was found that the quasi-normal mode spectrum of such objects could be used to discern between them and Kerr black holes. Of course, such tests still require the detection of ringdown gravitational waves with the right properties, such that more than one mode can be discerned and extracted from the signal (see Section 5.4.3).

6 Musings About the Future

Gravitational waves hold the key to testing Einstein’s theory of general relativity (GR) to new exciting levels in the previously unexplored strong-field regime. Depending on the type of wave that is detected, e.g., compact binary inspirals, mergers, ringdowns, continuous sources, supernovae, etc, different tests will be possible. Irrespective of the type of wave detected, two research trends seem currently to be arising: direct tests and generic tests. These trends aim at answering different questions. With direct tests, one wishes to determine whether a certain modified theory is consistent with the data. Generic tests, on the other hand, ask whether the data is statistically consistent with our canonical beliefs. Or put another way: are there any statistically-significant deviations present in the data from what we expected to observe? This approach is currently used in cosmological observations, for example by the WMAP team, and it is particularly well-suited when one tries to remain agnostic as to which is the correct theory of nature. Given that we currently have no data in the strong-field, it might be too restrictive to assume GR is correct prior to verifying that this is the case.

What one would like to believe is that gravitational waves will be detected by the end of this decade, either through ground-based detectors or through pulsar timing arrays. Given this, there is a concrete effort to develop the proper formalism and implementation pipelines to test Einstein’s theory once data becomes available. Currently, the research groups separate into two distinct classes: theory and implementation. The theory part of the research load is being carried out at a variety of institutions without a given focal point. The implementation part is being done mostly within the LIGO Scientific Collaboration and the pulsar timing consortia. Cross-communication between the theory and implementation groups has recently flourished and one expects more interdisciplinary work in the future.

So many accomplishments have been made in the past 50 years that it is almost impossible to list them here. From the implementation side, perhaps one of the most important is the actual construction and operation of the initial LIGO instruments at design sensitivity in all of their frequency domain. This is a tremendously important engineering and physics challenge. Similarly, the construction of impressive pulsar timing arrays, and the timing of these pulses to nano-second precision is an instrumental and data analysis feat to be admired. Without these observatories no waves would be detectable in the future, and of course, no tests of Einstein’s theory would be feasible. On the theory side, perhaps the most important accomplishment has been the understanding of the inspiral phase to really-high post-Newtonian order and the merger phase with numerical simulations. The latter, in particular, had been an unsolved problem for over 50 years, until very recently. It is these accomplishments that then allow us to postulate modified inspiral template families, since we understand what the GR expectation is. This is particularly true if one is considering small deformations away from Einstein’s theory, as it would be impossible to perturb about an unknown solution.

The main questions that are currently at the forefront are the following. On the theory side of things, one would wish to understand the inspiral to high post-Newtonian order in certain strong-field modifications to GR, like dynamical Chern-Simons gravity or Einstein-Dilaton-Gauss-Bonnet theory. One would also like to investigate theories with preferred frames, such as Einstein-Aether theory or Hořava-Lifshitz gravity, which will lead to Lorentz violating observables. Understanding these theories to high post-Newtonian order is particularly important for those that predict dipolar gravitational emission, such as Einstein-Dilaton-Gauss-Bonnet theory. Such corrections dominate over Einstein’s quadrupole emission at sufficiently low velocities.

Of course, a full inspiral-merger-ringdown template is not complete unless we also understand the merger. This would require full numerical simulations, which are very taxing even within GR. Once one modifies the Einstein field equations, the characteristic structure of the evolution equations will also likely change, and it is unclear whether the standard evolution methods will continue to work. Moreover, when dealing with the merger phase, one is usually forced to treat the modified theory as exact, instead of as an effective theory. Without the latter, it is likely that certain modified theories will not have a well-posed initial value problem, which would force any numerical evolution to fail. Of course, one could order-reduce these equations and then use these to evolve black-hole spacetimes. Much work still remains to be done to understand whether this is feasible.

On the implementation side of things, there is also much work that remains to be done. Currently, efforts are only beginning on the implementation of Bayesian frameworks for hypothesis testing. This seems today like one of the most promising approaches to testing Einstein’s theory with gravitational waves. Current studies concentrate mostly on single-detectors, but by the beginning of the next decade we expect four or five detectors to be online, and thus, one would like to see these implementations extended. The use of multiple detectors also opens the door to the extraction of new information, such as multiple polarization modes, a precise location of the source in the sky, etc. Moreover, the evidence for a given model increases dramatically if the event is observed in several detectors. One therefore expects that the strongest tests of GR will come from leveraging the data from all detectors in a multiply-coincident event.

Ultimately, research is moving toward the construction of robust techniques to test Einstein’s theory. A general push is currently observed toward the testing of general principles that serve as foundations of GR. This allows one to answer general questions, such as: Does the graviton have a mass? Are compact objects represented by the Kerr metric and the no-hair theorems satisfied? Does the propagating metric perturbation possess only two transverse-traceless polarization modes? What is the rate of change of a binary’s binding energy? Do naked singularities exist in nature and are orbits chaotic? Is Lorentz-violation present in the propagation of gravitons? These are examples of questions that can be answered once gravitational waves are detected. The more questions of this type that are generated and the more robust the methods to answer them are, the more stringent the test of Einstein’s theories and the more information we will obtain about the gravitational interaction in a previously unexplored regime.