Future Prospects for Probing Scalar-Tensor Theories with Gravitational Waves from Mixed Binaries

TL;DR

This paper explores how gravitational wave observations from mixed black hole-neutron star binaries can test and constrain scalar-tensor theories beyond general relativity, with current bounds improvable and future detections promising significant advancements.

Contribution

It analyzes the potential of gravitational wave data to constrain scalar-tensor theories, including those with spontaneous scalarization and string theory motivations, highlighting prospects for future improvements.

Findings

Current bounds can be improved with existing detectors for small black hole masses.

Future detections can enhance constraints by several orders of magnitude.

Combining multiple events significantly boosts the testing power.

Abstract

The extreme-gravity collisions of binaries with one black hole and one neutron star provide for excellent tests of general relativity. We here study how well one can constrain theories beyond general relativity with additional scalar fields that allow for spontaneous scalarization of neutron stars, and those motivated from string theory. We find that existing bounds can be improved with current gravitational-wave detectors if the black hole mass is sufficiently small. Bounds will further improve by many orders of magnitude with future detections, especially by combining multiple events.

Figures

Tables

Table 1. Table 1: Upper and lower frequency limits used in the Fisher matrix integrations of Ref.  [ 65 ] . For ground-based detectors the upper cutoff frequency is determined to be f ISCO subscript 𝑓 ISCO f_{\text{ISCO}} for a BH-NS binary with masses 10 ​ M ⊙ 10 subscript 𝑀 direct-product 10M_{\odot} and 1.4 ​ M ⊙ 1.4 subscript 𝑀 direct-product 1.4M_{\odot} , while space-based detectors use lower cutoff frequencies determined to be one year prior to merger  [ 24 ] .

Detector	$f_{\text{low}}$ [Hz]	$f_{\text{high}}$ [Hz]
A+	10	386
Voyager	10	386
CE	5	386
ET	5	386
B-DECIGO	0.12	100
DECIGO	0.12	100

Table 2. Table 2: Pessimistic, mean, and optimistic 1 year detection rates for 10 ​  M ⊙ − 1.4 ​  M ⊙ 10 subscript  M direct-product 1.4 subscript  M direct-product 10\text{ M}_{\odot}-1.4\text{ M}_{\odot} BH-NS binaries assuming a local BH-NS coalescence rate of ℛ ∈ [ 0.6 , 610 ] ​ Gpc − 3 ​ yr − 1 ℛ 0.6 610 superscript Gpc 3 superscript yr 1 \mathcal{R}\in[0.6,610]\text{Gpc}^{-3}\text{yr}^{-1} , with a “realistic value” of 30 ​  Gpc − 3 ​ yr − 1 30 superscript  Gpc 3 superscript yr 1 30\text{ Gpc}^{-3}\text{yr}^{-1}   [ 45 , 12 ] .

Detector	Detection Rate
	Pessimistic	Realistic	Optimistic
A+	5	270	5,500
Voyager	72	3,600	74,000
CE	720	36,000	730,000
ET	510	25,000	520,000
B-DECIGO	43	2,200	44,000
DECIGO	730	37,000	730,000

Equations

$$N=\Delta T\underset{0}{\overset{z_h}{\int }}4\pi [a_{0}r(z){]}^2\mathcal{R}R(z)\frac{d\tau }{dz}dz,$$

$${\sigma }_{{\theta }^a}^{-2}=\Delta T\underset{0}{\overset{z_h}{\int }}4\pi [a_{0}r(z){]}^2\mathcal{R}R(z)\frac{d\tau }{dz}{\sigma }_{{\theta }^a}(z{)}^{-2}dz,$$

$${\beta }_{\text{mbox}DEF}=-\frac{5{\eta }^{2/5}(\Delta \alpha {)}^2}{7168},n=-1,$$

$$\frac{{\dot{P}}_b^D}{{\dot{P}}_{\text{GR}}}=\frac{5}{96}\frac{G}{1+{\alpha }_0^2}\frac{({\Omega }_b{m}_2{)}^{-2/3}}{(1+q{)}^{2/3}}\frac{1+e^2/2}{1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4}(\Delta \alpha {)}^2<{\delta }_{{\dot{P}}_b},$$

$$\delta \psi ={\beta }_{\text{mbox}DEF}{u}^{-7}\left(1+\frac{2623+2640\eta }{4280}{u}^2-6\pi u^3\right),$$

$${\sigma }_{{\theta }^a}^{-2}=\Delta T\underset{5{\text{ M}}_{\odot }}{\overset{40{\text{ M}}_{\odot }}{\int }}\underset{1{\text{ M}}_{\odot }}{\overset{2{\text{ M}}_{\odot }}{\int }}\underset{0}{\overset{z_h(m_{\text{mbox}NS},m_{\text{mbox}BH})}{\int }}4\pi [a_{0}r(z){]}^2\mathcal{R}R(z)$$

$$\times f(m_{\text{mbox}BH},m_{\text{mbox}NS})dzd{m}_{\text{mbox}NS}d{m}_{\text{mbox}BH},$$

$$f_{\mbox{\tiny BH}}(m_{\mbox{\tiny BH}})\propto\left\{\begin{array}[]{ll}{\left(\frac{m_{\mbox{\tiny BH}}}{\text{M}_{\odot}}\right)^{-\alpha}}&{\text{ if }m_{\min}\leq m_{\mbox{\tiny BH}}\leq m_{\max}}\\ {0}&{\text{ otherwise }}\end{array}\right\},\hskip 14.22636ptf_{\mbox{\tiny NS}}(m_{\mbox{\tiny NS}})\propto\mathcal{N}(\mu_{\mbox{\tiny NS}},\sigma_{\mbox{\tiny NS}}).$$

