Standard sirens with a running Planck mass

TL;DR

This paper investigates how a time-varying Planck mass affects gravitational wave propagation and standard siren measurements, showing that GW observations can constrain modifications to gravity and impact cosmological parameter estimates.

Contribution

It demonstrates that GW observations can jointly constrain the running of the Planck mass and the Hubble constant, providing constraints comparable to current cosmological methods.

Findings

GW observations constrain $c_M$ to within ±0.9 with 100 BNS events

Current GW data from GW170817 constrains $c_M$ to approximately -9^{+21}_{-28}

A nonzero $c_M$ can bias $H_0$ measurements if not properly accounted for

Abstract

We consider the effect of a time-varying Planck mass on the propagation of gravitational waves (GWs). A running Planck mass arises naturally in several modified gravity theories, and here we focus on those that carry an additional dark energy field responsible for the late-time accelerated expansion of the universe, yet--like general relativity (GR)--propagate only two GW polarizations, both traveling at the speed of light. Because a time-varying Planck mass affects the amplitude of the GWs and therefore the inferred distance to the source, standard siren measurements of $H_0$ are degenerate with the parameter $c_M$ characterizing the time-varying Planck mass, where $c_M=0$ corresponds to GR with a constant Planck mass. The effect of non-zero $c_M$ will have a noticeable impact on GWs emitted by binary neutron stars (BNSs) at the sensitivities and distances observable by ground-based GW detectors such as advanced LIGO and A+, implying that standard siren measurements can provide joint constraints on $H_0$ and $c_M$. Assuming a $\Lambda$CDM evolution of the universe and taking Planck's measurement of $H_0$ as a prior, we find that GW170817 constrains $c_M = -9^{+21}_{-28}$ ($68.3\%$ credibility). We also discuss forecasts, finding that if we assume $H_0$ is known independently, then 100 BNS events detected by advanced LIGO can constrain $c_M$ to within $\pm0.9$. This is comparable to the current best constraints from cosmology. Similarly, for 100 LIGO A+ BNS detections, it is possible to constrain $c_M$ to $\pm0.5$. When analyzing joint $H_0$ and $c_M$ constraints we find that $\sim 400$ LIGO A+ events are needed to constrain $H_0$ to $1\%$ accuracy. Finally, we discuss the possibility of a nonzero value of $c_M$ biasing standard siren $H_0$ measurements from 100 LIGO A+ detections, and find that $c_M=+1.35$ could bias $H_0$ by 3-4$\sigma$ too low if we incorrectly assume $c_M=0$.

Equations

$$d{s}^2=a(\tau {)}^2\left[-d{\tau }^2+d{\vec{x}}^2\right],$$

$$3{\mathcal{H}}^2{M}_P^2=a^2({\rho }_m+{\rho }_{\mathrm{DE}}),$$

$${\rho }_i^{\prime }+3\mathcal{H}({\rho }_i+P_i)=0,i=m,\mathrm{DE}.$$

$$g_{\mu \nu }={\bar{g}}_{\mu \nu }+h_{\mu \nu },|h|\ll |\bar{g}|,$$

$$S=\frac{1}{2}\int d^{3}xd\tau M_{\ast }^2{a}^2\left[h_A^{{}^{\prime }2}-c_T^2(\vec{\nabla }{h}_A{)}^2\right],$$

$$h_A^{{}^{\prime \prime }}+\left[2+{\alpha }_M(t)\right]\mathcal{H}{h}_A^{{}^{\prime }}+k^2{c}_T^2{h}_A=0,$$

$${\alpha }_M\equiv \frac{d\ln (M_{\ast }/M_P{)}^2}{d\ln a}=\frac{2}{\mathcal{H}}\frac{M_{\ast }^{{}^{\prime }}}{M_{\ast }}.$$

$$S_s=\int d^{4}x\sqrt{-g}\left[G_4(\varphi )R+K(X,\varphi )-G_3(X,\varphi )\square \varphi \right],$$

$$S_v=\int d^{4}x\sqrt{-g}\left[R+G_2(X,F,Y)+G_3(X){\nabla }_{\mu }{A}^{\mu }\right],$$

$$d{s}^2=-a(\tau {)}^2[(1+2\Phi )d{\tau }^2+(1-2\Psi )d{\vec{x}}^2],$$

$$\bar{\square }\Phi =4\pi G_{\text{eff}}{\rho }_m{\Delta }_m,$$

$$\gamma =\Psi /\Phi .$$

$${\hat{h}}_A^{{}^{\prime \prime }}+\left(k^2-\frac{a_{\text{GW}}^{{}^{\prime \prime }}}{a_{\text{GW}}}\right){\hat{h}}_A=0.$$

$$h_A^o=\frac{M_{\ast }(z)}{M_{\ast }(z=0)}{h}_{A,\mathrm{GR}}^o\propto \frac{M_{\ast }(z)}{M_{\ast }(z=0)}\frac{1}{d_L(z)}\equiv \frac{1}{d_{\mathrm{GW}}},$$

$$d_L=\frac{1+z}{H_0}{\int }_0^z\frac{d\tilde{z}}{\hat{H}(\tilde{z})}$$

$$\hat{H}(z)=\sqrt{{\Omega }_{m,0}(1+z{)}^3+{\rho }_{\mathrm{DE}}(z)/{\rho }_c}.$$

$${\alpha }_M(z)=c_M\frac{{\Omega }_{\mathrm{DE}}(z)}{{\Omega }_{\mathrm{DE},0}},$$

$$\frac{d_{\mathrm{GW}}}{d_L}=\frac{M_P}{M_{\ast }}=\exp \left\{\frac{1}{2}{\int }_0^z\frac{d{z}^{\prime }}{1+z^{\prime }}{\alpha }_M(z^{\prime })\right\},$$

$$\frac{d_{\mathrm{GW}}}{d_L}=\exp \left\{\frac{c_M}{2{\Omega }_{DE,0}}\ln \frac{1+z}{({\Omega }_{m,0}(1+z{)}^3+{\Omega }_{DE,0}{)}^{1/3}}\right\},$$

$$\frac{d_{\mathrm{GW}}}{d_L}\approx 1+\frac{1}{2}{c}_{M}z+\frac{1}{8}{c}_M\left(c_M-2-6{\Omega }_{m,0}\right)z^2+\mathcal{O}(z^3).$$

$$H_0=\frac{z}{d_{\mathrm{GW}}}+z^2\frac{\left[1-(3/4){\Omega }_{m,0}+(1/2)c_M\right]}{d_{\mathrm{GW}}}+\mathcal{O}(z^3),$$

$$\underset{z\to \infty }{\lim }\frac{|d_{\mathrm{GW}}-d_L|}{d_{\mathrm{GW}}}=\left|1-\exp \left\{-\frac{c_M}{2{\Omega }_{DE,0}}\ln {\Omega }_{m,0}^{-1/3}\right\}\right|.$$

$$d_{\mathrm{GW}}=d_L(z,H_0,{\Omega }_{m,0})R(z,c_M,{\Omega }_{m,0})$$

$$p(d_{\mathrm{GW}}^{\mathrm{obs}}\mid d_{\mathrm{GW}})=N_{[\mu =d_{\mathrm{GW}},\sigma =1.8/\rho ]}(d_{\mathrm{GW}}^{\mathrm{obs}})$$

$$p(x_{\mathrm{GW}}\mid d_{\mathrm{GW}},\omega )=\int p(x_{\mathrm{GW}}\mid d_{\mathrm{GW}},\omega ,\iota ,z,\xi )d\iota dzd\xi .$$

$$p(x_{\mathrm{GW}},x_{\mathrm{EM}}\mid c_M,H_0)=\frac{\int p(x_{\mathrm{GW}}\mid d_{\mathrm{GW}}={\hat{d}}_{\mathrm{GW}}(z,c_M,H_0,\Xi ),\omega )p(x_{\mathrm{EM}}\mid z,\omega )p(z,\omega )p(\Xi )dzd\omega d\Xi }{\beta (H_0,c_M)},$$

$$p(x_{\mathrm{EM}}\mid z,\omega )=\delta (z-z_{\mathrm{obs}})\delta (\omega -{\omega }_{\mathrm{obs}}).$$

$$p(x_{\mathrm{EM}}\mid z,\omega )=N_{[\mu =z,\sigma ={\sigma }_z]}(z_{\mathrm{obs}})\delta (\omega -{\omega }_{\mathrm{obs}}),$$

$$S_s=\int d^{4}x\sqrt{-g}\left[R+K(X,\varphi )-G_3(X,\varphi )\square \varphi \right].$$

$$\beta (c_M,H_0)={\int }_{x_{\mathrm{GW}}>x_{\mathrm{GW}}^{\mathrm{thresh}},x_{\mathrm{EM}}>x_{\mathrm{EM}}^{\mathrm{thresh}}}p(x_{\mathrm{GW}}\mid d_{\mathrm{GW}}={\hat{d}}_{\mathrm{GW}}(z,c_M,H_0,\Xi ),\omega )p(x_{\mathrm{EM}}\mid z,\omega )p(z,\omega )p(\Xi )dzd\omega d\Xi d{x}_{\mathrm{GW}}d{x}_{\mathrm{EM}},$$

Full text

I Introduction

