Modified gravitational-wave propagation and standard sirens

# Modified gravitational-wave propagation and standard sirens

Enis Belgacem Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ansermet, CH–1211 Genève 4, Switzerland    Yves Dirian Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ansermet, CH–1211 Genève 4, Switzerland    Stefano Foffa Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ansermet, CH–1211 Genève 4, Switzerland    Michele Maggiore Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ansermet, CH–1211 Genève 4, Switzerland
###### Abstract

Studies of dark energy at advanced gravitational-wave (GW) interferometers normally focus on the dark energy equation of state . However, modified gravity theories that predict a non-trivial dark energy equation of state generically also predict deviations from general relativity in the propagation of GWs across cosmological distances, even in theories where the speed of gravity is equal to . We find that, in generic modified gravity models, the effect of modified GW propagation dominates over that of , making modified GW propagation a crucial observable for dark energy studies with standard sirens. We present a convenient parametrization of the effect in terms of two parameters , analogue to the parametrization of the dark energy equation of state, and we give a limit from the LIGO/Virgo measurement of with the neutron star binary GW170817. We then estimate the sensitivity of the Einstein Telescope (ET) to , combining standard sirens with other cosmological datasets and performing a Markov Chain Monte Carlo analysis. We discuss the prediction of a specific nonlocal modification of gravity, recently developed by our group, and we show that its predictions are within the reach of ET. Modified GW propagation also affects the GW transfer function, and therefore modifies the present-day energy density of stochastic backgrounds of GWs, as well as the tensor contribution to the ISW effect.

## I Introduction

In the last few years the spectacular observations of the gravitational waves (GWs) from binary black-hole coalescences by the LIGO/Virgo collaboration Abbott et al. (2016, 2016, 2017, 2017, 2017), as well as the observations of the GWs from the binary neutron star merger GW170817 Abbott et al. (2017), of the associated -ray burst Goldstein et al. (2017); Savchenko et al. (2017); Abbott et al. (2017), and the follow-up studies of the electromagnetic counterpart Abbott et al. (2017) have opened the way for gravitational-wave astrophysics and cosmology.

It has long been recognized Schutz (1986) that the detection of GWs from coalescing compact binaries allows us to obtain an absolute measurement of their luminosity distance. Therefore coalescing compact binaries are the GW analogue of standard candles, or “standard sirens”, as they are usually called. The standard expression of the luminosity distance as a function of redshift, , is

 dL(z)=1+zH0∫z0d~zE(~z), (1)

where

 E(z)=√ΩM(1+z)3+ρDE(z)/ρ0, (2)

and, as usual, , is the DE density and is the present matter density fraction (and we have neglected the contribution of radiation, which is negligible at the redshifts relevant for standard sirens). In the limit we recover the Hubble law , so from a measurement at such redshifts we can only get information on . This is the case of GW170817, which is at . Indeed, from the observation of GW170817 has been extracted a value  Abbott et al. (2017), that rises to if one includes in the analysis a modeling of the broadband X-ray to radio emission to constrain the inclination of the source, as well as a different estimate of the peculiar velocity of the host galaxy Guidorzi et al. (2017). The cosmological significance of this measurement can be traced to the discrepancy between the local measurement Riess et al. (2016, 2018) and the value obtained from the Planck CMB data Ade et al. (2016), that are in tension at the level. While the local measurement is a direct measurement of , independent of the cosmological model, the CMB data can be translated into a measurement of only by assuming a cosmological model, and performing Bayesian parameter estimation for the parameters of the given model. The discrepancy occurs if one assumes a standard CDM model. Thus, this tension could be a signal of deviations from CDM. The current accuracy on from the measurement with the single standard siren GW170817 is not accurate enough to discriminate between the local measurements and the Planck value. However, each standard siren provides an independent measurement of , so with standard sirens with comparable signal-to-noise ratio the error scales approximately as . The analysis of Chen et al. (2017); Feeney et al. (2018) indicates that with about 50-100 standard sirens one could discriminate between the local measurement and the Planck/CDM value.

The next generation of GW interferometers, such as the space interferometer LISA Audley et al. (2017), which is expected to fly by 2034, as well as third-generation ground-based interferometers currently under study, such as the Einstein Telescope (ET) Sathyaprakash et al. (2012) in Europe and Cosmic Explorer in the US, will have the ability to detect standard sirens at much higher redshifts. The information that one could get is then potentially much richer, since the result is now in principle sensitive to the dark energy (DE) density or, equivalently, to the DE equation of state (EoS) . Several studies have been performed to investigate the accuracy that one could obtain in this way on the DE EoS  Dalal et al. (2006); MacLeod and Hogan (2008); Nissanke et al. (2010); Cutler and Holz (2009); Sathyaprakash et al. (2010); Zhao et al. (2011); Del Pozzo (2012); Nishizawa et al. (2012); Taylor and Gair (2012); Camera and Nishizawa (2013); Tamanini et al. (2016); Caprini and Tamanini (2016); Cai and Yang (2017).