Full text

1 Introduction

Einstein’s famous theory of general relativity (GR) describes the relationship between matter, and the curvature of spacetime through a single tensor, the metric $g_{\mu\nu}$. Over the last 100 years, several alternative theories of gravity have been proposed – yet countless observations and tests of GR in a variety of environments have proven to be consistent with Einstein’s theory: solar system [1], pulsar timing [2, 3] and cosmological [4, 5, 6, 7, 8] observations. With such a rigorous upholding to all of our tests, why should we bother to further prove or disprove the solid theory of GR? The answer lies within the many unanswered questions we still have today, obscuring our understanding of the universe. Some examples of these include the unification of GR and quantum mechanics, the inflationary period in the early universe, dark matter and its influence on the galactic rotation curves, and dark energy and its impact on the expansion of our universe. Alternative theories of gravity, if found to be true, could explain the missing links to our unanswered questions.

More recently, a new observational opportunity into the extreme-gravity (strong, non-linear, and highly-dynamical) regime [9, 10] has been unveiled with the recent gravitational wave (GW) detections of coalescing black holes (BHs) [11, 12] and neutron stars (NSs) [13]. To date, eleven GW observations have been made [12], and have further found no statistically significant deviations from GR [14, 15, 16].

In this analysis, we consider the present and future implications on constraining non-GR theories with an additional massless scalar field, known as scalar tensor theories (STTs), using both GWs and pulsar timing observations. The latter has been studied for binary pulsar systems [3, 17], pulsar-white dwarf (WD) systems [18, 3, 19, 17], and triple pulsar-WD-WD systems [20, 21]. Here, we consider the present and future constraints obtained from the GW detections of BH-NS coalescences. See also Refs. [22, 23] for constraints on STTs from GWs. In STTs, compact objects acquire scalar charges that source the scalar field. A scalar force acts between two scalarized objects, giving rise to a fifth force which depends on the internal structure of the massive objects and violates the strong equivalence principle (SEP). Binaries consisting of scalarized astrophysical objects further emit scalar dipole radiation (on top of gravitational quadrupolar radiation in GR), causing the binaries to evolve faster.

Such radiation becomes larger when the difference between the scalar charges of the binary constituents become larger, and thus a mixed binary consisting of one black hole and one neutron star system is ideal for probing such theories [24, 25, 26, 27]. Specifically, the increased mass difference $m_{1}^{2}-m_{2}^{2}$ and small total mass $M$ of a BH/NS system will minimize the allowable constraints on $\sqrt{\alpha}\sim M$ (for non-GR coupling parameter $\alpha$), while the offset acquired from the increased SNR of the alternative scenario of a large-mass binary BH system111In many cases, low-mass binary BH systems and BH/NS systems with slowly rotating BHs may be indistinguishable. However, in the Einstein-dilaton Gauss-Bonnet theory of gravity, the dipole-radiation slightly decreases when comparing a BH/NS system with the equivalent binary BH system with the NS replaced by a slowly rotating BH but the effect is insignificant, and thus the constraints are not significantly affected. only improves constraints by a small factor of $\sim\text{SNR}^{-1/4}$, while the SNR only itself increases by $\sim M^{5/6}$. Moreover, a smaller binary system has a lower relative velocity for a fixed frequency, which leads to an enhanced dipole radiation. It is then extremely advantageous to decrease the total mass of the system rather than maximize it. Such sources are particularly interesting and extremely timely to consider as two of the candidates in the O3 run by the LIGO/Virgo Collaborations, S190426c and S190814bv, are likely to be the merger of a black hole and a neutron star, if they are of astrophysical origin [28, 29]222We note that the LIGO/Virgo Collaborations categorize BH/NS candidate events as $m_{1}>5\text{ M}_{\odot}$ and $m_{2}<3\text{ M}_{\odot}$ [30], and therefore O3 candidates S190426c and S190814bv potentially could be binary BH mergers. In addition, S190426c currently has a 58% possibility of being terrestrial noise.. Because BHs have vanishing scalar charges in the quasi-Brans-Dicke theory of gravity considered in this analysis, the presence of a NS is required to place constraints on such a theory. Therefore, it is of vital importance that, not only the events be of astrophysical origin, but they must with high confidence also be a BH/NS system rather than a binary BH system.

A particularly interesting class of theories within STTs are quasi-Brans-Dicke theory and Einstein-dilaton Gauss-Bonnet (EdGB) gravity. The former was introduced first by Damour and Esposito-Far$\acute{e}$se (DEF) [31, 32] which induces a non-linear growth of the scalar charges onto neutron stars called *spontaneous scalarization * [31, 33, 32], while black holes remain hairless as in GR. The latter is motivated from string theory and the dilaton scalar field is coupled with a quadratic curvature term (Gauss-Bonnet invariant) in the gravitational action [34, 35]. In this theory, black holes have non-vanishing scalar charges [36, 37, 38, 39] while neutron stars do not if the scalar field coupling is linear [38, 40]. We here consider the single BH-NS detections with future GW detectors, as well as the multi-band detections between both space- and ground-based detectors [41, 42, 43, 44], and finally the combination of multiple observations [25, 45] made on future detectors with expanded horizons.