Einstein’s general theory of relativity (GR) is the foundation of gravity. On Solar System scales, not only do its predictions show remarkable agreement with astrophysical data, but precise measurements of phenomena such as the deflection of light around the Sun and the perihelion shift of Mercury rule out many modifications to GR Will (2014); Berti et al. (2015). Nonetheless, GR exhibits weaknesses at both the very high and the very low energy regimes. At high energies, unavoidable singularities arise during gravitational collapse, and the so-called renormalization problem limits our understanding of quantum gravity Joshi and Malafarina (2011); ’t Hooft and Veltman (1974); Deser (2000). In order to fit observational data on cosmological scales, GR relies on the presence of exotic unknown matter components, namely dark matter and dark energy, to make up $95\%$ of the total energy content of the Universe Aghanim et al. (2018). These issues show the current limitations in our understanding of how gravity behaves and interacts with matter in extreme energy regimes. As a result, a number of modified gravity theories have been proposed (see e.g. Clifton et al. (2012); Koyama (2016); Joyce et al. (2015) for summaries and reviews), and it is necessary to analyse their consistency, viability, and consequences.

Recently, multimessenger astrophysics has shown great potential for revealing new details of physical phenomena and testing gravity. In particular, combined gravitational wave (GW) and electromagnetic (EM) science is becoming a central topic, allowing us to test different aspects of cosmology and fundamental physics, including alternative dark energy candidates, spatial curvature, the expansion rate of the universe, strong lensing sources, and the graviton mass Wei (2018); Wei and Wu (2017); Baker and Trodden (2017); Collett and Bacon (2017), among others. For instance, the recent detection of GWs from a binary neutron star (BNS) merger, GW170817 et. al. (2017a), by the advanced Laser Interferometer Gravitational Wave Observatory (LIGO) LIGO Scientific Collaboration et al. (2015) and advanced Virgo Acernese et al. (2015), in conjunction with the detection of an EM counterpart et. al. (2017b, c) constrained the propagation speed of GWs, $c_{T}$, to be $|c_{T}/c-1|\lesssim 10^{-15}$ relative to the speed of light $c$. This high-precision constraint provides information about possible modifications to GR and strongly disfavors a number of gravity theories proposed in the literature Baker et al. (2017); Creminelli and Vernizzi (2017); Sakstein and Jain (2017); Ezquiaga and Zumalacarregui (2017); Wang et al. (2017); Lombriser and Lima (2017).

The detection of GW170817 and its EM counterpart enabled the first standard siren measurement Abbott et al. (2017) of the current rate of expansion of the Universe, $H_{0}$, assuming $\Lambda$CDM cosmology and general-relativistic GW propagation. Besides Abbott et al. (2017)’s result, independent methods have been used to constrain the Hubble constant, falling into two major categories: large-scale and local observations. From the cosmic microwave background (CMB), the Planck mission The Planck Collaboration (2006) found $H_{0}=67.4\pm 0.5$ km/s/Mpc Aghanim et al. (2018), whereas from local observations of Type-Ia supernovae, the SHoES survey Riess (2006) found $H_{0}=73.48\pm 1.66$ km/s/Mpc Riess et al. (2018). These two measurements disagree at the $3.5\sigma$ level. Consequently, obtaining independent GW constraints on $H_{0}$ with $\sim 1\%$ accuracy could substantially improve understanding of this tension and potentially reveal new physics or sources of systematic error. However, if gravity is described by a theory other than GR, an accurate inference of $H_{0}$ could be hampered by degeneracies with modified gravity effects. Accordingly, in this paper we constrain $H_{0}$ with combined GW and EM observational data in the context of alternative theories of gravity with modified GW propagation.

The value of $H_{0}$ can be obtained if we know the redshift of the source and its luminosity distance. The GW amplitude is inversely related to the luminosity distance, and the scaling constant can be directly inferred from the evolution of the gravitational waveform, which depends on the theory of gravity (Schutz, 1986). While the source’s redshift cannot be inferred directly from GWs, for events with EM counterparts it is possible to obtain the redshift from the EM spectrum. Sources of this kind are known as “standard sirens” Holz and Hughes (2005), after their EM analogues known as “standard candles”. In this context, the first BNS detection, GW170817, has already enabled the first standard siren measurement, $H_{0}=70.0^{+12.0}_{-8.0}$ km $\text{s}^{-1}$ $\text{Mpc}^{-1}$ Abbott et al. (2017), and Chen et al. (2017) suggested that a detection of $\sim 100$ BNS events could yield a measurement with percent-level accuracy (see also (Nissanke et al., 2013; Feeney et al., 2018)). This target may be attainable in as little as five years, as scheduled improvements bring the LIGO and Virgo detectors to their design sensitivity Chen et al. (2017). Within the next decade, the addition of two more detectors, KAGRA Aso et al. (2013) and LIGO-India Iyer et al. (2011), to the LIGO-Virgo network, as well as the planned upgrade of advanced LIGO to LIGO A+ McClelland et al. (2015), will further increase sensitivity to BNS mergers, which occur at an astrophysical rate of 110 Gpc*-3* yr*-1* or greater The LIGO Scientific Collaboration and the Virgo Collaboration (2018).

Inspired by modified gravity models that introduce new degrees of freedom as an alternative to dark energy, we analyze how modified propagation of GWs affects the standard siren estimation of $H_{0}$. In the context of cosmology, some alternative models have indeed been found to affect $H_{0}$ measurements and reduce the tension between the cosmological $\Lambda$CDM and local $H_{0}$ values (see Ezquiaga and Zumalacarregui (2018) for a recent review and discussion). Since the estimated luminosity distance depends on the gravitational theory, we expect it to be different in models where the GW propagation changes due to the presence of additional fields. The impact of modified GW propagation on the measurement of $H_{0}$ was recently considered in Belgacem et al. (2018), which studied a specific model where the amplitude of the GWs is damped due to modified gravity effects in addition to the damping caused by the expansion of the universe. Specifically, this work considered a non-local gravity model and analyzed how constraints on $H_{0}$ are affected by modifications to gravity and to the equation of state of dark energy.

In this paper, we stay agnostic regarding the specific underlying gravity theory and study $H_{0}$ in a model-independent manner, although for simplicity we assume the background universe to evolve according to the $\Lambda$CDM model, while allowing for general modifications for cosmological perturbations and gravitational waves. We compare constraints on $H_{0}$ with and without considering cosmological data, and quantify the impact of damping effects caused by modified gravity. In particular, we identify the extra damping effect with an effective running of the Planck mass, as it is typically found in common modified gravity theories, whose time evolution is given by the fractional dark energy density $\Omega_{DE}$. This time evolution will have one constant arbitrary parameter $c_{M}$, representing the Planck mass rate today, with $c_{M}=0$ corresponding to GR with a constant Planck mass. The observable consequences of a running Planck mass have been studied in a number of regimes (see e.g. Uzan (2011) for a review), and its effect on the propagation of gravitational waves has been considered in Pettorino and Amendola (2015); Lombriser and Taylor (2016); Nishizawa (2018); Arai and Nishizawa (2018); Amendola et al. (2017); Linder (2018). We find that external cosmological data from Planck, in conjunction with GW170817, constrain $c_{M}$ to be $-81<c_{M}<28$ at the $95\%$ credible level. These constraints are very weak compared to those from cosmological data alone when assuming a specific modified gravity theory, as $c_{M}$ also affects CMB and structure formation. For instance, for Horndeski theories, current cosmological constraints give $-0.62<c_{M}<+1.35$ at $95\%$ confidence level Noller and Nicola (2018).

In addition, we consider populations of events and discuss future forecasts for standard sirens with advanced LIGO and LIGO A+. We find that with A+, 100 BNSs with detected EM counterparts can lead to cosmology-independent constraints on $H_{0}$ with an accuracy of $\sim 3\%$, and on $c_{M}$ with $\sigma(c_{M})\sim 0.9$. From these results we estimate the need for $\sim 400$ events in order to obtain a $1\%$-accurate constraint on $H_{0}$ in the presence of a running Planck mass, thereby matching current local and cosmological constraints. Furthermore, we find that $H_{0}$ and $c_{M}$ are highly degenerate, which highlights the importance of testing for the parameter $c_{M}$ to avoid biasing the inferred value of $H_{0}$ by assuming GR (i.e. $c_{M}=0$) in a setting where the true gravitational wave physics is described by $c_{M}$ of $\mathcal{O}(1)$. In particular, we show that if we have a population of 100 events detected by LIGO A+ with $c_{M}=1.35$, then the inferred $H_{0}$, assuming $c_{M}=0$, will typically be $>3\sigma$ away from the true value. In this case, the true value of $H_{0}$ would be outside the posterior $99\%$ credible interval due to the misplaced assumption of $c_{M}=0$.

The paper is structured as follows. In Section II we describe how a time-dependent Planck mass modifies the propagation of GWs, as well as other local and cosmological observables. In Section III we discuss the use of BNS mergers as standard sirens in the context of a running Planck mass. We describe our inference of $H_{0}$ in Section IV, and in Section V we show how estimates of $H_{0}$ change compared to GR. Finally, in Section VI we summarize our results and give an outlook for future work. We set the speed of light to unity ($c=1$) throughout.

II Running Planck mass

