Testing massive-field modifications of gravity via gravitational waves

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.


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 [1]. Testing GR has been investigated by several authors using these event data and no significant deviation from GR has been reported [2][3][4][5]. Although one of the most interesting regime to investigate is the ringdown phase [6,7], the concrete templates of merger-ringdown phase waveform in modified gravity theories are not established yet [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][10][11][12][13], 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 [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 is 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 10 −27 eV for stability on cosmological scales [19], (2) m 10 −16 eV from the observations of a pulsarwhite dwarf binary [17], and (3) 10 −13 eV m 10 −11 eV [20,21], which relies on the measurements of high spins of stellar mass BHs [22].
In the presence of matter, dynamical scalarization in scalar-tensor theories has been investigated [23][24][25][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 Einsteindilaton-Gauss-Bonnet (EdGB) theory with a massive scalar field, in which a BH can have a scalar hair [3]. Very recently, Nair et al. discussed to constrain EdGB theory from the GW observations [29]. Their order-of-magnitude estimation on the constraint is √ α EdGB 1.0 km, where α EdGB is the coupling parameter of EdGB theory, while their Fisher-estimated and Bayesian constraints are roughly √ α EdGB 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 √ α EdGB 4.3 km [30]. Current constraint on EdGB theory is √ α EdGB 1.9 km, which is obtained by using low-mass X-ray binaries [31].
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 [32] 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 m 10 −13 eV, while they reduce to the massless limit for m 10 −15 eV. We also evaluate the measurement accuracy of α 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.

The ppE waveform
In this section, we briefly review the ppE waveform following Ref. [9,33]. 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ḟ . The frequency domain (FD) ppE waveform for the inspiral phase of compact binaries is expressed ash whereh GR is the waveform in GR and its phase is given by with the coalescence time t 0 and phase φ 0 . αu a and δΨ represent modifications to the amplitude and the phase, respectively. We introduce where M = M η 3/5 is the chirp mass with the total mass of the binary M = m 1 + m 2 and the symmetric mass ratio η = m 1 m 2 /M 2 . We can also model the phase correction as α, β, 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 [3,34]. 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ḟ . Then, let us denotė f asḟ where the GR termḟ GR (f ) is expressed aṡ Therefore, we obtain 3/17 at the leading order in δḟ . In the stationary phase approximation the gravitational waveform in the FD,h(f ), can be expressed as where A (f ) is the amplitude of time domain (TD) waveform evaluated att =t(f ). Substituting Eqs. (4) and (5) into Eq. (6), we find at the leading order, where

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) where ω = πf is the angular frequency of the dipole radiation, for which the dispersion relation is ω 2 = k 2 + m 2 with the wave number k and the mass of the field m, Θ 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 Since the correction to the binding energy of the system is higher order in the PN expansion, we obtain In the massless limit, this becomes Substituting Eq. (12) into Eq. (9), we obtain the modification to the GW phase as 4/17 where we define withf = πf /m. In our matched filtering analysis, G(f ) is irrelevant because 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 (f ) in the following. The asymptotic behavior of the integral, Eq. (15), becomes 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 in Eq. (7) with Eq. (18). Since we have assumed that the modification is sufficiently small, i.e. δA 1, and neglected the non-linear corrections, the modified gravitational waveform (7) becomes invalid when δ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 [35][36][37], in this paper we consider EdGB type coupling as an example. In EdGB theory, the energy flux of the dipole radiation is expressed as [38][39][40] where ζ EdGB ≡ 16πα 2 EdGB /M 4 with the coupling parameter α EdGB and is the spin-dependent factor of the BH scalar charges in the theory. Then, comparing Eqs. (11) and (21), we obtain the relation 4. Analysis

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 [32]. Using KAGRA Algorithmic Library (KAGALI) [41], we evaluate the likelihood, following the standard procedure of the matched filtering [42][43][44] (see also Appendix B). We adopt the published noise power spectrum for each event [32]. The minimum and the maximum frequencies, f min and f max , of the datasets used in the analysis are summarized 6/17 in Table 1. As the GR waveformh GR , we adopt IMRPhenomD [45,46], which is an up-todate 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 +/×-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 M, η, χ 1 , and χ 2 as well as t 0 and φ 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 α EdGB , we compute the Fisher information matrix defined by whereĥ is an appropriately normalized template and θ i are its parameters. For a large SNR, the root-mean-square error in θ i can be evaluated as We use a polynomial fit given in Eq. (C.1) in Ref. [3] as the noise power spectrum density of Advanced LIGO around GW151226 and GW170608. We take the lowest frequency to be f low =10 Hz. We use the restricted TaylorF2 PN aligned-spin model as the waveform model [46][47][48][49].   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 of each event, which is chosen to be min(0.5f peak , f max ) 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 M and η. 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 min(0.5f peak , f max ), where f 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 5 × 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. [28] 1 . 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, 1 In their paper, the magnitude of modification due to the dipole radiation is governed by γ, which is related with A as γ = 48 5η 2/5 A. For a binary NSs with SNR = 25 for the design sensitivity of LIGO, they estimated the constraint as γ 10 −4 , which corresponds to A 10 −6 . 8/17 even for a quite large coupling constant, and no one can constrain the modification to such a theory.

Results
By combining the results, we can obtain more stringent and reliable constraints. Figure 3 shows the 68% CL, 90% CL, 95% CL constraints on √ α 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 √ α EdGB is uniform, we derive the confidence intervals for √ α EdGB for each . However, recalling that spins are almost zero-consistent, i.e., s i 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, (δA) 2 , becomes large [(δA) 2 ∼ 0.01] and the frequency at which the modification occurs exceeds 0.5f 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 √ α EdGB is constrained as √ α EdGB 2.47 km for all mass below 6 × 10 −14 eV for the first time with 90% CL including √ α EdGB 1.85 km in the massless limit. This constraint in the massless limit is much stronger than the results of in Ref. [29,30], while it is accidentally almost the same with that by low-mass X-ray binaries [31].  In order to investigate this disagreement with Ref. [29,30], we calculated the Fisher matrix with respect to the source parameters {ln d L , t 0 , φ 0 , M, η, χ 1 , χ 2 , α 2 EdGB } where d L is the luminosity distance. In Table 2, we show the estimated constraints on α 2 EdGB , the measurement accuracy of η, and the correlation coefficient between η and α 2 EdGB 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 [1]. 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 χ 1 and χ 2 must necessarily be smaller than unity, following Ref. [50] (see also Ref. [51]). The presented values are for a single detector. When we use spin prior, constraints on √ α 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. [29,30], while the constraints with spin prior are in good agreement with our main analysis with strain data. GW170608-like signal gives stronger constraints on √ α EdGB than GW151226-like signal. This is because the measurement accuracy of η for GW151226-like signal is worse than that for GW170608-like signal, and the correlation between the symmetric mass-ratio and √ α 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 α EdGB is not so important as long as variations of the spins are limited in a relevant range as shown in Fig. C1. 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 and χ 2 must necessarily be smaller than unity. The measurement accuracy of symmetric mass-ratio and correlation coefficient between η and √ α EdGB are also shown.

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 M = 10.6 M in the detector frame, η = 0.24 (so that m 1 = 15 M and m 2 = 10 M ), and χ 1 = χ 2 = 0. Figure 4 shows the 90% CL constraints on A for each mock 10/17  Fig. 4 The 90% CL constraints on A for each mock data. This implies 90% CL constraints are approximately in inverse proportion to SNR. 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.
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 5 × 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 11/17 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. (B1)] 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 M, η, χ 1 and χ 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 α EdGB as √ α EdGB 2.47 km for any mass less than 6 × 10 −14 eV for the first time including √ α EdGB 1.85 km in the massless limit. In this analysis, we choose the ratios of the +/×-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 [52]. 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 [53].

