Testing Gravity with Gravitational Waves from Binary Black Hole Mergers: Contributions from Amplitude Corrections

TL;DR

This paper investigates the impact of amplitude corrections in gravitational wave signals from binary black hole mergers to improve tests of alternative gravity theories, revealing that amplitude effects can be as significant as phase effects.

Contribution

It introduces a method to incorporate amplitude corrections into gravitational wave tests of gravity, providing more comprehensive constraints on modified theories.

Findings

Amplitude corrections can be comparable to phase corrections for massive binaries.

Including amplitude corrections changes constraints by up to 4%.

First constraints on scalar field evolution using gravitational waves.

Abstract

The detection of gravitational waves has offered us the opportunity to explore the dynamical and strong-field regime of gravity. Because matched filtering is more sensitive to variations in the gravitational waveform phase than the amplitude, many tests of gravity with gravitational waves have been carried out using only the former. Such studies cannot probe the non-Einsteinian effects that may enter only in the amplitude. Besides, if not accommodated in the waveform template, a non-Einsteinian effect in the amplitude may induce systematic errors on other parameters such as the luminosity distance. In this paper, we derive constraints on a few modified theories of gravity (Einstein-dilaton-Gauss-Bonnet gravity, scalar-tensor theories, and varying-$G$ theories), incorporating both phase and amplitude corrections. We follow the model-independent approach of the parametrized post-Einsteinian formalism. We perform Fisher analyses with Monte-Carlo simulations using the LIGO/Virgo posterior samples. We find that the contributions from amplitude corrections can be comparable to the ones from the phase corrections in case of massive binaries like GW150914. Also, constraints derived by incorporating both phase and amplitude corrections differ from the ones with phase corrections only by 4% at most, which supports many of the previous studies that only considered corrections in the phase. We further derive reliable constraints on the time-evolution of a scalar field in a scalar-tensor theory for the first time with gravitational waves.

Figures

Tables

Table 1. Table 1: 90% credible constraints on representative parameters of various modified theories of gravity from GW150914 and GW151226. For each of the GW events, “phase” and “amplitude” correspond to the cases where we include non-GR corrections only to the GW phase and amplitude respectively, while “combined” is the case where we include both corrections in the waveform and reduce the two constraints to a single one according to Sec.  II.2 . α ¯ EdGB subscript ¯ 𝛼 EdGB \bar{\alpha}_{{\mbox{\tiny EdGB}}} is the EdGB coupling parameter which is related to the dimensionless coupling by ζ EdGB ≡ 16 ​ π ​ α ¯ EdGB 2 / m 4 subscript 𝜁 EdGB 16 𝜋 superscript subscript ¯ 𝛼 EdGB 2 superscript 𝑚 4 \zeta_{\mbox{\tiny EdGB}}\equiv 16\pi\bar{\alpha}_{\mbox{\tiny EdGB}}^{2}/m^{4} with m 𝑚 m being the total mass of the binary. m 1 ​ ϕ ˙ subscript 𝑚 1 ˙ italic-ϕ m_{1}\dot{\phi} corresponds to a dimensionless parameter in scalar-tensor theories where m 1 subscript 𝑚 1 m_{1} is the mass of the primary BH while ϕ italic-ϕ \phi is the scalar field. The bounds are derived by assuming subdominant non-GR corrections, which is realized whenever ζ EdGB ≪ 1 much-less-than subscript 𝜁 EdGB 1 \zeta_{\mbox{\tiny EdGB}}\ll 1 ( m 1 ​ ϕ ˙ ≪ 1 much-less-than subscript 𝑚 1 ˙ italic-ϕ 1 m_{1}\dot{\phi}\ll 1 ) in EdGB (scalar-tensor) gravity. Numbers inside brackets mean such criterion is violated and the constraints are unreliable. G 𝐺 G is the gravitational constant with the subscript 0 representing the time of coalescence. An overhead dot denotes a derivative with respect to time.

Theories	Repr. Parameter	Constraints
		GW150914			GW151226
		Phase	Amplitude	Combined	Phase	Amplitude	Combined
EdGB Yagi et al. (2012)	$\sqrt{\left|{\overline{\alpha }}_{\text{EdGB}}\right|}$ [km]	(50.5)	(76.3)	(51.5)	4.32	10.5	4.32
	${\zeta }_{\text{EdGB}}$	3.62	32.4	3.91	0.0207	0.709	0.0207
Scalar-Tensor Scharre and Will (2002); Berti et al. (2005)	$\left|\dot{\varphi }\right|$ [${10}^4/\sec$]	(3.64)	(7.30)	(3.77)	1.09	(5.60)	1.09
	$\left|m_1\dot{\varphi }\right|$	6.87	16.4	7.15	0.688	3.66	0.688
Varying-$G$ Tahura and Yagi (2018); Yunes et al. (2010b)	$\left|\dot{G_0}/G_0\right|$ [${10}^6/\mathrm{yr}$]	7.30	137	7.18	0.0224	0.382	0.0220

Equations

$$\tilde{h}(f)={\tilde{h}}_{\text{mbox}GR}(1+{\alpha }_{\text{mbox}PPE}{u}^a)e^{i\delta \Psi },$$

$$\delta \Psi ={\beta }_{\text{mbox}PPE}{u}^b.$$

$$r=r_{\text{mbox}GR}(1+{\gamma }_r{u}^{c_r}),$$

$$\dot{f}={\dot{f}}_{\text{mbox}GR}\left(1+{\gamma }_{\dot{f}}{u}^{c_{\dot{f}}}\right).$$

$$r_{\text{mbox}GR}={\left(\frac{m}{{\Omega }^2}\right)}^{1/3},{\dot{f}}_{\text{mbox}GR}=\frac{96}{5}{\pi }^{8/3}{\mathcal{M}}^{5/3}{f}^{11/3},$$

