Binary White Dwarfs as Laboratories for Extreme Gravity with LISA

TL;DR

This paper proposes using low-frequency gravitational waves from white dwarf binaries observed by LISA to test for deviations from General Relativity, specifically dipolar dissipation, and derives constraints on scalar-tensor theories.

Contribution

It introduces a model-independent null test for dipolar dissipation in gravitational waves from white dwarf binaries and applies it to scalar-tensor theories to set new constraints.

Findings

White dwarf binaries can constrain dipolar dissipation.

Projected constraints on scalarization are comparable or stronger than binary pulsars.

A framework for null tests of gravity with LISA data is developed.

Abstract

The observation of low-frequency gravitational waves with the Laser Interferometer Space Antenna will allow the study of new sources of gravitational radiation that are not accessible by ground-based instruments. Gravitational wave sources provide invaluable information both about their astrophysics, as well as the nature of the gravitational interaction in their neighborhoods. One low frequency source that has not received much attention regarding the latter are galactic binaries composed of two white dwarves or a white dwarf and a neutron star. We here show that, contrary to the common lore, such gravitational wave sources can indeed be used to constrain an important feature of the gravitational interaction: the absence of pre-Newtonian, dipolar dissipation. We propose a model-independent framework to carry out a null test for the presence of this feature in the data that is very much analogous to tests of General Relativity with radio-observations of binary pulsars. We then go one step further and specialize this test to scalar-tensor theories to derive projected constraints on spontaneous scalarization. We find that these constraints can be comparable to current bounds with binary pulsars, and in some optimistic cases, they can be even stronger.

Figures

Tables

Table 1. Table 1: Parameters of the ten WD-WD systems and the five WD-NS systems considered in our analysis. The parameters ( f 0 , f ˙ , f ¨ ) subscript 𝑓 0 ˙ 𝑓 ¨ 𝑓 (f_{0},\dot{f},\ddot{f}) describe the orbital evolution, while 𝒜 𝒜 \mathcal{A} is the dimensionless gravitational wave amplitude. The sky location is encoded in ( cos ⁡ θ , ϕ ) 𝜃 italic-ϕ (\cos\theta,\phi) where θ 𝜃 \theta is the co-latitude and ϕ italic-ϕ \phi is the longitude in ecliptic coordinates. The angular parameters that encode the binary’s orientation are ( cos ⁡ ι , ψ GW , φ 0 𝜄 subscript 𝜓 GW subscript 𝜑 0 \cos\iota,\psi_{\rm GW},\varphi_{0} ) where ι 𝜄 \iota is the inclination angle, ψ GW subscript 𝜓 GW \psi_{\rm GW} is the polarization angle, and φ 0 subscript 𝜑 0 \varphi_{0} is the orbital phase at t = 0 𝑡 0 t=0 . Quantities derived from the source parameters are ℳ ℳ \mathcal{M} and the signal to noise ratio S / N 𝑆 𝑁 S/N of the binary after one year of observing (which will grow roughly as T 𝑇 \sqrt{T} ). Dashees indicate the parameter value is identical to the row above. All sources are simulated to be consistent with GR, so f ¨ ¨ 𝑓 \ddot{f} can be derived from f 0 subscript 𝑓 0 f_{0} and f ˙ ˙ 𝑓 \dot{f} . The WDWD systems are drawn from a spatial distribution consistent with the Milky Way galaxy, with random orientations uniform on a sphere. The WDNS systems are chosen to be edge-on (therefore eclipsing, to maximize what can be learned from EM observations) simulated with 𝒜 𝒜 \mathcal{A} to be consistent with the source at a distance of 8 kpc.

WD-WD Parameters
Source	$f_0$ (s-1)	$\dot{f}$ (s-2)	$\ddot{f}$ (s-3)	$𝒜$	$\cos \theta$	$\varphi$ (rad)	$\cos \iota$	${\psi }_{\mathrm{GW}}$ (rad)	${\phi }_0$ (rad)	$ℳ$ ($M_{\odot }$)	S/N (1 yr)
0	0.015248	2.5834e-14	1.6049e-25	7.7243e-23	-0.0708	4.6600	-0.2016	0.0430	4.2766	0.39	61
1	0.008994	1.7922e-14	1.3094e-25	1.4888e-22	0.2779	4.8717	0.5001	1.2662	6.1245	0.99	186
2	0.012124	1.6886e-14	8.6230e-26	7.8345e-23	-0.5828	4.3243	0.5973	0.1514	3.1676	0.50	108
3	0.012163	1.5602e-14	7.3380e-26	1.0910e-22	0.0909	4.7667	0.2371	1.9003	2.5800	0.47	97
4	0.018509	4.4578e-14	3.9367e-25	6.4157e-23	-0.5394	4.3525	0.7049	0.4413	2.2722	0.35	95
5	0.007984	7.9289e-15	2.8873e-26	4.7292e-22	-0.7201	4.0816	-0.4880	0.8372	2.2434	0.79	560
6	0.015234	1.3389e-14	4.3146e-26	8.3676e-23	0.0719	4.6291	0.2790	2.1810	4.3866	0.26	77
7	0.015144	2.0886e-14	1.0562e-25	1.0816e-22	-0.7246	3.8508	0.4231	2.0689	3.1987	0.35	113
8	0.017695	1.9882e-13	8.1911e-24	3.6914e-22	0.7288	5.4111	-0.3325	1.9704	3.6658	0.95	329
9	0.012800	1.5497e-14	6.8793e-26	7.5127e-23	-0.3924	4.4753	0.4679	1.3264	0.4712	0.42	89
WD-NS Parameters
Source	$f_0$ (s-1)	$\dot{f}$ (s-2)	$\ddot{f}$ (s-3)	$𝒜$	$\cos \theta$	$\varphi$ (rad)	$\cos \iota$	${\psi }_{\mathrm{GW}}$ (rad)	${\phi }_0$ (rad)	$ℳ$ ($M_{\odot }$)	S/N (1 yr)
$1.4+0.3M_{\odot }$	0.003000	1.1422e-16	1.5945e-29	5.4222e-23	-0.1494	4.600	0	0	0	0.53	36
–	0.005000	7.4332e-16	4.0519e-28	7.6221e-23	–	–	–	–	–	–	88
–	0.010000	9.4396e-15	3.2672e-26	1.2099e-22	–	–	–	–	–	–	91
–	0.020000	1.1988e-13	2.6345e-24	1.9207e-22	–	–	–	–	–	–	214
$1.4+1.3M_{\odot }$	0.010000	3.5059e-14	4.5069e-25	4.4938e-22	–	–	–	–	–	1.17	753

Equations

$${\dot{f}}_{\text{mbox}GR}$$

$${\ddot{f}}_{\text{mbox}GR}$$

$$\dot{f}$$

$$\ddot{f}$$