A. Dipole Radiation
By separation of variables φ(t, r) = ϕ(t)ψ(r), the scalar field equation is decomposed to a set of equations as d dt 2 + ω 2 ϕ(t) = 0, (A1) where ω 2 = k 2 + m 2 . Here we assume the validity of the PN expansion and hence the source S(r) is a function of r. The Green's function of Eq. (A2) is expanded by spherical harmonics as where r = |r|, r = |r |, we suppose r > r , and j and h (1) are, respectively, the spherical Bessel function and the spherical Hankel function. Therefore, the solution of Eq. (A2) can be expressed as Thus, the k-dependence of the energy loss by the dipole radiation of the scalar can schematically be written as where T φ µν is the energy-momentum tensor of the scalar field, dΣ is the surface element, and we assume kr 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 where D and Q are the dipole and quadruple moments, respectively, · · · denotes a temporal average over several periods of GWs, and we have used ω = ω orbit ∝ r −3/2 with the orbital angular velocity ω orbit . Combining Eqs. (A5)-(A6), the energy loss by dipole radiation can be estimated as where A is an overall factor and we introduce the step function to represent the sudden activation of dipole radiation. 13/17

B. Matched Filtering
Define the scalar product between two real functions h(t) and g(t) by where Re denotes the real part,˜indicates the Fourier transformation of the TD functions, and S n (f ) is the noise power spectrum density as The Gaussian probability distribution for the noiseñ is We are assuming that the output of the detector satisfies the condition for claiming detection, i.e. it is of the form s(t) = ρh(t; θ) + n 0 (t), where n 0 is a specific realization of the noise, and h(t; θ) is the signal with the parameter θ and the SNR of the event ρ = (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 θ, is obtained plugging n 0 = s − ρh(θ) into the above equation, To investigate the effect of the spins on the constraints, we compare the 90% CL constraints on √ α 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 C1 shows the comparison between them. It is found that the effect of the spins on the constraints is small enough. 14/17

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 − ρ iĥi = ρĥ − ρ iĥi + n, we have a relation between the match and the chi-squared as where we assumeĥ i is chosen to be the best fit within a given template bank. Therefore, the difference of χ 2 between different template banks is 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 χ 2 -distribution, so that the value χ 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 |ρ 2 | 2 − |ρ 1 | 2 = 2.13934 by averaging all eleven GW events.

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 × 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 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. E1 implies that for the case of M = 10.6 M we may detect the modification due to the massive field if A 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. E1), 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.