The outline of the paper is as follows. We begin in Sec. 2 with a discussion on the gravitational waveform model, GW interferometers, Fisher analysis techniques used in our analysis, together with information on combining uncertainties from multiple events. In Sec. 3, we review the theoretical background on STTs and present our main results. Section 4 provides an analysis and description of the validity of the various approximations made in our investigation. Finally, we conclude and offer avenues of future direction in Sec. 5. Throughout this paper, we have adopted the geometric units of $G=1=c$.

2 Non-GR gravitational waveforms and data analysis

Let us begin by considering how to capture modifications to GR in the gravitational waveform from compact binary mergers. Typically, one strives to be agnostic towards the multitudes of alternative theories of gravity available, by modifying the phase of the GR waveform by $\Psi_{\text{GR}}+\Delta\Psi$, where $\Delta\Psi\equiv\beta u^{2n-5}$. Here, $u\equiv(\pi\mathcal{M}f)^{1/3}$ is the effective relative velocity of binary constituents with GW frequency $f$ and chirp mass $\mathcal{M}\equiv(m_{1}^{3}m_{2}^{3}/M)^{1/5}$ with individual masses $m_{A}$ and total mass $M=m_{1}+m_{2}$. The post-Newtonian (PN) parameter333An $n$PN-order term is proportional to $(u/c)^{2n}$ relative to the leading-order term in the waveform. $n$ categorizes the order at which the modifications affect the GW phase, and $\beta$ prescribes the magnitude of the modification. This parameterized post-Einsteinian (ppE) formalism [46] allows for one to effectively constrain any modified theory of gravity by knowing the ppE expression $\beta$ (which can be mapped to the coupling parameters of SEP-violating theories) and the PN order $n$ at which the leading correction enters the waveform, many of which are tabulated in e.g. [47].

Let us start with the detector sensitivity curves utilized in this analysis: aLIGO (representative of observing run O3) [48, 49] (O3), LIGO A+ [50, 51], LIGO Voyager [51, 52], Cosmic Explorer [51, 53, 52] (CE, standard configuration described in Fig. 1 of [53]), Einstein Telescope configuration D [51, 52, 53] (ET), B-DECIGO [54], and DECIGO [55]. Figure 1 shows the resulting sensitivity curves of each interferometer, as well as the characteristic amplitudes $2\sqrt{f}|\tilde{h}(f)|$ for $10\text{ M}_{\odot}-1.4\text{ M}_{\odot}$ BH-NS systems located $100$ Mpc and $1$ Gpc away. Here, it can be seen that only the former system can be detected by the O3 detector.

We implement the sky-averaged IMRPhenomD waveform model (such that the inclination angle and sky-location parameters are averaged over) found in Refs. [56, 57]. Because accurate BH/NS tidal waveform models do not yet exist to sufficient accuracy for the purposes of this investigation444See Ref. [58], where a phenomenological BH/NS waveform model was constructed where the phase exceeded the NR results by 30%. See also Ref. [59], which updates the model from Ref. [58] with a more accurate baseline binary BH model, and Ref. [60] for a BH/NS amplitude model, Ref. [61] for an effective-one-body model applicable to BH/NS systems, or Ref. [62] for BH/NS models computed with tidal splicing. Finally, refer to Ref. [63] for an analysis on the waveform systematic uncertainties present in such models., we modify the waveform with the simple $5$PN+$6$PN finite size tidal corrections for the NS [64]. However, it is unlikely that such tidal effects (first entering at 5PN) are vital in the constraint of STT parameters (first entering at -1PN) due to minuscule correlations between the two. Additionally, we do not consider the effects of spin precession. See Sec. 4.2 for a discussion on the inclusion of this effect. We also note that in this waveform model, there is no contribution from the finite-size effect to the spin-induced quadrupole moment, first entering at $2$PN order. However, this is likely to make a negligible effect in our results due to the small spin priors assumed on the NS, as described below.

To estimate constraints on $\beta$ from future GW observations, we utilize a Fisher analysis as described in Ref. [65]. Assuming sufficiently loud signals and Gaussian-distributed noise, Fisher-techniques allow one to predict the root-mean-square error on waveform template parameters $\theta^{a}$ to be $\Delta\theta^{a}\approx\sqrt{(\Gamma^{-1})^{aa}}$. Here, $\bm{\Gamma}$ is the Fisher information matrix detailed in Eq. (2.7) of Ref. [65], and depends on the detector sensitivity and waveform model. The waveform model used has template parameters of luminosity distance $D_{L}$, masses $m_{A}$ and spins $\chi_{A}$ of the BH and NS (See B for a demonstration of the importance of including spin-effects in the waveform), the time and phase at coalescence $t_{c}$ and $\phi_{c}$, NS tidal deformability parameter $\Lambda$555The tidal deformability of a non-rotating BH is zero [66, 67, 68, 69, 70]. See also Ref. [71] in which it was shown that the standard computation of the tidal deformability relies on the comparison to the BH one and the effective BH tidal deformability may have a small non-zero effect on the gravitational waveform., and the ppE parameter $\beta$. We assume Gaussian priors on the spins to be $|\chi_{\mbox{\tiny BH}}|<1$ and $|\chi_{\mbox{\tiny NS}}|<0.05$ [72], and on the tidal parameter of $0<\Lambda<5000$ [73]666Such priors are constructed by translating the (non-Gaussian) values quoted above into the 68% confidence intervals of Gaussian probability distributions.. The Fisher analysis technique described here is only used an approximation to the more comprehensive (and computationally expensive) Bayesian analysis. However Ref. [74] showed that bounds on ppE parameters sufficiently agree between the two methods for GW150914 with SNR$\sim 24$. Thus, we expect the events considered in this analysis to give sufficiently valid order-of-magnitude estimations.