Let us start by considering a perfectly homogeneous and isotropic spatially flat cosmological background. In this case, the metric line element in conformal time $\tau$ takes the form of the Friedmann-Robertson-Walker (FRW) solution:

[TABLE]

where $a(\tau)$ is the scale factor determining the expansion of the universe. We next suppose that gravity is described by a modified theory with a time-dependent Planck mass. As previously mentioned, alternative theories typically include new degrees of freedom, which can interact non-trivially with the metric to produce an effective running of the Planck mass. In such a case, there is an ambiguity when defining the stress-energy tensor: the new effects can be interpreted as modifying gravity (i.e. the Einstein tensor $G_{\mu\nu}$), or as modifying the matter content of the universe (i.e. the stress-energy tensor $T_{\mu\nu}$). However, without loss of generality, we can always adopt the latter perspective and write the background equations of motion in the standard form

[TABLE]

where $\mathcal{H}=a^{\prime}/a$ is the conformal Hubble rate (with ${}^{\prime}\equiv d/d\tau$), $M_{P}$ is the constant Planck mass, $\rho_{m}(\tau)$ is the energy density of cold dark matter (CDM), and $\rho_{\rm DE}(\tau)$ is the energy density of dark energy, which is to encapsulate all modified gravity effects. According to this definition, each of the fluid energy-density components is separately conserved:

[TABLE]

We note that the choice to write the background equation as in eq. (2) is arbitrary, and we could have instead defined it with an effective time-dependent mass $M_{*}(\tau)$ instead of $M_{P}$ on the left-hand side and a different $\rho_{\rm DE}$ on the right-hand side. However, in that case, the energy-density components would not have been separately conserved, and thus for simplicity we avoid this choice.

Here we have assumed that there is an additional degree of freedom that couples to the metric in a non-trivial way, which leads to an arbitrary $\rho_{\rm DE}$. In contrast, standard matter components such as baryons and photons have been assumed to be minimally coupled to the metric, as usual, and therefore the way they contribute to the background equations is unchanged. Furthermore, their propagation and evolution are determined by standard geodesics in the given metric background. In particular, for perturbations about FRW, light still propagates at speed $c$.

Next, let us consider small cosmological perturbations around this background universe, and write the total metric as

[TABLE]

where $\bar{g}_{\mu\nu}$ is the background FRW metric and $h_{\mu\nu}$ is a linear perturbation whose transverse and traceless part $h^{TT}_{\mu\nu}$ encodes the GW amplitude. The other components of $h_{\mu\nu}$ determine the amplitude of matter energy-density perturbations, which will be ignored for our purposes.

GR is a single metric theory for a massless spin-2 particle, and hence $h^{TT}_{\mu\nu}$ propagates two physical degrees of freedom corresponding to two tensor polarizations. We consider modified gravity theories that do not propagate additional tensor modes, and therefore $h^{TT}_{\mu\nu}$ alone carries all the information on the evolution of the metric polarizations. Generically, such theories have the following quadratic action

[TABLE]

where we have expanded $h^{TT}_{\mu\nu}$ into two independent polarization components $h_{A}$ with $A=+,\times$. The equation of motion for GWs will thus be given by

[TABLE]

where we have transformed $h_{A}$ to the spatial 3D Fourier space; its amplitude therefore implicitly depends on time $\tau$ and wavenumber $k$. Here, $\alpha_{M}$ has been defined as

[TABLE]

We recover GR when $M_{*}=M_{P}$ (i.e. $\alpha_{M}=0$) and $c_{T}=1$. However, in modified gravity theories, additional gravitational fields such as scalars or vectors will generically have a time-dependent solution in cosmological backgrounds, which can induce a running of the Planck mass $M_{*}(\tau)$ due to conformal couplings, even in a local frame. Additional fields can also modify the propagation speed of GWs $c_{T}(\tau)$ due to non-minimal and non-conformal couplings. We remark that, in this background, both $M_{*}$ and $c_{T}$ are functions of time only, and are thus isotropic and polarization-independent.

Note that we have not added any non-derivative terms to eq. (5). This is a consequence of our assumption that the graviton remains massless in the modified theory, and propagates only two polarizations, like in GR. Non-derivative terms like $m^{2}h_{A}^{2}$ can appear in massive gravity theories de Rham et al. (2011); de Rham and Gabadadze (2010), but their GWs are described by five different polarization modes, instead of the two modes $h_{A}$. More complicated models have also been developed, such as ones in which GWs are directly coupled to additional fields; the perturbations of these fields would appear explicitly in eq. (5). This is the case of bigravity theories Hassan and Rosen (2012), which propagate one massless and one massive graviton. In these models, GWs can oscillate between the two gravitons (by analogy with neutrinos) and lead to distinctive signals in the waveform Max et al. (2017, 2018). A similar phenomenon is present in multi-vector theories with internal SU(2) symmetry Beltran Jimenez and Heisenberg (2018). Nonetheless, we will not consider such cases in this paper, focusing exclusively on gravity theories that lead to the action of eq. (5).

Due to the aforementioned constraints on the speed of GWs from GW170817, we will focus only on running Planck mass models where $c_{T}=1$. Well-known modified gravity theories that lead to eq. (5) with $c_{T}=1$ are scalar-tensor theories in the Horndeski Horndeski (1974); Deffayet et al. (2011) and Beyond Horndeski Gleyzes et al. (2015); Zumalacarregui and Garcia-Bellido (2014) families that have $c_{T}=1$ Baker et al. (2017). Their action takes the form

[TABLE]

where $\phi$ is an additional scalar gravitational field responsible for dark energy, $X=-\nabla_{\mu}\phi\nabla^{\mu}\phi/2$ is the kinetic term of the scalar field, and $G_{4}$, $K$ and $G_{3}$ are arbitrary functions. In this case, $\alpha_{M}=G_{4}(\phi)$. Specific models belonging to this category are quintessence, $f(R)$ gravity, kinetic gravity braiding and Jordan-Brans-Dicke theory Baker et al. (2017). Ignoring the matter sector, generalizations of the gravitational action in eq. (8) can be obtained by performing a disformal transformation of the metric (see e.g. Ezquiaga and Zumalacarregui (2017)) to obtain theories belonging to the Degenerate Higher Order Scalar-Tensor theory (DHOST) family Langlois and Noui (2016); Crisostomi et al. (2016); Achour et al. (2016), which retain the same structure for the GW action. See e.g. Kase and Tsujikawa (2018); Langlois (2018); Kobayashi (2019) for reviews on the status of scalar-tensor gravity theories after GW170817.

We also mention that vector-tensor theories belonging to the Generalized Proca Heisenberg (2014) family with $c_{T}=1$ have an action of the form Baker et al. (2017)

[TABLE]

where $A^{\mu}$ is an additional vector gravitational field responsible for dark energy; $G_{2}$ and $G_{3}$ are arbitrary functions of $X=-A^{\mu}A_{\mu}/2$, $F=-F_{\mu\nu}F^{\mu\nu}/4$, or $Y=A^{\mu}A^{\nu}F_{\mu}^{\alpha}F_{\nu\alpha}$, where $F_{\mu\nu}=\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}$. Since there are no conformal couplings between the vector field and the metric, it is straightforward to see that these models will also have constant $M_{*}$, and hence no modification will be seen for tensor polarizations (although the evolution of energy-density matter perturbations will be modified) Amendola et al. (2017).

The same happens for Lorentz-breaking vector-tensor theories of the Einstein-Aether family Jacobson and Mattingly (2001); Zlosnik et al. (2006). In this case, the vector field is time-like, and hence defines a preferred frame of reference, and the subclass satisfying $c_{T}=1$ also has constant $M_{*}$. We emphasize, however, that theories involving multiple vector fields do allow for non-trivial derivative couplings while still maintaining $c_{T}=1$ Beltran Jimenez and Heisenberg (2018). However, this case is not encompassed by eq. (5), as such models have an additional field explicitly coupled to $h_{A}$. Extended scalar-vector-tensor theories can also have $c_{T}=1$ and a non-trivial $M_{*}(t)$; these are encompassed by eq. (5) Heisenberg et al. (2018). Generalizations to action (9) including higher derivative interactions have been also studied in Kreisch and Komatsu (2018), where theories with $c_{T}=1$ and non-trivial predictions for GW astronomy have been found.

Having discussed examples of dark energy theories that can lead to an effective running Planck mass, we now turn to the observable consequences of such a modification. The time variation of fundamental constants has previously been studied in different regimes. Below, we survey its effect on various observables (see e.g. Uzan (2011) for a detailed review).

Cosmology: The effects of a running Planck mass on cosmology have been considered in Ade et al. (2015); Pettorino and Amendola (2015); Ade et al. (2016); Huang (2016) among others. The time-dependence of the Planck mass arises due to the presence of additional fields, which can modify the background evolution of the universe through a time-dependent energy density $\rho_{\rm DE}$, in addition to modifying the evolution of perturbations propagating in the background. Early-time background modifications are well constrained by Big Bang nucleosynthesis (see e.g. Copi et al. (2004)), and the late-time background modifications are usually constrained using simple parametrizations for the equation of state of dark energy, such as $w_{\rm DE}=w_{0}+(1-a)w_{a}$, with $w_{0}$ and $w_{a}$ two free constants. We recover the standard $\Lambda$CDM model with $w_{0}=-1$ and $w_{a}=0$. The tightest constraints come from combined CMB, supernova, and baryon acoustic oscillation (BAO) measurements yielding, $w_{0}=-0.961\pm 0.077$ and $w_{a}=-0.28^{+0.31}_{-0.27}$ Aghanim et al. (2018). Since only small non-vanishing $w_{a}$ are allowed, from now on we will assume that the background evolution is exactly $\Lambda$CDM. This has the added benefit of disentangling effects coming from the modified background and those coming from the evolution of the perturbations.

