Testing massive-field modifications of gravity via gravitational waves

Abstract

1. Introduction

The opening of gravitational wave (GW) astronomy/astrophysics allows us to approach a new test of general relativity (GR) in a 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 [1]. Testing GR has been investigated by several authors using these event data and no significant deviation from GR has been reported [2–5]. Although one of the most interesting regimes to investigate is the ringdown phase [6,7], concrete templates for merger ringdown phase waveforms in modified gravity theories are not yet established [8]. 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 [9–13], a waveform in this framework does not cover all viable extensions of gravity. One such extension 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 [14,15].

Spontaneous scalarization of compact objects for massless fields was first investigated by Damour and Esposito-Fareèse [16]. 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 [17]. 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 has been discussed by several authors (e.g. see Refs. [18,19]). 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 [19], (2) $m\lesssim {10}^{-16}$$$eV from the observations of a pulsar–white dwarf binary [17], and (3) ${10}^{-13}$$$eV $\lesssim m\lesssim {10}^{-11}$$$eV [20,21], which relies on measurements of high spins of stellar-mass BHs [22].

In the presence of matter, dynamical scalarization in scalar–tensor theories has been investigated [23–26]. 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 [27,28]. 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 scalar hair [3]. Very recently, Nair et al. discussed constraining EdGB theory from GW observations [29]. Their order-of-magnitude estimation on the constraint was $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 1.0$$$km, where ${\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$ is the coupling parameter of EdGB theory, while their Fisher-estimated and Bayesian constraints were roughly $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 6.0$$$km. In addition, Tahura et al. independently derived constraints on a few modified theories of gravity, including EdGB theory, and their constraint on EdGB theory was roughly $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 4.3$$$km [30]. The current constraint on EdGB theory is $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 1.9$$$km, which was obtained using low-mass X-ray binaries [31].

In this paper we discuss activation of dipole radiation due to massive fields and modifications to the gravitational waveform. We analyze the GW events in the GWTC-1 [32] catalog 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 LIGO, Virgo, and KAGRA ground-based detectors are sensitive at frequencies around 10–1000$$Hz, gravitational wave signals allow us to investigate the effect of 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 }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$ 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 modified gravitational wave templates and Fisher analysis on the determination error of the coupling parameter. Section 5 is devoted to a summary and discussion. Throughout this paper we use geometric units in which $G=1=c$⁠.

2. The ppE waveform

3. Massive field

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

Fig. 1.

Top panel: The exact integral of Eq. (15) (red solid line) and the fitting function in Eq. (18) (cyan dashed line). Bottom: The difference of the fitting function from Eq. (15).

Open in new tabDownload slide

4. Analysis

4.1. Setup

We employ the waveform in 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 [32]. Using the KAGRA Algorithmic Library (KAGALI) [41] we evaluate the likelihood, following the standard procedure of matched filtering [42–44] (see also Appendix B). We adopt the published noise power spectrum for each event [32]. The minimum and the maximum frequencies, $f_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and $f_{\mathrm{m}\mathrm{a}\mathrm{x}}$⁠, of the datasets used in the analysis are summarized in Table 1. As the GR waveform ${\tilde{h}}_{\mathrm{G}\mathrm{R}}$ we adopt IMRPhenomD [45,46], which is an up-to-date version of the 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 so that the angular parameters of the waveform template are set to the best-fit values in an analytical way.

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 ${\varphi }_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 both variations in our calculations. Therefore, we fix the spins to reduce the computational costs (see also Appendix C).

We use a polynomial fit given in Eq. (C.1) of Ref. [3] as the noise power spectrum density of Advanced LIGO around GW151226 and GW170608. We take the lowest frequency to be $f_{\mathrm{l}\mathrm{o}\mathrm{w}}=10$$$Hz. We use the restricted TaylorF2 PN aligned-spin model as the waveform model [46–49].

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. (B.2), of the best-fit template for each point of the $(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 the 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 for each event, which is chosen to be $\mathrm{m}\mathrm{i}\mathrm{n}(0.5f_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}},f_{\mathrm{m}\mathrm{a}\mathrm{x}})$⁠, where $f_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ is the peak frequency of GWs. Because of the existence of the noise, the constraints seem to oscillate at a 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 a $-1$PN 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. [28].¹ 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 quite a large coupling constant, and no one can constrain the modification to such a theory.

Fig. 2.

The 90% CL constraints on $A$ for six events possessing relatively high SNR. The calculations are terminated when the mass scale reaches a threshold frequency for each event, which is chosen to be $\mathrm{m}\mathrm{i}\mathrm{n}(0.5f_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}},f_{\mathrm{m}\mathrm{a}\mathrm{x}})$⁠.

Open in new tabDownload slide

