Testing massive-field modifications of gravity via gravitational waves

TL;DR

This paper uses gravitational wave data from LIGO-Virgo to test for massive-field modifications of gravity, especially dipole radiation, and places new constraints on Einstein-dilaton-Gauss-Bonnet coupling parameters.

Contribution

It extends the parametrized post-Einsteinian framework to include massive fields and provides the first constraints on EdGB coupling from gravitational wave observations.

Findings

Constraints on dipole radiation activation in the mass range $10^{-14}$ to $10^{-13}$ eV.

First 90\% confidence level bounds on EdGB coupling parameter from gravitational wave data.

Upper limit of $oxed{ ext{~}2.47 ext{ km}}$ on $oxed{ ext{~} oot{ ext{~}4} ext{~} ext{of } ext{~} ext{ } ext{ extit{ extbf{ ext{alpha}}}_{ ext{EdGB}}}}$ for masses less than $6 imes 10^{-14}$ eV.

Abstract

The direct detection of gravitational waves now provides a new channel of testing gravity theories. Despite that the parametrized post-Einsteinian framework is a powerful tool to quantitatively investigate effects of modification of gravity theory, the gravitational waveform in this framework is still extendable. One of such extensions is to take into account the gradual activation of dipole radiation due to massive fields, which are still only very weakly constrained if their mass $m$ is greater than $10^{-16}$ eV from pulsar observations. Ground-based gravitational-wave detectors, LIGO, Virgo, and KAGRA, are sensitive to this activation in the mass range, $10^{-14}$ eV $\lesssim m \lesssim 10^{-13}$ eV. Hence, we discuss a dedicated test for dipole radiation due to a massive field using the LIGO-Virgo collaboration's open data. In addition, assuming Einstein-dilaton-Gauss-Bonnet (EdGB) type coupling, we combine the results of the analysis of the binary black hole events to obtain the 90\% confidence level constraints on the coupling parameter $\alpha_{\rm EdGB}$ as $\sqrt{\alpha_{\rm EdGB}} \lesssim 2.47$ km for any mass less than $6 \times 10^{-14}$ eV for the first time, including $\sqrt{\alpha_{\rm EdGB}} \lesssim 1.85$ km in the massless limit.

Tables

Table 1. Table 1: The GR “best-fit parameters” in the detector frame

event	$(f_{\min },f_{\max })$/Hz	$ℳ/M_{\odot }$	$\eta$	${\chi }_1$	${\chi }_2$	SNR
GW150914	(20, 1024)	31.2	0.249	0.79	-1.00	24.4
GW151012	(20, 1024)	18.0	0.249	-0.36	0.22	9.00
GW151226	(20, 1024)	9.70	0.211	0.50	-0.46	12.0
GW170104	(20, 1024)	24.9	0.243	-0.87	1.00	13.2
GW170608	(20, 1024)	8.48	0.250	0.47	-0.49	15.4
GW170729	(20, 1024)	48.4	0.201	0.74	-1.00	10.4
GW170809	(20, 1024)	30.5	0.228	0.72	-1.00	12.0
GW170814	(20, 1024)	26.9	0.246	0.85	-1.00	16.2
GW170817	(23, 2048)	1.20	0.231	-0.41	0.96	32.2
GW170818	(16, 1024)	32.9	0.250	0.98	-1.00	10.5
GW170823	(20, 1024)	39.3	0.249	0.99	-1.00	11.5

Table 2. Table 2: Fisher-estimated constraints on EdGB gravity model from GW151226-like and GW170608-like signals. The lower layer shows the case when all parameters are unconstrained, while the upper layer takes into account the prior information that χ 1 subscript 𝜒 1 \chi_{1} and χ 2 subscript 𝜒 2 \chi_{2} must necessarily be smaller than unity. The measurement accuracy of symmetric mass-ratio and correlation coefficient between η 𝜂 \eta and α EdGB subscript 𝛼 EdGB \sqrt{\alpha_{\rm EdGB}} are also shown.

signal	$\sqrt{{\alpha }_{\mathrm{EdGB}}}$ [km]	$\Delta \eta /\eta$	$c_{\eta ,\sqrt{{\alpha }_{\mathrm{EdGB}}}}$
with spin prior
GW151226-like	3.35	361%	-0.708
GW170608-like	2.96	120%	-0.245
without spin prior
GW151226-like	5.80	11600%	0.968
GW170608-like	5.23	7110%	0.957

Equations

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

$${\Psi }_{\mathrm{GR}}=2\pi f{t}_0-{\varphi }_0-\frac{\pi }{4}+\frac{3}{128}{u}^{-5}+\cdots$$

$$u\equiv (\pi \mathcal{M}f{)}^{1/3},$$

$$\delta \Psi =-\beta u^b.$$

$$\dot{f}(f)$$

$${\dot{f}}_{\mathrm{GR}}(f)=\frac{96}{5\pi {\mathcal{M}}^2}{u}^{11}+\cdots .$$

$$\varphi (f)$$

$$={\int }^f\frac{2\pi f^{\prime }}{{\dot{f}}_{\mathrm{GR}}(f^{\prime })}\left(1-\frac{\delta {\dot{f}}^{\prime }}{{\dot{f}}_{\mathrm{GR}}(f^{\prime })}\right)d{f}^{\prime },$$