Nonetheless, even if the background evolution of the universe is unmodified, the evolution of linear cosmological perturbations can still differ from $\Lambda$CDM. It is customary to encode the modifications determining the evolution of the energy-density matter perturbations in two parameters, which affect the evolution of the two Newtonian potentials $\Phi$ and $\Psi$ in the metric

[TABLE]

which describes the modifications relative to the FRW solution. The first parameter is an effective Newton’s strength $G_{\text{eff}}$ such that the Poisson equation becomes

[TABLE]

where $\bar{\Box}$ is the d’Alembert operator using the covariant derivatives of the background metric (1), $\Delta_{m}=\delta\rho_{m}/\rho_{m}+3Hv_{m}$ is the comoving gauge-invariant matter perturbation, $\delta\rho_{m}$ is the energy-density perturbation and $v_{m}$ is its velocity potential. The effective Newton’s strength depends on time only for sub-horizon perturbations with wavenumber $k>aH$.

The second parameter, $\gamma$, describes an effective anisotropic stress which makes the two metric potentials differ:

[TABLE]

The parameter $\gamma$ depends on time only for sub-horizon scales as well. Matter and metric perturbations behave in the same way as general relativity when $\gamma=1$ and $G_{\text{eff}}=G_{N}$, where $G_{N}$ is Newton’s constant. In scalar-tensor theories, it can be seen that the running of the Planck mass causes dark energy to cluster, leading to $G_{\text{eff}}\not=G_{N}$, and also induces anisotropic stress, making $\gamma\not=1$ (see e.g. Bellini and Sawicki (2014); Lagos et al. (2018)).

These two parameters are independent of the background evolution, and can differ from GR even when $w_{\rm DE}=-1$. In this unmodified background, we can therefore still obtain different predictions, for instance, for the CMB temperature anisotropies, or large scale galaxy distributions. If we knew exactly the gravitational theory leading to these modifications, we could calculate exactly how $\alpha_{M}$ affects matter perturbations and use EM cosmological data to constrain the running Planck mass (e.g. Bellini et al. (2016); Noller and Nicola (2018) for scalar-tensor theories). However, in this paper we will remain agnostic about the underlying gravity theory, and instead constrain $\alpha_{M}$ using gravitational wave data.

Solar System: The effects of a running Planck mass at Solar system and laboratory scales have been studied in Williams et al. (2004); Hellings et al. (1983); Bonanno and Froehlich (2017). The relevant modification is a time varying Newton’s constant, which affects, for instance, the period and radius of planetary orbits. However, to remain consistent with the remarkable agreement between GR and observations in this regime, modified gravity theories come equipped with a “screening mechanism” that hides all modified gravity effects in dense regions, where one recovers $M_{*}=M_{P}$ (see e.g. Jain and Khoury (2010); Joyce et al. (2015) for reviews). The mechanism can be due to the additional field acquiring a large mass in dense environments, and effectively mediating an undetectable short-range force (chameleon mechanism); due to changing its coupling with matter in this regime, becoming negligible (symmetron mechanism); or having dominant non-linear kinetic terms effectively produce a negligible coupling to matter (Vainshtein mechanism), among other possibilities. These screening mechanisms are expected to act in the Solar System as well as clusters of galaxies. Currently, observations of the Solar System and laboratory experiments constrain variations of the Planck mass to be of order $10^{-3}$ Bonanno and Froehlich (2017); Williams et al. (2004).

Gravitational Waves: The presence of an additional dark energy field may affect the original emission of the waveform from binary mergers as well as GW propagation. In this paper, we will assume that the additional gravitational degrees of freedom will not modify the emitted waveform, and hence we can use GR to predict the emission of the waveform. Consistency checks must be performed in the future for specific dark energy models to make sure that this is the case, as it has been previously shown that even when stationary black hole solutions may not be affected by the dark energy field, dynamical situations may excite it and leave an imprint Tattersall et al. (2018). Under our assumptions, we can still use compact binaries with EM counterparts as standard sirens.

Regarding the propagation of gravitational waves, modifications can occur even if the original emitted waveform is the same as in GR. In particular, as shown in eq. (6) and considered previously in Pettorino and Amendola (2015); Lombriser and Taylor (2016); Nishizawa (2018); Arai and Nishizawa (2018); Amendola et al. (2017); Linder (2018), a running Planck mass can affect GW propagation. In this paper, we will assume that the propagation can be modified inside and outside galaxies, where screening will not be active for GWs and hence the effects of the dark energy field become relevant. However, we will assume that today, in our galaxy, the effective Planck mass is given by $M_{P}$ (although the present-day value of $\alpha_{M}$ need not vanish). Multi-messenger observations allow us to place independent constraints on the running of the Planck mass Lombriser and Taylor (2016); Arai and Nishizawa (2018); Amendola et al. (2017); Linder (2018) as well as constraints on the present Hubble constant $H_{0}$ (see Belgacem et al. (2018) for a specific model).

III Standard sirens

In this section, we show explicitly how a running Planck mass affects standard siren measurements. Since we take the running of the Planck mass to be the result of a new dark energy field, we expect $\alpha_{M}$ to change over cosmological time scales and only affect the late-time universe. We will consider a specific parametrization for $M_{*}(t)$ satisfying these requirements, and discuss in detail its consequences.

Following the approach in Belgacem et al. (2018); Amendola et al. (2017), we start by canonically normalizing the field $h_{A}$ in eq. (5), and obtain the following equation of motion for the propagation of the two GW polarizations:

[TABLE]

Here, $a_{\text{\rm GW}}$ is an effective scale factor given by $a_{\text{\rm GW}}(z)=a(z)(M_{*}(z)/M_{P})$, and $\hat{h}_{A}=a_{\text{\rm GW}}h_{A}$ is the canonically normalized amplitude of GWs. For very small wavelengths, namely for $k^{2}\gg a_{\rm GW}^{{}^{\prime\prime}}/a_{\rm GW}$ with $a_{\rm GW}^{{}^{\prime\prime}}/a_{\rm GW}=[2\alpha_{M}^{{}^{\prime}}\mathcal{H}+2\mathcal{H}^{\prime}(1+\alpha_{M})+\mathcal{H}^{2}(2+\alpha_{M})^{2}]/4$,111Note that this condition limits how large $\alpha_{M}$ can be, and how quickly it can evolve in time. Typically, we will consider evolution over cosmological times, and $\alpha_{M}\ll k^{2}/\mathcal{H}^{2}$. the solution for $\hat{h}_{A}$ simply is a plane wave with constant amplitude. This means that the original metric perturbation $h_{A}$ is also a plane wave with a decreasing amplitude due to the factor of $1/a_{\rm GW}$. In GR, the amplitude would simply decay as $1/a$, so we can interpret $a_{\rm GW}$ as the effective scale factor felt by GWs due to the combined effects of the background metric and the background’s additional degree of freedom. In this case, the present-day observed amplitude $h_{A}^{o}$ is given by

[TABLE]

where $z$ is the redshift of the source and $h_{A,{\rm GR}}^{o}$ is the expression for the observed amplitude in GR; $h_{A,{\rm GR}}^{o}$ decays as $1/d_{L}$, with

[TABLE]

the luminosity distance. $\hat{H}(z)=H(z)/H_{0}$ is the normalized Hubble rate, which is explicitly

[TABLE]

We have introduced the fractional energy-density parameter $\Omega=\rho/(3H^{2}M^{2})$, and the critical energy density $\rho_{c}=3H_{0}^{2}M_{P}^{2}$. We have also used the fact that $a=(1+z)^{-1}$, with $a=1$ today.

The omitted proportionality factor in eq. (14) characterizes the emitted waveform (which is a function of the GW frequency, chirp mass, the effective Newton’s constant felt by the compact binary, and the equation of state for neutron stars), and we assume that it is the same as in GR.

We note that eq. (14) depends only on the value of $M_{*}$ at the source and observation points, and not on the intermediate evolution. This is because we can only measure the cumulative change in the amplitude compared to GR, which depends exclusively on the initial and final values of $M_{*}$. In keeping with the previous section’s discussion, we assume that $M_{*}(z=0)=M_{P}$, and that $M_{*}$ depends only on time. As a consequence, the result depends only on the distance to the source galaxy, not on the properties of the host galaxy itself. Furthermore, we assume that the effective Planck mass evolves in the same way as the cosmological one, and hence it does not depend on whether we are inside a galaxy or not.

In order to make a concrete estimation of the effect of the running Planck mass, we need to assume a specific functional form for $M_{*}$. Since we are interested in modified-gravity models of dark energy, a common time parametrization for $\alpha_{M}$ is Bellini et al. (2016); Alonso et al. (2017)

[TABLE]