In CDM is a constant, while in a generic modified gravity model it will be a non-trivial function of . The evolution of the DE density is determined by its EoS function through the conservation equation

 ˙ρDE+3H(1+wDE)ρDE=0, (3)

which implies

 ρDE(z)/ρ0=ΩDEexp{3∫z0d~z1+~z[1+wDE(~z)]}, (4)

where . Studies of standard sirens usually assume a simple phenomenological parametrization of , given just by a constant , resulting in the CDM model, or use the parametrization Chevallier and Polarski (2001); Linder (2003)

 wDE(a)=w0+(1−a)wa, (5)

where is the scale factor, and then provide forecasts on , or on  Dalal et al. (2006); MacLeod and Hogan (2008); Nissanke et al. (2010); Cutler and Holz (2009); Sathyaprakash et al. (2010); Zhao et al. (2011); Del Pozzo (2012); Nishizawa et al. (2012); Taylor and Gair (2012); Camera and Nishizawa (2013); Tamanini et al. (2016); Caprini and Tamanini (2016), or else try to reconstruct the whole function  Cai and Yang (2017).

In this paper, elaborating on results presented in Belgacem et al. (2018) (see also Belgacem et al. (2018)) we perform a more complete analysis of the predictions of generic modified gravity models for standard sirens. In general, if the dark energy sector of a theory differs from a simple cosmological constant, this affects both the background evolution and the cosmological perturbations. The change in the background evolution is expressed by a non-trivial DE EoS . Cosmological perturbations will also be affected, both in the scalar and in the tensor sector (vector perturbations usually only have decaying modes and are irrelevant, in GR as well as in typical modified gravity models). In particular, modifications in the tensor sector can be very important for standard sirens and, as we will see, their effect on the luminosity distance can be more easily observable than that due to a non-trivial DE EoS.

The paper is organized as follows. In Section II we discuss the structure of cosmological perturbations in modified gravity theories, studying in particular the tensor sector. We will use as an example an explicit modified gravity model, the RR model, proposed and developed in the last few years by our group (based on the addition of nonlocal terms in the quantum effective action), which is very predictive and fits very well the current cosmological observations, but the arguments that we will present are more general. This explicit example will also allow us to propose a simple parametrization of the effect of modified propagation on the luminosity distance, in terms of a single parameter [or at most a pair of parameters ], that complements the pair that parametrizes the modification of the background evolution. We will see that, among these four parameters, can be the most important for observational purposes with standard sirens, so in a theory with modified GW propagation a minimal truncation of this parameter space should be to the pair . In Section III we show that standard sirens at low redshift are sensitive to the combination , and we find that the standard siren measurement of of LIGO/Virgo from GW170817 already gives a limit on this quantity. In Section IV, using the Markov Chain Monte Carlo (MCMC) method, we study the accuracy with which we can measure , and , in different combinations, using the estimated sensitivity of ET and combining it with Planck  CMB data, SNe and BAO to reduce the degeneracies between these parameters and and . In Section V we turn to a concrete model, rather than just a phenomenological parametrization, studying the perspectives for discriminating the a nonlocal modification of gravity from CDM, as a function of the number of standard sirens observed. In Section VI we will study further observable effects related to modified GW propagation, due to the modification of the transfer function that connects a primordial GW spectrum to that observed at later epochs. Section VII contains our conclusions.

## Ii Cosmological perturbations and GW propagation in modified gravity

Models such as CDM, or the parametrization (5), are not fundamental theories of dark energy but just phenomenological parametrizations, assumed to catch the effects coming from some fundamental theory that modifies general relativity (GR) at cosmological scales. Once a fundamental theory is specified, however, its consequences are much richer. First of all, the theory will modify the cosmological evolution at the background level. This will usually be equivalent to introducing an extra form of energy density , or equivalently a DE EoS , and this might (or might not) be caught by a simple expression such as the parametrization. On top of this, a specific modified gravity model will also generate cosmological perturbations that differ from those in CDM.

In GR, scalar perturbations are expressed in terms of the two gauge-invariant Bardeen potentials and , or, equivalently, in terms of their Fourier modes and . In a modified gravity theory the equations obeyed by these potential are modified, and, much as one introduces at the background level, one can define some functions that describe the deviation of the evolution of scalar perturbations from that in GR. However, at the level of perturbations, these functions will depend both on time (or, equivalently, on redshift ) and on , where is the wavenumber of the mode. For instance, a commonly used parametrization is Daniel et al. (2010); Amendola et al. (2008)

 Ψk(z) = [1+μ(z;k)]ΨGRk(z), (6) Ψk(z)−Φk(z) = [1+Σ(z;k)][Ψk(z)−Φk(z)]GR, (7)