$$\bar{t}(f)$$

$$={\int }^f\frac{1}{{\dot{f}}_{\mathrm{GR}}(f^{\prime })}\left(1-\frac{\delta \dot{f}(f^{\prime })}{{\dot{f}}_{\mathrm{GR}}(f^{\prime })}\right)d{f}^{\prime },$$

$$\tilde{h}(f)\simeq \frac{\mathcal{A}(f)}{2\sqrt{\dot{f}}}{e}^{i\left[\varphi (\bar{t})-2\pi f\bar{t}\right]},$$

$$\tilde{h}(f)={\tilde{h}}_{\mathrm{GR}}(1+\delta \mathcal{A})e^{i\delta \Psi }$$

$$\delta \mathcal{A}$$

$$\delta \Psi$$

$$\frac{(dE/dt{)}_D}{(dE/dt{)}_Q}=A\frac{({\omega }^2-m^2{)}^{3/2}}{{\omega }^3}\Theta ({\omega }^2-m^2)u^{-2},$$

$$\frac{(dE/dt{)}_D}{(dE/dt{)}_Q}=A{u}^{-2}.$$

$$\frac{\delta \dot{f}}{{\dot{f}}_{\mathrm{GR}}}\simeq \frac{(dE/dt{)}_D}{(dE/dt{)}_Q}=A\frac{({\omega }^2-m^2{)}^{3/2}}{{\omega }^3}\Theta ({\omega }^2-m^2)u^{-2}.$$

$$\frac{\delta \dot{f}}{{\dot{f}}_{\mathrm{GR}}}\simeq A{u}^{-2}.$$

$$\delta \Psi$$

$$\simeq \frac{5}{48}A{\left(\frac{1}{m\mathcal{M}}\right)}^{7/3}\left(F(\hat{f})+G(\hat{f})\right),$$

$$F(\hat{f})$$

$$G(\hat{f})$$

$$G(\hat{f})\simeq 0.0808004\hat{f}-0.153688$$