where $c_{M}$ is a free constant parameter and $\Omega_{\rm DE}$ is the fractional dark energy density. So far the background has been kept arbitrary, but we now assume it to be given by a fiducial $\Lambda$CDM expansion history, and thus from now on we assume that $\rho_{\rm DE}$ is constant, and set $\Omega_{{\rm DE},0}$ to the best-fit Planck value Aghanim et al. (2018). Note that the parametrization in eq. (17) is not well-suited to $f(R)$ models Linder (2017), and hence different time evolutions may also be of interest (see for instance Kreisch and Komatsu (2018)). In any case, simple parametrizations with a few free constants are sufficient for the time being, as observational data does not have the constraining power to test more complicated functions Gleyzes (2017).

This parametrization assumes that there are no modified gravity effects at early times, but at late times new degrees of freedom come into play and modify the evolution of the universe through $\Omega_{\rm DE}(z)$. Other parametrizations with early-time modifications can also be considered (see for instance Alonso et al. (2017)).

For the specific case of scalar-tensor theories, it is known how $\alpha_{M}$ affects the effective Newton constant $G_{\text{eff}}$ and the effective anisotropic stress $\gamma$, and thus cosmological data has been used to constrain $c_{M}$ to the range $-0.62<c_{M}<+1.35$ at $95\%$ CL Noller and Nicola (2018), which means $c_{M}$ is allowed to be of order unity. Thus, in contrast to constraints on the propagation speed of GWs, this kind of modified gravity effect is not yet ruled out. Forecasts for this same parametrization were studied in Alonso et al. (2017), where it was found that the 1$\sigma$ uncertainties on the parameter $c_{M}$ may improve by a factor of 5 when taking into consideration future photometric redshift surveys such as LSST Abell et al. (2009) and SKA Carilli and Rawlings (2004), as well as the Stage IV CMB experiment Abazajian et al. (2016).

Regardless of the underlying modified gravity theory causing $\alpha_{M}$, from standard sirens we can also place constraints on $c_{M}$ by measuring the difference between the luminosity distance and the GW distance. As shown in Belgacem et al. (2018), according to eq. (14) we have that:

[TABLE]

and for the specific parametrization (17), we explicitly obtain

[TABLE]

where we have used the assumption that $\rho_{\rm DE}$ is constant. This ratio describes the cumulative difference in the GW amplitude in GR compared to that in modified gravity, if the wave had the same emitted amplitude. In order to illustrate the effect of $c_{M}$ on this ratio , we Taylor expand this expression for low redshifts $z\ll 1$ to find

[TABLE]

From here we explicitly see that if $c_{M}=0$, then $d_{\rm GW}/d_{L}=1$, and the leading-order correction is proportional to $z$, hence no considerable modifications are expected for low-redshift events.

When estimating $H_{0}$ from standard sirens, we first obtain $d_{\rm GW}$ from the waveform and the redshift $z$ from the EM counterpart. Then, we combine eqs. (19) and (15) to obtain $H_{0}$. At low redshifts, we have that

[TABLE]

and thus the estimates of $H_{0}$ in GR and modified gravity would differ by a factor $z^{2}c_{M}/(2d_{\rm GW})$ for a given value of $d_{\rm GW}$. We note that the sign of $c_{M}$ determines whether the measured $H_{0}$ in the modified gravity model will be larger or smaller than the GR value. In particular, if $c_{M}<0$, then $c_{M}$ contributes with a negative term to $H_{0}$, yielding a smaller $H_{0}$, and vice versa.

Going beyond the low-redshift limit, Figure 1 shows the evolution of the fractional difference $(d_{\rm GW}-d_{L})/d_{\rm GW}$ as a function of $d_{\rm GW}$ for different values of $c_{M}$ (blue, green and orange lines), up to a GW distance of $1.5$ Gpc. For comparison, the GW distances for BNSs detected by aLIGO at design sensitivity and A+ are expected to follow the distributions in Figure 2 under the assumptions discussed in Section IV. The ratio $(d_{\rm GW}-d_{L})/d_{\rm GW}$ illustrates the fractional cumulative difference in the GW amplitude that would be detected in GR compared to that detected in modified gravity, for a source at a given distance with the same emitted amplitude. As expected, the larger the $|c_{M}|$, the larger the fractional difference between the GW distance and $d_{L}$, indicating a larger deviation from GR. We note that for values of $c_{M}\gtrsim-3$, the ratio flattens out at large distances due to the fact that our $\alpha_{M}$ in eq. (17) decreases with redshift (or, equivalently, with distance), and hence the accumulated difference between GR and modified gravity becomes negligible at large distances. The distance at which a given curve starts flattening out depends on the value of $c_{M}$ and the cosmological parameters. Explicitly, from eq. (19), we find that the maximum fractional difference is given by:

[TABLE]

Using best-fit Planck cosmological parameters, we find that if $c_{M}=1$, for instance, then the maximum fractional difference will be about $32\%$. However, for sources at aLIGO or LIGO A+ horizons (dotted vertical lines, corresponding to the largest detectable distances in Figure 1), the fractional difference will be about $4\%$ or $7\%$, respectively, for $c_{M}=1$.

Due to the aforementioned behavior, we conclude that more distant events have a priori more constraining power on $c_{M}$ than nearby events; however, this effect starts to diminish for events at distances much greater than a few Gpc. We note that the measurement uncertainties in $d_{\rm GW}$ tend to grow with distance, as the signal to noise ratio (SNR) of a GW signal scales inversely with the source distance. The range of expected $1\sigma$ uncertainties in the measured GW distance is shown as solid colored regions in Fig. 1: yellow for aLIGO and pink for LIGO A+, under the simplified assumptions discussed in Section IV. The boundary closest to zero of each colored region corresponds to the best sources, which produce the optimal SNR and the smallest distance uncertainties, and thus the tightest constraints on $c_{M}$. In the above discussion, we assumed that the uncertainty in the EM distance, $d_{L}$, is negligible, and that all the uncertainty in the fraction $(d_{\rm GW}-d_{L})/d_{\rm GW}$ comes from $d_{\rm GW}$. Indeed, for a fixed cosmology, the uncertainty on $d_{L}$ comes exclusively from the peculiar velocity of the source, which is typically around $150$–$250$ km/s at all distances, and therefore for a fixed background cosmology the fractional error of the luminosity distance $\sigma_{d_{L}}/d_{L}$ decreases with distance.

From Figure 1, we see that in order to distinguish GR from a nonzero $c_{M}$, we need to measure $d_{\mathrm{\rm GW}}$ to a precision better than the deviation caused by the nonzero value of $c_{M}$. For a single event, this is only possible for the most extreme values of $c_{M}$, because we will rarely get an event with a $1\sigma$ distance uncertainty that is comparable to the deviation caused by $-1<c_{M}<1$. However, by combining a population of events, it will be possible to place tight constraints on $c_{M}$, as discussed in Section V.2. A population of BNS events detected by LIGO A+ is especially promising, because the typical source detected by A+ at a given distance has a much higher SNR, and therefore a smaller relative distance uncertainty, than the typical source detected by aLIGO. However, it becomes less useful to extend the detection horizon much beyond A+, because the effect of nonzero $c_{M}\sim\mathcal{O}(1)$ starts to saturate at $\sim 1.5$ Gpc (additionally, it is increasingly difficult to find counterparts and identify host galaxies for events at high redshift).

Finally, we recall that we have a schematic relationship between $d_{\rm GW}$ and redshift $z$ of the form

[TABLE]

where $R$ is a function corresponding to the ratio of $d_{\rm GW}/d_{L}$ given on the righthand side of eq. (19). It is clear that, given cosmological parameters ($H_{0}$, $\Omega_{m,0}$), we can constrain $c_{M}$ by measuring $d_{\rm GW}$ and $z$. We emphasize that if the cosmological parameters are fixed by some cosmological data, then the resulting constraint on $c_{M}$ is not independent of these cosmological data sets. However, if the cosmological parameters $H_{0}$ and $\Omega_{m,0}$ are taken to be free constants, then eq. (23) can be used to find joint constraints on $c_{M}$, $H_{0}$ and $\Omega_{m,0}$ which are completely independent from the cosmological data sets. Note that the constraints would still depend on the cosmological model (assumed to be $\Lambda$CDM here), but not on the best-fit values from external data. If $c_{M}=0$, then $R=1$, and we can use this relationship to find independent constraints on the cosmological parameters, as has been done for $H_{0}$ in Abbott et al. (2017).

While for low-redshift sources $\Omega_{m,0}$ has a negligible effect and the only relevant cosmological parameter is $H_{0}$, for high-redshift sources $\Omega_{m,0}$ does become relevant and must be taken into consideration. For simplicity, in the rest of the paper we will always fix $\Omega_{m,0}=0.315$, the best-fit value from Planck 2018 Aghanim et al. (2018), while keeping $H_{0}$ free. In this sense, the constraints that we quote are not fully cosmology-independent, although we have checked that if we free $\Omega_{m,0}$ and use a flat prior with a $6\sigma$ width around the best-fit Planck value, i.e.  $0.330<\Omega_{m,0}<0.372$, our results are unaffected, as the uncertainty on $\Omega_{m,0}$ is subdominant, affecting the distance-redshift relation to less than $1\%$ over the redshift range of interest ($z\lesssim 0.1$ for aLIGO and $z\lesssim 0.2$ for LIGO A+) for BNSs. For this reason, in the rest of the paper we will focus on the joint constraints for $c_{M}$ and $H_{0}$ only, and these will be referred to as cosmology-independent constraints. In particular, we will analyze the data from GW170817 as well as forecasts for advanced LIGO and LIGO A+.