$$p(\vec{\lambda }|d)\approx \frac{1}{(2\pi \mathbf{\Sigma }{)}^{1/2}}{e}^{-\frac{1}{2}(\vec{\lambda }-\vec{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\vec{\lambda }-\vec{\mu })}$$

$$S$$

$$+S_M\left[{\chi }_M,A^2(\psi )g_{\mu \nu }^{\ast }\right],$$

$$g_{\mu \nu }=A^2(\psi )g_{\mu \nu }^{\ast }.$$

$$A(\psi )=\exp \left[{\alpha }_0\left(\psi -{\psi }_0\right)+\frac{1}{2}{\beta }_0{\left(\psi -{\psi }_0\right)}^2\right],$$

$$G_{\mu \nu }^{\ast }$$

$${\square }_{\ast }\psi$$

$$T_{\ast }^{\psi }{}_{\mu \nu }$$

$$\left|{\gamma }_{\text{mbox}ppN}-1\right|$$

$$\left|{\beta }_{\text{mbox}ppN}-1\right|$$

$${\alpha }_A^{\text{mbox}WD}={\alpha }_0\left[1-2C_A+\mathcal{O}(C_A^2)\right],$$

$$\ddot{f}$$

$$b_{\text{mbox}ppE}$$

$${\beta }_{\text{mbox}ppE}$$

$$\kappa$$

$$+\frac{1}{6}\left({\alpha }_1-{\alpha }_2\right)\left({\alpha }_1\frac{m_1}{m}+{\alpha }_2\frac{m_2}{m}\right)\left(\frac{m_1}{m}-\frac{m_2}{m}\right)$$

$$+\frac{5}{48}\frac{{\alpha }_1-{\alpha }_2}{1+{\alpha }_1{\alpha }_2}\left({\beta }_2{\alpha }_1\frac{m_1}{m}-{\beta }_1{\alpha }_2\frac{m_2}{m}\right)+d_2{\left({\alpha }_1-{\alpha }_2\right)}^2.$$

$${\mathcal{M}}_{\ast }=\frac{{\kappa }^{3/5}}{{\left(1+{\alpha }_1{\alpha }_2\right)}^{2/5}}\mathcal{M},$$

Full text

I Introduction

Sooner than all but the most hard-line optimists would have dared to imagine, gravitational waves (GWs) were observed by the advanced LIGO detectors in August 2015 Abbott et al. (2016a); Abramovici et al. (1992); Abbott et al. (2009); Harry (2010); Abramovici et al. (1992), forever changing physics and astronomy. That momentous discovery single-handedly confirmed the reality of GWs and their direct interaction with matter, an issue that although taken for granted by the relativity community, was still contentious in some astrophysical circles. And at that point, the flood gates opened. Discovery after discovery was heralded with great enthusiasm, leading to the first triple-coincident discovery when the advanced Virgo detector joined the network Caron et al. (1997); Giazotto (1990); Acernese et al. (2015); vir ; Abbott et al. (2017a), and later the first coincident discovery when light from across the electromagnetic spectrum, from gamma ray telescopes in orbit to radio telescopes on the ground, was detected in an unprecedented, coordinated observing campaign Abbott et al. (2017b, c); Tanvir et al. (2017).

The era of GW astronomy has thus begun, and with it, great discoveries and some surprises have trickled in from the interpretation of the observations made. These discoveries can be grossly classified into two separate categories: astrophysics and fundamental physics. On the astrophysics side, we have learned that the most common sources of GWs are the coalescence of black holes in quasi-circular inspirals, and perhaps the biggest surprise is the black hole masses detected. Before the GW era, the heaviest stellar-mass black holes observed were less than ${\sim}16M_{\odot}$ Orosz et al. (2007), but several of the GW observations are consistent with relatively heavy, stellar-mass black holes merging into objects of above $40M_{\odot}$ Abbott et al. (2016a, b). Naturally, this surprise had important implications on the astrophysics of black hole formation, restricting certain models and supporting others Abbott et al. (2016c).

On the fundamental physics side, the GW era has done much more than simply confirm the existence of GWs: it has allowed us to explore, for the first time, the extreme gravity regime, where the gravitational interaction is phenomenally large and dynamically changing, matter transcends into new forms not found on Earth, and physics is intrinsically non-linear and dynamical Yunes and Siemens (2013). The latest GW observation of merging NSs has allowed for the first stringent constraints on the mass and radius of NSs, and thus on the equation of state of supranuclear matter Annala et al. (2018); Radice et al. (2018); Abbott et al. (2018); Chatziioannou et al. (2018); Yagi and Yunes (2016, 2013a, 2013b). These constraints, in turn, have important implications on the theory of nuclear interactions at high densities, providing support for certain approximate methods to solve the many body Schrödinger equation Chatziioannou et al. (2015). Moreover, this same observation, together with electromagnetic counterparts, has allowed for the first constraints on the speed of GWs, and thus, the speed of gravity Abbott et al. (2017b, c); Yunes et al. (2016). These constraints, in turn, have dramatic implications on theories that modify General Relativity (GR) to explain late-time cosmological observations or to unify GR with quantum mechanics Yunes et al. (2016).

These first advanced LIGO/Virgo discoveries, however, are only the tip of the iceberg. Current ground-based detectors are expected to be upgraded to third-generation designs, allowing for detections from effectively everywhere in the universe where stellar-mass or intermediate-mass black hole collisions occur. But even with this technology, such detectors will be limited to observing GWs with frequencies above a few Hz due to inescapable Newtonian noise. The observation of low frequency GWs will therefore require different measurement techniques. Access to the low-frequency band between ${\sim}0.1$ and ${\sim 10}$ mHz, where the richest population of sources is expected, requires a space-based GW observatory: enter the Laser Interferometer Space Antenna Prince (2003); Audley et al. (2017), or LISA for short.

By opening up the low-frequency band, LISA will allow for a plethora of new astrophysical and fundamental physics discoveries because it will give us access to entirely different classes of sources: supermassive black hole coalescences, extreme mass-ratio inspirals, and ultra compact binaries, typically comprised of white-dwarf/white-dwarf or neutron-star/white-dwarf pairs, in our galaxy. The observation of supermassive black holes will reveal invaluable information about structure formation in the Universe, including how galaxies collide and form, as well as how black hole remnants relax to a stationary configuration after a major collision with another black hole Gair et al. (2013). GWs produced as a small compact object zooms and whirls into a supermassive black hole, a so-called extreme mass-ratio inspiral, will allow for a geodesic mapping of the spacetime created by the supermassive object, allowing for tests of the Kerr hypothesis Gair et al. (2013); Vigeland et al. (2011); Gair and Yunes (2011).

But what about GWs emitted by galactic binaries? Binaries composed of two white dwarves (WDs) or a WD and a NS in the Milky Way are known to exist because they are observed electromagnetically. Several of these known sources are guaranteed to be detectable for LISA, and are thus referred to as verification binaries Stroeer and Vecchio (2006); Korol et al. (2017); Kupfer et al. (2018). Moreover, given such observations and using population synthesis simulations, it is predicted that tens of millions of such binaries are emitting GWs in the LISA band. Tens of thousands of these will be individually detectable by LISA, while the rest blend together to form a confusion-limited foreground that is the dominant source of noise in the LISA band between $0.4$ mHz and $2$ mHz Benacquista (2006); Cornish and Robson (2017). The catalog of resolvable binaries, as well as the spectral shape of the confusion noise, will contain a trove of new information about the formation and evolution of compact binaries in the galaxy, and the complicated physical interactions between binaries including mass transfer, common envelop interaction, and the internal physics of compact stellar remnants Takahashi and Seto (2002); Crowder and Cornish (2004); Barack and Cutler (2004); Umstatter et al. (2005a, b); Cornish and Porter (2005); Rubbo et al. (2006); Kopparapu and Tohline (2007); Benacquista (2006); Benacquista et al. (2007); Littenberg (2011); Cornish and Robson (2017).

But could such sources also be used to learn something new about gravity itself? The expected answer to this question is a rotund no and the lore is cemented on (at least) the following three arguments.

Theoretical Physics Argument: GW-emitting galactic binaries detectable by LISA will be well-separated, with distances of roughly $10^{5}$ km for typical WD masses at GW frequencies of $(10^{-2},10^{-3})$ Hz. Such large separations imply orbital velocities smaller than $10^{-2}c$, which is certainly in the regime of velocities accurately probed by binary pulsar observations. Because binary pulsars have verified GR to amazing accuracy, modifications to gravity are expected to only be important (if they arise at all) at much higher velocities or smaller orbital separations. If so, then there is no point in using galactic binary signals to test GR since these will be probing a regime of gravity that has already been tested to a great degree by binary pulsar observations.

Data Analysis Argument: Such small separations imply orbital periods of tens of minutes, which in turn imply an extremely slow chirping rate of $\dot{f}<10^{-15}\;{\rm{Hz}}^{2}$. This is why galactic binaries are sometimes said to be monochromatic sources of GWs, emitting mostly at a single frequency during the observation period. Thus, such GW signals contain less accessible information than the familiar merger signals from ground-based observatories, since the former cannot be used to measure higher frequency derivatives, which encode information about the mass ratio and the spin of the binary components. The limited information contained in galactic binary signals implies that they would lead to suboptimal tests of GR.

Astrophysical Argument: Consider the effect of astrophysical processes on the orbital evolution of galactic binaries. WDs are considerably more deformable than NSs due to their weaker gravitational compactness. As such, WDs in a binary can easily be deformed if the orbital separation is small enough. What is worse, if the separation is even smaller, a WD in a binary will transfer its mass to a more compact companion, once its Roche Lobe is filled, or there could be mass transfer through winds. All of these astrophysical effects can contaminate or even control the orbital motion, and thus, the GW emission.

Putting these three arguments together (theoretical, data analysis and astrophysical) has led the field to conclude that galactic binaries would not be good laboratories for extreme gravity experiments, but is this really the case? Let us show through a simple two-fold argument that there is important information contained in the GWs emitted by galactic binaries that could be used to gain new knowledge about gravity.

Our first counter-argument is related to astrophysics. It is definitely true that WDs can be deformed, and their orbital evolution can be contaminated by mass transfer and/or tidal effects, as presented in the astrophysical argument above. However, higher mass WD binaries are more compact, and they can reach smaller orbital separations before interactions between the binary components become dominant in the evolution. High mass binaries are also among the loudest LISA sources anywhere in the galaxy, making them excellent candidates for detailed parameter estimation studies. While a small subset of the total number of detectable binaries in the galaxy, the sheer volume of sources suggests that there will be a few exceptional systems in the high mass tail of the population that can be exploited as laboratories to test gravity.

Our second counter-argument is related to theoretical physics. Again, it is definitely true that extreme gravity modifications to GR are expected to be dominant in the late inspiral, and thus, be suppressed in the low velocity regime, as presented in the theoretical physics argument above. However, there are some (very popular) modifications to GR that do exactly the opposite: they are dominant at low velocities. Such modifications are sometimes called “pre-Newtonian,” because they enter as negative powers of the velocity relative to the leading-order GR expression. A popular example is dipole radiation through the activation of a scalar or vector field, as in the case of scalar-tensor theories Damour and Esposito-Farese (1992, 1993), quadratic gravities theories, like dynamical Chern-Simons gravity Alexander and Yunes (2009); Jackiw and Pi (2003); Campbell et al. (1991, 1990); Yunes and Pretorius (2009a); Yagi et al. (2012, 2012, 2013); Stein and Yagi (2014) or Einstein-dilaton-Gauss-Bonnet gravity Campbell et al. (1992); Kanti et al. (1996); Alexeev and Pomazanov (1997); Torii et al. (1997); Yunes and Stein (2011); Ayzenberg and Yunes (2014); Stein and Yagi (2014); Yagi et al. (2016), and Lorentz-violating theories, like Einstein-Æther theory Jacobson and Mattingly (2001); Eling et al. (2004, 2007); Yagi et al. (2014a, b) and Horava gravity Horava (2009); Blas et al. (2010); Yagi et al. (2014a, b).

The emission of dipole radiation through the activation of new degrees of freedom forces the binary to decay, and therefore the GWs to chirp, faster than in GR. For example, the rate of change of the GW frequency in GR scales as $v^{-11}\mu/m^{3}$, where $v$ is the orbital velocity, and $\mu$ and $m$ are the reduced and total masses of the binary. But in theories with dipole radiation, the chirping rate acquires a correction that scales as $v^{-13}\mu/m^{3}$, which is thus dominant when $v\ll 1$ and invalidates the theoretical physics argument above. The observational implication for galactic binaries is thus clear: if pre-Newtonian dissipation is present, the GW chirping rate will be at least four orders of magnitude larger than the GR prediction, leading to non-monochromatic galactic binary signals in the LISA band. The absence of such an effect should thus allow us to place stringent constraints on the existence of pre-Newtonian modifications to GR.

With the theoretical physics and the astrophysical arguments counter-acted, we are then left only with the data analysis argument, which we recall states that there is simply not enough information in a galactic binary GW signal to carry out a test of GR. Whether this is true or not requires a data analysis study, as it is hard to assess otherwise. We will here show through such a study that, in fact, the data analysis argument is unfounded and tests of GR are indeed possible. Such a study allows us to develop both a model-independent and a model-dependent framework to test GR.

The model-independent framework is a null test of GR and, in fact, it is very much analogous to complementary tests that have been carried out with radio-observations of binary pulsars. The main idea is to extract the GW signal with a model parameterized only by its overall amplitude, frequency, first and second frequency derivatives, while remaining agnostic about the particular source that produced it. In particular, this means we impose no prior between the frequency derivatives and the mass of the system or frequency of the signal when we do parameter estimation. The best fit model parameters can then be used to draw curves in the mass-frequency plane assuming GR is correct, and the test is passed if all curves intersect at the same location.

For such a model-independent test to work, one must first detect at least one galactic binary that is sufficiently bound so that its second derivative is measurable, while at the same time not so bound that its evolution is controlled by mass transfer. From a galaxy simulation consistent with those used for past LISA performance studies and mock data challenges Nelemans (2003), we select the ten highest frequency detached (i.e. non-interacting) binaries to study the LISA measurement of orbital evolution, and how those results can be used to study gravity theories.

Figure 1 presents an example of this null test in the chirp mass–frequency plane, where the former is a particular combination of the component masses of the system. The measured GW frequency is a vertical curve in this plane (black in the figure), whose thickness represents its $90\%$ credible interval. Assuming the null hypothesis that GR is correct, the inferred value of the first frequency derivative can be used to draw another curve in this plane (orange in the figure), where again its thickness is its $90\%$ credible region. The intersection of these two curves gives the inferred chirp mass of the system assuming the null hypothesis. Any additional information is thus redundant and can be used to test the null hypothesis. For example, the measured value of the second frequency derivative can be used to draw another curve in this plane (green in the figure), and if the null hypothesis is correct, then all curves intersect in the same region.

In addition to GR consistency tests, galactic binaries can also be used to test specific models that predict the existence of pre-Newtonian dissipation. One such model is massless scalar-tensor theory, perhaps one of the simplest modifications to GR Brans and Dicke (1961); Damour and Esposito-Farese (1992, 1993, 1996, 1998). In this theory, a massless scalar field couples to the metric tensor, thus affecting the structure of massive objects and their motion in strong gravitational fields. In one sub-class of such models Damour and Esposito-Farese (1992, 1993, 1996, 1998), Solar System observations can be easily passed because the scalar field is shielded and does not activate in weak gravitational fields. In strong gravitational fields, such as those produced by NSs and WDs, the scalar field does activate, producing pre-Newtonian dissipation effects Will and Zaglauer (1989); Damour and Esposito-Farese (1998). Therefore, constraints on pre-Newtonian dissipation could be used to place bounds on the two coupling constants of this class of theories $(\alpha_{0},\beta_{0})$.

We here study the strength of such bounds using GWs emitted by galactic binaries of the mixed WD-NS type. The pre-Newtonian dissipation in scalar-tensor theories does not only depend on $(\alpha_{0},\beta_{0})$, but also on the response of the star to the presence of a scalar field, quantified in terms of certain sensitivities, which in turn depend on the compactness of the stars. Such tests, therefore, require first the independent measurement of the masses of the stars, which could be obtained electromagnetically if the binaries are eclipsing Parsons et al. (2017). These masses can then be used to infer the radius of the stars, and thus their compactness, using knowledge of the stellar equation of state, which for neutron stars we assume will be known by the time LISA flies through a combination of advanced LIGO/Virgo observations Annala et al. (2018); Radice et al. (2018); Abbott et al. (2018); Chatziioannou et al. (2018); Yagi and Yunes (2016, 2013a, 2013b) and X-ray observations Miller and Lamb (2015); Ozel et al. (2016).

With these assumptions, a GW measurement of $(f_{0},\dot{f},\ddot{f})$ implies a constraint of an open region in the $(\alpha_{0},\beta_{0})$ plane. Such constraints would be very much analogous to those placed with radio-observations of binary pulsars Stairs (2003); Freire et al. (2012). We indeed find that the constraints LISA could place in the future will be complementary, comparable and sometimes better than what binary pulsars can do today. Figure 2 shows the constraints that LISA could place assuming different accuracies for the measurement of the masses and radii of the component stars. Observe that the $(f_{0},\dot{f},\ddot{f})$ constraints are comparable to those that can be placed with radio-observations of binary pulsars, and they could be better by as much as an order of magnitude if the masses and radii were known sufficiently accurately.

The study presented here is by no means intended as a final word on tests of GR with LISA observations of GWs emitted by galactic binaries; rather, it hopefully opens the door to future studies of the extraction of theoretical physics from a guaranteed LISA source that has not received much attention in the GW experimental relativity community. Both the model-independent and the model-dependent frameworks ought to be more thoroughly investigated to account for systematics that could arise from astrophysical effects, like mass transfer and winds. The latter should also be extended to other modified theories of gravity that also predict the existence of pre-Newtonian dissipation effects. Nonetheless, the work presented here will hopefully act as a igniter that will stimulate such future studies.

The remainder of this paper presents the details of the calculations that lead to the results described above. Section II describes how to model WD binaries outside of GR. Section III discusses the data analysis techniques we will use to determine the magnitude of projected constraints with LISA. Section IV presents the projected constraints we can place on model-independent deformations of the GR model prediction. Section V focuses on scalar-tensor theories and presents model-dependent constraints. Section VI briefly discusses caveats that could modify the conclusions presented earlier quantitatively but not qualitatively. Section VII summarizes the main results of this work and concludes by discussing possible future work. Throughout this paper, we use the following conventions: $G=c=1$, $(\mu,\nu,\rho,...)$ are spacetime indices.

II Galactic Binary Modeling with pre-Newtonian Dissipation

This section describes how to model GWs emitted by galactic binaries appropriate for LISA science. We then conclude with a discussion of astrophysical populations of galactic binaries

II.1 Modeling Individual Sources

We model the LISA response to GWs from galactic binaries using the fast-slow decomposition described in Cornish and Littenberg (2007). The waveform model is extended to include the second time derivative of the frequency, modifying the GW phase $\Psi(t)=2\pi f_{0}t+\pi\dot{f}t^{2}+\frac{1}{2}\pi\ddot{f}^{3}+\varphi_{0}$ where $f_{0}$ and $\varphi_{0}$ are the GW frequency and phase measured at some fiducial time, such as when observations begin, and the over dots denote time derivatives evaluated at the fiducial time.

In the waveform model, the parameters $\vec{\lambda}=\{f_{0},\dot{f},\ddot{f},\ldots\}$ are independent, but in any given theory of gravity, there are inter-relations between them. For example, in GR, we have that

[TABLE]

to leading order in an expansion about ${\cal{M}}f\ll 1$, where ${\cal{M}}=\eta^{3/5}m$ is the chirp mass, $\eta=m_{1}m_{2}/m^{2}$ is the symmetric mass ratio and $m=m_{1}+m_{2}$ is the total mass. Therefore, assuming the system dynamics are dominated by GW emission, and assuming the null hypothesis that GR is correct, a measurement of $(f_{0},\dot{f})$ suffices to extract the chirp mass ${\cal{M}}$. The measurement of any higher derivative can then serve as a consistency check of the null hypothesis.

When pre-Newtonian dissipation effects are present, the relation between the components of $\vec{\lambda}$ changes. Generically, one can show that the chirping rate now becomes Yunes and Pretorius (2009b)

[TABLE]

which then implies that

[TABLE]

where $\beta_{{\mbox{\tiny ppE}}}\in\Re$ and $b_{{\mbox{\tiny ppE}}}\in\Re$ are parameterized post-Einsteinian amplitude and exponent parameters Yunes and Pretorius (2009b). We have here neglected terms proportional to $\beta_{{\mbox{\tiny ppE}}}^{2}(\pi{\cal{M}}f)^{2b_{{\mbox{\tiny ppE}}}}$ since we assume that $\beta_{{\mbox{\tiny ppE}}}\ll(\pi{\cal{M}}f)^{-b_{{\mbox{\tiny ppE}}}}$. The amplitude parameter $\beta_{{\mbox{\tiny ppE}}}$ can in principle depend on other dimensionless system parameters, like the symmetric mass ratio. The exponent coefficient $b_{{\mbox{\tiny ppE}}}<0$ in order for the effect to be pre-Newtonian, with $b_{{\mbox{\tiny ppE}}}=-2/3$ representing the dipole effect we will focus on in this paper.

II.2 Astrophysical Population

Of the tens of thousands of galactic binaries expected to be detectable in the LISA band, population synthesis simulations suggest that several such systems will be exceptionally strong sources. These binaries will be detectable within the first few weeks of mission operations, and over years of observing will accumulate signal to noise ratios (S/N) in the hundreds. For testing GR with galactic binaries, many of the loudest LISA sources are the best laboratories not only because of their high S/N, but also because these tend to be higher-mass WD binaries (with $\mathcal{M}\sim 0.4\,\text{M}_{\odot}$, corresponding to component masses of $\sim 0.5\,\text{M}_{\odot}$ for equal-mass binaries). As a consequence, the individual stars will be more compact, suppressing contributions to the orbital evolution of the binary through gravitational interactions between the stars (e.g., tidal interactions or mass transfer).

For our investigations we will select “loud,” non-interacting binaries from a galaxy realization consistent with those used for the LISA Data Challenges111https://gitlab.in2p3.fr/stas/MLDC constructed from simulations by Toonen et al. (2012). From the galaxy simulation, binaries are ranked by S/N and the 10 loudest systems with $f_{0}>4$ mHz are selected for study. There are high S/N binaries at lower frequencies, but they will more likely be contending with overlapping signals, time-varying galactic confusion noise, and, for lower mass systems, dynamical influences other than gravity due to mass transfer, tidal interactions, etc. In general, exploiting galactic binaries for testing orbital dynamics will be most straightforward using isolated systems at high frequency. The binary parameters for the systems we will study are presented in Table 1.

In addition to the astrophysical population of WD binaries, we will also investigate what can be achieved with LISA should we observe systems composed of a WD and a NS. In our simulations, the NS-WD binaries are placed near the galactic center (where the source density of galactic binaries is likely the highest), and we study five canonical systems, varying the observing time, masses, and the orbital frequencies of the binary. The binary parameters for these systems are also included in Table 1.

III Data Analysis

III.1 LISA

LISA will be a constellation of three free-flying spacecraft, each on independent heliocentric orbits approximately one AU from the Sun. Each satellite is equipped with an optical metrology system to precisely monitor the distances between free-falling test masses housed within the spacecraft. From the inter-spacecraft ranging measurements, two Michelson-like interferometry signals are synthesized digitally to simultaneously measure the two GW polarizations, while a third Sagnac channel is insensitive to GWs (at low frequencies) and will be used for detector characterization. The simulated detector performance used in this study is consistent with the instrument noise levels quoted in the LISA L3 Mission Proposal Audley et al. (2017). This analysis takes advantage of the noise-orthogonal $A$ and $E$ data streams constructed from linear combinations of the Michelson-like time delay interferometry channels Tinto et al. (2002). The $A$ and $E$ data streams assume identical noise in each of the six interferometer links of the spacecraft, and so are an idealized case that will not be realized during mission operations, but are suitable for this proof-of-concept study.

III.2 Methodology (Bayes/MCMC)

The simulated galactic binary sources are analyzed using a parallel tempered Markov Chain Monte Carlo (MCMC) pipeline adapted from Littenberg (2011) specifically designed for galactic binary searches in LISA data. The MCMC pipeline has been modified to include the second time derivative of the frequency as a parameter. We assume uniform priors on the frequency parameters. Mapping the gravitational-wave observables to the modified theory parameters is done through post-processing of the posterior samples from the Markov chain, the specifics of which will be described in later sections.

III.3 Example Parameter Estimation Results

In this work we will be mapping the parameters used by the MCMC sampler to some other physical quantities and evaluating the posterior distribution function at arbitrary parameter values. The latter is generally a challenging problem as the MCMC yields a discrete set of samples drawn from the distribution, rather than the continuous distribution itself. Here we are spared much of this difficulty because our investigations involve only extraordinarily high S/N sources, and the high S/N limit is where the errors are well approximated by Gaussians distributions.

To demonstrate this, we show 1D histograms and 2D scatter plots of the samples from the MCMC for a standard analysis of Source 0 (i.e., not including $\ddot{f}$ as a parameter), assuming a three year observation time. The parameters shown in Fig. 3 are the GW amplitude $\mathcal{A}$, the GW frequency $f_{0}$, and its first time derivative $\dot{f}$. The posterior distribution function $p(\vec{\lambda}|d)$ is then approximated by a multivariate Gaussian

[TABLE]

centered at the mean of the parameters $\vec{\mu}$ with covariance matrix ${\bf\Sigma}$ calculated directly from the chain samples, $\Sigma_{ij}=\langle\lambda_{i}\lambda_{j}\rangle$. The green curves in Fig. 3 show the Gaussian representation of the 1D posteriors and the $2\sigma$ contours of the covariance matrix computed from the MCMC samples.

Figure 4 shows an example of how the MCMC samples can be mapped to different physical quantities. In particular, we here show the canonical example of mapping $(\mathcal{A},f_{0},\dot{f})$ to the luminosity distance and chirp mass $(D_{L},\mathcal{M})$. At every value of $(\mathcal{A},f_{0},\dot{f})$ sampled by the chains, we compute the corresponding $(D_{L},\mathcal{M})$, and then construct a histogram of the result (gray in the figure). We then use the Gaussian approximation to approximate the posterior distributions for $(D_{L},\mathcal{M})$ using the covariance matrix of the reparameterized samples (orange in the figure). Once more, observe how the Gaussian estimate is an excellent approximation to the sampled posteriors. Beyond confirming by eye the agreement, we test that the p-values of the MCMC samples evaluated with the approximate posterior are consistent with being uniformly distributed between $[0,1]$.

The most important parameter for our study is the second time derivative of the frequency, $\ddot{f}$, as later in this work it will either be used to place an independent constraint on the chirp mass, or to break degeneracies with the pre-Newtonian dipole term. The $\ddot{f}$ parameter is often neglected in GW studies of WD binaries because it is seldom measurable. Moreover, this parameter provides redundant (though less constraining) information about $\mathcal{M}$ and $D_{L}$, which are the two parameters typically most sought after when thinking about exploiting the astrophsical predictions of the frequency evolution in galactic binaries.

For this work, the story of how well WD binaries can be used to test gravity is the story of how well $\ddot{f}$ can be measured. To that end, we focus on golden galactic binaries, the best-of-the-best in a simulated galactic population: high mass, short orbital period, systems which will undergo the most significant frequency evolution during the LISA observations, and will be easily detectable anywhere in the galaxy. Even with these extraordinary sources, long integration times are needed for the observations to constrain $\ddot{f}$. Figure 5 shows the $2\sigma$ contours in the $\dot{f}\text{--}\ddot{f}$ plane from posterior distribution functions represented by the MCMCs covariance matrix. On display here is the same source as shown in Figs. 3 and  4, now including $\ddot{f}$ as one of the model parameters. This figure demonstrates how $\ddot{f}$ improves with observing time. During the first 3 years of observing the second derivative of $f_{0}$ is basically unmeasured. After 5 years there are hints of a detection, as the parameter is constrained to be non-zero. For this simulated source, a 10 year observation is required to unambiguously measure a non-zero $\ddot{f}$.

IV Projected Model-independent Constraints

In this section we present the details of the model-independent test. We begin by describing the injection and the template that we used. We then proceed with the details of the test and present results from our data analysis studies.

IV.1 GR injection and dipole template

As a first demonstration of using WD binaries to test GR we take a model-agnostic approach, instead constraining the magnitude of the dipole effect on the waveforms. To do so, we analyze the ten binaries selected from the simulated galaxy as the most promising sources for these tests using the MCMC pipeline with $\ddot{f}$ included as a parameter. The binaries are simulated with $\beta_{\rm ppE}=0$, i.e. consistent with GR, and are then recovered using $f_{0},\dot{f}$ and $\ddot{f}$ as independent parameters, instead of being constrained to follow the GR relationship. For the remainder of this section, we will drop the ppE subscript on $\beta$ for convenience. The marginalized posterior $p(f_{0},\dot{f},\ddot{f}|d)$ is then mapped to $p(\mathcal{M},\beta|d)$ by inverting Eqs. (3) and (4).

Figure 6 shows cumulative distribution functions of the marginalized posterior on $\beta$, $p(\beta|d)$ for each source as a function of observation length. Horizontal dashed lines denote the $5\%$ and $95\%$ percentiles, enclosing the $90\%$ credible interval of the measurement. Recall that the $\beta$ parameter can in principle depend on the symmetric mass ratio $\eta$, and therefore, if it does, then $\beta$ is not the same for all binaries in our study. As a result, we are not combining results to establish a joint posterior on $\beta$, instead treating the $\beta$ measurement of each source independently. The gray shaded region at $\beta\gtrsim 10^{-5}$ is roughly where the approximation in Eqs. (3) and (4) breaks down, with $\beta(\pi\mathcal{M}f_{0})^{-2/3}\gtrsim 1$. The exact place where the $\beta$ correction becomes too large depends on the system’s mass and orbital period, which is comparable for all of the binaries we consider here, but not identical.

From these results we see that uncertainties in the measurement of $\beta$ are typically of order $10^{-4}$ for observation times of a few years, improving to ${\lesssim}10^{-5}$ after monitoring the sources for ten years. For the binaries to be measured well enough to definitively constrain $\beta$ within the region where the ppE formalism is valid, ${\sim}$decade long monitoring campaigns are needed.

IV.2 Redundancy Test and Fundamental Bias

Let us now consider the effect on WD binary parameter estimation if we assume GR is correct in the analysis, but include the possibility of non-zero dipole radiation in the injected source simulations. Treating each marginalized posterior of the orbit parameters $p(f_{0}|d),p(\dot{f}|d),p(\ddot{f}|d)$ independently, we can compute three independent constraints on the frequency-dependent chirp mass function $\mathcal{M}(f_{0})$, derived from the inversion of Eq. (1), or Eqs. (3) and (4). Doing so ignores the covariances between the different frequency parameters, which are shown to be non-negligible in Figs. 3 and  5, leading to a conservative estimate.

If the sources’ frequency evolution are consistent with GR, the three constraints will intersect in $\mathcal{M}-f_{0}$ space. Figs. 1 and  7 show $90\%$ credible intervals for $\mathcal{M}(f)$ from two of the detached WD binaries used in this study. The vertical line is the constraint from $p(f_{0}|d)$, which spans the prior range on $\mathcal{M}$. The bands sweeping from top left to bottom right are the constraints from $p(\dot{f}|d)$ (thin) and $p(\ddot{f}|d)$ (thick). The horizontal dashed line indicates the true value of the chirp mass from the simulated signal.

Similar to Fig. 1, Fig. 7 shows constraints for different injected sources with a five year (left panel) and ten year (right panel) observing time. When the injected source is simulated with a GR frequency evolution (green bands), one obtains a redundancy test like the one shown also in Fig. 1. Observe again that as the observing time increases, the redundancy test becomes stronger. When the injected source is simulated with a non-GR frequency evolution (orange and purple bands), one sees how the failure of the redundancy test could indicate a departure from GR.

The non-GR detections were simulated as follows. We used the frequency evolution in Eqs. (3) and (4) with $\beta\neq 0$, but still within the regime where $\beta(\pi\mathcal{M}f_{0})^{-2/3}\ll 1$. We use Source 4 as the test candidate, as it yielded some of the best constraints on $\beta$. The non-zero $\beta$’s were selected to be at the 0.99 and 0.999 percentile of the $p(\beta|d)$ from the 10 year observation in Fig. 6. The orange bands show results for the injection with $p$-value of $10^{-2}$, while the purple bands are for the source with non-zero $\beta$ at $p$-value of $10^{-3}$.

This demonstration leads to two notable conclusions. First, the left panel of Fig. 7, coming from a simulated 5 year observation, shows that the uncertainties in the $\ddot{f}$ measurement are such that signals whose underlying evolution is not consistent with GR can still pass the $\mathcal{M}$ consistency tests insofar as the three different constraints intersect in the same region of $\mathcal{M}-f$ space. However, the altered dynamics lead to a systematic or stealth bias Yunes and Pretorius (2009b); Vallisneri and Yunes (2013) in the inferred chirp mass that is larger than the statistical error (roughly the width of the intersection between the $p(f_{0}|d)$ and $p(\dot{f}|d)$ constraints). In other words, beyond-relativistic frequency evolution can meaningfully affect our inferences about the source before it is large enough to be directly detected. Without an accurate, independent measurement of $\mathcal{M}$ from EM follow-up observations, this effect would go unnoticed.

The right panel of Fig. 7 panel shows how this stealth bias is exposed by continued monitoring of the source, showing the extreme case of a 10 year observation. Here the constraints on the frequency evolution are tight enough that the bands are clearly not intersecting, proving that the assumption of only GR-driven evolution is not supported by the data.

The notion of non-zero $\beta$ fundamentally biasing chirp mass measurements is further illustrated in Fig. 8, now showing the 2D marginalized posterior $p(\mathcal{M}(f_{0},\dot{f}),\mathcal{M}(f_{0},\ddot{f})|d)$ (and therefore taking advantage of the information encoded in the parameter correlations). Here the covariance matrix of the MCMC samples has been used to approximate the posterior as a multivariate Gaussian, and we are plotting the 1, 2, and 3$\sigma$ contours in green. The line $y=x$ (purple, dashed) denotes where the two chirp mass constraints are self-consistent, i.e. $\mathcal{M}(f_{0},\dot{f})$ and $\mathcal{M}(f_{0},\ddot{f})$ return the same value. The solid orange lines mark the true value of the simulated source. The left panel is from a 5 year observation of a source with $\beta=0$, and we find that the posterior encompasses $y=x$ and encloses the true value.

The right hand panel shows results from the same duration observation and same source, but with $\beta$ injected at the $p$-value of $10^{-2}$ of the 10 year constraint. Note that the posterior includes the line $y=x$ so, without independent information about the true $\mathcal{M}$, this source would appear to be consistent with GR-only frequency evolution. However, the true value (orange line) is several $\sigma$ away from the median of the distribution, resulting in a ${\sim}15\%$ bias in the inferred $\mathcal{M}$ (and closer to $30\%$ for the $p\sim 10^{-3}$ case, not shown here).

V Projected Model-Dependent Constraints

In this section we specialize the model-independent test of the previous section to scalar-tensor theories. We begin by describing the basics of this theory and establishing some notation. We proceed by presenting projected constraints assuming GW observations can be used to extract $(f_{0},\dot{f},\ddot{f})$ in connection with electromagnetic observations of the same source, which can provide the masses and radii of the binary components. Throughout this section, we mostly follow the notation of Damour and Esposito-Farese (1998) for scalar-tensor theories.

V.1 Basics of Scalar-Tensor Theories and the Modeling of Individual Sources

The class of scalar-tensor theories we will consider can be defined by the Einstein-frame action

[TABLE]

where $\kappa_{*}=8\pi G_{*}$, with $G_{*}$ the bare gravitational constant in the Einstein frame, $\psi$ is a massless scalar field in the Einstein frame (not to be confused with the gravitational wave polarization angle $\psi_{\rm GW}$ from above), $R_{*}$ is the Ricci scalar in the Einstein frame associated with the Einstein frame metric $g_{\mu\nu}^{*}$, with $g_{*}$ its determinant, $\chi_{M}$ are additional matter degrees of freedom and $A(\psi)$ is a coupling function.

The action above describes a theory in which matter follows geodesics of the physical (Jordan-frame) metric

[TABLE]

We are here interested in a particular subclass of scalar-tensor theories defined by the conformal coupling

[TABLE]

where $\psi_{0}$ is the “vacuum-expectation value” of the field at spatial infinity. One could conformally transform the Einstein-frame action to the Jordan frame, but the resulting action is more difficult to work with, and so, we here choose to present equations in the Einstein frame only. The physical world in which experiments are performed, however, exists in the Jordan frame, so one must always remember to re-cast any Einstein frame results to the Jordan frame before comparing against experiments.

The field equations in the Einstein frame take the form

[TABLE]

where $G_{\mu\nu}^{*}$ is the Einstein tensor in the Einstein frame, $T_{*}^{M}{\!\!\!}_{\mu\nu}$ is the matter stress-energy tensor in the Einstein frame, and the stress-energy tensor for the Einstein-frame scalar field is

[TABLE]

and $T^{M}_{*}$ is the Einstein-frame trace of the matter stress-energy tensor $T_{*}^{M}{\!\!\!}_{\mu\nu}$ in the Einstein frame. Taking the divergence of the field equations for the metric, one finds an additional compatibility condition, which guarantees the physical metric satisfies the weak-equivalence principle.

In the Solar System, this theory can reduce closely to GR. To show this, let us say the Einstein frame field $\psi\to\psi_{0}$ in the Solar System. Then, at 1 PN order in the Solar System, Damour and Esposito-Farése have shown that the parametrized post-Newtonian (ppN) parameters Nordtvedt (1968); Will and Nordtvedt (1972) take the form Damour and Esposito-Farese (1992, 1993, 1996, 1998)

[TABLE]

The observation of the Shapiro time delay through tracking of the Cassini spacecraft can be satisfied by requiring that $\alpha_{0}^{2}\lesssim 10^{-3}$, without imposing any condition on $\beta_{0}$ Will (2006); Bertotti et al. (2003).

Even with a small value of $\alpha_{0}$, large modifications to GR can still occur in systems that are strongly gravitating. Damour and Esposito-Farése have shown that an effect akin to ferromagnetism occurs in the vicinity of sufficiently compact and massive stars, which leads to a non-linear activation of the scalar field, an effect called spontaneous scalarization Damour and Esposito-Farese (1992, 1993, 1996, 1998). In numerical simulations, such scalarization has also been seen to occur dynamically in NS binaries when the binary’s binding energy exceeds a certain threshold Barausse et al. (2013); Palenzuela et al. (2013). Scalarization can only occur when $\beta_{0}<0$, and in fact, Damour and Nordtvedt have shown that GR is an attractor in theory phase space when $\beta_{0}>0$. Recently, Anderson et al have shown that when $\beta_{0}<0$, GR is a repeller in theory phase space, leading to a maximal violation of Solar System experiments, unless one fine-tunes the field’s initial conditions at the time of Big Bang nucleosynthesis Sampson et al. (2014); Anderson et al. (2016); Anderson and Yunes (2017). We will here ignore such problems, as done cavalierly in the literature.

The degree of scalarization for isolated compact object $A$ is controlled by its scalar charge $\alpha_{A}$, which in general is a complicated function of $\alpha_{0}$, $\beta_{0}$ and the object’s compactness. For WDs, one can approximate the scalar charge via

[TABLE]

where $C_{A}$ is the star’s compactness. The scalar charge is here independent of $\beta_{0}$ because these objects cannot become compact without collapsing into a black hole for the non-linear effects associated with $\beta_{0}$ to become relevant. For NSs, the scalar charge can only be calculated numerically by solving the equations of structure in scalar-tensor theories; we here use the tabulated expression of Anderson .

Although Solar System observations cannot constrain $\beta_{0}$ since weakly gravitating objects have tiny scalar charges, the radio observation of binary pulsars can Stairs (2003); Freire et al. (2012). These observations have accurately measured the rate at which the orbital period of certain binary pulsars decays, finding it to be in perfect agreement with GR. But when $\beta_{0}$ is sufficiently negative, then not only do NSs scalarize, but this activation of the scalar field also sources dipole radiation, which in turn forces the orbital period to decay more rapidly than in GR. This is precisely an example of the pre-Newtonian dissipation we have discussed previously. Given that such rapid decay is not observed, radio observations have placed constraints in the $(\alpha_{0},\beta_{0})$ plane. Due to degeneracies between these parameters in the timing model, only a certain combination of $(\alpha_{0},\beta_{0})$ can be bounded for any given observation. Constraints on the $(\alpha_{0},\beta_{0})$ plane from radio observations have already been presented in Fig. 2.

We use the same LISA response function as described in Sec. II.1. In principle, there is an additional scalar GW mode that ought to be included, but this is suppressed by the scalar charges, and we will thus ignore it here. The mapping between components of the model parameter vector $\vec{\lambda}$ is still given by Eq. (3) and second-order version of Eq. (4), namely

[TABLE]

with ${\cal{M}}\to{\cal{M}}_{*}$ and the particular choices

[TABLE]

where $\alpha_{1,2}$ are the scalar charges of the objects, while

[TABLE]

The quantity $\beta_{1,2}$ are the derivatives of the scalar charges, which we will ignore in this work for simplicity. The quantity $d_{2}$ is a numerical factor that has not yet been calculated, but that can be ignored because it is subdominant. Notice that the quantity ${\cal{M}}_{*}$ in these equations is not quite the chirp mass, but rather it is the tensor chirp mass, namely

[TABLE]

where we have set $G_{*}=1$.

V.2 Constraints on Spontaneous Scalarization with $(\dot{f},\ddot{f})$ Measurements

Consider a NS-WD system in the galaxy with GW frequency in the LISA band. While much rarer than the WD-WD systems, a NS-WD binary is a clearly identifiable source in the LISA data stream due to its high mass, and will therefore be well localized ($d\Omega\ll 1$ sq. deg. after several years of observing) in the galaxy and yield a precise determination of the frequency evolution parameters Cornish and Robson (2017). If electromagnetic follow-up observations constrain the masses and radii of the binary constituents, they would be used to predict the orbital evolution of the binary in scalar-tensor theories of gravity. Those predictions, compared to the frequency parameters observed by LISA, then serve as a test of the theory.

Exploring to what extent such a test yields novel insight into our understanding of gravity, we simulate a binary comprised of a 1.4 $M_{\odot}$ NS and a 0.3 $M_{\odot}$ WD with a GW frequency $f_{0}=10$ mHz near the galactic center, at a distance of 8 kpc (see also Table 1). Using the MCMC pipeline, we characterize how well the source parameters are measured. For electromagnetic follow up we would likely require the source be closer to Earth, and not towards the galactic center where EM observations will have to contend with extinction and source confusion. However, placing the sources near the galactic center is a conservative choice for the GW analysis: anything further away is difficult to observe electromagnetically, and anything closer results in better GW parameter estimation due to the increased S/N.

The marginalized posteriors from the MCMC analysis are approximated as multivariate Gaussians, as in Eq. (5) with parameters $\vec{\lambda}=(f_{0},\dot{f},\ddot{f})$. It is important to consider the full covariance matrix, rather than just the variance of the 1D marginalized posteriors for $\dot{f}$ and $\ddot{f}$, because the two parameters are strongly correlated (see Fig. 5). While the perils of blindly assuming Gaussian distributions for posteriors in GW astronomy are well documented, the sources we are analyzing have high S/N ratios where this approximation is sufficient for proof of concept studies. Furthermore, we again confirm that the Gaussian approximation is adequate by testing that the distribution of p-values for the Markov chain samples evaluated using the Gaussian approximation to the posterior are uniformly distributed between 0 and 1, as we did in Sec. III.3.

As a stand in for the electromagnetic observations, we compute the WD and NS radii using mass-radius relations appropriate for the different compact objects. For the WD we use a simple Newtonian relation, and for the NS we test two different equations of state, shown in Fig. 9. Although currently the NS equation of state is not known accurately, by the 2030s when LISA flies, additional LIGO/Virgo observations like the recent Abbott et al. (2017b), combined with observations of the X-ray pulse profile of bursting NSs by NICER Ozel et al. (2016), will have measured it much more accurately. We then assume Gaussian errors for the masses of the objects and the radius, centered at the true value, testing assumed standard deviations of $10\%$, $1\%$, $0.1\%$, and $0.01\%$ of the true value. The NS radius is determined from the mass using the assumed equation of state.

In lieu of a closed-form expression for mapping the EM and GW parameters to scalar-tensor theory parameters, we test whether a set of discrete values for $\alpha_{0}$ and $\beta_{0}$ are jointly consistent with $p(\dot{f},\ddot{f}|d)$ approximated from the MCMC, marginalizing over the inferred $f_{0}$ measurement from the GW data and the assumed mass and radius uncertainties from EM observations. We begin by randomly drawing $10^{6}$ combinations of $f_{0}$, masses, and radii for the NS and WD consistent with $p(f_{0}|d)$ and the assumed EM measurement capabilities, respectively. For each set of masses and radii, we compute the scalar charge of the WD using Eq. (15). The NS scalar charge is obtained from a look up table of numerical simulations on a discretized grid of mass, $\alpha_{0}$, and $\beta_{0}$ Anderson . Then, using Eqs. (17)-(20), and substituting into Eqs. (3) and (16), we calculate the predicted $\dot{f}^{*}$ and $\ddot{f}^{*}$ for that combination of parameters. If the $p$-value of the marginalized posterior $p(\dot{f},\ddot{f}|d)$ evaluated at the predicted $\dot{f}^{*}$ and $\ddot{f}^{*}$ is greater than $0.1$ the scalar tensor parameters are considered consistent with the GW and EM measurements, from which upper limits are placed on $\alpha_{0}$ as a function of $\beta_{0}$.

The bound that can be set using joint EM+GW observations is heavily dependent on how well the masses and radii of the binary components can be independently measured by EM observations. Upper bounds at 90% confidence for a ten year observation of our canonical NS-WD binary system, assuming mass and radius errors ranging from $10\%$ to (an admittedly optimistic) $0.01\%$, are compared to existing bounds in Fig. 2. For LISA observations to improve on existing bounds of scalar-tensor theories, long observing times and independent precision measurement of the mass and radii of the binary are necessary.

To see if this conclusion is generically true, or if it depends on properties of the source, we repeat the analysis for different observing times, orbital periods of the binary, WD mass, mass ratio, and NS equation of state. Below, we present the $1\%$ and $0.01\%$ accuracy cases for EM-derived mass and radius constraints as representative of plausible and best case outcomes, respectively.

Figure 10 shows how the LISA constraints improve over the mission lifetime. For systems with mass and radius constraints in the $1\%$ range (left panel), the uncertainty in the properties of the binary constituents dominates the constraint on $\alpha_{0}$, and we see minimal improvements even as the orbital evolution parameters become increasingly well measured with longer observing times. Adopting the optimistic mass and radius constraints (right panel), the bounds become competitive, and ultimately exceed, existing bounds from Cassini and observations of binary pulsars. Because we are assuming joint electromagnetic observations, these systems could be monitored well beyond the LISA mission lifetime, perhaps with follow-on GW observatories, to further improve the bounds.

We also investigate how the orbital period of the binary affects the bounds that can be inferred, which presents an interesting trade space to explore. Dipole radiation plays a more prominent role in the orbital evolution at lower orbital velocities, i.e. low GW frequency. However, LISA will have less sensitivity at lower frequencies due to the presumed levels of confusion noise from unresolved binaries in the galaxy. Furthermore, measuring frequency derivatives is more difficult at low frequency because the overall change to the GW phase during the observing time is smaller.

We test which of these effects wins out by comparing bounds derived for the $1.4+0.3M_{\odot}$ NS-WD system at different GW frequencies ($f_{0}=3,5,10$ and $20\ {\rm mHz}$) in Fig. 11. Observe how we again find, for $1\%$ uncertainties, that the mass/radius measurement dominates the bound on scalar-tensor parameters. If EM constraints could reach the $.01\%$ level, we see separation of the bounds depending on the frequency of the binary, with the higher frequency systems providing tighter bounds owing to the better measurement of $\dot{f}$ and $\ddot{f}$. Interestingly, the bound does not monotonically improve with frequency, as the $3\ {\rm mHz}$ source better constrains $\alpha_{0}$ than the $5\ {\rm mHz}$ system.

The final result presented here explores how the mass of the WD companion affects the bounds on the scalar-tensor parameters. The reason for this study is that a higher mass WD would be more compact, and therefore would have larger scalar charge, which, in turn, would be a stronger source of dipole radiation. Alas, we find little difference in the bounds derived from the GW measurement comparing $1.4+0.3M_{\odot}$ and $1.4+1.3M_{\odot}$ NS-WD binaries at $f_{0}=20\ {\rm mHz}$, as we can see in Fig. 12. We also investigated the effect of the NS composition, comparing bounds using the AP3 Akmal et al. (1998) and SLy4 Douchin and Haensel (2001) equations of state, and found no measurable difference between the derived bounds.

VI Caveats

In this section we present several caveats on our analysis that could be investigated further in the future. We begin by discussing the dependence of our results on the galactic binary population. We end with a short discussion of the effects of tidal and/or rotational deformations.

VI.1 Dependence on the Population

For the WD binary sources, we are using a small sample of the shortest period binaries in a simulated galaxy. These have shorter periods than any known binaries at the time of this study, and so there existence is a prediction of population simulations. From those simulations, there will be a large amount of variance in this high mass, short period sub-population, and the sample of such binaries in the observed galaxy may not be as kind. That being said, current observations of binaries are limited to nearby, bright sources, whereas sources with these masses and GW frequencies will be observable throughout the galaxy so, if they exist, they will be excellent LISA sources. A similar argument exists for NS-WD binaries in the LISA band. Further population synthesis studies, and electromagnetic surveys searching for short period binaries, may strengthen or weaken the case for these extraordinary systems. It may also be the case that we will just have to wait to see what LISA observations have in store.

VI.2 Tidal and/or Rotational Deformations

In this work we have intentionally isolated a particular modification to the binary evolution (dipole radiation) to determine how well it can be constrained using LISA observations. The effect of dipole radiation will be competing with other, similarly difficult to measure, processes that will modify the phase evolution of the binary from the point-particle-in-GR case. Mass transfer is an obvious culprit, though such binaries are expected to be a subset of the total observed population, particularly at the high end of the $\mathcal{M}$ distribution, where our simulated sources are found.

Tidal interactions between the stars, however, may prove to play an important role generically, especially for the high frequency binaries we consider here. Fortunately for us, tidal contributions enter at much higher post-Newtonian order. More specifically, the tidal deformability enters at $5$ post-Newtonian order, and thus, its effects would scale as $(\pi{\cal{M}}f_{0})^{21/3}$ in Eq. (1) and $(\pi{\cal{M}}f_{0})^{29/3}$ in Eq. (2). These modifications are multiplied by a tidal deformability parameter that for WDs is very large, as it scales as $C^{-5}$ where $C$ is the gravitational compactness, and that is enhanced by non-zero orbital eccentricity Willems et al. (2008, 2010); Piro (2011); Valsecchi et al. (2012). However, an order of magnitude estimate suggests that the 5PN suppression is at least 5 orders of magnitude stronger than the tidal deformability parameter enhancement, leading to an overall effect that should be very small. Ultimately, an holistic approach to modeling the phase evolution of galactic binaries in LISA data including these different effects will be needed to avoid biasing inferences about any one potential source of orbital evolution.

VII Conclusions

We studied the potential of using galactic binary GW observations with LISA to carry out tests of GR. In spite of the current lore, there are a large class of modified theories of gravity that predict large modifications to the temporal or spectral evolution of such GWs. We thus developed a model-independent framework to carry out null tests of GR with such GWs, which resemble those carried out with radio observations of binary pulsars. We performed several MCMC studies to determine the strength of such null tests, and found that long observing times (of order $5$–$10$ years) would be most probably required.

We then investigated model-dependent test of GR with future LISA observations of galactic binaries. In particular, we focused on tests of scalar tensor theories that predict dipole emission in compact binaries. We performed several MCMC studies to determine the strength of such model-dependent tests, and found that shorter observing times (of ${\cal{O}}(3$–$5)$ years) would suffice to place better than current constraints. The caveats here though are that such model-dependent tests require an observation of a mixed NS-WD galactic binary and knowledge of the radius of the NS, which could be obtained with coincident electromagnetic observations if the binary is eclipsing.

We also investigated whether modifications to GR in the signal could affect the astrophysical inferences one would obtain from LISA observations of galactic binaries analyzed with GR templates. We found that indeed fundamental bias can be present, inducing a systematic error in the extracted GR parameters of ${\cal{O}}(10\%)$, even before a deviation is measurable. Such stealth bias eventually disappears as the observing time increases and the deviations begin to induce a failure of the null tests we developed. This implies that accurate chirp mass inferences with GR templates may require marginalizing over these non-GR effects to within existing bounds.

We hope the analysis presented here catalyzes an explosion of future work. One issue that is particularly important is the impact of astrophysical systematics on the model-independent and model-dependent tests considered here. Perhaps the most important astrophysical effect are tidal deformations on WDs, which would induce a correction on the waveform template at high post-Newtonian order. These corrections should be taken into account to determine whether they can spoil the tests discussed here.

Another possible avenue for future work is the study of follow-on missions beyond LISA. Even after LISA flies and is de-commissioned, electromagnetic observing campaigns of the Milky Way will continue. Long enough observations could then allow for measurements of the radii of eclipsing NS-WD galactic binaries, even of these radii measurements are not available when LISA flies. Once this is done, one could then go back to the LISA data to carry out the model-dependent tests described here. Moreover, follow-on space-based GW missions, like the Big Bang observer, will observe galactic binaries later in their orbital evolution, and at much higher S/N, leading to even stronger constraints. It would be interesting to investigate the level to which such tests could be carried out.

The analysis described above, however, requires knowledge of the equation of state. By the time LISA flies in the 2030s, one expects the equation of state to have been constrained significantly with future observations by LIGO/Virgo and NICER. However, still some level of uncertainty will remain, which would have to be marginalized over to carry out the model-dependent tests described here. It would be interesting to carry out such a marginalization analysis to determine how the strength of the model-dependent tests would be affected.

Acknowledgements.

We would like to thank Neil Cornish for useful comments and discussion. We are also most appreciative to David Anderson, who provided thorough tables for the scalar charges of NSs as a function of compactness and the coupling constants $(\alpha_{0},\beta_{0})$. N. Y. acknowledges support from NASA grants 4W5883, NNX16AB98G and 80NSSC17M0041. T. L. was supported by NASA grant NNH15ZDA001N and the NASA LISA Study Office.