where and are, a priori, functions of and , and the superscript “GR” denotes the same quantities computed in GR, assuming a CDM model with the same value of as the modified gravity model. This parametrization is convenient because it separates the modifications to the motion of non-relativistic particles, which is described by , from the modification to light propagation, which is encoded in . Therefore affects cosmological structure formation and affects lensing. Alternatively, one can use an effective Newton constant , defined so that, for modes well inside the horizon, the Poisson equation for becomes formally the same as in GR, with replaced by , together with a second indicator such as that, in GR, in the absence of anisotropic stresses, vanishes, but can be non-vanishing in modified gravity theories  Amendola et al. (2008); Kunz (2012). For modes well inside the horizon, for typical modified gravity models the functions , or are actually independent of . Indeed, in the absence of another explicit length scale in the cosmological model, for dimensional reasons all dependence on in the recent cosmological epoch will be through the ratio (where ), which is extremely small. Thus, any dependence of and on can be expanded in powers of [or, in fact, more typically in powers of ], and for modes well inside the horizon we can stop to the zeroth-order term.

When studying standard sirens we are rather interested in the modification of the perturbation equations in the tensor sector, i.e. in the modification of the propagation equation of GWs over the cosmological background. Let us recall that, in GR, tensor perturbations over a Friedmann-Robertson-Walker (FRW) background satisfy (we use units )

 ~h′′A+2H~h′A+k2~hA=16πGa2~σA, (8)

where are the Fourier modes of the GW amplitude, labels the two polarizations, the prime denotes the derivative with respect to cosmic time , is the scale factor, , and the source term is related to the helicity-2 part of the anisotropic stress tensor (see e.g. Maggiore (2018)). In a generic modified gravity model each term in this equation could a priori be modified by a function of redshift and wavenumber. A modification to the source term would induce a change in the production mechanism and therefore in the phase of an inspiralling binary. To understand the effect of changes in the terms or let us consider first the free propagation in GR (we follow the discussion in Belgacem et al. (2018, 2018)). We then set and we introduce a field from

 ~hA(η,k)=1a(η)~χA(η,k). (9)

Then eq. (8) becomes

 ~χ′′A+(k2−a′′a)~χA=0. (10)

For modes well inside the horizon, such as the GWs targeted by ground-based and space-born detectors, the term is totally negligible with respect to ,111More precisely, for GWs from astrophysical sources with frequencies in the range of ground-based interferometers, corresponds to an effective change of the propagation speed [with a weak time dependence: since, in MD, , changes by a factor in the propagation from the source at redshift to us]. Even if this gives rise to an integrated effect over the propagation time, still this is totally negligible compared to the bound from GW170817, which of course also comes from an effect integrated over the propagation time. For wavelengths comparable to the horizon size, for which the term is not so small, one can use a WKB approximation, as in Nishizawa (2017). and we get a standard wave equation that tells us that GWs propagate at the speed of light.

In contrast, the factor in eq. (9), that was inserted to get rid of the “friction” term proportional to in eq. (8), tells us how the GW amplitude decreases in the propagation across cosmological distances, from the source to the observer. In particular, for inspiraling binaries this factor combines with other factors coming from the transformation of masses and frequency from the source frame to the detector frame, to produce the usual dependence of the GW amplitude from the luminosity distance,

 ~hA(η,k)∝1dL(z), (11)

see e.g. Section 4.1.4 of Maggiore (2007). From this discussion we see that tampering with the coefficient of the term in eq. (8) is very dangerous, since this would modify the speed of propagation of GWs. This is by now excluded, at the level , by the observation of GW170817/GRB 170817A Abbott et al. (2017), and indeed this has ruled out a large class of scalar-tensor and vector-tensor modifications of GR Creminelli and Vernizzi (2017); Sakstein and Jain (2017); Ezquiaga and Zumalacárregui (2017); Baker et al. (2017). We next study the effect of modifying the coefficient of the friction term . We then consider the propagation equation

 ~h′′A+2H[1−δ(η)]~h′A+k2~hA=0, (12)

where is a function that parametrizes the deviation from GR, and that we have taken to be independent of the wavenumber. In this case, to eliminate the friction term, we must introduce from

 ~hA(η,k)=1~a(η)~χA(η,k), (13)

where

 ~a′~a=H[1−δ(η)]. (14)

Then we get

 ~χ′′A+(k2−~a′′~a)~χA=0. (15)

Once again, inside the horizon the term is totally negligible, so GWs propagate at the speed of light. However, now the amplitude of is proportional to rather than . As a result, rather than being just proportional to , the GW amplitude observed today, after the propagation from the source to the observer, will have decreased by a factor

 ~aem~aobs≡~a(z)~a(0) (16)

instead of a factor , where the labels refer to the emission time (at redshift ) and the observation time, at redshift zero, respectively. Therefore

 ~hA∝~a(z)~a(0)a(0)a(z)1dL(z)=~a(z)a(z)1dL(z), (17)

where is the usual notion of luminosity distance appropriate for electromagnetic signals and, since only the ratios and enter, without loss of generality we can choose the normalizations . Then, we see that in such a modified gravity model we must in general distinguish between the usual luminosity distance appropriate for electromagnetic signal, , which is given by eqs. (1) and (2), and a GW luminosity distance , with

 dgwL(z)=a(z)~a(z)demL(z). (18)