IV Method

In this section, we describe the method used for our standard siren inference of $H_{0}$ in the context of a running Planck mass. Throughout the analysis, we fit only for $c_{M}$ and $H_{0}$, assuming that the other cosmological parameters in the $\Lambda$CDM background (namely, $\Omega_{k}$, $\Omega_{m}$ and $\Omega_{\Lambda}$), collectively denoted by $\Xi$, are known to a few percent. We note that $c_{M}$ is not expected to correlate with these other cosmological parameters in their impact on CMB observables (Huang, 2016), and it is therefore self-consistent to fix the background cosmology to the Planck 2018 values while measuring $c_{M}$ with standard sirens. Allowing these parameters to vary by up to 10% from their best-fit Planck 2018 values has a $\lesssim 1\%$ effect on the distance-redshift relation over the detectable redshift range, and so even if the current measurement errors were several times larger, marginalizing over the extra uncertainty would have a negligible impact on our results. Within the current Planck 2018 uncertainties, the distance-redshift relation only varies by $<0.05\%$ over the detectable redshift range. In the future, however, especially if standard sirens with counterparts are detectable to much higher redshifts $z>1$ by e.g. LISA, one could carry out a joint CMB-standard siren analysis that would incorporate all cosmological parameters.

GW Measurements: We assume that the GW measurement uncertainty of distance scales as $1.8/\rho$ (at $1\sigma$) where $\rho$ is the single-detector SNR of the source. We assume a threshold SNR of $\rho_{\mathrm{th}}=8$ for detection. From Fisher matrix arguments, we expect the GW distance to scale inversely with the SNR, with the proportionality factor $>1$ because of the distance-inclination degeneracy Cutler and Flanagan (1994). We choose the $1\sigma$ uncertainty of $1.8/\rho$ to match the expected $H_{0}$ convergence rate of $(13\%$–$15\%)$/$\sqrt{N}$, where $N$ is the number of GW detections Chen et al. (2017). For simplicity, when simulating a mock population of sources, we assume that the GW distance likelihood is approximated by a Gaussian distribution:

[TABLE]

where $N_{[\mu,\sigma]}$ denotes the standard normal distribution with mean $\mu$ and standard deviation $\sigma$. While this is not a realistic approximation for individual sources, when combining tens to hundreds of detections, it yields the expected convergence rate for the recovered cosmological parameters.

Likelihood: We denote the GW data by $x_{\mathrm{\rm GW}}$ and the EM data by $x_{\mathrm{EM}}$. The likelihood of the GW data depends on the source’s GW distance $d_{\mathrm{\rm GW}}$, sky position $\omega$, inclination $\iota$, and all other parameters of the signal, including its redshift $z$ (which affects the frequency of the waveform), and the source-frame masses, spins, etc., which we collectively denote by $\xi$. The likelihood given the extrinsic parameters (sky localization and inclination) is largely independent of the intrinsic parameters Singer and Price (2016); Chen and Holz (2017), and we marginalize over these other parameters to get the GW likelihood given its distance and sky position:

[TABLE]

Meanwhile, the likelihood of the EM data depends on the host galaxy’s sky position and luminosity distance, which is related to its cosmological redshift (the redshift it would have if it were in the Hubble flow, corrected for any peculiar velocities) by the standard $\Lambda$CDM relation.

We put a prior $p(z,\omega)$ on the redshifts and sky positions of the host galaxies. We choose this prior to match a merger rate density that is isotropic and roughly follows the star-formation rate (see Appendix). To avoid a biased measurement of the cosmological parameters, this prior must match the true redshift distribution of the host galaxies; however, we find that any likely deviation from a uniform merger rate density, including a merger rate that traces the star-formation rate, is too small to cause a noticeable bias. Similarly, deviations from an isotropic distribution of sources on the sky (for example, due to large-scale structure) are largely irrelevant to this analysis unless there are significant correlations between the underlying distribution of sources on the sky and the antenna power patterns of the detectors. Moreover, a significant deviation from the assumed merger rate density will be easily detectable with hundreds of host galaxies with well-measured redshifts and sky positions, and can be used to update our prior accordingly.

The relationship between $d_{L}$ and $d_{\mathrm{\rm GW}}$ is given in eq. (19), where $d_{L}$ is given by $z$, $H_{0}$, and the standard $\Lambda$CDM parameters $\Xi$. We denote the function that returns $d_{\mathrm{\rm GW}}$ given $z$, $H_{0}$, $c_{M}$ and $\Xi$ by $\hat{d}_{\mathrm{\rm GW}}$.

The likelihood for the data given $c_{M}$ and $H_{0}$ is

[TABLE]

where $\beta(H_{0},c_{M})$ ensures that the likelihood integrates to unity over detectable data sets, and accounts for selection effects in the GW detection and measurement process (see Appendix).

We take $p(\Xi)$ to be the posterior on these parameters from Planck 2018, although, as noted earlier, we can approximate these as being measured exactly and given by a $\delta$-function centered on their mean values. We also assume that the sky position is measured exactly, and that the redshift uncertainty is small (i.e. spectroscopic redshifts, so that the only significant source of uncertainty in the cosmological redshift is in the peculiar velocity correction, typically around 200 km/s (Wu and Huterer, 2017)). For a population where the majority of detected sources are at redshifts $z\gtrsim 0.05$ and the redshift uncertainty is subdominant to the $d_{\mathrm{\rm GW}}$ uncertainty, we can therefore approximate the EM likelihood term, $p(x_{\mathrm{EM}}\mid z,\omega)$ by a $\delta$-function centered at the true redshift and sky position:

[TABLE]

When analyzing GW170817, a very nearby event at redshift $z\sim 0.01$, we include the peculiar velocity uncertainty in the calculations, taking the EM likelihood to be:

[TABLE]

where $z_{\mathrm{obs}}$ is the observed, peculiar-velocity corrected redshift, $\sigma_{z}$ is the uncertainty, and $\omega_{\mathrm{obs}}$ is the observed sky position. For this analysis, we also choose priors that match the default priors in Abbott et al. (2017): $p(d_{\mathrm{\rm GW}})\propto d_{\mathrm{\rm GW}}^{2}$ and $p(H_{0})\propto 1/H_{0}$.

V Results

In this section we analyze the one BNS detection so far, GW170817, as well as a simulated population of BNSs detected by aLIGO at design and A+ sensitivities. In each case, we present joint constraints on $c_{M}$ and $H_{0}$ and highlight the correlations between them. Section V.1 discusses constraints from the single event GW170817, detected at $d_{\mathrm{\rm GW}}\sim 40$ Mpc. Then, in Section V.2 we discuss forecast constraints from populations of standard sirens detected with design-sensitivity LIGO and A+.

V.1 Single Event

In this section we study joint constraints on $c_{M}$ and $H_{0}$ from a single multimessenger signal. We start by considering GW170817, which has a measured222We use the publicly available posterior samples released with Abbott et al. (2018) and available at dcc.ligo.org/LIGO-P1800061/public. GW distance of $d_{\mathrm{\rm GW}}=41^{+4}_{-7}$ Mpc and a “Hubble” velocity of $v_{H}=3017\pm 166$ km/s Abbott et al. (2017). Combining the GW distance samples with the redshift of the host galaxy, we can calculate the joint posterior on $H_{0}$ and $c_{M}$. If we assume an $H_{0}$ prior given by Planck (2018) The Planck Collaboration (2006), we find $c_{M}=-9^{+21}_{-28}$ (68.3% credible interval); alternatively, taking a $H_{0}$ prior given by SHoES (2018) Riess (2006) gives $c_{M}=8^{+21}_{-30}$. The $c_{M}$ posterior under each assumption is shown in Figure 3. The 95% credible interval under the Planck (SHoES) $H_{0}$ prior is $c_{M}=-9^{+37}_{-72}$ ($c_{M}=8^{+39}_{-74}$).

We see that even though this event had a high SNR of 32.4, since it was very close-by ($\sim 40$ Mpc), the constraints on $c_{M}$ are very broad. As Figure 1 shows, at 40 Mpc, a large range of $c_{M}$ values would produce GW distances that are consistent with the $\sim 15\%$ distance uncertainty from GW170817. As a comparison, we mention that constraints on $c_{M}$ have been obtained for scalar-tensor theories from cosmological data, where it was found that $-0.62<c_{M}<+1.35$ at $95\%$ CL Noller and Nicola (2018). Therefore, current GW constraints allow for $c_{M}$ of $\mathcal{O}(10)$, whereas cosmological data require $c_{M}$ of $\mathcal{O}(1)$.

However, future detector networks with improved sensitivity, such as A+, will detect events out to much higher distances with comparable measurement uncertainties. For example, a single event detected by A+ at $d_{\mathrm{\rm GW}}=400$ Mpc with a $1\sigma$ distance uncertainty of 15% would constrain $-4<c_{M}<4$ for a true $c_{M}=0$, assuming the Planck $H_{0}$ measurement as a prior.