$$\tilde{\mathcal{A}}(f)={\tilde{\mathcal{A}}}_{\text{mbox}GR}\left(1+2{\gamma }_r{u}^{c_r}-\frac{1}{2}{\gamma }_{\dot{f}}{u}^{c_{\dot{f}}}\right),$$

$$\Psi ={\Psi }_{\text{mbox}GR}-\frac{15\text{γf˙}}{16(\text{cf˙}-8)(\text{cf˙}-5)}{u}^{c_{\dot{f}}-5},$$

$${\alpha }_{\text{mbox}PPE}=\frac{8}{15}(a-8)(a-5){\beta }_{\text{mbox}PPE},$$

$${\alpha }_{\text{mbox}PPE}=\frac{8}{15}\frac{(8-a)(5-a)(a^2-4a-6)}{a^2-2a-6}{\beta }_{\text{mbox}PPE}.$$

$$b=a-5.$$

$$s(t)=h(t)+n(t),$$

$$\left(A|B\right)=4\Re {\int }_0^{\infty }df\frac{{\tilde{A}}^{\ast }(f)\tilde{B}(f)}{S_n\left(f\right)}.$$

$$P\left(n=n_0(t)\right)\propto \text{exp}\left[-\left(n_0|n_0\right)\right],$$

$$\rho \equiv \sqrt{\left(h|h\right)}.$$

$$P\left({\theta }^a|s\right)\propto p^{(0)}\left({\theta }^a\right)\text{exp}\left[-\frac{1}{2}{\Gamma }_{ab}\Delta {\theta }^a\Delta {\theta }^b\right],$$

$${\Gamma }_{ab}=\left({\partial }_{a}h|{\partial }_{b}h\right),$$

$${\Sigma }^{ab}={\left({\tilde{\Gamma }}^{-1}\right)}^{ab}=⟨\Delta {\theta }^a\Delta {\theta }^b⟩,$$

$$p^{(0)}\left({\theta }^a\right)\text{exp}\left[-\frac{1}{2}{\Gamma }_{ab}\Delta {\theta }^a\Delta {\theta }^b\right]=\text{exp}\left[-\frac{1}{2}{\tilde{\Gamma }}_{ab}\Delta {\theta }^a\Delta {\theta }^b\right].$$

$${\theta }^a\equiv \left(\ln {\mathcal{M}}_z,\ln \eta ,\chi ,\ln D_L,\ln t_0,{\varphi }_0,\alpha ,\delta ,\psi ,\iota ,{\theta }_{\text{mbox}PPE}\right),$$

$$P\left(\xi \right)=\int P\left(\xi |{\sigma }_{\xi }\right)P\left({\sigma }_{\xi }\right)d{\sigma }_{\xi },$$

$${\beta }_{\text{mbox}MG}=\frac{{\pi }^2}{{\lambda }_g^2}\frac{\mathcal{M}}{1+z}D,b=3,$$

$$D=\frac{z}{H_0\sqrt{{\Omega }_M+{\Omega }_{\Lambda }}}\left[1-\frac{z}{4}\left(\frac{3{\Omega }_M}{{\Omega }_M+{\Omega }_{\Lambda }}\right)+\mathcal{O}(z^2)\right].$$

$${\beta }_{\text{mbox}EdGB}=-\frac{5}{7168}{\zeta }_{\text{mbox}EdGB}\frac{(m_1^2{\tilde{s}}_2^{\text{mbox}EdGB}-m_2^2{\tilde{s}}_1^{\text{mbox}EdGB}{)}^2}{m^4{\eta }^{18/5}},$$

$${\alpha }_{\text{mbox}EdGB}=-\frac{5}{192}{\zeta }_{\text{mbox}EdGB}\frac{(m_1^2{\tilde{s}}_2^{\text{mbox}EdGB}-m_2^2{\tilde{s}}_1^{\text{mbox}EdGB}{)}^2}{m^4{\eta }^{18/5}},$$

$${\beta }_{\text{mbox}ST}=-\frac{5}{7168}{\eta }^{2/5}({\alpha }_1-{\alpha }_2{)}^2$$

$${\alpha }_{\text{mbox}ST}=-\frac{5}{192}{\eta }^{2/5}({\alpha }_1-{\alpha }_2{)}^2$$

$${\alpha }_A=2m_A\dot{\varphi }[1+(1-{\chi }_A^2{)}^{1/2}],$$

$$m_A(t)\approx m_{A,0}+{\dot{m}}_{A,0}(t-t_0),$$

$$G(t)\approx G_0+{\dot{G}}_0(t-t_0).$$

$${\beta }_{\dot{G}}$$

$${\alpha }_{\dot{G}}$$

Full text

I Introduction

So far, general relativity (GR) is the most successful theory of gravitation. This century-old theory which exquisitely describes gravity as the curvature of spacetime has passed numerous tests with high precision Will (2014). Nonetheless, GR is not expected to be a complete description of gravity. The inconsistencies in galaxy rotation curves Bosma (1981a, b); Begeman et al. (1991); Rubin and Ford (1970); Rubin et al. (1980); Ostriker and Peebles (1973); Ostriker (1993) and the accelerated expansion of the universe Abbott (1988); Copeland et al. (2006); Perlmutter et al. (1999); Riess et al. (1998, 2004); Weinberg (1989); van Albada et al. (1985); Weinberg et al. (2013) are difficult to explain within the formulation of GR without introducing dark matter and dark energy. Moreover, a new theory is required to reconcile quantum mechanics with classical gravity Adler (2010); Ng (2003). Hence, one needs to continue testing GR through various experiments and observations. Gravitational wave (GW) observations are one of the most recent additions to this venture Abbott et al. (2016a); Yunes et al. (2016); Abbott et al. (2019, 2017a, 2018a) which have enabled us to probe the formerly inaccessible strong, highly non-linear and dynamical regime of gravity. Since the strong-field regime is precisely the place to look for evidence of beyond-GR phenomena due to quantum gravity corrections Giddings et al. (2019); Carballo-Rubio et al. (2018), it is important to extract as much physics as possible from the available GW data.