Standard sirens measure rather than . Equation (14) can be rewritten as

 (loga/~a)′=δ(η)H(η), (19)

whose integration gives Belgacem et al. (2018, 2018)

 dgwL(z)=demL(z)exp{−∫z0dz′1+z′δ(z′)}. (20)

To sum up, in modified gravity all terms in eq. (8) can in principle be different from GR. A modification of the source term affects the phase of the binary waveforms; the recent BH–BH observations, in particular of GW150914 and GW151226, have set some limit on such modifications, although for the moment not very stringent Abbott et al. (2016). A modification of the term changes the speed of gravity, and is now basically excluded. A modification of the term changes the amplitude of the GW signal received from a source at cosmological distance. This is particularly interesting because it implies that the luminosity distance measured with standard sirens is in principle different from that measured with standard candles or other electromagnetic probes such as CMB or BAO, and this could provide a “smoking gun” signature of modified gravity.

One might wonder how likely it is that a modified gravity theory predicts a non-vanishing , while still leaving untouched the term . Actually, as observed in Belgacem et al. (2018, 2018), this is just what happens in an explicit model in which gravity is modified by the addition of a nonlocal term proportional to , the so-called RR model Maggiore (2014); Maggiore and Mancarella (2014). The field theoretical motivations for this model, as well as its cosmological predictions, have been reviewed in detail in Maggiore (2017); Belgacem et al. (2018), whom we refer the interested reader. We will come back to this model in Section V. For the present discussion, the relevant point is that the model features a nonlocal term (which is assumed to emerge, at the level of the quantum effective action, because of infrared quantum effects), that behaves effectively as a dark energy, with a phantom equation of state shown in Fig. 1. The model is quite predictive and (in its simpler version, with a parameter set to zero, see Section V) has the same number of parameters as CDM, with the cosmological constant replaced by the mass parameter that appears in the term ; assuming flatness, both the parameter in the RR model and the cosmological constant in CDM can be taken as derived parameters, so the free parameters in the two models are eventually the same. The cosmological perturbations of the model have been worked out Dirian et al. (2014) and turn out to be stable, and quite close to those of CDM. The perturbations have then been inserted in a modified Einstein-Boltzmann code and compared to CMB, BAO, SNe and structure formation data Dirian et al. (2015, 2016); Dirian (2017); Belgacem et al. (2018). The result is that the model fits the cosmological observations at a level statistically indistinguishable from CDM. Furthermore, parameter estimation gives a larger value of the Hubble parameter, that significantly reduces the tension with local measurements, and provides a prediction for the sum of the neutrino masses consistent with oscillation experiments.

In the RR model the equation governing the tensor perturbations turns out to be Dirian et al. (2016); Belgacem et al. (2018, 2018)

 ~h′′A+2H[1−δ(η)]~h′A+k2~hA=16πGa2~σA, (21)

where the source term is the same as in GR to very high accuracy.222This is due to the fact that in this model GR is recovered smoothly in the limit in which the mass parameter (contrary, e.g., to massive gravity where there is the vDVZ discontinuity). The parameter is determined to be of order by the comparison with the data, and at a distance scale the corrections to GR are of order Kehagias and Maggiore (2014); Maggiore and Mancarella (2014). For of the order of the size of an astrophysical binary this is utterly negligible. We see that the speed of gravity is unchanged. The function is predicted explicitly by the model (see Belgacem et al. (2018, 2018)), and is shown in Fig. 2 as a function of redshift. At large redshifts it goes to zero because, in this model, the modifications from GR only appear close to the recent cosmological epoch. The corresponding ratio of the gravitational and electromagnetic luminosity distances is shown in Fig. 3. The fact that goes to zero at large implies that saturates to a constant value, while because there can be no effect from modified propagation when the redshift of the source goes to zero.

In Fig. 4 we show the same ratio as a function of the scale factor , and we compare it with the simple fitting function

 dgwL(a)demL(a)=Ξ0+an(1−Ξ0), (22)

with and . Observe that the parametrization (22) is such that, for small , i.e. at large redshift, goes to the constant value , while, at , . Note that this simple parametrization reproduces the exact function extremely well [indeed, in this model it performs much better than the parametrization for the equation of state, compare with Fig. 1]. In terms of redshift,

 dgwL(z)demL(z)=Ξ0+1−Ξ0(1+z)n. (23)

One could equivalently find a parametrization for the function . From eq. (20),

 δ(z)=−(1+z)ddzlog(dgwL(z)demL(z)). (24)

Then, in terms of , the parametrization (22) reads

 δ(z)=n(1−Ξ0)1−Ξ0+Ξ0(1+z)n. (25)

Note that

 δ(z=0)=n(1−Ξ0). (26)

We will however use eq. (22) or eq. (23), which are somewhat simpler, and also have the advantage of referring directly to the observable quantity .