We can also obtain constraints on $c_{M}$ independently of other data sets by using an uninformative prior on $H_{0}$. In this case, we find a strong positive correlation between $c_{M}$ and $H_{0}$, as shown in Fig. 4.

This correlation arises because a given pair ($z$, $d_{\rm GW}$) can also be achieved in a universe where $H_{0}$ is larger (and hence the source is closer) but the friction term $c_{M}$ is correspondingly larger too (and hence the amount of amplitude decay during its travel is larger). Equivalently, the same data could also be fitted by a universe with a smaller $H_{0}$ and a smaller $c_{M}$.

The posterior probability of $H_{0}$, with $c_{M}$ marginalized over a flat prior in the range $c_{M}\in[-150,150]$ and a flat-in-log prior on $H_{0}$, is shown in Fig. 5. In this case, the constraints become $H_{0}=76^{+53}_{-28}$ km/s/Mpc ($68.3\%$ highest density posterior interval). As a comparison, we also show the constraint on $H_{0}$ found by marginalizing over a narrow $c_{M}$ prior, $-2<c_{M}<2$, and assuming $c_{M}=0$. In both of these cases, we find $H_{0}=70^{+13}_{-7}$ km/s/Mpc, in agreement with Abbott et al. (2018).

In this case, we find that for this one event the uncertainties on $H_{0}$ grow by a factor of $>2$ when marginalizing over a very broad prior on $c_{M}$, but are unaffected by a more reasonable prior $-2<c_{M}<2$.

V.2 Population

We now consider a population of BNS events detected by aLIGO at design and A+ sensitivities, with EM counterparts that allow us to identify unique redshifts of the host galaxies. As the GW network sensitivity reaches design sensitivity for aLIGO and later upgrades to A+, a GW event at a given distance will be detected with higher SNR and yield a better-constrained $d_{\mathrm{\rm GW}}$ measurement, meaning that although GW170817 is only sensitive to $c_{M}\sim\mathcal{O}(50)$, a single event detected by aLIGO at design (A+) sensitivity will typically constrain $c_{M}$ to a $1\sigma$ width of $\lesssim 10$ ($\lesssim 5$). In Fig. 6 we show the joint and marginalized constraints on $c_{M}$ and $H_{0}$ for a simulated population of 250 BNSs detected by aLIGO, where the injected values are $c_{M}=0$ and $H_{0}=67.4$ km/s/Mpc (indicated by solid black lines). We assume flat priors on both $c_{M}$ and $H_{0}$.

The bottom left panel shows the joint posterior probability on $c_{M}$ and $H_{0}$, with the contours indicating 90% and 50% credibility levels. As expected, there is a positive correlation between $H_{0}$ and $c_{M}$, which leads to a broader recovered posterior on $H_{0}$ when marginalizing over $c_{M}$, as opposed to fixing $c_{M}=0$ (the correct value in this case). This is shown in the top left panel. When fixing $c_{M}=0$, the $1\sigma$ constraints on $H_{0}$ scale roughly as $(13$–$15\%)/\sqrt{N}$ (for aLIGO and A+), giving a $1\sigma$ interval of $0.6$ for 250 events (green dotted line); however, when marginalizing over a completely uninformative prior on $c_{M}$, the same number of events yields a $1\sigma$ interval that is twice as broad (solid blue line). This implies that with no external knowledge of $c_{M}$, it would take four times as many events to reach the same precision in the standard siren $H_{0}$ measurement, or around 200 events to reach $2\%$ as opposed to only $\sim 50$ events if fixing $c_{M}=0$ Chen et al. (2017).

Meanwhile, the bottom right panel shows the posterior probability of $c_{M}$ when marginalizing over $H_{0}$ with a flat, broad prior (solid blue line), as well as when fixing the $H_{0}$ prior to the Planck (2018) posterior (green dotted line). In the former case, taking the Planck $H_{0}$ measurement as a prior, this realization gives $c_{M}=0.16^{+0.58}_{-0.60}$ ($68.3\%$ credible interval), whereas the latter case, which assumes no external measurement of $H_{0}$, gives $c_{M}=0.35^{+1.08}_{-1.10}$. Including an external constraint on $H_{0}$ reduces the uncertainties on $c_{M}$ by almost a factor of 2.

Although Fig. 6 shows only a single realization of 250 simulated BNSs detected by design-sensitivity aLIGO, we find that the expected constraints on $c_{M}$ and $H_{0}$ and the $1/\sqrt{N}$ scalings are typical across many realizations for aLIGO. Generically, for aLIGO, we find constraints on $c_{M}$ with a $1\sigma$ width that scale roughly as $\sim 9.3/\sqrt{N}$ for an informative $H_{0}$ prior, and $\sim 16/\sqrt{N}$ for a flat $H_{0}$ prior. In particular, with 100 BNS events detected by aLIGO, assuming that $H_{0}$ is obtained from external information (such as from the cosmic microwave background), we would find that $|c_{M}|\lesssim 0.9$. This number is comparable to the current constraints for scalar-tensor theories obtained from cosmological observations Noller and Nicola (2018). Similarly, with A+ sensitivity, we find that the same number of events yields constraints on $c_{M}$ that are tighter by a factor of 2, constraining $c_{M}$ with a $1\sigma$ width that scales roughly as $\sim 4.7/\sqrt{N}$ for an informative $H_{0}$ prior, and $\sim 9.5/\sqrt{N}$ for a flat $H_{0}$ prior. In this case, we find that 100 BNS events would allow us to get a limit $|c_{M}|\lesssim 0.5$ when assuming $H_{0}$ is known.

For the previous example, we chose the true $c_{M}=0$ (i.e. assumed that GWs propagate as predicted by general relativity), and so by fitting a model with an uninformative prior for $c_{M}$ and $H_{0}$, we recovered the true values with a larger uncertainty than if we had assumed GR and fixed $c_{M}=0$. However, if the true $c_{M}\neq 0$ but GR is assumed in the usual standard siren analysis, we will recover a biased $H_{0}$ measurement. Due to the previously shown positive correlation between $H_{0}$ and $c_{M}$, if the true $c_{M}>0$, the $H_{0}$ measurement will be biased to low values if falsely assuming $c_{M}=0$, and if $c_{M}<0$, the $H_{0}$ measurement will be biased to large values.

As an example, Fig. 7 shows the joint and marginalized constraints on $c_{M}$ and $H_{0}$ for a simulation in which we injected $c_{M}=1.35$. We emphasize that this value of $c_{M}$ is currently allowed within the 95% confidence interval inferred by Noller and Nicola (2018) for scalar-tensor theories of gravity. For this simulation, we take 100 mock BNS events detected by LIGO A+.

As in the previous figure, in the bottom left panel we show the joint posterior probability, where the contours indicate 50% and 90% levels. The bottom right panel shows the posterior probability of $c_{M}$ marginalizing over an uninformative prior on $H_{0}$ (solid blue line), as well as when adopting the Planck (2018) measurement as the prior on $H_{0}$ (green dotted line). In the former case, it is found that $c_{M}=1.3\pm 0.8$ ($68.3\%$ credible interval), whereas in the latter case $c_{M}=1.5\pm 0.5$.

In the top left panel, we see that the posterior probability of $H_{0}$ is significantly biased away from its true value if we falsely assume $c_{M}=0$ (green dotted line). If we properly marginalize over a flat $c_{M}$ prior, we find $H_{0}=66.9^{+1.6}_{-1.5}$ km/s/Mpc, whereas fixing $c_{M}=0$ gives $H_{0}=64.8^{+0.8}_{-0.8}$ km/s/Mpc. In the latter case, the true $H_{0}$ is outside the $99\%$ credible interval. In other words, incorrectly assuming $c_{M}=0$ yields an $H_{0}$ measurement that is biased by more than $3\sigma$ with only 100 events detected by A+. In general, if the true value is positive ($c_{M}>0$) and we fit a model with $c_{M}=0$, we bias $H_{0}$ towards lower values than the true one. Conversely, if the true value is negative ($c_{M}<0$), we bias $H_{0}$ towards higher values.

VI Discussion

In this paper we describe a modification to GR which impacts the propagation of GWs. This extension corresponds to a possible running of the Planck mass, which we describe with one free parameter $c_{M}$, where $c_{M}=0$ represents GR with a constant Planck mass. This modification affects the friction of GW amplitudes when they propagate through a homogeneous and isotropic universe, and also affects the evolution of matter perturbations in different ways depending on the theory of interest. Here we focus on GW standard siren measurements, studying how a $c_{M}\neq 0$ modification is degenerate with the value of the local Hubble expansion rate, $H_{0}$. We also explore the ability to constrain both these quantities with future standard siren events detected by LIGO at design and A+ sensitivities.

Studying the event GW170817, we find that if we include external cosmological data, namely the Planck 2018 $H_{0}$ posterior, then we find $-81<c_{M}<28$ at $95\%$ credibility. This constraint is very weak compared to existing constraints from cosmological data when a specific modified gravity theory is considered. Since $c_{M}$ affects matter perturbations, it leaves potentially detectable imprints on CMB and structure formation. For Horndeski theories, current cosmological constraints give $-0.62<c_{M}<+1.35$ at $95\%$ CL Noller and Nicola (2018). Hence, current GW constraints allow for $c_{M}$ of $\mathcal{O}(10)$, whereas cosmological data require $c_{M}$ of $\mathcal{O}(1)$. On the other hand, for cosmology-independent results, we find that $c_{M}$ and $H_{0}$ are highly degenerate, and constraints on $H_{0}$ are degraded from $H_{0}=70^{+12}_{-8}$ km/s/Mpc (when $c_{M}=0$, that is, when GR is assumed to be correct) to $H_{0}=76^{+53}_{-28}$ km/s/Mpc (when marginalizing over a very broad $c_{M}$ prior, $-150<c_{M}<150$).