The upper and lower frequency integration limits used for each detector in the Fisher analyses are tabulated in Table 1. For ground-based detectors, the upper cutoff frequency of $f_{\text{ISCO}}=(6^{3/2}\pi M_{z})^{-1}$ with redshifted mass $M_{z}\equiv(1+z)M$ is used for each ground-based detector, with lower cutoff frequencies of 10 Hz, 10 Hz, 5 Hz, and 5 Hz for A+, Voyager, CE and ET respectively777Reference [53] mentions a starting frequency of 10 Hz for CE, yet the noise curve extends down to 5 Hz and is thus used here for completeness (such an extension is not made for A+ and Voyager, which also extend to $5$ Hz.). Additionally, very little SNR is accumulated between 5-10 Hz and thus the signal between $5-10$ Hz is expected to have a negligible impact on our results. This was indeed tested, and it was found that the difference was negligible to the accuracy of our results.. The redshift is computed in a cosmology such that $D_{L}=\frac{1+z}{H_{0}}\int^{z}_{0}\frac{dz^{\prime}}{(1-\Omega_{\Lambda})(1+z^{\prime})^{3}+\Omega_{\Lambda}}$ with Hubble constant $H_{0}=70$ km s*-1* Mpc*-1* and vacuum energy density $\Omega_{\Lambda}=0.67$. Following Ref. [24], the upper cutoff frequencies for space-based detectors is taken to be 100 Hz, with lower cutoff frequencies corresponding to one year prior to merger. We note that we use the redshifted chirp mass in the computation of such frequencies.

Additionally shown in Table 2 are the one-year BH-NS detection rates used in our analysis to combine uncertainties, for each detector considered. Such rates are computed following Ref. [25] as

[TABLE]

where $\Delta T=1$ yr is the observing time, $z_{h}$ is the horizon distance888The Horizon distance is defined as the luminosity distance such that an threshold SNR of $8$ is observed. redshift of each detector, $\mathcal{R}\in[0.6,610]\text{Gpc}^{-3}\text{yr}^{-1}$ (with a “realistic value” of $30Gpc^{-3}yr^{-1}$) is the local BH-NS coalescence rate [12, 45], and $a_{0}r(z)$, $d\tau/dz$, and $R(z)$ can all be found in Ref. [25]. Finally, the combined uncertainty on a parameter $\theta^{a}$ can be computed as

[TABLE]

where $\sigma_{\theta^{a}}(z)$ is the root-mean-square error on $\theta^{a}$ as a function of redshift, computed via Fisher analyses for increasing luminosity distances $D_{L}$. See Sec. 4.3 for an alternative analysis in which the BH and NS masses are varied for a more comprehensive combined uncertainty calculation.

3 Scalar-tensor theories of gravity

In this section, we discuss the two primary STTs of gravity presented in the following analysis. These include the quasi-Brans-Dicke theory DEF, as well as the EdGB theories of gravity.

3.1 Quasi-Brans-Dicke theories

Let us first focus on a class of scalar-tensor theories, quasi-Brans-Dicke theory. In this theory, matter fields couple to the scalar field $\varphi$ through the effective metric $A^{2}(\varphi)g_{\mu\nu}$ [19, 33, 32, 75]. One can then define the gradient and curvature of the conformal potential $\ln A(\varphi)$ to be $\alpha(\varphi)\equiv d\ln{A(\varphi)}/d\varphi$, and $\beta(\varphi)\equiv d\alpha/d\varphi$. In particular, we focus on the Damour and Esposito-Far$\acute{e}$se (DEF) [32, 31] theory999See A for a comparison with the similar Mendes-Ortiz (MO) [76] model., where the coupling function can be written in one of its simplest forms as $A(\varphi)=\exp{(\beta_{0}\varphi^{2}/2)}$. Such a theory can be completely characterized by the two weak-field parameters $\alpha_{0}=\alpha(\varphi_{0})=\beta_{0}\varphi_{0}$ and $\beta_{0}=\beta(\varphi_{0})$, where $\varphi_{0}=\alpha_{0}/\beta_{0}$ is the asymptotic value of the scalar field $\varphi$ at infinity 101010$\beta_{0}<0$ leads to cosmological runaway evolution of the scalar field that violates the current solar system bounds [77, 78, 79], unless one introduces a mass to the scalar field either directly or effectively by e.g. coupling the scalar field to an inflaton [80]..

Similarly in the strong-field case, NSs with mass $m_{A}$ couple to the scalar field with an effective coupling $\alpha_{A}=\partial\ln m_{A}/\partial\phi_{0}$, known as the (dimensionless) scalar charge 111111Scalar charges depend on the NSs underlying equation-of-state (EoS). In this analysis we assume the APR4 EoS, consistent with the binary NS observation GW170817 [73, 13]. See B for a comparison between results found with different EoSs. (the scalar charge for BHs is 0 [31]). Such scalar charges induce scalar dipole radiation in a compact binary, which enters at $-1$PN order relative to the leading GR quadrupole radiation and makes the binary evolve faster. Following Ref. [81, 19, 82], one can derive the corresponding ppE correction to the waveform to be

[TABLE]

where $\eta\equiv m_{1}m_{2}/M^{2}$ is the symmetric mass ratio, and $\Delta\alpha\equiv(\alpha_{1}-\alpha_{2})$ is the difference in scalar charges between orbiting compact objects. Additionally, see Sec. 4.1 for a discussion and comparison on the inclusion of higher-order PN corrections to the waveform phase, as well as to the amplitude. See also Ref. [83] for constraints from GW170817, and predictions for future binary NS detections.