Modified propagation equation of the form (12) have also been previously found in other modified gravity models. To the best of our knowledge this was first observed in Deffayet and Menou (2007) within the DGP model Dvali et al. (2000) (which, in its more interesting self-accelerated branch, is now ruled out by instabilities at the level of cosmological perturbations Luty et al. (2003); Nicolis and Rattazzi (2004); Gorbunov et al. (2006); Charmousis et al. (2006)). In this case the effect is due to the fact at, at large scales, gravity leaks into extra dimensions, and this affects the behavior of the amplitude of a gravitational signal. A more recent discussion of the modification to the GW luminosity distance induced by the leaking of gravity into extra dimension is given in Pardo et al. (2018), where this effect is used to put constraints on the number of extra dimensions or on the associated screening scale.333The fact that, in brane models, a gravitational signal can travel along geodesics in the extra dimensions, while electromagnetic signals are confined to the (3+1)-dimensional brane, can also lead to delays between the arrival time of a GW and the associated electromagnetic signal Chung and Freese (2000); Caldwell and Langlois (2001). Bounds on this effect from GW170817 are discussed in Visinelli et al. (2018). Modification of the propagation equation has also been found in Einstein-Aether models and in scalar-tensor theories of the Horndeski class in Saltas et al. (2014); Lombriser and Taylor (2016); Arai and Nishizawa (2017); Amendola et al. (2017); Linder (2018). This indicates that a modified propagation equation for the tensor modes of the form (12) is quite generic in alternatives to CDM; see also Gleyzes et al. (2014) for a discussion within the effective field theory approach to dark energy, and Nishizawa (2017) for a general formalism for testing gravity with GW propagation.

It is also interesting to observe that, in some theories, the function that parametrizes the deviations from CDM in the tensor sector can be related to the effective Newton constant that parametrizes the modification of the growth of structure. Indeed, it has been observed in Linder (2018) that, in a subclass of Horndeski theories called “no slip gravity”, eq. (20) can be rewritten as

 dgwL(z)=demL(z)√Geff(z)Geff(0). (27)

Quite remarkably, this relation holds also in the RR model [for modes well inside the horizon, where the dependence on wavenumber of disappears], as found in Belgacem et al. (2018). In turn, eq. (23) provides a simple fitting formula for .

Observe also that, in eq. (12), we could have inserted a function . The RR model actually predicts a function independent of wavenumber, for modes well inside the horizon. More generally, the GWs relevant for ground-based or space-borne GW interferometers have a wavelength so small, compared to the horizon size, that no significant dependence on momentum should be expected in typical modified gravity models, for the same reason discussed above for the functions and that parametrize the modifications in the scalar sector.

Observe also that, when at all redshifts, as in the RR model, (or more generally, when ) we have . Since the GW amplitude is proportional to , this means that a GW source is actually magnified, with respect to the GR prediction, and can therefore be seen to larger distances.

## Iii Limits on modified GW propagation from GW170817

We first compare the above results with the limits on modified GW propagation that can already be obtained by the single standard siren provided by the coincident detection of the GWs from the neutron star binary GW170817 Abbott et al. (2017) with the -ray burst GRB170817A Goldstein et al. (2017); Savchenko et al. (2017); Abbott et al. (2017). From the identification of the electromagnetic counterpart Abbott et al. (2017), the redshift of the source is , so in this case we are at very small redshift. After correcting for the peculiar velocity of the host galaxy, one finds a cosmological redshift  Hjorth et al. (2017). To first order in , eq. (20) becomes

 dgwL(z)demL(z)=1−zδ(0)+O(z2), (28)

so in this limit we are actually sensitive to . In the parametrization , but in fact the analysis for sources at such a low redshift can be carried out independently of any parametrization of the function . Note, however, that the deviation of from 1 is also proportional to . It is therefore clear that, with a single standard siren at a redshift , we cannot get a stringent limit on modified GW propagation. In any case, it is methodologically interesting to carry out this exercise more quantitatively. To this purpose, we assume that the correct value of is the one obtained from local electromagnetic measurements Riess et al. (2018)

 H0=73.48±1.66, (29)

(here and in the rest of the paper is given in units of ), that updates the value found in Riess et al. (2016).444We use the local measurement of , rather than the Planck value, since we want to compare standard sirens to standard candles at comparable redshifts. In this logic, the discrepancy between the local measurement of and the Planck value would be due to the fact that local (electromagnetic) measurements are independent of the cosmological model, while to translate the CMB observations into a value of we need a cosmological model. A discrepancy takes place if one assumes CDM, while in modified gravity model, in particular with a phantom De EoS, this tension is reduced or eliminated. In particular, in the (minimal) RR nonlocal model  Belgacem et al. (2018), so the tension with the updated value (29) is reduced to the level. The value of obtained from GW170817/GRB170817A in Abbott et al. (2017), assuming no modification in the GW propagation, is

 Hgw0=70.0+12.0−8.0, (30)

where we have added the superscript “gw” to stress that the measurement is obtained with standard sirens. As we mentioned in the Introduction, this value rises to

 Hgw0=75.5+11.6−9.6, (31)