By combining the results, we can obtain more stringent and reliable constraints. Figure 3 shows the 68% CL, 90% CL, and 95% CL constraints on $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ after combining five binary BH events possessing relatively high SNR (GW150914, GW170814, GW170608, GW170104, and GW151226), assuming EdGB type coupling as in Eq. (23). Similar to Fig. 2, assuming the prior distribution of $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ is uniform, we derive the confidence intervals for $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ 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 $\le \sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\le 5$$$km. In the case of EdGB type coupling the dipole modification vanishes if $m_1^2{s}_2^{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}=m_2^2{s}_1^{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$⁠. Therefore, in this case the leading signature of modification vanishes. In our calculations, we find no event which satisfies $m_1^2{s}_2^{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}=m_2^2{s}_1^{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$⁠. However, recalling that spins are almost zero-consistent, i.e. $s_i\simeq 1$⁠, one can find it very difficult to constrain the modification from nearly equal-mass binaries, such as GW170818. The regions where the error due to nonlinearity, $(\delta \mathcal{A}{)}^2$⁠, becomes large [$(\delta \mathcal{A}{)}^2\sim 0.01$] and the frequency at which the modification occurs exceeds $0.5f_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ 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 }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ is constrained as $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 2.47$$$km for all mass below $6\times {10}^{-14}$$$eV for the first time with 90% CL including $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 1.85$$$km in the massless limit. This constraint in the massless limit is much stronger than the results in Refs. [29,30], while it is accidentally almost the same as that for low-mass X-ray binaries [31].

Fig. 3.

The 68% CL, 90% CL, and 95% CL constraints on $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ after combining five binary BH events possessing relatively high SNR (GW150914, GW170814, GW170608, GW170104, and GW151226), assuming EdGB type coupling as in Eq. (23).

Open in new tabDownload slide

In order to investigate this disagreement with Refs. [29,30], we calculated the Fisher matrix with respect to the source parameters $\{\ln d_L,t_0,{\varphi }_0,\mathcal{M},\eta ,{\chi }_1,{\chi }_2,{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}^2\}$⁠, where $d_L$ is the luminosity distance. Table 2 shows the estimated constraints on ${\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}^2$⁠, the measurement accuracy of $\eta$⁠, and the correlation coefficient between $\eta$ and ${\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}^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 the GWTC-1 catalog [1]. The lower part of Table 2 shows the case when all parameters are unconstrained a priori, while the upper part takes into account the prior information that ${\chi }_1$ and ${\chi }_2$ must necessarily be smaller than unity, following Ref. [50] (see also Ref. [51]). The values presented are for a single detector. When we use the spin prior, the constraints on $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ are about 40% stronger than when we do not, for both GW151226-like and GW170608-like signals.

The constraints without the spin prior are consistent with those presented previously in Refs. [29,30], while the constraints with the spin prior are in good agreement with our main analysis with strain data. GW170608-like signals give stronger constraints on $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ than GW151226-like signals. This is because the measurement accuracy of $\eta$ for GW151226-like signals is worse than that for GW170608-like signals, and the correlation between the symmetric mass ratio and $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ for GW151226-like signals is larger than for GW170608-like signals. Finally, we should note that while the results of the Fisher analysis imply that incorporating the spin prior affects the constraints by a factor, the degeneracy between the spins and ${\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$ is not so important as long as variations of the spins are limited in a relevant range, as shown in Fig. C.1.

Fig. C.1.

Comparison of the 90% CL constraints on $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}$ between Case 1 and Case 2, in which we implement a grid search with/without varying the spins, respectively.

Open in new tabDownload slide

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 of about 10, 20, and 30, respectively. The parameters of the waveform are fixed to $\mathcal{M}=10.6M_{\odot }$ in the detector frame, $\eta =0.24$ (so that $m_1=15M_{\odot }$ and $m_2=10M_{\odot }$⁠), and ${\chi }_1={\chi }_2=0$⁠. Figure 4 shows the 90% CL constraints on $A$ for each mock dataset. Here, assuming a 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 that the 90% CL constraints are approximately in inverse proportion to the SNR for lighter $m$⁠. Similarly to Fig. 2, the constraints seem to oscillate at a 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 signal magnitudes, 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 still be appropriate even if the constraints are close to the border of the 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. Investigating the statistical features by examining a number of realizations is left for future work.

Fig. 4.

The 90% CL constraints on $A$ for each mock dataset. This implies that the 90% CL constraints are approximately in inverse proportion to the SNR.

Open in new tabDownload slide

Figure 5 shows the matches, Eq. (B.1), 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 ground-based detectors.

Fig. 5.

Left panel: Match between the modified template due to a massive field and the GR limit. Right panel: The modified template due to a massive field and the massless limit.

Open in new tabDownload slide

We discussed a dedicated test of the massive-field modification using the LIGO/Virgo open data from 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 }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}$ as $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 2.47$$$km for any mass less than $6\times {10}^{-14}$$$eV, for the first time including $\sqrt{{\alpha }_{\mathrm{E}\mathrm{d}\mathrm{G}\mathrm{B}}}\lesssim 1.85$$$km in the massless limit. In this analysis, we choose the ratios of the $+/\times$-mode amplitudes of the respective detectors so that the sky position and the inclination of the source are set to the best-fit values in an analytical way. Moreover, since the number of events with a lighter chirp mass, such as a 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, NS–BH binaries are also interesting [52] because, in addition to their light chirp mass, $A$ can be large in the theory that only either the NS or the BH has a charge, for example EdGB theory. Moreover, future multiband observations will improve the current upper bounds of various modified theories [53].

Acknowledgements

The authors are grateful to the developers of the 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 useful conversations and comments. K.Y. thanks Nicolás Yunes and his group for their hospitality and nice conversations. K.Y. also wishes 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 Japan Society for the Promotion of Science (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

Appendix B. Matched filtering