1. Introduction
Since the first direct detection of gravitational wave (GW) event GW150914 [1, 2] by the Laser Interferometer Gravitational-Wave Observatory (LIGO) Scientific Collaboration [3, 4] and Virgo Collaboration [5], so far there have been more than 90 confirmed GW detections [6–10]. Most of the GW signals were from the coalescence of binary black holes, except for two from binary neutron stars [11, 12] and two events, GW200115 and GW200105, from neutron star-black hole mergers [13]. These GWs provide useful information about the gravitational interaction in the strong-field and nonlinear regime [14–18] and are possible to be used to test General Relativity (GR) and understand the nature of gravity. Although GR has passed solar system and binary pulsar tests perfectly, the strong-field regime tests are still lacking [19, 20]. For example, observations have not yet been able to confirm GR's no-hair theorems or the existence of extra fields additional to gravitational fields, such as scalars, vectors or tensors in modified gravity theories. It would be very useful to assess how accurately such GW observations could actually constrain alternative theories of gravity.
One of the simple alternative theories of gravity is the Brans–Dicke (BD) theory [21–23]. In the BD theory, a scalar field φ is postulated in addition to the space time metric gμν, with an effective strength inversely proportional to the coupling parameter ωBD. In generalized scalar-tensor theories, which arise as byproducts of string theory and other unification schemes, ωBD is treated as a function of φ. The scalar field φ not only takes the role of 1/G but also mediates gravity and excites the scalar breathing mode in GWs. As ωBD → ∞ , the effect of the scalar field tends to zero and BD theory reduces to GR. The most stringent constraint ωBD > 40 000 [24] was obtained in the Cassini measurements of the Shapiro time delay in the solar system [25]. For a binary system composed of compact objects, the orbital period of the system will decrease due to the loss of energy by the emission of GWs. In the BD theory, the extra emission channel of GWs can further decrease the orbital period of a binary system [23, 26], which provides another way to constrain the BD theory [23, 26–33] via measuring the change of the orbital period of a binary system. The measurement of the orbital decay from the pulsar-white dwarf binary PSR J1738 + 0333 gives the constraint ωBD > 25 000 [34]. The above laboratory experiments and astrophysical observations set strong constraints on the parameter ωBD in the weak-field region, the validity for the BD theory in the strong-field and the nonlinear regime remains tested. Theoretically, using the simple approximate waveform template with the dipolar correction in the phase and the Fisher matrix method, it was estimated that the observation of a 0.7M⊙ neutron star (NS) on a quasicircular inspiral into a 3M⊙ black hole (BH) with a signal-to-noise (SNR) of 10 by LIGO/Virgo detectors can give the constraint ωBD ≳ 2000 [35]. In the extreme mass ratio inspiral (EMRI) system, the constraint is expected to be enhanced to ωBD > 105 [36]. However, the Fisher matrix method has limitations and the estimates are usually too optimistic because the Fisher matrix method models the likelihood as a covariant Gaussian [37, 38]. Furthermore, the Fisher matrix method does not provide actual posterior distributions for physical parameters [39]. In reality, the constraint from observed GW events for the the BD parameter ωBD in the strong-field region is much worse than the expectation. For example, with the Bayesian inference method, the GW event GW190426_152155 of a possible 1.5M⊙ NS/5.7M⊙ BH binary merge only yields the constraint ωBD ≳ 10 [40].
Previous works constructed modified gravitational waves in the BD theory considering only the dipolar emission from the circular orbits [40]. However, the orbit in reality usually has eccentricity for many reasons such as the binary emerges, captures, or their initial stages [41]. In addition to the dipolar emission, monopolar radiation is also present if the eccentricity is nonzero [27, 42], such additional emission channels in the BD theory for eccentric binaries can be helpful to distinguish GR from the BD theory. In a usual binary system, orbits circularize before the coalescence, so that the eccentricity is negligible and undetectable by the time they enter the band of ground-based GW detectors at 10 Hz [43]. Dynamic formation of binaries through strong GW radiation could make the binary eccentricity large at coalescence without enough orbits circularizations [44–53]. The detection of eccentric binaries and the measurement of eccentricity can be useful to infer the binary formation mechanism [39, 51, 53–60], moreover it can be helpful to constrain alternative theories of gravity because the eccentricity can induce additional radiation. In binary BHs, the orbital evolution and gravitational radiation from binary BHs in the BD theory are identical to those in GR, after a mass scaling, so we cannot distinguish GR from the BD theory in the BH-BH binary system. [27, 61, 62]. Henceforth, we examine the NS-BH binary system, carry out Bayesian parameter estimation for NS-BH events GW200105 and GW200115 by taking the eccentricity into account. Although there exists waveform approximations that describe the full inspiral, merger and ringdown of eccentric, nonspinning binary BHs in GR [63, 64], the waveform approximation that includes spin and higher order modes in GR or the BD theory is not yet available. In the discussion of measuring eccentricity in binaries with GWs, the inspiral-only waveform approximations such as EccentricFD and TalylorF2ECC for nonspinning eccentric binaries were usually used [39, 57, 58, 65, 66]. For the spinning eccentric binaries, the first eccentric binary waveform model SEOBNRE (Spinning Effective-One-Body-Numerical-Relativity model for Eccentric binary) based on effective-one-body-numerical-relativity the framework was developed in [67]. Based on SEOBNRE, the new waveform model SEOBNREHM including higher-multipole (l, ∣m∣) = (2, 2), (2, 1), (3, 3), (4, 4) was constructed in [68]. Furthermore, the improved waveform model to describe eccentric spin-precessing binary black hole coalescence was recently developed in [69].
Since we are more interested in the effect of the eccentricity on the constraint of the BD coupling parameter ωBD, we ignore the effect of eccentricity on the amplitude of GWs and implement the eccentricity by simply adding its effect to the phase of a waveform template with no eccentricity in GR by using the post-Newtonian results in the BD theory. With the recognition of the limitation of this procedure, we derive an analytical expression for the phase modification and use the modified waveform template to perform a Bayesian parameter estimation.
This paper is organized as follows. In section 2 , we derive the phase deviation caused by monopole and dipole gravitational radiation for eccentric compact binaries. In section 3 , we describe our data analysis set-up and apply our methodology to the selected events GW200105 and GW200115 with Bayesian inference. The conclusion is provided in section 4 . Throughout this paper, we use the units G = c = 1.
2. Binary gravitational waveforms in the Brans–Dicke theory
In the BD theory [21], the locally measured gravitational constant G varies with the scalar field φ and the gravitation is also mediated by the scalar field in addition to tensor fields. As the gravitational constant varies, the binding energy of a self-gravitating body changes and its inertial mass depends on the scalar gravitational field φ as m(φ) [23]. The variation in the inertial mass caused by the varying of G results in different accelerations for massive bodies with different structures and the equivalence principle is violated. The action for the BD theory with such massive bodies reads [21, 23]5 ) is9 ), the two terms in the square brackets represent the combined effects of quadrupole and monopole gravitational radiation, and the last term is the effect of dipole radiation. Carrying out the average over one orbit using Keplerian orbital formulae13 )17 ) can be simplified as19 ) with time, we get the time evolution of gravitational radiations19 ), (25 ), (26 ) and (27 ), we derive the phase of GWs to the leading order as
$\begin{eqnarray}\begin{array}{l}I=\displaystyle \frac{1}{16\pi }\displaystyle \int {{\rm{d}}}^{4}x\sqrt{-g}\left(\phi R-\displaystyle \frac{{\omega }_{{BD}}}{\phi }{\phi }_{,\alpha }{\phi }^{,\alpha }\right)\\ \quad -\displaystyle \sum _{a}\displaystyle \int {\rm{d}}{\tau }_{a}{m}_{a}(\phi ),\end{array}\end{eqnarray}$
where R is the Ricci scalar, g is the determinant of the metric gμν, and τa is the proper time along the trajectory of body a. The resulting field equations are $\begin{eqnarray}\begin{array}{l}{R}_{\mu \nu }-\displaystyle \frac{1}{2}{g}_{\mu \nu }R=\displaystyle \frac{8\pi }{\phi }{T}_{\mu \nu }\\ \quad +\displaystyle \frac{{\omega }_{{BD}}}{{\phi }^{2}}\left({\phi }_{,\mu }{\phi }_{,\nu }-\displaystyle \frac{1}{2}{g}_{\mu \nu }{\phi }_{,\lambda }{\phi }^{,\lambda }\right)\\ \quad +\displaystyle \frac{1}{\phi }\left({\phi }_{;\mu \nu }-{g}_{\mu \nu }{\square }_{g}\phi \right),\end{array}\end{eqnarray}$
$\begin{eqnarray}{\square }_{g}\phi =\displaystyle \frac{8\pi }{3+2{\omega }_{{BD}}}\left(T-2\phi \displaystyle \frac{\partial T}{\partial \phi }\right),\end{eqnarray}$
where ${\square }_{g}={\partial }_{\mu }(\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu })/\sqrt{-g}$, T = Tμνgμν, and $\begin{eqnarray}{T}^{\mu \nu }=\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \sum _{a}\displaystyle \frac{{m}_{a}(\phi ){u}^{\mu }{u}^{\nu }}{{u}^{0}}{\delta }^{3}({\boldsymbol{x}}-{{\boldsymbol{x}}}_{{\boldsymbol{a}}}).\end{eqnarray}$
The energy conservation takes the form $\begin{eqnarray}{T}_{;\nu }^{\mu \nu }-\displaystyle \frac{\partial T}{\partial \phi }{\phi }^{,\mu }=0.\end{eqnarray}$
If the equivalence principle holds, the inertial mass ma of body a is independent of φ and its sensitivity ${s}_{a}=[\mathrm{dln}{m}_{a}(\phi )/\mathrm{dln}\phi ]{| }_{{\phi }_{0}}=0$, then ∂T/∂φ = 0 and the conservation of the energy-momentum tensor ${T}_{;\nu }^{\mu \nu }=0$ is recovered. We expand φ about its asymptotic, cosmological value φ0, φ = φ0 + φ, and write the variations in ma with φ as $\begin{eqnarray}\begin{array}{l}{m}_{a}(\phi )={m}_{a}\left[1+{s}_{a}\displaystyle \frac{\varphi }{{\phi }_{0}}-\displaystyle \frac{{s}_{a}^{{\prime} }-{s}_{a}^{2}+{s}_{a}}{2}{\left(\displaystyle \frac{\varphi }{{\phi }_{0}}\right)}^{2}\right.\\ \quad \left.+\,O{\left(\displaystyle \frac{\varphi }{{\phi }_{0}}\right)}^{3}\right],\end{array}\end{eqnarray}$
where ${s}_{a}^{{\prime} }=[{{\rm{d}}}^{2}(\mathrm{ln}{m}_{a})/{\rm{d}}{\left(\mathrm{ln}\phi \right)}^{2}]{| }_{{\phi }_{0}}$. For a two-body system, the Newtonian limit of the equation of motion ( $\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}{\boldsymbol{x}}}{{\rm{d}}{t}^{2}}=-\displaystyle \frac{G{\mathscr{G}}m}{{r}^{2}}{\boldsymbol{n}},\end{eqnarray}$
where x = x2 − x1, r = ∣x∣, n = x/r, m = m1 + m2, ${\mathscr{G}}\,=1-{\left(2+{\omega }_{{BD}}\right)}^{-1}({s}_{1}+{s}_{2}-2{s}_{1}{s}_{2})$ and the measured gravitational constant G is $\begin{eqnarray}G=\displaystyle \frac{4+2{\omega }_{{BD}}}{{\phi }_{0}(3+2{\omega }_{{BD}})}.\end{eqnarray}$
The rate of energy loss of the binary system is [27] $\begin{eqnarray}\displaystyle \frac{{\rm{d}}E}{{\rm{d}}t}=-\left\langle \displaystyle \frac{{\mu }^{2}{m}^{2}}{{r}^{4}}\left[\displaystyle \frac{8}{15}\left({\kappa }_{1}{v}^{2}-{\kappa }_{2}{\dot{r}}^{2}\right)\right]+\displaystyle \frac{1}{3}{\kappa }_{D}{{\mathscr{S}}}^{2}\right\rangle ,\end{eqnarray}$
where the angular brackets denote an orbital average, $\begin{eqnarray}\begin{array}{rcl}{\kappa }_{1} & = & {{\mathscr{G}}}^{2}\left[12\left(1-\displaystyle \frac{1}{2}\xi \right)+\xi {{\rm{\Gamma }}}^{2}\right],\\ {\kappa }_{2} & = & {{\mathscr{G}}}^{2}\left[11\left(1-\displaystyle \frac{1}{2}\xi \right)+\displaystyle \frac{1}{2}\xi \left({{\rm{\Gamma }}}^{2}-\displaystyle \frac{5{\rm{\Gamma }}{\rm{\Lambda }}}{{\mathscr{G}}}-\displaystyle \frac{15{{\rm{\Lambda }}}^{2}}{2{{\mathscr{G}}}^{2}}\right)\right],\\ {\kappa }_{D} & = & 2{{\mathscr{G}}}^{2}\xi ,\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{rcl}\xi & = & {\left(2+{\omega }_{{BD}}\right)}^{-1},\qquad \qquad \quad {\mathscr{S}}={s}_{2}-{s}_{1},\\ {\rm{\Gamma }} & = & 1-\displaystyle \frac{2({m}_{1}{s}_{2}+{m}_{2}{s}_{1})}{{m}_{1}+{m}_{2}},\qquad \mu =\displaystyle \frac{{m}_{1}{m}_{2}}{{m}_{1}+{m}_{2}},\\ {\rm{\Lambda }} & = & {\mathscr{G}}(1-{s}_{1}-{s}_{2})-\xi \left[(1-{s}_{1})s{{\prime} }_{2}+(1-{s}_{2})s{{\prime} }_{1}\right].\end{array}\end{eqnarray}$
In equation ( $\begin{eqnarray}E=-\displaystyle \frac{{\mathscr{G}}\mu m}{2a},\qquad f={\pi }^{-1}{\left({\mathscr{G}}m\right)}^{1/2}{a}^{-3/2},\end{eqnarray}$
we get the rate of change of orbital period $\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\dot{P}}{P}=-\displaystyle \frac{\dot{f}}{f}=-\displaystyle \frac{2\dot{E}}{3E}=-\displaystyle \frac{96}{5}\displaystyle \frac{\mu {m}^{2}}{{a}^{4}}F(e)\\ \quad -\displaystyle \frac{2\mu m}{{a}^{3}}{\mathscr{G}}\xi {{\mathscr{S}}}^{2}G(e),\end{array}\end{eqnarray}$
where f is the frequency of emitted gravitational waves, $\begin{eqnarray}\begin{array}{l}F(e)=\displaystyle \frac{1}{12}{\left(1-{e}^{2}\right)}^{-7/2}\left[{\kappa }_{1}\left(1+\displaystyle \frac{7}{2}{e}^{2}+\displaystyle \frac{1}{2}{e}^{4}\right)\right.\\ \quad \left.-{\kappa }_{2}\left(\displaystyle \frac{1}{2}{e}^{2}+\displaystyle \frac{1}{8}{e}^{4}\right)\right],\end{array}\end{eqnarray}$
and $\begin{eqnarray}G(e)={\left(1-{e}^{2}\right)}^{-5/2}\left(1+\displaystyle \frac{1}{2}{e}^{2}\right).\end{eqnarray}$
For small eccentricities we can get an analytic expression for e as a function of f [70] $\begin{eqnarray}e={e}_{0}{\left(\displaystyle \frac{f}{{f}_{0}}\right)}^{-19/18}\left[1+b\left(1-{\left(\displaystyle \frac{f}{{f}_{0}}\right)}^{-2/3}\right)+O({e}_{0}^{2})\right],\end{eqnarray}$
where $b=\tfrac{25}{288}\xi {{\mathscr{S}}}^{2}{\left(\pi {f}_{0}{\mathscr{G}}m\right)}^{-2/3}$ and e0 is the initial eccentricity at the initial frequency f0. For simplicity, we define the renormalized mass ${m}_{r}={\mathscr{G}}m$. We can get the time and frequency evolutions of gravitational radiations from equation ( $\begin{eqnarray}\begin{array}{l}f\displaystyle \frac{{\rm{d}}t}{{\rm{d}}f}=\displaystyle \frac{5{a}^{4}}{2\mu \left(5{{am}}_{r}\xi {{\mathscr{S}}}^{2}+4{\lambda }_{1}{m}_{r}^{2}\right)}\\ \quad -\displaystyle \frac{5{e}^{2}{a}^{4}\left(15a\xi {{\mathscr{S}}}^{2}+28{\lambda }_{1}{m}_{r}-2{\lambda }_{2}{m}_{r}\right)}{2\mu {m}_{r}{\left(5a\xi {{\mathscr{S}}}^{2}+4{\lambda }_{1}{m}_{r}\right)}^{2}},\\ \quad =\displaystyle \frac{5\left(4{\lambda }_{1}{\left({{fm}}_{r}\right)}^{2/3}+5\xi {{\mathscr{S}}}^{2}\right)}{2{f}^{2}\mu {\left(4{\lambda }_{1}{\left({{fm}}_{r}\right)}^{2/3}+5\xi {{\mathscr{S}}}^{2}\right)}^{2}}\\ \quad -\displaystyle \frac{5{e}^{2}\left(28{\lambda }_{1}{\left({{fm}}_{r}\right)}^{2/3}-2{\lambda }_{2}{\left({{fm}}_{r}\right)}^{2/3}+15\xi {{\mathscr{S}}}^{2}\right)}{2{f}^{2}\mu {\left(4{\lambda }_{1}{\left({{fm}}_{r}\right)}^{2/3}+5\xi {{\mathscr{S}}}^{2}\right)}^{2}},\end{array}\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{rcl}{\lambda }_{1} & = & 12\left(1-\displaystyle \frac{1}{2}\xi \right)+\xi {{\rm{\Gamma }}}^{2},\\ {\lambda }_{2} & = & 11\left(1-\displaystyle \frac{1}{2}\xi \right)+\displaystyle \frac{1}{2}\xi \left({{\rm{\Gamma }}}^{2}-\displaystyle \frac{5{\rm{\Gamma }}{\rm{\Lambda }}}{{\mathscr{G}}}-\displaystyle \frac{15{{\rm{\Lambda }}}^{2}}{2{{\mathscr{G}}}^{2}}\right).\end{array}\end{eqnarray}$
For large ωBD, ξ ≪ 1 and equation ( $\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{\rm{d}}t}{{\rm{d}}f}=\displaystyle \frac{5\left(2{\lambda }_{1}-14{e}^{2}{\lambda }_{1}+{e}^{2}{\lambda }_{2}\right)}{16{\pi }^{8/3}{f}^{11/3}{\lambda }_{1}^{2}\mu {m}_{r}^{2/3}}\\ \quad +\displaystyle \frac{25\xi {{\mathscr{S}}}^{2}\left(11{e}^{2}{\lambda }_{1}-{e}^{2}{\lambda }_{2}-{\lambda }_{1}\right)}{32{\pi }^{10/3}{f}^{13/3}{\lambda }_{1}^{3}\mu {m}_{r}^{4/3}}+O(\xi ).\end{array}\end{eqnarray}$
In the limit ωBD → + ∞ , we get $\begin{eqnarray}\begin{array}{l}\xi =0,\qquad {\lambda }_{1}=12,\qquad {\lambda }_{2}=11,\\ \quad a={\pi }^{-2/3}{m}^{1/3}{f}^{-2/3},\end{array}\end{eqnarray}$
and the rate at which the frequency changes because of the back-reaction due to the emission of gravitational waves reduces back to the result of GR, $\begin{eqnarray}\displaystyle \frac{{\rm{d}}t}{{\rm{d}}f}=\displaystyle \frac{5}{96{f}^{11/3}\mu {m}^{2/3}}\left(1-\displaystyle \frac{157{e}^{2}}{24}\right).\end{eqnarray}$
The Fourier transformation of a signal h(t) is $\begin{eqnarray}\tilde{h}(f)={\int }_{-\infty }^{+\infty }{\rm{d}}{{te}}^{2\pi {\rm{i}}{ft}}h(t).\end{eqnarray}$
Under the stationary phase approximation, $\tilde{h}(f)$ becomes $\begin{eqnarray}\tilde{h}(f)={ \mathcal A }{f}^{-7/6}{e}^{{\rm{i}}{\rm{\Psi }}(f)},\end{eqnarray}$
where the amplitude ${ \mathcal A }$ and the phase $Psi$(f) are given by $\begin{eqnarray}{ \mathcal A }=\displaystyle \frac{1}{\sqrt{30}{\pi }^{2/3}}\displaystyle \frac{{{ \mathcal M }}^{5/6}}{{D}_{L}},\end{eqnarray}$
$\begin{eqnarray}{\rm{\Psi }}(f)=2\pi {ft}(f)-\phi (f)-\pi /4,\end{eqnarray}$
${ \mathcal M }={\left({m}_{1}{m}_{2}\right)}^{3/5}{\left({m}_{1}+{m}_{2}\right)}^{-1/5}$ is the chirp mass and DL is the luminosity distance. By integrating equation ( $\begin{eqnarray}t(f)={\int }_{{f}_{a}}^{f}\displaystyle \frac{{\rm{d}}t}{{\rm{d}}f}{\rm{d}}f.\end{eqnarray}$
From dφ/dt = 2πf, we get the phase evolution of gravitational radiations $\begin{eqnarray}\phi (f)=2\pi {\int }_{{f}_{a}}^{f}f\displaystyle \frac{{\rm{d}}t}{{\rm{d}}f}{\rm{d}}f.\end{eqnarray}$
Combining equations ( $\begin{eqnarray}\begin{array}{l}{{\rm{\Psi }}}_{{BD}}={{\rm{\Psi }}}_{{GR}}^{e=0}+\displaystyle \frac{9}{32{\pi }^{5/3}{f}^{5/3}\mu {m}_{r}^{2/3}}\left(\displaystyle \frac{1}{{\lambda }_{1}}-\displaystyle \frac{1}{12}\right)\\ \quad -\displaystyle \frac{45\xi {{\mathscr{S}}}^{2}}{224{\pi }^{7/3}{f}^{7/3}{\lambda }_{1}^{2}\mu {m}_{r}^{4/3}}\\ \quad -\displaystyle \frac{405{e}_{0}^{2}{f}_{0}^{19/9}\left(392{\lambda }_{1}{\left(\pi {{fm}}_{r}\right)}^{2/3}(14{\lambda }_{1}-{\lambda }_{2})+731\xi {{\mathscr{S}}}^{2}({\lambda }_{2}-11{\lambda }_{1})\right)}{4584832{\pi }^{7/3}{f}^{40/9}{\lambda }_{1}^{3}\mu {m}_{r}^{4/3}}\\ \quad +\xi {{\mathscr{S}}}^{2}{e}_{0}^{2}{f}_{0}^{19/9}\displaystyle \frac{225(14{\lambda }_{1}-{\lambda }_{2})\left(731-980{\left(f/{f}_{0}\right)}^{2/3}\right)}{36\,678\,656{\pi }^{7/3}{f}^{40/9}{\lambda }_{1}^{2}\mu {m}_{r}^{4/3}},\end{array}\end{eqnarray}$
where ${{\rm{\Psi }}}_{{GR}}^{e=0}$ is the result in GR without eccentricity. The contribution of the eccentricity, monopole and quadrupole gravitational radiations to the phase of GWs in the BD theory is $\begin{eqnarray}\begin{array}{l}{{\rm{\Psi }}}_{{em}}(f)=\displaystyle \frac{9}{32{\pi }^{5/3}{f}^{5/3}\mu {m}_{r}^{2/3}}\left(\displaystyle \frac{1}{{\lambda }_{1}}-\displaystyle \frac{1}{12}\right)\\ \quad -\displaystyle \frac{405{e}_{0}^{2}{f}_{0}^{19/9}\left(392{\lambda }_{1}{\left(\pi {{fm}}_{r}\right)}^{2/3}(14{\lambda }_{1}-{\lambda }_{2})+731\xi {{\mathscr{S}}}^{2}({\lambda }_{2}-11{\lambda }_{1})\right)}{4584832{\pi }^{7/3}{f}^{40/9}{\lambda }_{1}^{3}\mu {m}_{r}^{4/3}}\\ \quad +\xi {{\mathscr{S}}}^{2}{e}_{0}^{2}{f}_{0}^{19/9}\displaystyle \frac{225(14{\lambda }_{1}-{\lambda }_{2})\left(731-980{\left(f/{f}_{0}\right)}^{2/3}\right)}{36\,678\,656{\pi }^{7/3}{f}^{40/9}{\lambda }_{1}^{2}\mu {m}_{r}^{4/3}}.\end{array}\end{eqnarray}$
The contribution of dipole gravitational radiation to the phase of GWs in the BD theory is $\begin{eqnarray}{{\rm{\Psi }}}_{d}(f)=-\displaystyle \frac{45\xi {{\mathscr{S}}}^{2}}{224{\left(\pi f\right)}^{7/3}{\lambda }_{1}^{2}\mu {m}_{r}^{4/3}}.\end{eqnarray}$
For GWs coming from a binary system composed of compact objects in the BD theory, we only consider the post-Newtonian corrections to the phase of GWs while keeping the amplitude in the form of GR because the correction in the phase is dominant.3. Analysis of the selected neutron star-black hole signals
Our baseline waveform model is the frequency domain inspiral-merger-ringdown approximant IMRPhenomXPHM [71–74] with the corrections on the phase from equations (29 ) and (30 ). The waveform template IMRPhenomXPHM includes the effects of spin-induced orbital precession and higher order multipole GW moments. Given a waveform approximate $\tilde{h}(f,\bar{\theta })$ with the modified phase (28 ) and a detected GW strain data d, the posterior density function (PDF) for the parameters $\bar{\theta }$ is derived from Bayesian theorem
$\begin{eqnarray}p(\bar{\theta }| d)=\displaystyle \frac{p(d| \bar{\theta })p(\bar{\theta })}{p(d)},\end{eqnarray}$
where $p(\bar{\theta })$ is the prior density and p(d) is the evidence. The likelihood $p(d| \bar{\theta })$ is given by $\begin{eqnarray}p(d| \bar{\theta })\propto \exp \left[-\displaystyle \frac{\langle \tilde{d}-\tilde{h}(f,\bar{\theta })| \tilde{d}-\tilde{h}(f,\bar{\theta })\rangle }{2}\right],\end{eqnarray}$
where the noise-weighted inner product ⟨·∣ · ⟩ is defined as $\begin{eqnarray}\langle a| b\rangle =\mathrm{Re}\left[4{\int }_{{f}_{{\rm{low}}}}^{{f}_{{\rm{high}}}}{\rm{d}}f\displaystyle \frac{{\tilde{a}}^{* }(f)\tilde{b}(f)}{{S}_{n}(f)}\right],\end{eqnarray}$
flow is the detector's lower cut-off frequency, fhigh is the frequency at which a given signal ends, and Sn(f) is the power spectral density of the noise as shown in figure 1. In this paper we set flow = 20 Hz for LIGO Hanford (H1) and Virgo (V1), flow = 25 Hz for LIGO Livingston (L1), and fhigh = 2048 Hz for H1, L1 and V1.
Figure 1. The amplitude spectral density $\sqrt{{S}_{n}}$ of H1 with noise level (black line), L1 (blue line) and V1 (green line). |
Since the sensitivity s = 0.5 for a black hole, so binary black hole systems are not suitable for our analysis because ${\mathscr{S}}=0$. As the consequence of the no-hair theorem for black holes [75, 76], the BD scalar field is radiated away, so the motion for a binary black hole system in the BD theory is identical to that in GR after a mass rescaling [27, 62]. Similarly, binary neutron star systems are also not effective in distinguishing the BD theory from GR. This is because masses of neutron stars tend to be around the Chandrasekhar limit 1.4M⊙, and the sensitivity of neutron stars is not a strong function of mass for a given equation of state. Therefore, we use binaries consisting of a black hole and a neutron star to distinguish the BD theory from GR. Currently, there are four detected GW events that are possible neutron star-black hole binaries, GW190426_152155, GW190814, GW200105 and GW200115. In this paper, we consider the two events GW200105 and GW200115 only. Although GW200105 was considered as effectively a single-detector event in L1 because H1 was not operational and the signal-to-noise ratio (SNR) in V1 is small, we still use L1 and V1 data. The signal of GW200115 in V1 is not loud enough to further improve the significance of the coincident detection by the LIGO detectors, so we use the H1 and L1 data only.
The likelihood (32 ) is sampled using the Dynesty algorithm, and the waveform model described in the previous section is implemented as an extension to the approximant IMRPhenomXPHM [74] in Bilby [77]. The intrinsic parameters being sampled over are the mass m1 and spin of the black hole, the mass m2 and spin of the neutron star, and the scalar-tensor parameter $\xi ={\left(2+{\omega }_{{BD}}\right)}^{-1}$. In this paper, we take s1 = 0.5 and s2 = 0.3, and the sampled parameters are uniformly distributed. In particular, we take m1 ∈ [1, 20] M⊙, m2 ∈ [0.5, 3] M⊙, and ξ ∈ [0, 0.2], the amplitudes of the component spins a1 and a2 are uniform in the range [0, 0.99], and the directions are uniform on the unit sphere. For the extrinsic parameters, the sky location and orientation of the binary are uniform on the unit sphere.
The posterior probabilities for the parameters in GR for the events GW200105 and GW200115 are shown in figures 2 and 3. For GW200105, the mass m1 is between [8.09, 11.92] M⊙ and m2 is between [1.72, 2.29] M⊙ in table. 1. For GW200115, the mass m1 is between [6.80, 13.61] M⊙ and m2 is between [0.88, 1.44] M⊙ in table. 2. These results, which are consistent with the released results [13], show that our data and analysis methods are correct and credible.
Figure 2. Posterior probability distributions for parameters in GR theory for the event GW200105. The dashed vertical lines indicate 16% and 84% percentiles for the event. m1 is the black hole mass and m2 is the neutron star mass, a1,2 are their dimensionless spin magnitudes, θ1,2 are the tilt angle between their spins and the orbital angular momentum, and the two spin vectors describing the azimuthal angle separating the spin vectors Δφ and the cone of precession about the system's angular momentum φJL. The extrinsic parameters are the luminosity distance dL, the right ascension RA and declination DEC, the inclination angle between the observer's line of sight and the orbital angular momentum θJN, the polarization angle $Psi$ and the coalescence time tc. |
Figure 3. Posterior probability distributions for parameters in GR theory for the event GW200115. The parameters are the same as in figure 2. |
Table 1. Results of parameters with 90% confidence level for the event GW200105. |
| GR | BD | ||||||
|---|---|---|---|---|---|---|---|
| Parameters | Published | e = 0 | e ≠ 0 | EccentricFD | e0 = 0.001 | e0 = 0.01 | e0 = 0.1 |
| m1 | ${8.9}_{-1.5}^{+1.2}$ | ${9.78}_{-1.69}^{+2.14}$ | ${10.57}_{-1.81}^{+1.94}$ | ${9.68}_{-0.14}^{+0.17}$ | ${9.37}_{-1.85}^{+2.33}$ | ${9.91}_{-2.02}^{+5.87}$ | ${11.7}_{-2.33}^{+2.34}$ |
| m2 | ${1.9}_{-0.2}^{+0.3}$ | ${1.98}_{-0.26}^{+0.31}$ | ${1.89}_{-0.22}^{+0.28}$ | ${2.0}_{-0.04}^{+0.02}$ | ${1.98}_{-0.33}^{+0.32}$ | ${1.90}_{-0.59}^{+0.35}$ | ${1.76}_{-0.22}^{+0.31}$ |
| dL | ${280}_{-110}^{+110}$ | ${341.49}_{-91.33}^{+79.71}$ | ${320.31}_{-78.44}^{+80.92}$ | ${315.97}_{-128.05}^{+112.58}$ | ${342.95}_{-82.24}^{+76.65}$ | ${341.01}_{-85.40}^{+78.69}$ | ${322.59}_{-83.67}^{+80.48}$ |
| a1 | ${0.08}_{-0.08}^{+0.22}$ | ${0.14}_{-0.12}^{+0.23}$ | ${0.23}_{-0.18}^{+0.15}$ | × | ${0.22}_{-0.18}^{+0.24}$ | ${0.24}_{-0.20}^{+0.28}$ | ${0.30}_{-0.14}^{+0.12}$ |
| χeff | $-{0.01}_{-0.15}^{+0.11}$ | ${0.01}_{-0.17}^{+0.16}$ | ${0.09}_{-0.15}^{+0.13}$ | × | ${0.0}_{-0.17}^{+0.19}$ | ${0.05}_{-0.19}^{+0.33}$ | ${0.20}_{-0.18}^{+0.12}$ |
| e0 | × | × | ${0.07}_{-0.03}^{+0.02}$ | ≤0.2 | × | × | × |
| ξ | × | × | × | × | ${0.09}_{-0.08}^{+0.09}$ | ${0.08}_{-0.07}^{+0.11}$ | ≤0.01 |
Table 2. Results of parameters with 90% confidence level for the event GW200115. |
| GR | BD | ||||||
|---|---|---|---|---|---|---|---|
| Parameters | Published | e = 0 | e ≠ 0 | EccentricFD | e0 = 0.001 | e0 = 0.01 | e0 = 0.1 |
| m1 | ${5.7}_{-2.1}^{+1.8}$ | ${7.47}_{-0.67}^{+6.14}$ | ${7.44}_{-0.60}^{+0.83}$ | ${6.54}_{-0.17}^{+0.14}$ | ${7.64}_{-0.68}^{+0.86}$ | ${7.68}_{-0.69}^{+1.16}$ | ${7.43}_{-0.5}^{+0.50}$ |
| m2 | ${1.5}_{-0.3}^{+0.7}$ | ${1.34}_{-0.46}^{+0.10}$ | ${1.35}_{-0.10}^{+0.09}$ | ${1.5}_{-0.03}^{+0.03}$ | ${1.29}_{-0.10}^{+0.10}$ | ${1.29}_{-0.16}^{+0.11}$ | ${1.35}_{-0.07}^{+0.07}$ |
| dL | ${300}_{-100}^{+150}$ | ${427.99}_{-121.62}^{+112.56}$ | ${435.50}_{-131.69}^{+114.0}$ | ${345.59}_{-148.32}^{+220.21}$ | ${446.58}_{-132.48}^{+112.49}$ | ${428.67}_{-142.81}^{+114.48}$ | ${454.28}_{-137.65}^{+112.90}$ |
| a1 | ${0.33}_{-0.29}^{+0.48}$ | ${0.10}_{-0.09}^{+0.52}$ | ${0.09}_{-0.08}^{+0.20}$ | × | ${0.11}_{-0.10}^{+0.21}$ | ${0.13}_{-0.12}^{+0.25}$ | ${0.07}_{-0.06}^{+0.12}$ |
| χeff | $-{0.19}_{-0.35}^{+0.23}$ | $-{0.01}_{-0.09}^{+0.37}$ | ${0.01}_{-0.09}^{+0.10}$ | × | ${0.05}_{-0.09}^{+0.09}$ | ${0.05}_{-0.11}^{+0.12}$ | ${0.06}_{-0.08}^{+0.07}$ |
| e0 | × | × | ${0.06}_{-0.05}^{+0.05}$ | ≤0.04 | × | × | × |
| ξ | × | × | × | × | ${0.04}_{-0.03}^{+0.15}$ | ${0.03}_{-0.03}^{+0.12}$ | ${0.02}_{-0.02}^{+0.05}$ |
Before performing Bayesian analysis to constrain parameter ωBD in the BD theory, it is necessary to verify whether the events GW200105 and GW200115 have eccentricity. We can get the modified waveform with small eccentricity in GR theory by setting the scalar-tensor parameter ξ = 0. We take e0 uniformly distributed in [0, 0.2] and match the modified waveform with the events GW200105 and GW200115 to get Bayesian inference of the parameters and the results are shown in figures 4 and 5. From figures 4 and 5, we can see that the events GW200105 and GW200115 probably have eccentric orbit and their eccentricities are close to 0.1. By assuming different initial eccentricity e0 at f0 = 20 Hz, we match the modified waveform in the BD theory with the events GW200105 and GW200115 to get Bayesian inference of the parameters and the results are shown in tables 1 and 2. For GW200105, we get the constraints ωBD > {9.1, 10.5, 98} for different initial eccentricity e0 = {0.001, 0.01, 0.1}. For GW200115, we get the constraints ωBD > {23, 31.3, 48} for different initial eccentricity e0 = {0.001, 0.01, 0.1}. The larger initial eccentricity helps constrain the BD parameter ωBD. As we can see, the eccentricity can excite additional energy flux channels from monopole mode and quadrupolar mode for the scalar field [78]. So the dephasing caused by the scalar field for the eccentric orbit will be larger than the circular orbit, thus constraining the BD parameter ωBD better. Additionally, the events GW200115 and GW200105 may have eccentric orbits in reality and the modified GW matches the data well to give the better constraint on the BD parameter ωBD.
Figure 4. Posterior probability distributions for parameters in GR theory for the event GW200105. The parameters are the same as in figure 2 except the eccentric parameter e0. |
Figure 5. Posterior probability distributions for parameters in GR theory for the event GW200115. The parameters are the same as in figure 2 except the eccentric parameter e0. |
4. Conclusion
Previous works constructed modified gravitational waves in the BD theory considering only the dipolar emission from the circular orbits. However, in reality, the orbit usually has eccentricity. In this paper, we analyze the influence of gravitational waves from eccentricity in the BD theory and use the modified gravitational waveform model to match the observed GW events such as GW200105 and GW200115 to constrain the BD parameter ωBD. We first give the full phase deviation from small eccentric compact binaries in the BD theory caused by the combined effects of monopole, dipole and quadrupole gravitational radiation. We construct the waveform model from the system of small eccentric neutron star-black hole coalescing in the BD theory and use this model to conduct a Bayesian analysis based on the actual signals from the event GW200115 and GW200105 to constrain the BD theory. The contribution of the eccentricity to the phase of GWs in the BD theory is important for the reason that monopole and quadrupole gravitational radiation can be influenced by the eccentric orbit. We can see that the phase deviation between GR and the BD theory for eccentric trajectory is more observable than that for circular orbit. We give bounds on the BD coupling parameter ωBD > {9.1, 10.5, 98} for the event GW200105 and ωBD > {23, 31.3, 48} for the event GW200115 with different initial eccentricity e0 = {0.001, 0.01, 0.1} at f0 = 20 Hz. The constraint ωBD from the event GW200105 is tighter than that from the event GW200115. The reason is that the signal from the event GW200105 is stronger and more eccentric than that from the event GW200115. So we can constrain the BD coupling parameter ωBD > 98 from the event GW200105 and ωBD > 48 from the event GW200115 with small eccentricity e0 = 0.1. Including the eccentricity can help us constrain the BD parameter ωBD better in the strong-field region.