if one includes in the analysis a modeling of the broadband X-ray to radio emission to constrain the inclination of the source, as well as a different estimate of the peculiar velocity of the host galaxy Guidorzi et al. (2017). Equation (29) is obtained from electromagnetic probes, and then using . In contrast, eqs. (30) and (31) are obtained from the measurement of the GW luminosity distance of this source, and then evaluating the quantity . This is the same as only if there is no modified GW propagation, so that . If we rather take into account the possibility of modified GW propagation, the correct value of obtained from a standard siren at low redshift is rather

 H0 ≡ zdemL(z)+O(z2) (32) = [1−zδ(0)]zdgwL(z)+O(z2) = [1−zδ(0)]Hgw0+O(z2),

and therefore, for a source for which we can neglect the terms, such as GW170817,

 δ(0)=Hgw0−H0Hgw0z. (33)

Using in eq. (33) the value of given in eq. (29), the value of from eq. (30) and the cosmological redshift of the source , we get555According to eq. (33), the separate errors , and induce an error on given, respectively, by , and . We have computed the error on by adding the errors in quadrature, both for the upper limit and for the lower limit. In principle there could be a correlation between the and , since the uncertainty in the cosmological redshift has already been used in Abbott et al. (2017) when determining . However, numerically is actually negligible with respect to and , so the issue here is irrelevant.

 δ(0)=−5.1+18.5−12.5, (34)

while using the value in eq. (31) we get

 δ(0)=2.7+15.4−12.8. (35)

By comparison, we have seen above that the RR model predicts and , so (actually, this is the prediction of the “minimal” version of the RR model, where a parameter is set to zero, see Section V). As expected, the limits (34) or (35), that are obtained from a single standard siren at very low redshift, are not very stringent. However, with standard sirens with comparable accuracy the error scales as , so with standard sirens at advanced LIGO/Virgo, at a redshift comparable to that of GW170817, the error on would become . Furthermore, sources at higher redshift (that should be available at advanced LIGO/Virgo at design sensitivity) could significantly improve these limits since, as we saw in eq. (28), in the low- regime the deviation of from 1 is proportional to .

Note that, at the level of this analysis, there is no degeneracy between and the DE EoS . Indeed, the value (29) of is obtained from standard candles at sufficiently low redshift, so that it is basically independent of the cosmological model, and the comparison of the value (29) with (30) or (31), that are obtained from standard sirens at low redshift, only depends on the ratio of the gravitational and electromagnetic luminosity distances, and therefore only on , but not on .

## Iv Measuring w0, wa, Ξ0 with standard sirens at large redshift

We next discuss the prospects for measuring modified GW propagation and the dark energy EoS for sources at larger redshift, as appropriate for LISA and ET. While a large part of our discussion is general, the specific comparison with observations will be made using sources (binary neutron stars) and the sensitivity appropriate for ET.

### iv.1 Understanding the role of degeneracies

Before moving directly to the results of our MCMC runs, it is useful to understand in physical terms why the parameter (or, more generally, the function that describes modified GW propagation) can be more relevant than the DE EoS , for studies of DE with standard sirens. To this purpose, let us first start from a simple CDM model, with a fixed value of , and let us ask how a set of measurements of the luminosity distances with standard sirens could help in discriminating it from CDM. For CDM, is constant and eq. (4) gives

 ρDE(z)/ρ0=ΩDE(1+z)3(1+w0), (36)

where is fixed in terms of by the flatness condition, . Thus,

 dL(z;H0,ΩM,w0) = 1+zH0 ×∫z0d~z√ΩM(1+z)3+(1−ΩM)(1+z)3(1+w0),

where we have written explicitly the dependence of on the cosmological parameters. Let us first consider

 ΔdLdL≡dL(z;H0,ΩM,w0)−dΛCDML(z;H0,ΩM)dΛCDML(z;H0,ΩM) (38)

This is the relative difference between the luminosity distance in CDM with a given value of , and the luminosity distance in CDM (where ), at fixed . For this quantity is shown as the green, dot-dashed curve in Fig. 5, while for it is given by the green, dot-dashed curve in Fig. 6.

However, this is not the quantity relevant to observations. The only way of obtaining reasonably accurate values of both and , currently, is to use the Planck CMB data, in combination with other datasets such as baryon acoustic oscillations (BAO) and supernovae, assume a cosmological model, and determine and through Bayesian parameter estimation in that model.666At least for , one might argue that one can use the value from local measurements, which is independent of the cosmology. In any case, no sufficiently accurate measurement of can be currently obtained without using the CMB data, and therefore without assuming a cosmological model. In other words, if we want to compute the prediction of CDM (with a given ) for the luminosity distance, we must not only use the chosen value of in eq. (IV.1), but we must also use the predictions of CDM for and , that are obtained by comparing CDM (with the chosen value of ) to CMB+BAO+SNe and performing the corresponding parameter estimation. This must be compared with the prediction of CDM, which is of course also obtained using in eq. (IV.1) the values of and obtained by fitting the same dataset to CDM. Thus the relevant quantity, when comparing the predictions of CDM with a given to the predictions of CDM, is not the relative difference given in eq. (38), but rather

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(ΔdLdL)≡dL(z;Hw00,Ωw0M,w0)−dΛCDML(z;H0,ΩM)dΛCDML(z;H0,ΩM), (39)