One can adopt either a model-independent or a theory-specific method for testing gravity, although the former is more efficient if one wishes to achieve constraints on multiple theories with GW observations. One of the first works on the theory-agnostic approach was taken in Arun et al. (2006a, b); Mishra et al. (2010) where each post-Newtonian (PN) term in the GR waveform phase were treated independent and the authors proposed to study the consistency among them. One drawback of such an approach is that it cannot capture the non-GR effects entering at PN orders that are absent in GR (like $-1$PN order common in scalar-tensor theories). To overcome this, Yunes and Pretorius proposed a new framework called parametrized post-Einsteinian (PPE) formalism by introducing generic corrections at any PN order to both the phase and the amplitude Yunes and Pretorius (2009); Chatziioannou et al. (2012). A theory-agnostic data analysis pipeline named TIGER has been developed Agathos et al. (2014); Meidam et al. (2014) and the LIGO and VIRGO Scientific Collaboration (LVC) recently employed the generalized IMRPhenom (gIMR) waveform model which has a one-to-one mapping with the PPE formalism in the inspiral part of the waveform phase. With such waveforms, tests of gravity with the GW phase have been carried out in Abbott et al. (2017a, 2016b); Yunes et al. (2016); Abbott et al. (2017b, 2019, 2016a).

Many of the previous studies on tests of GR with GWs focused only on the phase corrections, though scenarios where amplitude corrections bear importance are not uncommon. In some parity-violating theories, one of the circularly-polarized modes is amplified while the other one is suppressed, an effect called amplitude birefrigence Alexander et al. (2008); Yunes and Finn (2009); Yunes et al. (2010a); Yagi and Yang (2018). Such an effect enters only in the GW amplitude of circularly-polarized modes. Probing amplitude corrections is also important in constraining gravitational theories with GW stochastic backgrounds Maselli et al. (2016). Furthermore, theories with flat extra dimensions Cardoso et al. (2003), Horndeski gravity Saltas et al. (2014), and $f(R)$ gravity Hwang and Noh (1996) may predict amplitude damping that scales with the cosmological distance. Such phenomena have been studied in Ref. Nishizawa (2018) in terms of a generalized GW propagation framework. Possible bounds on the PPE amplitude parameters at various PN orders were studied in Cornish et al. (2011) while both amplitude and phase corrections were included in Arun (2012) for generic theories with scalar dipole radiation.

We here study how much impact the amplitude corrections may bring to tests of GR with GWs and provide justifications for previous studies that only considered the phase corrections. PPE amplitude corrections due to generation mechanisms in various example theories have been derived analytically in Ref. Tahura and Yagi (2018). We compute the constraints on some of those theories from both the phase and amplitude, focusing on leading PN corrections to the tensorial modes only. We choose theories where the leading correction enters at a negative PN order, and the sensitivities of black holes (BHs) are known. Such criteria lead us to choose Einstein-dilaton-Gauss-Bonnet (EdGB) gravity, scalar-tensor theories, and varying-$G$ theories. We carry out Fisher analyses with Monte-Carlo simulations utilizing the parameter posterior samples of GW151226 and GW150914 released by LVC lig 111We choose GW151226 and GW150914 as representatives of low-mass and massive binaries respectively, following Yunes et al. (2016).. Such analyses with actual posterior samples produce more reliable results compared to the ones with sky-averaged waveforms. In fact, when implementing such samples, we can determine the credibility of the small coupling approximation in scalar-tensor theories, which allows us to place reliable bounds on the time-evolution of the scalar field from GW observations for the first time.

We find that the constraints derived from the phase and the amplitude can be comparable in case of massive binary systems like GW150914. Whereas for less massive binaries with a larger number of GW cycles, the phase always yields stronger constraints. Moreover, inclusion of an amplitude correction to the waveform impacts the bound on the phase correction as well since the former can easily be related to the latter provided the dissipative and conservative corrections do not enter at the same order. The amount and direction of such effects vary with the PN order of the corrections. All such constraints in the theories under consideration are summarized in Table 1.

The rest of the paper is organized as follows. Section II.1 briefly reviews PPE formalism while Sec. II.2 summarizes the data analysis techniques. Section III.1 is devoted to justifying our formalism against the one by LVC in massive gravity Abbott et al. (2016a) while we derive constraints on EdGB, scalar-tensor, and varying-$G$ theories in Secs. III.2- III.4. Section IV presents a summary of our work while discussing the effects of an amplitude correction on that of phase. Appendix A compares the PhenomB and PhenomD waveforms for constraining PPE parameters.

II Methodology

In this section, we explain how we perform our analysis. We first explain the PPE formalism and the non-GR waveform template. We then describe the Fisher analysis and how we construct probability distributions of non-GR parameters.

II.1 PPE Waveform

We begin by reviewing the PPE formalism briefly. PPE gravitational waveform for a compact binary inspiral in the frequency domain is given by Yunes and Pretorius (2009)

[TABLE]