$$F(\hat{f})\simeq \{\begin{array}{ll}\text{vspace}0.2cm\frac{9}{70}{\hat{f}}^{-7/3}\left(1-\frac{105}{208}{\hat{f}}^{-2}+\frac{105}{1672}{\hat{f}}^{-4}+\frac{1}{160}{\hat{f}}^{-6}+\frac{105}{67456}{\hat{f}}^{-8}\right),& (\hat{f}\gtrsim 1.36)\\ {\left(\frac{3.12}{(\hat{f}-1{)}^{7/2}}+\frac{9.2}{(\hat{f}-1{)}^{5/2}}+\frac{8.4}{(\hat{f}-1{)}^2}+\frac{1}{0.08(\hat{f}-1)}\right)}^{-1}-G(\hat{f}).& (\hat{f}\lesssim 1.36)\end{array}$$

$$\delta \mathcal{A}$$

$$\delta \Psi$$

$$\frac{(dE/dt{)}_D}{(dE/dt{)}_Q}=\frac{5}{96}{\zeta }_{\mathrm{EdGB}}\frac{(m_1^2{s}_2^{\mathrm{EdGB}}-m_2^2{s}_1^{\mathrm{EdGB}}{)}^2}{M^4{\eta }^{18/5}}{u}^{-2}.$$

$$s_i^{\mathrm{EdGB}}$$

$$A=\frac{5\pi }{6}\frac{{\alpha }_{\mathrm{EdGB}}^2}{M^4}\frac{(m_1^2{s}_2^{\mathrm{EdGB}}-m_2^2{s}_1^{\mathrm{EdGB}}{)}^2}{M^4{\eta }^{18/5}}.$$

$${\Gamma }_{ij}\equiv \left(\frac{\partial \hat{h}}{\partial {\theta }^i}|\frac{\partial \hat{h}}{\partial {\theta }^j}\right),$$

$$\Delta {\theta }^i\equiv \sqrt{⟨({\theta }^i{)}^2⟩}=\sqrt{({\Gamma }^{-1}{)}^{ii}}.$$

$$\gamma =\frac{48}{5{\eta }^{2/5}}A.$$

$$\left(\frac{d}{d{t}^2}+{\omega }^2\right)\phi (t)$$

$$\left({\nabla }^2+k^2\right)\psi (\mathbf{r})$$

Full text

1 Introduction

The opening of gravitational wave (GW) astronomy/astrophysics allows us to approach a new test of general relativity (GR) in strong gravity regime. Regarding O1 and O2 by LIGO-Virgo collaboration (LVC), ten binary black hole (BH) mergers and one neutron star (NS) binary merger are summarized in the catalog GWTC-1 LIGOScientific:2018mvr . Testing GR has been investigated by several authors using these event data and no significant deviation from GR has been reported TheLIGOScientific:2016src ; Yunes:2016jcc ; Abbott:2018lct ; LIGOScientific:2019fpa . Although one of the most interesting regime to investigate is the ringdown phase Brito:2018rfr ; Nakano:2018vay , the concrete templates of merger-ringdown phase waveform in modified gravity theories are not established yet Berti:2018vdi . On the other hand, the inspiral phase can be studied by using the post-Newtonian (PN) approximation.

Although a possible test is the one in the parametrized post-Einsteinian (ppE) framework, in which the so-called ppE parameters are introduced to describe modifications of gravity theory in a generic way Yunes:2009ke ; Chatziioannou:2012rf ; Sampson:2013lpa ; Loutrel:2014vja ; Huwyler:2014gaa , a waveform in this framework does not cover the whole viable extensions of gravity. One of such extensions is to consider a massive field which is coupled to a compact object through an additional charge exciting dipole radiation in the frequency range above the mass scale of the field. Such an effect can let the modification of the gravitational waveform turn on as the frequency increases during the inspiral Sampson:2013jpa ; Sampson:2014qqa .

Spontaneous scalarization of compact objects for massless fields was first investigated by Damour and Esposito-Fareèse Damour:1993hw . They showed that even scalar-tensor theories which pass the present weak-field gravity tests exhibit nonperturbative strong-field deviations from GR in systems involving NSs, whereas such massless scalar fields have been strongly constrained by the recent observation of a pulsar-white dwarf binary Antoniadis:2013pzd . A possible extension is to add a mass term to the scalar field since the observational bounds can be avoided if the mass of the scalar field is larger than the binary orbital frequency. Such an extension to the spontaneous scalarization is discussed by several authors (e.g. see Refs. Chen:2015zmx ; Ramazanoglu:2016kul ). In the following mass ranges of the field, the coupling is constrained by several observations: (1) $m\lesssim 10^{-27}$ eV for stability on cosmological scales Ramazanoglu:2016kul , (2) $m\lesssim 10^{-16}$ eV from the observations of a pulsar-white dwarf binary Antoniadis:2013pzd , and (3) $10^{-13}$ eV $\lesssim m\lesssim 10^{-11}$ eV Arvanitaki:2014wva ; Brito:2017zvb , which relies on the measurements of high spins of stellar mass BHs Narayan:2013gca .

In the presence of matter, dynamical scalarization in scalar-tensor theories has been investigated Barausse:2012da ; Palenzuela:2013hsa ; Shibata:2013pra ; Taniguchi:2014fqa . In this case, dipole radiation abruptly turns on at a given threshold, where the scalar charge of a NS in a binary grows suddenly. Possible constraints on hidden sectors, such as scalar-tensor gravity or dark matter, which induce a Yukawa-type modification to the gravitational potential are discussed by several authors for future NS binary merger observations Kopp:2018jom ; Alexander:2018qzg . However, such activation of dipole radiation in vacuum spacetimes, i.e., in the case of binary BHs has not been discussed. It may be interesting to consider Einstein-dilaton-Gauss-Bonnet (EdGB) theory with a massive scalar field, in which a BH can have a scalar hair Yunes:2016jcc . Very recently, Nair et al. discussed to constrain EdGB theory from the GW observations Nair:2019iur . Their order-of-magnitude estimation on the constraint is $\sqrt{\alpha_{\rm EdGB}}\lesssim 1.0$ km, where $\alpha_{\rm EdGB}$ is the coupling parameter of EdGB theory, while their Fisher-estimated and Bayesian constraints are roughly $\sqrt{\alpha_{\rm EdGB}}\lesssim 6.0$ km. In addition, Tahura et al. have independently derived constraints on a few modified theories of gravity including EdGB theory and their constraint on EdGB theory is roughly $\sqrt{\alpha_{\rm EdGB}}\lesssim 4.3$ km Tahura:2019dgr . Current constraint on EdGB theory is $\sqrt{\alpha_{\rm EdGB}}\lesssim 1.9$ km, which is obtained by using low-mass X-ray binaries Yagi:2012gp .

In this paper, we discuss activation of dipole radiation due to massive fields and modifications to gravitational waveform. We analyze the GW events in the catalog GWTC-1 GWOSCweb to constrain the magnitude of the modification. In particular, assuming EdGB theory with a massive scalar field, we investigate the bound for the coupling parameter. Since the ground based detectors, LIGO, Virgo, and KAGRA, are sensitive at frequencies around 10-1000 Hz, gravitational-wave signals allow us to investigate the effect of the mass in the range $10^{-14}$ eV $\lesssim m\lesssim 10^{-13}$ eV, while they reduce to the massless limit for $m\lesssim 10^{-15}$ eV. We also evaluate the measurement accuracy of $\alpha_{\rm EdGB}$ by using the Fisher information matrix.

This paper is organized as follows. Section 2 briefly reviews the gravitational waveform of the inspiral phase in the ppE framework. Section 3 presents how the gravitational waveform is modified by a massive field. Section 4 shows the results of our analysis of LVC open data of the gravitational wave catalog using the modified gravitational-wave templates and the Fisher analysis on the determination error of the coupling parameter. Section 5 is devoted to summary and discussion. Throughout this paper we use geometric units in which $G=1=c$.

2 The ppE waveform

In this section, we briefly review the ppE waveform following Ref. Yunes:2009ke ; Tahura:2018zuq . In the ppE framework, the waveform is formulated by considering modifications to the binding energy and GW luminosity, which correspond to conservative and dissipative corrections, respectively. Here, we employ another approach, in which we focus on the correction to the GW frequency evolution $\dot{f}$. The frequency domain (FD) ppE waveform for the inspiral phase of compact binaries is expressed as

[TABLE]

where $\tilde{h}_{\rm GR}$ is the waveform in GR and its phase is given by

[TABLE]

with the coalescence time $t_{0}$ and phase $\phi_{0}$. $\alpha u^{a}$ and $\delta\Psi$ represent modifications to the amplitude and the phase, respectively. We introduce

[TABLE]

where $\mathcal{M}=M\,\eta^{3/5}$ is the chirp mass with the total mass of the binary $M=m_{1}+m_{2}$ and the symmetric mass ratio $\eta=m_{1}m_{2}/M^{2}$. We can also model the phase correction as

[TABLE]

$\alpha$, $\beta$, $a$, and $b$ are the ppE parameters.

In this paper, we focus only on the dominant correction to the energy flux due to the dipole radiation, which appears in -1PN order as is known well, while the correction to the binding energy of the system is the Newtonian order Stein:2013wza ; Yunes:2016jcc . Since only the correction to the phase accumulates through the binary evolution, we neglect the correction to the amplitude and incorporate only the correction to the evolution of the frequency $\dot{f}$. Then, let us denote $\dot{f}$ as

[TABLE]

where the GR term $\dot{f}_{\rm GR}(f)$ is expressed as

[TABLE]

Therefore, we obtain

[TABLE]

at the leading order in $\delta\!\!\>\dot{f}$. In the stationary phase approximation the gravitational waveform in the FD, $\tilde{h}(f)$, can be expressed as

[TABLE]

where $\mathscr{A}(f)$ is the amplitude of time domain (TD) waveform evaluated at $\bar{t}=\bar{t}(f)$. Substituting Eqs. (2) and (2) into Eq. (6), we find

[TABLE]

at the leading order, where

[TABLE]

3 Massive field

At the leading order in the PN expansion, we obtain the ratio between the energy losses due to the dipole radiation and the quadrupole one as (see Appendix A)

[TABLE]

where $\omega=\pi f$ is the angular frequency of the dipole radiation, for which the dispersion relation is $\omega^{2}=k^{2}+m^{2}$ with the wave number $k$ and the mass of the field $m$, $\Theta$ is the Heaviside step function, and $A$ is a parameter which denotes the relative amplitude of the dipole radiation. In the massless limit, we have

[TABLE]

Since the correction to the binding energy of the system is higher order in the PN expansion, we obtain

[TABLE]

In the massless limit, this becomes

[TABLE]

Substituting Eq. (12) into Eq. (9), we obtain the modification to the GW phase as

[TABLE]

where we define

[TABLE]

with $\hat{f}=\pi f/m$. In our matched filtering analysis, $G(\hat{f})$ is irrelevant because

[TABLE]

can be absorbed by the shifts of the merger time and phase (see the first and the second terms in the left hand side of Eq. (1)). Therefore, we consider only the first term $F(\hat{f})$ in the following. The asymptotic behavior of the integral, Eq. (15), becomes

[TABLE]

In practice, we use a fitting function Eq. (18) instead of Eq. (15) in order to save the computational cost. Figure 1 shows the difference between Eq. (15) and the fitting function (18).

Substituting Eq. (12) into Eq. (7), we obtain the gravitational waveform in FD modified by the dipole radiation of a massive field, which is summarized as

[TABLE]

in Eq. (7) with Eq. (18). Since we have assumed that the modification is sufficiently small, i.e. $\delta\mathcal{A}\ll 1$, and neglected the non-linear corrections, the modified gravitational waveform (7) becomes invalid when $\delta\mathcal{A}$ becomes large. Moreover, since we only focus on the inspiral phase in this paper, the modified waveform is not valid if the mass of the field becomes too heavy for the modifications to occur mainly in the inspiral phase. It is worth mentioning that the modified waveform, Eq. (7) with Eqs. (19) and (20), can model not only a scalar field modification but also the vector field one.

In this parametrization, we have two free parameters, i.e., the mass of the field $m$ and the relative amplitude of dipole radiation $A$. In general, $A$ depends on the parameters of the system such as the masses and spins of compact objects, as well as on the modified gravity theory that one considers. Applying our modified template to the actual LVC events, we will obtain the likelihood in the two parameter space $(m,A)$ for each event. To combine the results of 10 binary BH events, one should assume an extended theory which allows hairy BHs and rewrite the relative amplitude $A$ by the model parameter common to all events, such as the coupling constant of the theory. Since EdGB theory has been studied in detail and known to allow hairy BHs Zwiebach:1985uq ; Moura:2006pz ; Pani:2009wy , in this paper we consider EdGB type coupling as an example. In EdGB theory, the energy flux of the dipole radiation is expressed as Yagi:2011xp ; Berti:2018cxi ; Prabhu:2018aun

[TABLE]

where $\zeta_{\rm EdGB}\equiv 16\pi\alpha_{\rm EdGB}^{2}/M^{4}$ with the coupling parameter $\alpha_{\rm EdGB}$ and

[TABLE]

is the spin-dependent factor of the BH scalar charges in the theory. Then, comparing Eqs. (11) and (21), we obtain the relation

[TABLE]

4 Analysis

4.1 Setup

We employ the waveform Eq. (7) with Eqs. (19) and (20) as templates to be matched with the strain data of Hanford and Livingston taken from the confident detections cataloged in GWTC-1 GWOSCweb . Using KAGRA Algorithmic Library (KAGALI) Oohara:2017aix , we evaluate the likelihood, following the standard procedure of the matched filtering maggiore2008gravitational ; 2009agwd.book…..J ; Creighton:2011zz (see also Appendix B). We adopt the published noise power spectrum for each event GWOSCweb . The minimum and the maximum frequencies, $f_{\rm min}$ and $f_{\rm max}$, of the datasets used in the analysis are summarized in Table 1. As the GR waveform $\tilde{h}_{\rm GR}$, we adopt IMRPhenomD Husa:2015iqa ; Khan:2015jqa , which is an up-to-date version of inspiral-merger-ringdown (IMR) phenomenological waveform for binary BHs with aligned spins. In this paper, since we are not interested in estimating the sky position and the inclination of the source, we choose the ratios of the $+/\times$-mode amplitudes of the respective detectors, i.e. these parameters are analytically optimized.

First, we implement a grid search to find the “best-fit parameters” of GR templates for each event varying the parameters $\mathcal{M}$, $\eta$, $\chi_{1}$, and $\chi_{2}$ as well as $t_{0}$ and $\phi_{0}$. The results in the detector frame are summarized in Table 1. Next, we calculate the likelihood for the modified templates around the GR best-fit. As is well known, there is an approximate degeneracy among the mass ratio and the spins in the inspiral phase waveform. In fact, we find that it is unnecessary to take into account the both variations in our calculations. Therefore, here we fix the spins to save the computational costs (see also Appendix C).

Moreover, to quantify the measurement accuracy of $\alpha_{\rm EdGB}$, we compute the Fisher information matrix defined by

[TABLE]

where $\hat{h}$ is an appropriately normalized template and $\theta^{i}$ are its parameters. For a large SNR, the root-mean-square error in $\theta^{i}$ can be evaluated as

[TABLE]

We use a polynomial fit given in Eq. (C.1) in Ref. Yunes:2016jcc as the noise power spectrum density of Advanced LIGO around GW151226 and GW170608. We take the lowest frequency to be $f_{\rm low}=$10 Hz. We use the restricted TaylorF2 PN aligned-spin model as the waveform model Buonanno:2009zt ; Arun:2008kb ; Mikoczi:2005dn ; Khan:2015jqa .

4.2 Results

Figure 2 shows the 90% confidence level (CL) constraints on $A$ for each event. To obtain the constraints, we calculate the likelihood, Eq. (32), of the best fit template for each point of $(A,m)$ plane. The likelihood can be regarded as the unnormalized posterior distribution when the prior of $A$ is uniform as well as those of $\mathcal{M}$ and $\eta$. Thus, we integrate it with respect to $A$ for each $m$ to evaluate 90% CL constraints. Here, we show only six events possessing relatively high SNR because the constraints on $A$ obtained from the other events are much larger and they are in the region where the error due to linearization becomes too large. Since our modified template becomes no longer valid if the mass of the field is too heavy, we terminate the calculations when the mass scale reaches a threshold frequency of each event, which is chosen to be ${\rm min}(0.5f_{\rm peak},f_{\rm max})$, where $f_{\rm peak}$ is the peak frequency of GWs. Because of the existence of the noise, the constraints seem to oscillate at the frequency corresponding to the region $m\gtrsim 5\times 10^{-14}$ eV, where the detectors are sensitive to the mass effects. For all events, larger values of $A$ remain unconstrained for heavier masses because the frequency that the dipole radiation turns on is too high for the detectors to be sensitive. This figure indicates that the constraints tend to be stronger for lighter chirp mass events. This is reasonable because events with a lighter chirp mass have longer inspiral phase at higher frequencies and the dipole radiation is basically -1PN correction which becomes more efficient in the case of a lighter chirp mass for fixed $A$ and $f$ as one can see in Eq. (10). Our constraint on $A$ for GW170817 is consistent with the estimated constraint in Ref. Alexander:2018qzg 111In their paper, the magnitude of modification due to the dipole radiation is governed by $\gamma$, which is related with $A$ as

$\gamma=\frac{48}{5\eta^{2/5}}A.$

For a binary NSs with SNR = 25 for the design sensitivity of LIGO, they estimated the constraint as $\gamma\lesssim 10^{-4}$ , which corresponds to $A\lesssim 10^{-6}$.. We should notice that in general a strong constraint on $A$ indicates that the modification to the waveform is strongly constrained, but it does not always imply a strong constraint on the modification to theory. This is because in general in scalar-tensor theory $A$ is proportional to the squared difference of the scalar charges of constituents of the binary as shown in Eq. (23). Therefore, if the difference vanishes, the modification to the waveform vanishes, i.e., $A=0$, even for a quite large coupling constant, and no one can constrain the modification to such a theory.

By combining the results, we can obtain more stringent and reliable constraints. Figure 3 shows the 68% CL, 90% CL, 95% CL constraints on $\sqrt{\alpha_{\rm EdGB}}$ after combining five binary BH events possesing relatively high SNR, i.e. GW150914, GW170814, GW170608, GW170104, GW151226, by assuming the EdGB type coupling (23). Similar to Fig. 2, assuming the prior distribution of $\sqrt{\alpha_{\rm EdGB}}$ is uniform, we derive the confidence intervals for $\sqrt{\alpha_{\rm EdGB}}$ for each $m$ from the integral of the likelihood. We plot the regions $10^{-15}$ eV $\lesssim m\lesssim$ $10^{-13}$ eV and $0.1$ km $\leq\sqrt{\alpha_{\rm EdGB}}\leq$ $5$ km. In the case of EdGB type of coupling the dipole modification vanishes if $m_{1}^{2}s_{2}^{\rm EdGB}=m_{2}^{2}s_{1}^{\rm EdGB}$. Therefore, in this case the leading signature of modification vanishes. In our calculations, we find no event which satisfies $m_{1}^{2}s_{2}^{\rm EdGB}=m_{2}^{2}s_{1}^{\rm EdGB}$. However, recalling that spins are almost zero-consistent, i.e., $s_{i}\simeq 1$, one can find it is very difficult to constrain the modification from nearly-equal mass binaries, such as GW170818. The regions where error due to nonlinearity, $\left(\delta\mathcal{A}\right)^{2}$, becomes large [$\left(\delta\mathcal{A}\right)^{2}\sim 0.01$] and the frequency at which the modification occurs exceeds $0.5f_{\rm peak}$ are hatched by black dashed lines. Here, we stacked only five events. This is because the region where our approximations are valid becomes smaller by including heavier chirp mass events, such as GW170729, which do not improve the constraints. Strictly speaking, the parameters in the hatched region are not excluded. However, the likelihood rapidly decreases before the hatched region and it seems quite unnatural that the likelihood increases in the hatched region to give the maximum even if we use templates which are valid there. Therefore, we believe that our constraints here are appropriate. We find that $\sqrt{\alpha_{\rm EdGB}}$ is constrained as $\sqrt{\alpha_{\rm EdGB}}\lesssim 2.47$ km for all mass below $6\times 10^{-14}$ eV for the first time with 90% CL including $\sqrt{\alpha_{\rm EdGB}}\lesssim 1.85$ km in the massless limit. This constraint in the massless limit is much stronger than the results of in Ref. Nair:2019iur ; Tahura:2019dgr , while it is accidentally almost the same with that by low-mass X-ray binaries Yagi:2012gp .

In order to investigate this disagreement with Ref. Nair:2019iur ; Tahura:2019dgr , we calculated the Fisher matrix with respect to the source parameters $\{\ln d_{L},t_{0},\phi_{0},{\cal M},\eta,\chi_{1},\chi_{2},\alpha_{\rm EdGB}^{2}\}$ where $d_{L}$ is the luminosity distance. In Table 2, we show the estimated constraints on $\alpha_{\rm EdGB}^{2}$, the measurement accuracy of $\eta$, and the correlation coefficient between $\eta$ and $\alpha_{\rm EdGB}^{2}$ in the cases of GW151226-like and GW170608-like signals, which are modeled by the inspiral PN waveform TaylorF2. The Fisher matrix is evaluated at the median value of the posterior distributions of GWTC-1 catalog LIGOScientific:2018mvr . The lower layer in Table 2 shows the case when all parameters are unconstrained a priori, while the upper layer takes into account the prior information that $\chi_{1}$ and $\chi_{2}$ must necessarily be smaller than unity, following Ref. Cutler:1994ys (see also Ref. Poisson:1995ef ). The presented values are for a single detector. When we use spin prior, constraints on $\sqrt{\alpha_{\rm EdGB}}$ are about 40% stronger than that when we do not, for both GW151226-like and GW170608-like signals.

The constraints without spin prior are consistent with those presented previously in Ref. Nair:2019iur ; Tahura:2019dgr , while the constraints with spin prior are in good agreement with our main analysis with strain data. GW170608-like signal gives stronger constraints on $\sqrt{\alpha_{\rm EdGB}}$ than GW151226-like signal. This is because the measurement accuracy of $\eta$ for GW151226-like signal is worse than that for GW170608-like signal, and the correlation between the symmetric mass-ratio and $\sqrt{\alpha_{\rm EdGB}}$ for GW151226-like signal is larger than that of GW170608-like signal. Finally, we should note that while the results of Fisher analysis imply that incorporating the spin prior affects the constraints by a factor, but the degeneracy between the spins and $\alpha_{\rm EdGB}$ is not so important as long as variations of the spins are limited in a relevant range as shown in Fig. 6.

5 Summary and discussion

In order to investigate the role of the noise in our calculations in the parameter space $(m,A)$, we prepare three mock datasets in which an artificial GR IMRPhenomD waveform is injected into the Gaussian noise with SNR about 10, 20, and 30, respectively. The parameters of the waveform are fixed to $\mathcal{M}=10.6\,M_{\odot}$ in the detector frame, $\eta=0.24$ (so that $m_{1}=15\,M_{\odot}$ and $m_{2}=10\,M_{\odot}$), and $\chi_{1}=\chi_{2}=0$. Figure 4 shows the 90% CL constraints on $A$ for each mock data. Here, assuming the uniform prior distribution of $A$, we integrate the likelihood with respect to $A$ for each $m$ as we did in Figs. 2 and 3. The regions where our approximation breaks down are hatched by black dashed lines. This indicates 90% CL constraints are approximately in inverse proportion to SNR for lighter $m$. Similarly to Fig. 2, the constraints seem to oscillate at the frequency corresponding to the region $m\gtrsim 5\times 10^{-14}$ eV because of the existence of the noise. The 90% CL boundaries for different SNR cases cross with each other in the range with heavier mass of the field, because we use a different noise realization in each case. If we use the same noise but with different magnitudes of signals, then the constraints are just scaled. As mentioned above, it seems quite unnatural that the likelihood increases in the hatched region to give the maximum even if we use templates which are valid there, although logically this possibility cannot be completely excluded. Therefore, 90% CL constraints would be still appropriate even if the constraints are close to the border of validity region (see also Appendix E). Since we use only one realization of injection for mock data, the result may depend on the noise realization. To investigate statistical feature by examining a number of realizations is left for a future work.

Figure 5 shows the matches [Eq. (31)] between the template modified by the effect of a massive field and the GR limit (left panel) and between the former and the massless limit for each value of $A$ (right panel), where we use the same values for $\mathcal{M}$, $\eta$, $\chi_{1}$ and $\chi_{2}$. The regions where our approximation is invalid are hatched by white dashed and black dashed lines, respectively. From the match with the GR limit, larger values of $A$ tend to be allowed for heavier masses even if there is no deviation from GR because the threshold frequency beyond which the dipole radiation is excited is too high and the detectors are less sensitive there. This is consistent with Fig. 2. On the other hand, the right panel shows that as expected it becomes difficult to distinguish a massive field from a massless field if the mass of the field gets smaller, while the match gets drastically smaller if the field is massive enough. Therefore, we cannot substitute the massless templates for the massive one if we try to constrain the magnitude of modification starting with the frequency band of the ground base detectors.

We discussed a dedicated test of the massive-field modification using the LIGO/Virgo’s open data of GWTC-1. In addition, assuming EdGB type coupling, we stacked the results of five events possessing relatively high SNR, and found the 90% CL constraints on the coupling parameter $\alpha_{\rm EdGB}$ as $\sqrt{\alpha_{\rm EdGB}}\lesssim 2.47$ km for any mass less than $6\times 10^{-14}$ eV for the first time including $\sqrt{\alpha_{\rm EdGB}}\lesssim 1.85$ km in the massless limit. In this analysis, we choose the ratios of the $+/\times$-mode amplitudes of the respective detectors, i.e. these parameters are analytically optimized. Moreover, since the number of events with a lighter chirp mass, such as binary NS, must increase in the near future, it is interesting to stack those events by assuming various theories in which charged NSs are allowed consistently. Furthermore, the NS-BH binary is also interesting Carson:2019fxr . Because, in addition to their light chirp mass, $A$ can be large in the theory that only either NS or BH has a charge, for example EdGB theory. Moreover, future multiband observations will improve the current upper bounds of various modified theories Gnocchi:2019jzp .

Acknowledgments

The authors are grateful to the developers of KAGRA Algorithmic Library (KAGALI), especially to Hideyuki Tagoshi, Nami Uchikata, and Hirotaka Yuzurihara for their support. The authors would like to thank Hiroyuki Nakano for the useful conversations and comments. K.Y. thanks to Nicolás Yunes and his group for their hospitality and nice conversations. K.Y. also wish to acknowledge the Yukawa Institute for Theoretical Physics at Kyoto University, where this work was developed during “The 3rd Workshop on Gravity and Cosmology by Young Researchers” (YITP-W-18-15). This work was supported by JSPS KAKENHI Grant Number JP17H06358 (and also JP17H06357), A01: Testing gravity theories using gravitational waves, as a part of the innovative research area, “Gravitational wave physics and astronomy: Genesis”. T.N.’s work was also supported in part by a Grant-in-Aid for JSPS Research Fellows. T. T. acknowledges support from JSPS KAKENHI Grant No. JP15H02087.

Appendix A Dipole Radiation

By separation of variables $\phi(t,\boldsymbol{r})=\varphi(t)\psi(\boldsymbol{r})$, the scalar field equation is decomposed to a set of equations as

[TABLE]

where $\omega^{2}=k^{2}+m^{2}$. Here we assume the validity of the PN expansion and hence the source $S(\boldsymbol{r})$ is a function of $\boldsymbol{r}$. The Green’s function of Eq. (25) is expanded by spherical harmonics as

[TABLE]

where $r=|\boldsymbol{r}|$, $r^{\prime}=|\boldsymbol{r}^{\prime}|$, we suppose $r>r^{\prime}$, and $j_{\ell}$ and $h^{(1)}_{\ell}$ are, respectively, the spherical Bessel function and the spherical Hankel function. Therefore, the solution of Eq. (25) can be expressed as

[TABLE]

Thus, the $k$-dependence of the energy loss by the dipole radiation of the scalar can schematically be written as

[TABLE]

where $T^{\phi}_{\mu\nu}$ is the energy-momentum tensor of the scalar field, $d\Sigma$ is the surface element, and we assume $kr^{\prime}\ll 1$ and expand the spherical Bessel function. The dot and dash denote the derivatives with respect to time and radial coordinates, respectively. The order of magnitude of dipole radiation is larger than that of quadruple because the energy flux of the dipole and quadrupole radiation can be estimated as

[TABLE]

where $D$ and $Q$ are the dipole and quadruple moments, respectively, $\langle\cdots\rangle$ denotes a temporal average over several periods of GWs, and we have used $\omega=\omega_{\rm orbit}\propto r^{-3/2}$ with the orbital angular velocity $\omega_{\rm orbit}$. Combining Eqs. (28)–(29), the energy loss by dipole radiation can be estimated as

[TABLE]

where $A$ is an overall factor and we introduce the step function to represent the sudden activation of dipole radiation.

Appendix B Matched Filtering

Define the scalar product between two real functions $h(t)$ and $g(t)$ by

[TABLE]

where $\operatorname{Re}$ denotes the real part, $\tilde{}$ indicates the Fourier transformation of the TD functions, and $S_{n}(f)$ is the noise power spectrum density as

[TABLE]

The Gaussian probability distribution for the noise $\tilde{n}$ is

[TABLE]

We are assuming that the output of the detector satisfies the condition for claiming detection, i.e. it is of the form $s(t)=\rho h(t;\theta)+n_{0}(t)$, where $n_{0}$ is a specific realization of the noise, and $h(t;\theta)$ is the signal with the parameter $\theta$ and the SNR of the event $\rho=(s|h)$ is sufficiently high. The likelihood function for the observed output $s(t)$, given that there is a GW signal corresponding to the parameters $\theta$, is obtained plugging $n_{0}=s-\rho h(\theta)$ into the above equation,

[TABLE]

Appendix C Approximate degeneracy among the mass ratio and the spins

To investigate the effect of the spins on the constraints, we compare the 90% CL constraints on $\sqrt{\alpha_{\rm EdGB}}$ in two cases, after combining five binary BH events possesing relatively high SNR. Case 1: the same with the 90% CL constraint in Fig. 3. Case 2: we additionary vary the spins and estimate the 90% CL constraint. In Case 2, we reduce the number of sample points to save the computational costs. Figure 6 shows the comparison between them. It is found that the effect of the spins on the constraints is small enough.

Appendix D Effects on SNR by increasing the number of parameters

How much SNR can become larger if the number of parameters increases? To see this, let us consider the number of degrees of freedom in the parameter search. Let $(x_{i}|x_{j})$ denote the match between the datasets $x_{i}$ and $x_{j}$. Substituting $x_{i}=s-\rho_{i}\hat{h}_{i}=\rho\hat{h}-\rho_{i}\hat{h}_{i}+n$, we have a relation between the match and the chi-squared as

[TABLE]

where we assume $\hat{h}_{i}$ is chosen to be the best fit within a given template bank. Therefore, the difference of $\chi^{2}$ between different template banks is

[TABLE]

Let $i=1,2$ represent the GR and modified theory, respectively. In general, increasing the number of parameters causes decreasing the number of degrees of freedom of $\chi^{2}$-distribution, so that the value $\chi^{2}$ decreases. Therefore, the mean value of the left-hand side of the above equation becomes 2. On the other hand, the right-hand side of the above equation is the deference of SNR. For example, we obtain $|\rho_{2}|^{2}-|\rho_{1}|^{2}=2.13934$ by averaging all eleven GW events.

Appendix E Injection test of detectability of massive field modification

In order to demonstrate the detectability of the modification caused by a massive field, we prepare nine mock datasets for which artificial modified waveforms are injected into the Gaussian noise with $m=5\times 10^{-14}$ eV, $A=10^{-2},10^{-3}$, and $10^{-4}$ and SNR about 10, 20, and 30. The other parameters of the waveform are fixed to $\mathcal{M}=10.6\,M_{\odot}$ in the detector frame, $\eta=0.24$ (so that $m_{1}=15\,M_{\odot}$ and $m_{2}=10\,M_{\odot}$), and $\chi_{1}=\chi_{2}=0$. Figure 7 shows the likelihood normalized by the GR values in the $(A,m)$ plane for each mock data. The region hatched by white dashed lines indicates the outside of the range of validity. The red circle shows the injected values of $A$ and $m$. In general, as we discussed in Appendix D, increasing the number of parameters of templates leads to increasing SNR. In our case, since the number of the additional parameters is two, the likelihood normalized by the GR values can reach to $e\simeq 2.7$ regardless of whether the signal is GR or not. Therefore, the likelihood needs to be much larger than this value for a confident detection of the modification. From this point of view, Fig. 7 implies that for the case of $\mathcal{M}=10.6\,M_{\odot}$ we may detect the modification due to the massive field if $A\gtrsim 10^{-3}$ and SNR is larger than 20. For the other cases it is possible only to give an upper bound to $A$. However,in the case that $A=10^{-2}$ with SNR is 10 (Right-Top panel in Fig. 7), the likelihood remains large even at the boundary of the hatched region and hence a meaningful constraint will not be obtained. Since we use only one realization of injection for mock data, the result may depend on the noise realization. To investigate statistical feature by a number of realizations is left for a future work.