Now let us discuss how one can constrain STTs with pulsar timing measurements. The first way to do this is through orbital decay rate $\dot{P}_{b}$ measurement. The dominant correction to orbital decay rate in STTs is from the dipole radiation of the scalar field $\dot{P}_{b}^{D}$. Thus, we constrain $\dot{P}_{b}^{D}/\dot{P}_{\text{GR}}$ by the fractional measurement accuracy of the orbital decay rate $\delta_{\dot{P}_{b}}$. In STTs, the expression for $\dot{P}_{b}^{D}/\dot{P}_{\text{GR}}$ is

[TABLE]

where $\Omega_{b}$ is the orbital frequency, $m_{2}$ is the pulsar’s companion mass, $q$ is the mass ratio $m_{1}/m_{2}$, and $e$ is the orbital eccentricity. Using the upper bound on dipole radiation in Eq. (4), we can place constraints on the theory parameters. The second way of testing STTs with pulsar timing is through constraints on the SEP. Recently, PSR J0337+1715 has placed the most stringent bounds on SEP violation [21]. This SEP violation measurement $\Delta$ constrains the theory parameter of STTs with the inequality $|\alpha_{\mbox{\tiny WD,out}}(\alpha_{\mbox{\tiny PSR}}-\alpha_{\mbox{\tiny WD,in}})|<\Delta$ where $\alpha_{\mbox{\tiny PSR}}$, $\alpha_{\mbox{\tiny WD,out}}$, and $\alpha_{\mbox{\tiny WD,in}}$ are the scalar charges of the pulsar, outer WD, and inner WD respectively. We update the constraints on STTs from Ref. [21] by using a softer EOS, because stiff EoSs are inconsistent with the recent GW observations [13, 84, 73].

We first discuss the present considerations of DEF constraints using GW and pulsar timing observations. Figure 2 presents the estimated constraints in the DEF theory parameter $\alpha_{0}-\beta_{0}$ plane for the various observations considered in this analysis. Observe that the combination of Cassini and pulsar timing measurements from PSR J0337 and PSR J1738 places the strongest constraints on DEF gravity. Moreover, even if the O3 candidates S190426c and S190814bv were BH-NS merger events [28], they struggle to place competitive bounds on DEF theory. Thus, this motivates why we must consider future bounds on DEF from GW measurements.

We conclude with an expedition into the future of GW astronomy. We consider the BH-NS system described previously, with fixed BH and NS masses of $10\text{ M}_{\odot}$ and $1.4\text{ M}_{\odot}$, respectively. We assume detections on the future GW interferometers A+, Voyager, CE, ET, B-DECIGO, and DECIGO, and following the spirit of Ref. [25] we combine the bounds on $\Delta\alpha$ from $N$ BH-NS detections falling within the horizon of each detector over one observing year, as described in Sec. 2. Further, we consider the multi-band observations [90, 42, 91] of such binaries between both ground-based detector ET and space-based detectors DECIGO/B-DECIGO. Unlike the multi-band case with space-based detector LISA, in which the possibly large SNR threshold of $\sim 9$ [92] would prevent one from obtaining event rates larger than $\mathcal{O}(1)$, multi-band detections between ET and DECIGO/B-DECIGO will instead be limited by the ET event rate because of the large SNRs obtained by the space-based telescopes (this reasoning primarily applies to DECIGO, with event rates approximately equivalent or greater than those on ET and CE, rather than B-DECIGO with significantly smaller rates). Such event rates are still significantly large at $\sim(500-500,000)$ (see Table 2).

Figure 3 presents the bounds in the DEF theory placed for the above-mentioned procedures. Observe how all of the current constraints can be improved upon with the optimistic number of detections on the A+ detector, while CE and ET begin to approach the same point with only the pessimistic number of detections. Further, all predicted bounds placed with DECIGO/B-DECIGO (single-event, multiple-event, and multi-band) improve the current constraints by several orders of magnitude. Of course, existing bounds from solar system experiments and binary pulsar observations will also improve in future. For example, bounds on $\alpha_{0}$ from Gaia will improve those from Cassini by a factor of a few [89], while the bounds from the pulsar triple system PSR J0337 will improve by a factor of $\sim 10$ with SKA [20]. Future GW bounds with 3rd generation detectors (ET/CE) and space-based detectors (B-DECIGO/DECIGO) are likely to be even stronger than them. We also note that the bounds for Brans-Dicke theory with $\beta_{0}=0$ for ET and DECIGO are consistent with those in [93, 25].

3.2 Einstein-dilaton Gauss-Bonnet gravity