where we have denoted by the values obtained from parameter estimation in CDM with the given , and by the values obtained in CDM (more precisely, one should use the relative priors in the two models; in this discussion we will use the best-fit values for making the presentation simpler). One might think that this is a point of details, and that it will not change the order of magnitudes involved. However, this is not true, and the effect of parameter estimation is quite significant. This can be understood by observing that, when we perform Bayesian parameter estimation, we are basically comparing the predictions of a model to a set of fixed distance indicators, such as the scales given by the peaks of the CMB, or the BAO scale, or the observed luminosity distances of type Ia SNe. Thus, the best-fit values of the cosmological parameters in a modified gravity model change, compared to their values in CDM, just in the direction necessary to compensate for the change in luminosity distance (or in comoving distance, or in angular diameter distance) induced by the non-trivial DE EoS, so that in the end the luminosity distance at large redshift retains basically the same value.

In order to check and quantify this statement, we have run a series of Markov Chain Monte Carlo (MCMC), fitting both CDM and CDM (with and with ) to the same dataset of cosmological observations. In particular, we use the CMB+BAO+SNe dataset described in detail in section 3.3.1 of Belgacem et al. (2018), that includes Planck CMB data for temperature and polarization, a compilation of isotropic and anisotropic BAO measurements, and the JLA supernovae dataset. For such datasets, Bayesian parameter estimation for CDM gives the best-fit parameters

 H0=67.64,ΩM=0.3087. (40)

In contrast, for CDM with we get

 H0=70.096,ΩM=0.2908, (41)

while, for CDM with , we find

 H0=65.658,ΩM=0.32406. (42)

The magenta dashed curves in Figs. 5 and 6 show the relative difference in luminosity distance (39), obtained using for each model its own best-fit values of and . We see two important effects.

1. At redshifts the relative difference of luminosity distances becomes much smaller (in absolute value) than that obtained by keeping and fixed (and given by the green dot-dashed curves), and this suppression is of about one order of magnitude. For instance, for , keeping fixed and , the relative difference of luminosity distances at is , while, once parameter estimation in the respective models is taken into account, this becomes , with a drop in absolute value of about a factor of 10.

2. As , the green curves in Figs. 5 and 6 go to zero. This is of course a consequence of the fact that, for , eqs. (1) and (2) reduce to , and to compute the green curves we have used the same value of in the two models. In contrast, the magenta curves do not go to zero, since for each model we are using its own value of . Observe that the fact that the relative difference (39) does not go to zero at is precisely the reason that allows the LIGO/Virgo measurement of to have potentially interesting cosmological consequences. Bayesian parameter estimation to the CMB data in different cosmological models predict different values of , and therefore a local measurement of , whether with standard candles or with standard sirens, can discriminate among different cosmological models.

If GW propagation is the same as in GR, the “signal” that we must detect with standard sirens, to distinguish the modified gravity model from CDM, is then given by the magenta lines in Fig. 5 and 6 and is much smaller, in absolute value, than the signal that would be obtained if and were externally fixed quantities, determined independently of the cosmological model. Let us assume now that, in the modified gravity theory, on top of a modified DE EoS, there is also a modified GW propagation. Then, for standard sirens, the relevant quantity is the GW luminosity distance given by eq. (20). In particular, for models (such as the RR model) where the parametrization (22) is valid, at the redshifts relevant for LISA and ET, in the modified gravity model basically differs from by the factor . In contrast, in CDM the luminosity distance for standard sirens is the same as the standard electromagnetic luminosity distance. Thus, the relevant quantity for discriminating a modified gravity model from CDM is now

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(ΔdLdL)gw≡dm,gwL(z;Hm0,ΩmM)−dΛCDML(z;H0,ΩM)dΛCDML(z;H0,ΩM). (43)

where the superscript “” denotes the quantities relative to the modified gravity model. Writing , we get

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(ΔdLdL)gw ≃(Ξ0−1)dm,emL(z;Hm0,ΩmM)dΛCDML(z;H0,ΩM) (44) xx+dm,emL(z;Hm0,ΩmM)−dΛCDML(z;H0,ΩM)dΛCDML(z;H0,ΩM) ≃(Ξ0−1)+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(ΔdLdL).

The last term is the relative difference in the electromagnetic luminosity distances introduced in eq. (39). We have seen in the example of CDM with or that, even if differ from by , eventually, because of the compensating effect of and , this term in absolute value is only of order . Thus, if differs from one by more than this, it will dominate the signal.