In addition, we consider populations of events and discuss future forecasts for standard sirens with aLIGO and LIGO A+. We find that 100 BNSs detected by A+ with identified EM counterparts can lead to cosmology-independent constraints on $H_{0}$ with an accuracy of $\sim 3\%$, and on $c_{M}$ with $\sigma(c_{M})\sim 0.9$ (or $\sigma(c_{M})\sim 1.6$ with 100 aLIGO detections). From these results we estimate the need for $400$ detections in order to obtain a $\sim 1\%$ constraint on $H_{0}$, or four times what is required if $c_{M}$ is known exactly to be zero. Furthermore, we find that $H_{0}$ and $c_{M}$ are highly degenerate, which highlights the importance of testing for the parameter $c_{M}$ to avoid biasing the inferred value of $H_{0}$ by wrongly assuming GR. In particular, we show that if we have a population of 100 events with $c_{M}=1.35$, then the inferred $H_{0}$, assuming $c_{M}=0$, will be $>3\sigma$ below the true value. In this case, the actual $H_{0}$ of the population may be ruled out at more than $99\%$ confidence due to the incorrect assumption that $c_{M}=0$. This result emphasizes the importance of testing the minimal assumptions of one’s models. Finding a bias of this magnitude could help arbitrate the current discrepancy between local and cosmological $H_{0}$ constraints.

It is important to discuss some caveats of our results and calculations. In general, all population results depend on the time-evolution of the background, assumed to be $\Lambda$CDM here. They also depend on the specific values of the cosmological parameters considered, and in some cases we also assumed $H_{0}$ to be fixed. One possible extension to the analysis made in this paper could involve a change in the background, for example changing $\Lambda$ ($w=-1$) to a more general form of dark energy, $w_{\rm DE}=w_{0}+w_{a}(1-a)$.

In addition, the numbers quoted here also depend on the parametrization adopted for $\alpha_{M}(t)$, which means that even if these constraints are found to be in tension with future measurements, we cannot conclude that all models with nonzero $\alpha_{T}$ are disfavored, but rather that the specific time behaviour assumed here, $\alpha_{M}(t)=c_{M}\frac{\Omega_{\rm DE}(z)}{\Omega_{{\rm DE},0}}$, is disfavored. This is why alternative parametrizations, such as the ones considered in Noller and Nicola (2018), must also be tested.

Furthermore, we have made the crucial assumption that the GW emission of the BNS source is the same in GR as in the modified gravity theory. One typically justifies this by arguing that some compact object solutions, such as those for black holes, are exactly the same as the solutions in GR. However, as shown in Tattersall et al. (2018), even if the stationary solutions in the strong-field regime of modified gravity look the same as in GR, the emission of GWs in dynamical environments can still differ. Another argument is that the extra gravitational field may be suppressed due to a screening mechanism in the intermediate and high-energy regimes. However, even if this is true, one should check that the suppressed effects are negligible given LIGO’s current measurement uncertainty. To date, GR waveforms have been found to describe all detected compact binary GWs. Nonetheless, it would be interesting to have analytical or numerical calculations (for instance for Horndeski theories) that allow us to estimate the size of modified gravity effects in the waveform, and determine whether they can be seen with LIGO A+ or next-generation detectors.

Finally, we highlight that given the current constraints on the propagation speed of gravitational waves $|c_{T}/c-1|\lesssim 10^{-15}$, a constraint on $\alpha_{M}$ would have a large impact on modified gravity theories, as this is the only possible effect that could be detected with GW data for second-order-derivative scalar-tensor theories, and the only remaining non-trivial extension to GR that can be achieved in this case. For the model in eq. (8), a constraint pointing to $\alpha_{M}\approx 0$ (i.e. no modification to GW propagation) would disallow all scalar-tensor interactions but those given by

[TABLE]

In the sub-horizon regime, these models always give $\gamma=1$, as well as $G_{\text{eff}}=G_{N}$ for structure formation when the coupling term $G_{3}$ is sufficiently small, and hence modified gravity would only affect scales near the cosmological horizon.

Acknowledgments

We would like to thank Ignacy Sawicki for useful discussions during the development of this work. ML, MF, PL, and DEH are supported by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation and its founder Fred Kavli. MF was also supported by the NSF Graduate Research Fellowship Program under grant DGE-1746045. PL was supported in part by the Natural Sciences and Engineering Research Council of Canada, and by NSF grant PHY-1505124. MF, PL, and DEH were supported by NSF grant PHY-1708081.

Appendix A Statistical Formalism

In this appendix, we give additional details regarding the statistical formalism summarized in Sec. IV. In order to arrive at the likelihood (eq. 26) for a given EM and GW dataset given values of $c_{M}$ and $H_{0}$, we follow the formalism of (Abbott et al., 2017; Chen et al., 2017; Pardo et al., 2018; Mandel et al., 2018) for incorporating measurement uncertainty and selection effects.

In the likelihood eq. 26, we include a term $\beta(H_{0},c_{M})$ to account for selection effects in the measurement process. This term is given by the integral of the numerator over all detectable EM and GW datasets (Mandel et al., 2018):

[TABLE]

where we assume that a given EM or GW dataset is detected if and only if it is above certain threshold. In reality, the detectability of an EM counterpart to a GW event may depend on details such as the inclination, masses, and apparent magnitude of the source, but for simplicity we assume that all BNS mergers detected in GWs will have an observed EM counterpart and an identified host galaxy. If such counterparts are similar to the kilonova associated with GW170817, their detection is certainly feasible with current telescopes for aLIGO sources (which will be at redshifts $z\lesssim 0.1$), and with future telescopes such as LSST for higher-redshift A+ sources Chen et al. (2017). We therefore assume that the integral over $x_{\mathrm{EM}}$ is independent of the other terms, and ignore it. However, if it becomes the case with future detections that only a subset of GW BNS events have identified host galaxies,this term must be modeled and incorporated into the likelihood.

For the GW selection effects, we assume that a BNS is detected if it produces a single-detector SNR $\rho>8$. When a real population of BNSs is detected in GWs, this assumption can be easily modified to consider the network SNR, or calibrated to injection campaigns in real data (The LIGO Scientific Collaboration and The Virgo Collaboration, 2018). As in (Chen et al., 2017), we define:

[TABLE]

(Recall that our assumed prior $p(z,\omega)$ is separable, $p(z,\omega)=p(z)p(\omega)$.) We evaluate the term $P_{\mathrm{det}}(d_{\mathrm{\rm GW}})$ with the procedure described in (Chen et al., 2017). In particular, we make the simplifying assumptions that the detectability of a GW waveform is independent of its redshift, which affects the observed frequency and therefore the SNR of the source, but only by a negligible amount for the redshifts $z\lesssim 0.2$ considered here. We assume that all BNS sources are nonspinning (the dimensionless spin is expected to be very small, $a<0.05$, for BNS sources) and 1.4–1.4 $M_{\odot}$ in mass. The mass distribution will, in general, affect the term $P_{\mathrm{det}}(d_{\mathrm{\rm GW}})$, since the SNR of a GW source is a strong function of the binary’s mass. However, to leading order the SNR depends only on the chirp mass, and so the term $P_{\mathrm{det}}(d_{\mathrm{\rm GW}})$ depends only on the underlying distribution of BNS chirp masses. The chirp mass of each source is measured extremely well, and so with $\mathcal{O}(100)$ sources, the distribution of chirp masses will be accurately determined and can be used to update the function $P_{\mathrm{det}}(d_{\mathrm{\rm GW}})$ used in standard siren analyses. Note that we also assume that any running of the Planck mass (nonzero $c_{M}$) only affects the amplitude of the signal and not the frequency evolution of the waveform; otherwise, the recovered masses would be affected.

We therefore have that the $\beta(c_{M},H_{0})$ term in the likelihood is given by:

[TABLE]

The prior on the redshift of the source, $p(z)$, enters into this equation. In general, to avoid a biased measurement, the prior $p(z)$ must match the true redshift distribution. In our simulations, we assume that the underlying redshift distribution matches a merger rate that roughly traces the low-redshift star-formation rate:

[TABLE]

where $V_{c}$ is the comoving volume, $(1+z)^{2.7}$ approximates the Madau-Dickinson star-formation rate (Madau and Dickinson, 2014) at low redshift and the factor of $\frac{1}{1+z}$ accounts for difference in clocks between the source-frame and the detector-frame. In reality, for the redshifts considered here, $z\lesssim 0.2$, any reasonable redshift distribution, including one that traces the star-formation rate, is a very small deviation from the uniform-in-comoving volume and source-frame time redshift distribution, and is unlikely to significantly affect the results. Furthermore, the true redshift distribution will be accurately measured given a precise redshift measurement of each identified host galaxy, and can be used to update the prior $p(z)$.