We now show how the GW observations of BH-NS binaries can be applied to constrain another alternative scalar-tensor theory: EdGB gravity. In this string-inspired theory, the dilaton scalar field $\varphi$ non-minimally couples to a quadratic curvature term. In particular, we consider a linear coupling, where the Einstein-Hilbert action is corrected with a term $\alpha_{\mbox{\tiny EdGB}}\varphi R^{2}_{\mbox{\tiny GB}}$ [34, 35, 94]121212See Refs. [95, 96, 97] for general couplings.. Here, $\alpha_{\mbox{\tiny EdGB}}$ is the coupling parameter of the theory while $R_{\mbox{\tiny GB}}$ is the Gauss-Bonnet invariant. In this SEP-violating theory of gravity, BHs can accumulate scalar charges [36, 37, 25, 39] while ordinary stars like NSs do not [38, 40]. Similar to quasi-Brans-Dicke theory, such charges will induce a scalar dipole radiation in a binary involving at least one BH, which accelerates the rate of inspiral between gravitating bodies. Such an effect modifies the gravitational waveform phase at $-1$PN order, and is proportional to the coupling parameter $\alpha_{\mbox{\tiny EdGB}}$ of the theory, as well as the masses $m_{A}$ and sensitivities $s_{A}$ of the compact bodies [38] (see Ref [47] for the appropriate ppE expression). The spin-dependent sensitivities are non-zero for BHs only, and are taken to be 1 (spinless) for this analysis. The current constraints on EdGB gravity have been found to be $\sqrt{\alpha_{\mbox{\tiny EdGB}}}\lesssim 2$ km [98, 99, 100]131313We note that the small-coupling approximation $\zeta_{\mbox{\tiny EdGB}}\equiv 16\pi\alpha_{\mbox{\tiny EdGB}}^{2}/M^{4}\ll 1$ for binaries with total mass $M$ must be satisfied in order for valid constraints on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ to be placed.. We urge caution that the constraint $\sqrt{\alpha_{\mbox{\tiny EdGB}}}\lesssim 2$ found in Refs. [99, 100] take into account certain approximations which warrant such results to be a rough estimate. For example, the authors of Ref. [99] used posterior samples provided by the LIGO/Virgo Collaborations though data more finely-sampled around the GR value seem to be necessary to derive a more reliable posterior distribution on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$. In addition, the authors utilize $-1$PN constraints obtained by the LIGO-Virgo analysis of the binary black hole signals, where it is warned that they cannot be interpreted as dipole radiation constraints, due to higher-order non-negligible terms from dipole radiation. Also, Ref. [100] utilized a simplified grid-search method rather than a full stochastic sampling procedure.

We begin by discussing the current observational constraints on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$, had a BH-NS coalescence been observed by the current iteration of LIGO interferometers. Figure 4 projects the prospective constraints on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ for BH-NS binaries with $m_{\mbox{\tiny NS}}=1.4\text{ M}_{\odot}$ as a function of $m_{\mbox{\tiny BH}}$ for detection SNRs ranging between 8 and 20 on the aLIGO O3 detector. Observe how for BHs with mass less than $16.5\text{ M}_{\odot}$ ($19.5\text{ M}_{\odot}$), the current constraint in the literature of $2$ km can be improved upon for events with SNR $=8$ (20). Thus, if S190426c or S190814bv are NS-BH merger events with sufficiently low-mass BHs, such events would place a bound in EdGB gravity that is stronger than the existing bounds. Reference [101] estimated the properties of S190426c from the probability of the system being in specific categories, such as BH/NS. In particular, the BH mass is estimated as $\sim 6M_{\odot}$. If S190426c was indeed a BH/NS system (58% probability of terrestrial origin) and if this mass estimate was correct, we can place strong constraints on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ of $0.4$ ($0.5$) km for a 20 (8) SNR event – a factor of 4–5 improvement from the current observational constraint.

We follow this up with a discussion of future constraints placed on EdGB gravity. Similar to the previous section, we estimate the constraints placed on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ from a $10\text{ M}_{\odot}-1.4\text{ M}_{\odot}$ BH-NS merger event detected on each detector, in the single-event, multiple-event, and multi-band cases. Figure 5 displays the corresponding bounds for each scenario. Observe that the single-event rates can place constraints between $0.02-1$ km, all stronger than the current bound of $2$ km. Further, we see that the multi-band constraints do not offer much improvement from the single-band case, while the combined event bounds can reach down to $\sim 10^{-5}$ km with DECIGO, improving the current bounds by up to five orders-of-magnitude. These bounds with DECIGO are consistent with the rough estimate presented in [98] and a recent analysis of [91] for binary black holes with single events.

4 Validity of Approximations

In this section we explore the simplifying approximations made in the above analysis, and their effects on the presented results. We begin with a discussion on the number of higher-order PN corrections added to the gravitational waveform. In the main analysis, only the leading order $-1$PN correction term was taken into account in the gravitational waveform, and here we discuss the effect of including corrections up to $0.5$PN order. Next we consider the effect of rotating BHs, rather than the static ones considered in the analysis. Finally we conclude with an alternative method to combine multiple events with varying BH and NS masses, rather than the fixed masses considered in the main analysis.

4.1 Higher-order STT corrections to the waveform

Let us begin with a discussion of the higher-order PN corrections present in the gravitational waveform. We start by taking into account the $l=m=2$ dipolar contributions to the Fourier domain phase found in Eq. (81c) of Ref. [102]141414Note the non-dipolar phase corrections are not relative to the $-1$PN phase correction, proportional to $\Delta\alpha$, so were not included in this approximation for simplicity. up to $0.5$PN order. We make the approximation $\alpha_{\text{A}}\ll 1$ such that only terms proportional to $\Delta\alpha^{2}$ (and thus also proportional to $\beta_{\mbox{\tiny DEF}}$) remain, making it possible to keep correlations minimal by allowing only one non-GR parameter to remain in the waveform, $\beta_{\mbox{\tiny DEF}}$151515Note that the addition of new parameters to the gravitational waveform template may act to weaken the obtained constraints on $\Delta\alpha$ due to increased correlations between the parameters. A detailed analysis on the magnitude of this effect is beyond the scope of this analysis. The error due to linearization of the dipole radiation is expected to be minimal because $\beta_{\mbox{\tiny DEF}}\ll(\Delta\alpha)^{2}$.. The resulting DEF corrections to the gravitational waveform up to $0.5$PN order are given by

[TABLE]