where $\tilde{h}_{{\mbox{\tiny GR}}}$ is the gravitational waveform in GR. $\alpha_{{\mbox{\tiny PPE}}}\,u^{a}$ is a correction to the GW amplitude with $u\equiv(\pi\mathcal{M}f)^{1/3}$, $\mathcal{M}\equiv(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the chirp mass with component masses $m_{1}$ and $m_{2}$, and $f$ is the frequency of the GW. The constant $\alpha_{{\mbox{\tiny PPE}}}$ controls the overall magnitude of the correction, while the index $a$ specifies at which PN order the correction enters. One can write the non-GR phase correction $\delta\Psi$ in a similar manner as that of the amplitude as

[TABLE]

Together $\left(\alpha_{{\mbox{\tiny PPE}}},a\right)$ and $\left(\beta_{{\mbox{\tiny PPE}}},b\right)$ are called the PPE parameters.

PPE modifications in Eq. (1) can enter through the non-GR corrections to the binding energy and the GW luminosity Yunes and Pretorius (2009); Chatziioannou et al. (2012), or alternatively to the frequency evolution and the Kepler’s law Tahura and Yagi (2018). We will follow the latter approach and write the modified Kepler’s law as

[TABLE]

and the frequency evolution as

[TABLE]

Here, $\left(\gamma_{r},c_{r}\right)$ and $\left(\gamma_{\dot{f}},c_{\dot{f}}\right)$ parametrize the non-GR corrections to the binary separation $r$ and the frequency evolution $\dot{f}$ respectively. To leading PN order, the GR contribution is given by Cutler and Flanagan (1994); Blanchet et al. (1995)

[TABLE]

where $m$ represents the total mass of the binary and $\Omega=\pi f$ is the orbital angular frequency.

Utilizing the stationary phase approximation Damour et al. (2000); Yunes et al. (2009) and the quadrupole formula for the metric perturbation Blanchet (2002), one can easily derive the amplitude and phase of the dominant quadrupolar mode in Fourier space from Eqs. (3) and (4) as

[TABLE]

and

[TABLE]

respectively. Eq. (7) is already in the PPE format, while Eq. (6) can be reduced to such a form by keeping only the dominant correction222A detailed derivation can be found in Ref. Tahura and Yagi (2018).

In fact, the PPE phase and the amplitude parameters may be related as follows. If the dissipative correction (correction entering in the GW luminosity) dominates over the conservative one (correction entering in the binding energy and Kepler’s law), we find

[TABLE]

while for the conservative-dominated case, we obtain

[TABLE]

On the other hand, when the aforementioned corrections enter at the same PN order, no direct relation between $\alpha_{{\mbox{\tiny PPE}}}$ and $\beta_{{\mbox{\tiny PPE}}}$ exists. The exponents $a$ and $b$ in the correction terms are related by the following equation which is valid for all three cases:

[TABLE]

The above formalism needs to be slightly modified for theories containing time-varying gravitational constants. Variations in the gravitational constants cause the masses of the binary components to vary as well Nordtvedt (1990), and one needs to take this into account when deriving the PPE parameters Tahura and Yagi (2018).

II.2 Data Analysis Formalism

We adopt a Fisher analysis Cutler and Flanagan (1994) to estimate the statistical errors of the non-GR parameters in various theories. Such an analysis is valid for GW events with sufficiently large signal-to-noise (SNR) ratios. We make the assumptions that the detector noise is Gaussian and stationary. Let us write the detector output as

[TABLE]

where $h(t)$ and $n(t)$ are the GW signal and the noise respectively. Let us also define the inner product of two quantities $A(t)$ and $B(t)$ as

[TABLE]

Here $\tilde{A}(f)$ is the Fourier component of $A$, an asterisk ($*$) superscript means the complex conjugate and $S_{n}\left(f\right)$ is the noise spectral density. With the above definitions, the probability distribution of the noise can be written as

[TABLE]

and the SNR for a given signal $h(t)$ can be defined as

[TABLE]

Under the assumptions of Gaussian and stationary noise, the posterior probability distribution of binary parameters $\theta^{a}$ takes the following form:

[TABLE]

where $\Delta\theta^{a}=\hat{\theta}^{a}-\theta^{a}$ with $\hat{\theta}^{a}$ being the maximum likelihood values of $\theta^{a}$. $p^{(0)}\left(\theta^{a}\right)$ gives the probability distribution of the prior information, which we take to be in a Gaussian form for simplicity. $\Gamma_{ab}$ is called the Fisher information matrix which is defined as

[TABLE]

where $\partial_{b}\equiv\frac{\partial}{\partial\theta^{b}}$. One can estimate the root-mean-square of $\Delta\theta^{a}$ by taking the square root of the diagonal elements of the inverse Fisher matrix $\Sigma^{ab}$:

[TABLE]

where $\tilde{\Gamma}_{ab}$ is defined by

[TABLE]

To save computational time, we use IMRPhenomB waveform. Reference Yunes et al. (2016) showed that the difference in constraints on PPE phase parameters between IMRPhenomB and IMRPhenomD waveforms are negligible for propagation mechanisms at any PN order and for generation mechanisms at negative PN orders. In App. A, we perform a similar comparison for generation mechanism corrections in the amplitude using sky-averaged waveforms and show that the former is at least suitable for constraining generation mechanisms that enter at negative PN orders, which is what we will consider in Sec. III.

We choose the following parameters as our variables for the Fisher analysis:

[TABLE]

where $\mathcal{M}_{z}$ is the redshifted chirp mass, $\eta\equiv m_{1}m_{2}/(m_{1}+m_{2})^{2}$ is the symmetric mass ratio and $\chi$ is the effective spin parameter333The effective spin parameter is defined as $\chi\equiv\left(m_{1}\chi_{1}+m_{2}\chi_{2}\right)/\left(m_{1}+m_{2}\right)$, where $\chi_{A}$ with $A=(1,2)$ is the dimensionless spin of the $A$th body.. $\alpha,\delta,\psi$, and $\iota$ are the right ascension, declination, polarization and inclination angles respectively in the detector frame. The non-GR parameter is represented by $\theta_{\mbox{\tiny PPE}}=\alpha_{{\mbox{\tiny PPE}}}$ or $\beta_{{\mbox{\tiny PPE}}}$. We perform a Monte Carlo simulation by using each set of the posterior samples released by LIGO lig for $(\mathcal{M}_{z},\eta,D_{L},\chi,\alpha,\delta,\iota)$, while we randomly sample the polarization angle $\psi$ and the coalescence phase $\phi_{0}$ in $[0,\pi]$ and $[0,2\pi]$ respectively. We impose prior information such that $-1\leq\chi\leq 1$, $-\pi\leq(\phi_{0},\alpha,\psi)\leq\pi$, and $-\pi/2\leq(\delta,\iota)\leq\pi/2$.

We use the detector sensitivity of Advanced LIGO (aLIGO) O1 run Abbott et al. (2018b), and we consider the two detectors at Hanford and Livingston. For simplicity, we assume that the Livingston noise spectrum is identical to that of Hanford Yunes et al. (2009). For the Fisher integration, the minimum frequency is taken to be 20 Hz while the maximum frequency is same as the cutoff frequency above which the signal power is negligible Ajith et al. (2011).

Now we are going to discuss how we compute the probability distribution of a non-GR parameter from the output of a Fisher analysis with a Monte Carlo simulation. We set the fiducial value of any non-GR parameter to be zero for our analysis. We perform the following integration numerically to obtain the compound probability density function444If the distribution of a random variable $y$ depends on a parameter $x$, and if $x$ follows a certain distribution $P(x)$ (called the mixing or latent distribution), the marginal distribution of $y$ is called mixture distribution or compound probability distribution and is given by $P\left(y\right)=\int P\left(y|x\right)P\left(x\right)dx$ Röver and Friede (2016). of any parameter $\xi$ :

[TABLE]

where $P\left(\xi\right)$ is the marginal (unconditional) probability density function of $\xi$. $P\left(\xi|\sigma_{\xi}\right)\propto\exp[-(\xi-\bar{\xi})^{2}/2\sigma_{\xi}^{2}]$ is the conditional probability density function of $\xi$ which we assume to be a Gaussian distribution with a mean $\bar{\xi}$ and a standard deviation $\sigma_{\xi}$. $P(\sigma_{\xi})$ is the probability distribution of $\sigma_{\xi}$ computed from the Fisher analysis for the entire posterior distribution.

Let us finish this section by explaining how we can utilize both amplitude and phase corrections to derive constraints on some theory. One can include $\alpha_{{\mbox{\tiny PPE}}}$ or $\beta_{{\mbox{\tiny PPE}}}$ as variables to the Fisher analysis as in Eq. (II.2) and map them to a non-GR parameter of a theory to derive constraints from the phase and amplitude independently. We refer to such constraints as the “phase-only” and “amplitude-only” bounds respectively. How can we achieve a single constraint that accommodates both of them? Recall the relations between the PPE parameters in Sec. II.1. One can rewrite $\alpha_{{\mbox{\tiny PPE}}}$ in the waveform in terms of $\beta_{{\mbox{\tiny PPE}}}$ according to Eqs. (8) or (9) and eliminate the former variable from the analysis. We refer to such constraints as the “phase & amplitude combined” bounds555Alternatively, one can rewrite the PPE corrections in the phase and the amplitude in terms of non-GR parameters of a theory. Performing Fisher analyses with such parameters as variables lead to similar constraints as the ‘phase & amplitude combined” bounds, although such an approach is not theory-agnostic..

III Results in Example Theories

We now apply our analysis to some example theories. We begin by studying massive gravity that yields corrections in the phase through propagation mechanisms. We compare bounds from the Fisher analysis to those from LVC’s Bayesian analysis to justify the former. We next study EdGB gravity, scalar-tensor theories and varying-$G$ theories, which achieve the corrections through generation mechanisms entering at negative PN orders. Bounds on these theories are summarized in Table 1.

III.1 Validation of the Fisher Analysis: Massive Gravity

The idea of introducing the mass to gravitons is rather old Fierz and Pauli (1939), and many attempts have been made to construct a feasible theory that allows one to do so de Rham (2014). Such a theory may arise in higher-dimensional setups Hinterbichler (2012) and has the potential to solve the cosmic acceleration problem de Rham (2014). Although gravitons with non-vanishing masses may have additional polarizations as well de Paula et al. (2004), we here restrict our attention to the non-GR effects on the tensor modes due to a massive dispersion relation.

We will focus on the non-GR corrections specifically to the GW phase. Thus, the purpose of this section is simply to compare our Fisher analysis with the Bayesian one performed by the LVC. Gravitons with a non-vanishing mass travel at a speed smaller than the speed of light and the non-GR effects accumulate over the distance. Modified dispersion relation for such gravitons is given by $E^{2}=p^{2}c^{2}+m_{g}^{2}c^{4}$, where $m_{g}$ is the mass of the graviton while $E$ and $p$ are the energy and the momentum respectively. The PPE phase parameters are Will (1998)

[TABLE]

where

[TABLE]

Here, $\Omega_{M}$ and $\Omega_{\Lambda}$ are the energy density of matter and dark energy respectively. $H_{0}$ is the Hubble constant while $z$ is the redshift of the source. $\lambda_{g}$ is the Compton wavelength of the graviton that is related to $m_{g}$ as $\lambda_{g}\equiv h/\left(m_{g}c\right)$, where $h$ is Planck’s constant.

We compute the probability distribution of $\lambda_{g}$ from GW150914 according to the procedure outlined in Sec. II and compare with the one obtained by the LVC Abbott et al. (2016a) (Fig. 1). The Fisher analysis with Monte Carlo simulations yields $\lambda_{g}<1.2\times 10^{13}$ km at 90% CL, which is in a good agreement with the LVC bound of $1.0\times 10^{13}$ km and thus shows the validity of the former. The difference in the two cumulative distributions of $\lambda_{g}$ presented in Fig. 1 can be attributed to the fact that the LVC used a more accurate Bayesian analysis and imposed a uniform prior on the graviton mass. The GW bound has recently been updated by combining multiple events Abbott et al. (2019). The new bound is stronger than binary pulsar constraints Finn and Sutton (2002); Miao et al. (2019) but slightly weaker than the updated solar system bounds Will (2018). The bound is also weaker than the ones from the observations of galactic clusters Goldhaber and Nieto (1974); Gupta and Desai (2018); Desai (2018), gravitational lensing Choudhury et al. (2004), and the absence of superradiant instability in supermassive BHs Brito et al. (2013).

III.2 Einstein-Dilaton-Gauss-Bonnet Gravity

EdGB gravity endows one of the simplest high-energy modifications to GR Moura and Schiappa (2007); Pani and Cardoso (2009). Such a theory is motivated from low-energy effective string theories and also arises as a special case of Horndeski gravity Zhang et al. (2017); Berti et al. (2015). The EdGB action is given by introducing a quadratic-curvature correction (Gauss-Bonnet invariant) to the GR action which is non-minimally coupled to a scalar field (dilaton) with a coupling constant $\bar{\alpha}_{\mbox{\tiny EdGB}}$ Kanti et al. (1996)666We use barred quantities for coupling constants in order to distinguish them from the PPE parameters..

In EdGB gravity, BHs acquire scalar monopole charges which may generate scalar dipole radiation if they form binaries Yagi et al. (2012); Sotiriou and Zhou (2014); Berti et al. (2018); Prabhu and Stein (2018). Such radiation leads to an earlier coalescence of BH binaries compared to that of GR and modifies the gravitational waveform with the PPE parameters given by Yunes et al. (2016); Yagi et al. (2012)

[TABLE]

with $b=-7$ and

[TABLE]

with $a=-2$. Here, $\zeta_{\mbox{\tiny EdGB}}\equiv 16\pi\bar{\alpha}_{\mbox{\tiny EdGB}}^{2}/m^{4}$ is the dimensionless EdGB coupling parameter and $\tilde{s}_{A}^{\mbox{\tiny EdGB}}$ are the spin-dependent factors of the BH scalar charges given by $\tilde{s}_{A}^{\mbox{\tiny EdGB}}\equiv 2(\sqrt{1-{\chi_{A}}^{2}}-1+{\chi_{A}}^{2})/{\chi_{A}}^{2}~{}$ Berti et al. (2018); Prabhu and Stein (2018)777For ordinary stars like neutron stars $\tilde{s}_{A}^{\mbox{\tiny EdGB}}$ are zero Yagi et al. (2012, 2016)..

We now derive constraints on EdGB gravity from GW150914 and GW151226. First, we estimate how well these events satisfy the small coupling approximation $\zeta_{{\mbox{\tiny EdGB}}}<1$. To do so, we extract the 90% CL upper bound $\Delta\zeta_{{\mbox{\tiny EdGB}}}$ from each sample of the posterior distribution of a particular event. We then create histograms with all the samples (see Fig. 2) and calculate the fraction satisfying $\Delta\zeta_{{\mbox{\tiny EdGB}}}<1$. For GW150914, 72% (42%) of the samples satisfy the small coupling approximation if $\Delta\zeta_{{\mbox{\tiny EdGB}}}$ is derived from the phase (amplitude) correction only, while 71% of the posterior distribution satisfies such approximation if the phase and amplitude corrections are combined. A similar analysis with GW151226 gives 98% and 87% for the phase and amplitude corrections respectively while combining the two yields almost the same result as that of the phase-only case. Since the fraction of samples satisfying $\zeta_{{\mbox{\tiny EdGB}}}<1$ is much higher for GW151226 than GW150914 due to a larger number of GW cycles and slower relative velocity of the binary constituents, the former event places more reliable constraints on EdGB gravity compared to the latter one.

Figure 3 presents cumulative probability distributions of $\bar{\alpha}^{2}_{{\mbox{\tiny EdGB}}}$888We show the distribution of $\bar{\alpha}^{2}_{{\mbox{\tiny EdGB}}}$ instead of $\sqrt{\bar{\alpha}_{{\mbox{\tiny EdGB}}}}$ as it is the former that directly enters in the waveform. for GW150914 for three different cases with vertical lines representing the 90% CL of the corresponding distribution. We found the 90% CL constraints on $\sqrt{\bar{\alpha}_{{\mbox{\tiny EdGB}}}}$ from each of the phase and amplitude corrections as 50.5 km and 76.3 km respectively. Notice that these bounds have the same order of magnitude. On the other hand, combining the amplitude and phase corrections leads to an upper bound of 51.5 km, which is weaker than the phase-only constraint by 2% (to be discussed more in Sec. IV). Though above constraints may not be reliable as the 90% CL bounds on $\zeta_{\mbox{\tiny EdGB}}$ do not satisfy the small coupling approximation, which is shown in Table 1.

We now look at bounds on GW151226. We found that this event yields 4.32 km and 10.5 km respectively from the phase and amplitude corrections, while combining the two only changes the result from the phase-only case by 0.01%. These bounds are consistent with those in a recent paper Nair et al. (2019) that utilized the LVC posterior samples including the non-GR phase corrections at $-1$PN order while Ref. Yamada et al. (2019) found even stronger bounds by combining multiple GW events. These GW bounds are comparable to the one obtained from low-mass X-ray binaries Yagi (2012).

Although GW150914 leads to weaker constraints on EdGB gravity compared to GW151226, the effect of amplitude correction is more manifest for the former event. This is because GW150914 has a smaller number of GW cycles, and thus the amplitude contribution becomes relatively higher than GW151226.

III.3 Scalar-Tensor Theories

Scalar-Tensor theories of gravity emerge from the dimensional reduction of higher-dimensional theories such as the Kaluza-Klein theory Fujii and Maeda (2007); Overduin and Wesson (1997) and string theories Polchinski (1998a, b). In addition to the spacetime curvature, scalar fields mediate an additional force which is introduced through non-minimal couplings of scalar fields and gravity Berti et al. (2015); Chiba et al. (1997); Harrison (1972). Such theories can explain the accelerated expansion of the universe Brax et al. (2004); Kainulainen and Sunhede (2006); Baccigalupi et al. (2000); Riazuelo and Uzan (2002); Schimd et al. (2005), inflation Burd and Coley (1991); Barrow and Maeda (1990); Clifton et al. (2012), primordial nucleosynthesis Coc et al. (2006); Damour and Pichon (1999); Larena et al. (2007); Torres (1995), and the structure formation Brax et al. (2006).

Certain scalar-tensor theories predict scalarization of neutron stars Damour and Esposito-Farèse (1993); Barausse et al. (2013), which can also happen to BHs if the scalar field evolves with time cosmologically Jacobson (1999); Horbatsch and Burgess (2012). Compact binaries formed by such objects emit dipole radiation which modifies the GW phase with Freire et al. (2012); Wex (2014); Tahura and Yagi (2018)

[TABLE]

with $b=-7$, and the GW amplitude as Tahura and Yagi (2018)

[TABLE]

with $a=-2$. Here $\alpha_{A}$ represents the scalar charge of the $A$th binary component and depends on specific theories and the type of compact objects. If we consider a binary consisting of BHs in a theory where the scalar field $\phi$ obeys a massless Klein-Gordon equation, $\alpha_{A}$ is given by Horbatsch and Burgess (2012)

[TABLE]

where $\dot{\phi}$ is the rate of change of $\phi$ with time.

One can use Eqs. (25)–(27) and the numerical analysis described in Sec. II.2 to find constraints on $\dot{\phi}$ as long as the small coupling approximation $m_{A}\dot{\phi}<1$ is satisfied. In this regard, only 11.7% (13.5%) of the samples of GW150914 satisfies such approximation with 90% confidence level for the phase (combined) correction, while all of the samples fail to do so for the amplitude correction. Hence GW150914 cannot place any meaningful bound on scalar-tensor theories considered here. On the other hand, 90.4% of the samples from GW151226 meets the small coupling criterion for the phase-only and combined analyses, while the fraction is only 25% for the amplitude correction. Thus, we derive reliable constraints from GW151226 with the phase-only and combined analyses, with both leading to $\dot{\phi}<1.1\times 10^{4}$/sec999A previous analysis with the sky-averaged waveform in Ref. Yunes et al. (2016) could not place a reliable bound on scalar-tensor theories. Since the posterior distributions of the GW events were not available then, one could not determine how well those events satisfied the small coupling approximation from a simple Fisher analysis.. This constraint is 10 orders of magnitude weaker than the current most stringent bound obtained from the orbital decay rate of quasar OJ287 Horbatsch and Burgess (2012), though this is the first bound obtained in the strong/dynamical regime.

III.4 Varying- $G$ Theories

Many metric theories of gravity that violate the strong equivalence principle Di Casola et al. (2015); Will (2014); Bertotti and Grishchuk (1990) predict time variation in the gravitational coupling parameter $G$ Uzan (2011). Scalar-tensor theories are examples where $G$ varies as a function of the asymptotic scalar field Will (2006), which may vary over time. Any such time-dependence of $G$ leads to a variation in the effective masses of compact bodies, which in turn makes them experience anomalous cosmic acceleration Nordtvedt (1990). Such phenomena alter the gravitational waveform through the modifications of the binary orbital evolution and the energy balance law Tahura and Yagi (2018).

We now show the PPE modifications due to a time variation in the gravitational constant. In fact, the amount of gravitational coupling that appears in different sectors of a gravitational theory may not be unique. Einstein-Æther theory Yagi et al. (2014) and Brans-Dicke theory with a cosmologically evolving scalar field Will (2006) are examples of such theories in which various gravitational constants exist. Reference Tahura and Yagi (2018) studied a generic case with two distinct gravitational constants in Kepler’s law (conservative sector) and GW luminosity (dissipative sector). Here we place constraints on the special case where these two constants coincide with each other. Let the masses and the Newton’s constant vary according to the following equations:

[TABLE]

Here, a subscript [math] denotes that the quantity is measured at the time $t=t_{0}$, and an overhead dot means a derivative with respect to time. Eqs. (28) and (29) modify the GW phase and amplitude as Tahura and Yagi (2018)

[TABLE]

with $b=-13$, and

[TABLE]

with $a=-8$ respectively. Here $s_{A}$ is the sensitivity of the $A$th binary component defined as

[TABLE]

Employing Eqs. (30) and (31), GW150914 (GW151226) imposes constraints on $|\dot{G_{0}}/G_{0}|$ from the phase-only and amplitude-only analyses as $7.30\times 10^{6}\,\mathrm{yr}^{-1}$ ($2.24\times 10^{4}\,\mathrm{yr}^{-1}$) and $1.37\times 10^{8}\,\mathrm{yr}^{-1}$ ($3.82\times 10^{5}\,\mathrm{yr}^{-1}$) respectively, with the combined analyses yielding slight improvements over the phase-only results. Unlike the EdGB and scalar-tensor cases, the amplitude-only analyses yield much worse bound than that from the phase-only cases even with GW150914. Notice also that for varying-$G$ theories, the combined bound is slightly stronger than the phase-only bound (to be discussed more in the subsequent section). These bounds are much less stringent compared to the other contemporary constraints Will (2006). However, future space-borne detectors such as LISA Seoane et al. (2013); Audley et al. (2017) will be able to obtain constraints up to 13 orders of magnitude stronger compared to the aLIGO ones Yunes et al. (2010b); Chamberlain and Yunes (2017).

IV Conclusion and Discussion

In this analysis, we have derived constraints on scalar-tensor, EdGB and varying-$G$ theories from GW150914 and GW151226. To do so, we performed Fisher analyses with Monte-Carlo simulations using the posterior samples constructed by LVC. In particular, we derived reliable constraints on the time-evolution of the scalar field in scalar-tensor theories from GW observations for the first time.

We explored how amplitude corrections contribute to the constraints on such theories. We derived three sets of bounds on each theory: phase-only, amplitude-only, and from both phase and amplitude combined. We found that for binaries with large masses such as GW150914, where we have less number of cycles, the bounds from the amplitude and phase can be comparable to each other. On the other hand, combined analyses yield constraints that differ from the phase-only case at most by 3.6% for the theories under consideration. Hence, at least in theories where the leading corrections enter at negative PN orders, the phase-only analyses as done in previous literature Yunes et al. (2016); Chamberlain and Yunes (2017); Nair et al. (2019); Yamada et al. (2019); Carson et al. (2019) can produce sufficiently accurate constraints.

Depending on the prior information and the PN order of the non-GR correction, a combined analysis can yield stronger or weaker constraint compared to a phase-only one. With the priors mentioned in Sec. II.2, the fractional difference between $\beta_{{\mbox{\tiny PPE}}}$ for the two cases is presented in Fig. 4. From $-4$ PN to $-2.5$ PN correction, the combined analyses give rise to slight improvements over the phase-only constraints, while for other cases, the former is weaker with a maximum deterioration of 8.5% at $-1$ PN order. Nonetheless, it would be safer to include both phase and amplitude corrections in the analysis as a lack of the former in the waveform may cause systematic errors on GR parameters such as luminosity distance if non-GR corrections exist in nature.

In this paper, we considered only the leading PN corrections in the inspiral part of the waveform, but how important are higher-PN corrections and modifications in the merger-ringdown portion? Reference Yunes et al. (2016) partially addressed this question by taking Brans-Dicke theory as an example whose leading correction enters at $-1$PN order, similar to EdGB gravity and scalar-tensor theories considered here. Appendix B of  Yunes et al. (2016) shows that including higher-PN corrections only affects the bound from the leading PN correction by 10% at most for GW150914. Moreover, for EdGB gravity, including the correction to the black hole ringdown frequency and damping time only affects the bound from the leading PN corrections in the inspiral by 4.5% for GW150914 Carson and Yagi . Thus, it is likely that the bounds presented here are valid as order-of-magnitude estimates.

A possible avenue for future work includes repeating the calculation presented here but with a Bayesian analysis using a more accurate waveform such as PhenomD, PhenomPv2 or effective-one-body ones. In particular, it would be interesting to investigate whether the amplitude correction contribution entering at positive PN orders is negligible like the negative PN cases reported here. It is also interesting to repeat the analysis here to all the other events in GWTC-1 Abbott et al. (2018b) and study how much the bounds on each theory improve by combining these events. Another possibility is to take into account non-tensorial polarization modes following e.g. Chatziioannou et al. (2012). For future detectors with improved sensitivities, the ability to measure amplitude corrections may be dominated by calibration errors101010The calibration error on the amplitude for the O2 run was 3.8% Abbott et al. (2018b)..

Acknowledgements.

We thank Carl-Johan Haster for providing valuable comments on the LIGO posterior samples. K.Y. and Z.C. acknowledge support from NSF Award PHY-1806776. K.Y. would like to also acknowledge support by the COST Action GWverse CA16104 and JSPS KAKENHI Grants No. JP17H06358.

Appendix A Comparison of Bounds on PPE parameters with PhenomB and PhenomD Waveforms

Even though the PhenomD waveform produces more accurate results, we utilized the PhenomB one in this paper throughout because the latter is simpler and saves computational time when performing Monte Carlo simulations. In this appendix, we compare constraints on the PPE parameter $\alpha_{\mbox{\tiny PPE}}$ from both waveforms to justify our method.

Let us discuss the distinct features of the two waveforms first. Both PhenomB and PhenomD waveforms are spin-aligned (non-precessing) frequency-domain phenomenological models of gravitational waveforms Ajith et al. (2011); Khan et al. (2016). The PhenomB waveform is calibrated for mass ratios up to $m_{1}/m_{2}=4$ and spin components of $\chi_{i}\in[-0.85,0.85]$ are unified into a single effective spin. On the other hand, the PhenomD waveform covers a larger region of the parameter space with mass ratios upto 18 and spins of $\chi_{i}\in[-0.95,0.95]$, with both spins introduced independently. The waveform contains a much higher order in PN terms in the inspiral than the PhenomB waveform and further introduces an intermediate phase connecting the inspiral and merger-ringdown portions, which make such waveforms more reliable than the PhenomB ones.

We now estimate the constraints on the PPE amplitude modification from the two waveforms. Since modifications to propagation mechanisms used for massive gravity in Sec. III.1 do not give rise to amplitude corrections, we here focus on modifications to generation mechanisms. We performed Fisher analyses with sky-averaged PhenomB and PhenomD waveforms and derived upper bounds on $\alpha_{{\mbox{\tiny PPE}}}$ at different PN orders. As shown in Fig. 5, the results from the two waveforms agree very well at negative PN corrections but deviate from each other at the positive ones. On the other hand, truncating the Fisher analyses at the end of the inspiral phase show significant agreement between the two waveforms at positive PN orders (Fig. 6), suggesting that the deviation in Fig. 5 originates mainly from the intermediate/merger-ringdown portion.

The example theories considered in this paper acquire leading non-GR corrections either from propagation effects or from the generation effects with the latter entering in negative PN orders. A comparison between the PhenomB and PhenomD results for constraining $\beta_{{\mbox{\tiny PPE}}}$ presented in Ref. Yunes et al. (2016) reveals consistency on constraining modifications to propagation mechanisms at both positive and negative PN orders, while the two waveforms show agreement only at negative PN orders for constraining modifications to generation mechanisms. Together with the results on amplitude corrections discussed above confirms that the results of this paper should not change significantly if one utilizes the PhenomD waveform instead. On the other hand, for constraining theories like dynamical Chern-Simons or noncommutative gravity where the leading correction enters at a positive PN order, the PhenomB waveform is not expected to produce reliable results.