This is indeed what happens in the RR model, where differs from by about , similarly to the CDM model with , and , as we see from Fig. 3. The situation for the relative difference of luminosity distances in the RR model is illustrated in Fig. 7, where the green and magenta curves are obtained as in Fig. 5 and 6, while the blue solid curve is the relative difference of luminosity distances (43) where, for the RR model, we use the GW luminosity distance (while, of course, for CDM the GW luminosity distance is the same as the electromagnetic one). We see that the “signal” that allows us to distinguish the RR model from CDM, represented by the blue curve, is much larger, in absolute value, than that obtained neglecting modified GW propagation, represented by the dashed magenta curve.

These results indicate that, in a generic modified gravity theory where GW propagation differs from GR, the measurement of luminosity distances of standard sirens should in general be more sensitive to modified propagation than to the DE EoS. It also implies that the prospects for detecting deviations from CDM are better than previously expected.777Of course, in a given specific modified gravity model, the function could simply be zero, or anyhow such that , in which case the main effect would come from . What our argument shows is that, in a generic modified gravity model where the deviation of from zero and the deviation of from are of the same order, the effect of dominates; see also the discussion in Belgacem et al. (2018).

### iv.2 Standard sirens and modified gravity with the Einstein Telescope

We now wish to determine more quantitatively the prospects for studying dark energy and modified gravity with future GW experiments, using the MCMC method for standard sirens combined with CMB+BAO+SNe data. For definiteness we will focus on ET. At its projected sensitivity, ET could have access to binary neutrons star (BNS) mergers up to redshifts , corresponding to events per year Sathyaprakash et al. (2010). However, only a fraction of the GW events will have an observed associated -ray burst. Estimates of the probability of observing the -ray burst are quite uncertain, depending on the opening angle of the jet (typically estimated between and ) and of the efficiency of the network of existing and future -ray telescopes Regimbau et al. (2015). A typical working hypothesis is that ET might observe to BNS with electromagnetic counterpart over a three-year period Sathyaprakash et al. (2010); Zhao et al. (2011); Regimbau et al. (2015).888Information of the redshift could also be obtained statistically, by exploiting the narrowness of the neutron star mass function (see Taylor and Gair (2012) and references therein) or by using tidal effects in neutron stars Messenger and Read (2012). In this paper we restrict to sources with detected counterpart.

We then proceed as follows. We generate a catalog of BNS detections for ET, with source, all taken to have an electromagnetic counterpart. We choose a fiducial model, that we take to be CDM with and , and we generate our simulated catalog of events by assuming that, for a source at redshift , the actual luminosity distance will be . The measured value of the luminosity distance is then randomly extracted from a Gaussian distribution centered on , and with a width obtained from an estimate of the error on the luminosity distance at ET. For this error, we add in quadrature the instrumental error and the error due to lensing,

 ΔdL(z)dL(z)=⎡⎣(ΔdL(z)dL(z))2ET+(ΔdL(z)dL(z))2lensing⎤⎦1/2. (45)

For the instrumental error we assume the expression given in Zhao et al. (2011),

 (ΔdL(z)dL(z))ET≃0.1449z−0.0118z2+0.0012z3, (46)

while for the error due to lensing we use the estimate Sathyaprakash et al. (2010); Zhao et al. (2011)

 (ΔdL(z)dL(z))lensing≃0.05z. (47)

Observe that the error due to lensing is subdominant compared to the estimate of the instrumental error in eq. (46).

To generate our catalog of events we use the standard expression of the number density of the observed events between redshift and , which is given by , where

 f(z)=4πNr(z)d2L(z)H(z)(1+z)3, (48)

(see e.g. Zhao et al. (2011)) and is the coalescence rate at redshift .999Equation (48) is derived by observing that the comoving volume between comoving distances and is . One then uses , and . Thus the observed event distribution is proportional to where is the number of events per unit time in the observer frame. Time in the observer frame, , is related to time in the source frame, , by , which provides the extra factor of at the denominator. Thus is the number of event per unit time, with respect to the unit of time relevant at redshift . The normalization constant is determined by requiring that the total number of sources be given by

 Ns=∫zmaxzmindzf(z). (49)

We take , as in Zhao et al. (2011), to have a typical signal-to-noise ratio above 8, and we also use a lower cutoff to exclude sources for which a modelisation of the local Hubble flow is necessary, before including them in the analysis. This cutoff can be estimated by observing that typical uncertainties on the recessional velocities of galaxies are around  km/s Chen et al. (2017). In the small- regime, where and , the error on the Hubble parameter is given by

 (ΔH0H0)2≃(Δvv)2+(ΔdLdL)2. (50)

Setting  km/s, , and using eq. (46) (that, in the low- regime, can be approximated by the first term), we get

 (ΔH0H0)2≃(6.67×10−4z)2+(0.15z)2. (51)

Observe that goes as , dominating at low redshifts, while is basically proportional to (because it is proportional to the inverse of the signal-to-noise ratio). Therefore there is an optimal redshift where is minimum, whose value depends of course on the sensitivity of the detector. For the sensitivity given in eq. (46), we get , corresponding to a distance  Mpc, where is defined as usual by .

For we follow Cutler and Holz (2009); Zhao et al. (2011); Cai and Yang (2017) and we use the form