where $\beta_{\mbox{\tiny DEF}}$ is the $-1$PN ppE parameter given in Eq. (3). We now include these additional corrections into the gravitational waveform and recompute bounds on the DEF theory for an SNR 10 event with a BH mass of $10\text{ M}_{\odot}$ on the LIGO O3 detector. We find constraints on $\Delta\alpha$ to be $\sim 0.104$ with leading-order $-1$PN corrections included, and $\sim 0.105$ (negligible difference on a plot such as Fig. 2) with higher-order corrections to $0.5$PN included. Such results have a $\sim 0.4$% difference, and are likely to be negligible in the presented analysis. A similar result was found in [74] for Brans-Dicke theory. We also bring to attention Ref. [103], where it was discovered that the additional presence of non-GR amplitude corrections to the waveform only differs from those in phase by $\sim 4$% - another negligible difference in our analysis.

4.2 Rotating BHs

Now let us discuss the effect of considering rotating BHs in our analysis, rather than the static ones considered previously. To perform this simple comparison, we compute constraints on $\Delta\alpha$ in the DEF theory of gravity, and $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ in the EdGB theory of gravity for a SNR 10 event with a BH mass of $10\text{ M}_{\odot}$ on the LIGO O3 detector, with the assumptions of BH spins $\chi_{\mbox{\tiny BH}}=(0,0.5,1)$. The resulting constraints were found to be $\Delta\alpha<(0.105,0.104,0.103)$ and $\sqrt{\alpha_{\mbox{\tiny EdGB}}}<(1.80\text{ km},1.79\text{ km},1.78\text{ km})$ respectively, with only a maximal difference of $\sim 1.5$% found between them in either case. Thus, we conclude that the effect of rotating BHs in our analysis is sufficiently negligible.

We also point out that the inclusion of spin precession does not play a crucial role in the order-of-magnitude constraint of $-1$PN effects. This is demonstrated in Table IV of Ref. [104], where it was shown that spin precession strengthens constraints on the Brans-Dicke parameter by $\sim 43$%, due to small correlations between the $-1$PN ppE parameter and the spins entering at $1.5$PN. Such an effect is magnified for increasing mass-ratio systems, much larger than those considered in this analysis, as they tend to increase the effects of precession.

4.3 BH/NS mass populations

In this section we model an appropriate BH-NS mass distribution function and implement it into the procedure used to combine the uncertainty in non-GR parameters from $N$ events, now with variational masses. To do this, we modify the expression given in Eq. (2) by injecting a mass-distribution function $f(m_{\mbox{\tiny BH}},m_{\mbox{\tiny NS}})$ like so

[TABLE]

where $\sigma_{\theta^{a}}(z,m_{\mbox{\tiny BH}},m_{\mbox{\tiny NS}})$ and $z_{h}(m_{\mbox{\tiny NS}},m_{\mbox{\tiny BH}})$ are now interpolated functions that also depend on the individual binary masses. For simplicity, we assume that $f(m_{\mbox{\tiny BH}},m_{\mbox{\tiny NS}})$ for BH/NS is simply given by a product of the individual mass distributions $f_{\mbox{\tiny BH}}(m_{\mbox{\tiny BH}})$ and $f_{\mbox{\tiny NS}}(m_{\mbox{\tiny NS}})$ as $f(m_{\mbox{\tiny BH}},m_{\mbox{\tiny NS}})=\mathcal{C}f_{\mbox{\tiny BH}}(m_{\mbox{\tiny BH}})f_{\mbox{\tiny NS}}(m_{\mbox{\tiny NS}})$, where the constant $\mathcal{C}$ is determined by normalizing the function to be unity when being integrated over $m_{\mbox{\tiny BH}}$ and $m_{\mbox{\tiny NS}}$. We use $f_{\mbox{\tiny BH}}(m_{\mbox{\tiny BH}})$ as the mass distribution of primary black holes in stellar-mass BH binaries derived by the LIGO/Virgo Collaborations [105], while we adopt $f_{\mbox{\tiny NS}}(m_{\mbox{\tiny NS}})$ as a Gaussian distribution used e.g. in [106]. Namely, we have

[TABLE]

The relevant parameters $\alpha$, $m_{\text{min}}$, $m_{\text{max}}$, $\mu_{\mbox{\tiny NS}}$, and $\sigma_{\mbox{\tiny NS}}$ have been fit to be $0.4$, $5\text{ M}_{\odot}$, $41.6\text{ M}_{\odot}$, $1.34\text{ M}_{\odot}$, and $0.06\text{ M}_{\odot}$ respectively. We perform a grid search with 8 redshift values between 0 and 8, 10 NS masses between $1\text{ M}_{\odot}$ and $2\text{ M}_{\odot}$ and 10 BH masses between $5\text{ M}_{\odot}$ and $40\text{ M}_{\odot}$ to compute $\sigma_{\sqrt{\alpha_{\mbox{\tiny EdGB}}}}(z,m_{\mbox{\tiny BH}},m_{\mbox{\tiny NS}})$ for 800 mass/redshift samples, which is then interpolated.

We perform this example computation for the EdGB theory of gravity on the CE detector with the “realistic” number of events ($730,000$) and compare the resulting constraint on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ to the case of fixed-mass binaries presented in the main analysis. Under the described circumstances, we find a constraint on $\sqrt{\alpha_{\mbox{\tiny EdGB}}}$ of $0.003$ km with the new variational mass model. Compared to the static-mass model result of $0.004$ km, we find that the two methods agree to within 25%. Interestingly, we find the new constraints to be stronger than the old ones as the new analysis includes BH masses lower than $10M_{\odot}$ considered originally, and thus the results displayed in the main text can be presented as a conservative estimate. Because the relationship between $\alpha_{0}$ and $\beta_{0}$ in DEF theory of gravity themselves depend on the constituent masses, this method can not be used to reliably compute bounds in the $(\alpha_{0},\beta_{0})$ plane. However, we expect to find similar results to the EdGB case (where $\alpha_{\mbox{\tiny EdGB}}^{2}$ doesn’t depend on the masses). This is verified by instead estimating constraints on $\Delta\alpha$ only, which agrees to the static-mass model to within 25% as well.

5 Conclusion

In this analysis, we demonstrated the present and future considerations on constraining STTs which violate the SEP. We considered both the DEF and EdGB theories, which predict massless scalar fields $\varphi$ which couple to matter and alter the consequent trajectories of gravitating bodies. We investigate constraints placed on these theories’ coupling parameter spaces for the possible detection of BH-NS coalescences, both on the current iteration of LIGO interferometers, and with future GW detectors both on the ground and in space. In the DEF theory, we find that if such an event (such as the possible candidates S190426c or S190814bv in the O3 run) were to be observed with the present GW detection capabilities, competing bounds to those from pulsar timing observations can be presented. In EdGB theory, we find that with BH masses less than $19.5\text{ M}_{\odot}$, improvements to the current constraint on the coupling parameter $\sqrt{\alpha_{\mbox{\tiny EdGB}}}<2$ km can be made to the order of $\mathcal{O}(0.1)$ km. Such events detected on future GW detectors (single-event, multi-band observations, and multiple-event stacking) have been demonstrated to improve upon the current bounds by several orders of magnitude in many cases.

Future work in this direction can improve upon this analysis by considering a Bayesian approach to parameter estimation, rather than the Fisher one considered here. Further, more accurate BH-NS population simulations other than those found in Ref. [45] may be utilized in future analyses, together with different masses for different events. Finally, one could consider a more comprehensive list of STTs to study, rather than the select few examples investigated here: DEF, MO (see A for a comparison between the two), and EdGB. One could even consider theories other than STTs, such as those involving vector fields and/or additional tensor fields.

6 Acknowledgments

We thank Takahiro Tanaka for carefully reading over the work. We also thank Benjamin Lackey for generously providing us tabulated data for the WFF1, APR4, and MPA1 EoSs used in  [107]. Z.C. and K.Y. acknowledges support from NSF Award PHY-1806776. K.Y. would like to also acknowledge support by the COST Action GWverse CA16104 and JSPS KAKENHI Grants No. JP17H06358. B.C.S. and K.Y. acknowledge support from the Mead Endowment.

Appendix A Quasi-Brans-Dicke theory comparison

In this appendix, we compare results between the DEF [108, 109] and MO [76] quasi-Brans-Dicke theories of gravity. The latter theory is defined by the coupling $\alpha(\varphi)=\text{tanh}(\sqrt{3}\beta_{0}\varphi)/\sqrt{3}$, while the former relies upon only the first term in the above expansion about $\varphi_{0}$, namely $\alpha(\varphi)=\beta_{0}\varphi$. Figure 6 compares the results for the PSR J0337 [85, 21] system from the SEP-violation test, assuming an APR4 EoS161616The BH-NS system constraints were found to be indistinguishable between theories, due to their lack of the horn structure present only in pulsar-WD binaries.. Observe how the “horn” structure171717The horns arise in pulsar-WD systems because certain values of $\alpha_{0}$ and $\beta_{0}$ suppresses the dipole term and deteriorating the constraints [17]. is more pronounced in MO theory, and the drop-off in the lower left region is shifted. Otherwise, the two theories predict nearly-identical values. This finding is consistent with that in [75] for the orbital decay rate measurement of pulsar-WD binaries.

Appendix B Equation of state comparison and the effects of spin in the waveform

In this appendix, we present a comparison between the assumption of different NS EoSs, as well as an investigation into spin effects in the gravitational waveform. Figure 7 compares bounds in the $\alpha_{0}-\beta_{0}$ plane for $10\text{ M}_{\odot}-1.4\text{ M}_{\odot}$ BH-NS system detected by ET, assuming three different EoSs: WFF1 [110, 111, 107], APR4 [112, 107], and MPA1 [113, 107]. Such EoSs were chosen to be consistent with the GW observation of binary NSs, GW170817 [73, 13, 114], and with increasing degrees of stiffness. We observe that, while the constraints do not differ much, the softer EoSs produce stronger bounds for small values of $\beta_{0}\lessapprox-3$, while the stiffer EoSs give stronger results for large values of $\beta_{0}\gtrapprox-3$. Thus, for consistency, we present results in the main text for the APR4 EoS.

We now consider the advisement of including spin effects in the gravitational waveform when computing constraints on quasi-Brans-Dicke theories. Such bounds were computed for binary NS systems found in Ref. [19] with a waveform template not including any spin effects. In our analysis, we utilize the PhenomD [56, 57] gravitational waveform which does indeed include spin effects. Figure 8 compares the constraints formed from the ET observation of a $10\text{ M}_{\odot}-1.4\text{ M}_{\odot}$ BH-NS system, both with and without spin effects included in the PhenomD waveform. We see that the latter under-estimates bounds on $\Delta\alpha$ by a factor of 2, indicating the necessity to include spin effects in the waveform. These discrepancies arise from the correlations between spin, and the other parameters in the waveform, in particular the non-GR parameter, which ultimately increases the uncertainties in parameter estimation.