Testing modified gravity at cosmological distances with LISA standard sirens

TL;DR

This paper explores how gravitational-wave observations with LISA can test modifications to General Relativity by measuring deviations in GW propagation, using a parametrized model and combining multiple cosmological datasets.

Contribution

It introduces a generic parametrization for modified GW propagation and demonstrates its effectiveness in constraining modified gravity theories with LISA and other cosmological data.

Findings

LISA standard sirens can measure the GW propagation parameter  to percent accuracy.

Combining LISA with CMB, BAO, and SNe data improves constraints on modified gravity.

Modified GW propagation is a generic feature in many theories of modified gravity.

Abstract

Modifications of General Relativity leave their imprint both on the cosmic expansion history through a non-trivial dark energy equation of state, and on the evolution of cosmological perturbations in the scalar and in the tensor sectors. In particular, the modification in the tensor sector gives rise to a notion of gravitational-wave (GW) luminosity distance, different from the standard electromagnetic luminosity distance, that can be studied with standard sirens at GW detectors such as LISA or third-generation ground based experiments. We discuss the predictions for modified GW propagation from some of the best studied theories of modified gravity, such as Horndeski or the more general degenerate higher order scalar-tensor (DHOST) theories, non-local infrared modifications of gravity, bigravity theories and the corresponding phenomenon of GW oscillation, as well as theories with extra or varying dimensions. We show that modified GW propagation is a completely generic phenomenon in modified gravity. We then use a simple parametrization of the effect in terms of two parameters $(\Xi_0,n)$, that is shown to fit well the results from a large class of models, to study the prospects of observing modified GW propagation using supermassive black hole binaries as standard sirens with LISA. We construct mock source catalogs and perform detailed Markov Chain Monte Carlo studies of the likelihood obtained from LISA standard sirens alone, as well as by combining them with CMB, BAO and SNe data to reduce the degeneracies between cosmological parameters. We find that the combination of LISA with the other cosmological datasets allows one to measure the parameter $\Xi_0$ that characterizes modified GW propagation to the percent level accuracy, sufficient to test several modified gravity theories. [Abridged]

Tables

Table 1. Table 1 : Mapping of the parametrisation in Eq. ( 2.31 ) to a number of frequently studied, representative modified gravity models embedded in the Horndeski action ( 3.1 ) with luminal speed of gravitational waves. For simplicity, we have employed the approximations α M ​ 0 ≪ 1 much-less-than subscript 𝛼 𝑀 0 1 \alpha_{M0}\ll 1 (and n ∼ 1 similar-to 𝑛 1 n\sim 1 ).

Model	${\Xi }_0-1$	$n$	Refs.
HS $f\left(R\right)$ gravity	$\frac{1}{2}{f}_{R0}$	$\frac{3\left(\tilde{n}+1\right){\Omega }_m}{4-3{\Omega }_m}$	[68] $X$
Designer $f\left(R\right)$ gravity	$-0.24{\Omega }_m^{0.76}{B}_0$	$3.1{\Omega }_m^{0.24}$	[69] $X$
Jordan–Brans–Dicke	$\frac{1}{2}\delta {\varphi }_0$	$\frac{3\left(\tilde{n}+1\right){\Omega }_m}{4-3{\Omega }_m}$	[70] $X$
Galileon cosmology	$\frac{\beta {\varphi }_0}{2M_{\mathrm{Pl}}}$	$\frac{{\dot{\varphi }}_0}{H_0\varphi }$	[71] $X$
${\alpha }_M={\alpha }_{M0}{a}^{\tilde{n}}$	$\frac{{\alpha }_{M0}}{2\tilde{n}}$	$\tilde{n}$	[67] $X$
${\alpha }_M={\alpha }_{M0}\frac{{\Omega }_{\Lambda }\left(a\right)}{{\Omega }_{\Lambda }}$	$-\frac{{\alpha }_{M0}}{6{\Omega }_{\Lambda }}\ln {\Omega }_m$	$-\frac{3{\Omega }_{\Lambda }}{\ln {\Omega }_m}$	[72, 67] $X$
$\Omega =1+{\Omega }_{+}{a}^{\tilde{n}}$	$\frac{1}{2}{\Omega }_{+}$	$\tilde{n}$	[6] $X$
Minimal self-acceleration	$\lambda \left(\ln a_{acc}+\frac{C}{2}{\chi }_{acc}\right)$	$\frac{C/H_0-2}{\ln a_{acc}^2-C{\chi }_{acc}}$	[66] $X$

Table 2. Table 2 : Forecasts for the 1 σ 𝜎 \sigma errors on Ξ 0 subscript Ξ 0 \Xi_{0} and w 0 subscript 𝑤 0 w_{0} from CMB+BAO+SNe+LISA standard sirens, for the three seed scenarios. Left: realistic scenario for the error on the luminosity distance. Right: optimistic scenario.

seeds	$\Delta {\Xi }_0$	$\Delta w_0$
hnd	$0.023$	$0.045$
hQ3	$0.036$	$0.046$
popIII	$0.044$	$0.045$

Equations

$$h_A^{\prime \prime }+2\mathcal{H}{h}_A^{\prime }+k^2{h}_A={\Pi }_A,$$

$$h_A^{\prime \prime }+2\left[1-\delta (\eta )\right]\mathcal{H}{h}_A^{\prime }+\left[c_T^2(\eta )k^2+m_T^2(\eta )\right]h_A={\Pi }_A.$$

$$\delta (\eta )=-\frac{d\ln M_{\mathrm{eff}}}{d\ln a}.$$

$$\mathcal{F}\equiv \frac{\mathcal{L}}{4\pi d_L^2}.$$

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

$$H(z)=H_0\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_R(1+z{)}^4+{\Omega }_{\mathrm{DE}}(z)}.$$

$$h_A(\eta ,\mathbf{k})=\frac{1}{a(\eta )}{\chi }_A(\eta ,\mathbf{k}),$$

$${\chi }_A^{\prime \prime }+\left(k^2-\frac{a^{\prime \prime }}{a}\right){\chi }_A=0.$$

$${\chi }_A(\eta ,\mathbf{k})\simeq {\mathcal{A}}_A\sin (k\eta +{\phi }_A),$$

$$h_{+}(t)$$

$$h_{\times }(t)$$

$$h_A^{\prime \prime }+2\mathcal{H}[1-\delta (\eta )]h_A^{\prime }+k^2{h}_A=0,$$

$$h_A(\eta ,\mathbf{k})=\frac{1}{\tilde{a}(\eta )}{\chi }_A(\eta ,\mathbf{k}),$$

$$\frac{{\tilde{a}}^{\prime }}{\tilde{a}}=\mathcal{H}[1-\delta (\eta )].$$

$${\chi }_A^{\prime \prime }+\left(k^2-\frac{{\tilde{a}}^{\prime \prime }}{\tilde{a}}\right){\chi }_A=0.$$

$$d_L^{\mathrm{gw}}(z)=\frac{a(z)}{\tilde{a}(z)}{d}_L^{\mathrm{em}}(z).$$

$$d_L^{\mathrm{gw}}(z)=d_L^{\mathrm{em}}(z)\exp \left\{-{\int }_0^z\frac{d{z}^{\prime }}{1+z^{\prime }}\delta (z^{\prime })\right\}.$$

$$\ln \frac{M_{\mathrm{eff}}(z)}{M_{\mathrm{eff}}(0)}={\int }_0^z\frac{d{z}^{\prime }}{1+z^{\prime }}\delta (z^{\prime }),$$

$$\frac{d_L^{\mathrm{gw}}(z)}{d_L^{\mathrm{em}}(z)}=\frac{M_{\mathrm{eff}}(0)}{M_{\mathrm{eff}}(z)}.$$

$$\frac{d_L^{\mathrm{gw}}(z)}{d_L^{\mathrm{em}}(z)}=\sqrt{\frac{G_{\mathrm{eff}}(z)}{G_{\mathrm{eff}}(0)}}.$$

$${\tilde{h}}_A(\eta ,\mathbf{k})=\frac{1}{a(\eta )}{\mathcal{A}}_A\sin (k\eta +{\phi }_A).$$

$${\dot{h}}_A\simeq \frac{{\mathcal{A}}_{A}k\cos (k\eta +{\phi }_A)}{a^2},$$

$${\rho }_{\mathrm{gw}}=\frac{1}{16\pi G}\sum _{A=+,\times }⟨{\dot{h}}_A^2⟩,$$

$$\frac{1}{G_{\mathrm{eff}}(z)}\frac{1}{(a\tilde{a}{)}^2}=\frac{1}{G_{\mathrm{eff}}(z=0)}\frac{1}{a^4},$$

$$S=\frac{1}{8\pi G}\int d^{4}x\sqrt{-g}A(\varphi )R+\dots ,$$

$${\chi }_A^{\prime \prime }+\left(k^2{c}_T^2(\eta )-\frac{{\tilde{a}}^{\prime \prime }}{\tilde{a}}\right){\chi }_A=0.$$

$${\chi }_A(\eta ,\mathbf{k})\simeq {\mathcal{A}}_A\sqrt{\frac{c_T({\eta }_0)}{c_T(\eta )}}\sin \left[k{\int }^{\eta }d{\eta }^{\prime }{c}_T({\eta }^{\prime })+{\phi }_A\right],$$

$$d_L^{\mathrm{gw}}(z)=\sqrt{\frac{c_T(z)}{c_T(0)}}\frac{a(z)}{\tilde{a}(z)}(1+z)r_{\mathrm{com}}(z),$$

$$r_{\mathrm{com}}(z)={\int }_0^{z}d\tilde{z}\frac{c_T(\tilde{z})}{H(\tilde{z})}.$$

$$d_L^{\mathrm{gw}}(z)=\sqrt{\frac{c_T(z)}{c_T(0)}}\exp \left\{-{\int }_0^z\frac{d{z}^{\prime }}{1+z^{\prime }}\delta (z^{\prime })\right\}(1+z){\int }_0^{z}d\tilde{z}\frac{c_T(\tilde{z})}{H(\tilde{z})},$$

$$d_L^{\mathrm{em}}(z)=(1+z){\int }_0^{z}d\tilde{z}\frac{1}{H(\tilde{z})}.$$

Full text

1 Introduction

The physics of gravitational waves (GWs) provides important opportunities to improve our understanding of gravitational interactions, and to test theories of gravity alternative to General Relativity (GR). For example, the recent simultaneous detection of a gravitational and an electromagnetic signal from GW170817 [1, 2, 3, 4], and the follow-up studies of the electromagnetic counterpart [5] impose severe constraints on theories predicting a different speed for gravitational and electromagnetic waves [6, 7, 8, 9, 10].

The next generation of GW experiments, and in particular the LISA mission [11], has the potential of performing stringent tests of other aspects of modified gravity theories, by studying the propagation of GWs across cosmological distances. In GR, the linearised evolution equation for GWs traveling on an FRW background in four dimensional space-time is

[TABLE]

where the primes indicate derivatives with respect to conformal time $\eta$, related to the physical time through the usual relation $d\eta=dt/a(t)$, with $a(t)$ the scale factor, $A=+,\times$ labels the two polarizations, and $\Pi_{A}$ is the source term, related to the anisotropic stress tensor. In a generic modified gravity model the above equation is modified into

[TABLE]

More generally, $\delta(\eta)$ could be a function also of wavenumber $k$, and the expression $(c_{T}^{2}(\eta)\,{k^{2}}+m_{T}^{2}(\eta))$ could be replaced by an expression with a more complicated $k$ dependence, corresponding to a non-trivial dispersion relation (parity-violating theories could also a priori introduce a dependence of these quantities on the polarization index $A$). Each of the quantities $\delta(\eta),c_{T}(\eta)$ and $m_{T}(\eta)$ appearing in the equation (1.2) can in principle be tested with GW observations. With respect to GR, the GW evolution equation (1.2) contains several new ingredients:

The cosmological homogeneous background controlling the Hubble parameter ${\cal H}$ is in general distinct from $\Lambda$CDM, given the different background solutions in modified gravity. This, as we shall see, must be taken into due account when investigating the dynamics of gravitational waves propagating through cosmological distances. 2. 2.

The function $\delta(\eta)$ modifies the friction term in the propagation equation. As we will recall below, this affects the amplitude of a GW propagating across cosmological distances, giving rise to a notion of “gravitational-wave luminosity distance” [12, 13, 6, 14, 15, 16].333Notice that $\delta$ is also indicated with other names in the literature, such as $\nu$ and $\alpha_{M}$, related to $\delta$ by $\nu\,=\,\alpha_{M}\,=\,-2\delta$. See for example [17]. In several modified gravity theories the function $\delta(\eta)$ is related to the time dependence of the effective Planck mass $M_{\rm eff}(\eta)$ (that could arise because of non-minimal coupling of tensor modes with other fields), via the relation

[TABLE]

In GR, $\delta=0$. Notice however that (as discussed in [16] and as we will recall below), the relation (1.3) is not universally valid in modified gravity theories, and is not obeyed in some interesting models; thus, one should not a priori identify a non-vanishing $\delta(\eta)$ with a time-dependent effective Planck mass. 3. 3.

The tensor velocity $c_{T}$ can be in general time (and scale) dependent. In GR, $c_{T}=c$, with $c$ the light speed (that we have set to one). 4. 4.

While in GR the tensor modes are massless, in theories of modified gravity the tensor mode can be massive, with $m_{T}$ its mass. This can occur in theories such as massive gravity or bigravity. 5. 5.

In the presence of anisotropic stress, or in theories where tensors couple with additional fields already at linearised level (as in theories breaking spatial diffeomorphisms), the tensor evolution equation contains a “source term” $\Pi_{A}$ in the right hand side of eq. (1.2). In absence of anisotropic stress, and in cosmological scenarios where spatial diffeomorphisms are preserved, we have $\Pi_{A}=0$.

The physical consequences of these parameters have been discussed at length in the literature (see [18] for a review on their implications for GW astronomy). In this paper we investigate how they affect a specific observable, the GW luminosity distance, which can be probed by LISA standard sirens.

The space-based interferometer LISA can qualitatively and quantitatively improve our tests on the propagation of gravitational waves in theories of modified gravity. LISA can probe signals from standard sirens of supermassive black hole mergers (MBHs) at redshifts $z\sim{\cal O}(1-10)$, much larger than the redshifts $z\sim{\cal O}(10^{-1})$ of typical sources detectable from second-generation ground-based interferometers. This implies that LISA can test the possible time dependence of the parameters controlling deviations from GR or the standard $\Lambda$CDM model, since GWs travel large cosmological distances before reaching the observer. Moreover, as we will review in section 4, LISA can measure the luminosity distance to MBHs with remarkable precision, thereby reaching an accuracy not possible for second-generation ground-based detectors.

It is also interesting to observe that LISA can probe GWs in the frequency range in the milli-Hz regime (more precisely, in the interval $10^{-4}-10^{0}$ Hz), much smaller than the typical frequency interval of ground-based detectors, $10^{1}-10^{3}$ Hz. This is a theoretically interesting range to explore since several theories of modified gravity designed to explain dark energy, such as Horndeski, degenerate higher order scalar-tensor (DHOST) theories or massive gravity, have a low UV cutoff, typically of order $\Lambda_{\text{cutoff}}\sim\left(H_{0}^{2}\,M_{\text{Pl}}\right)^{1/3}\,\sim\,10^{2}\,{\text{Hz}}$. This cutoff is within the frequency regime probed by LIGO, making a comparison between modified gravity predictions and GW observations delicate [19]. The frequency range tested by LISA, instead, is well below this cutoff, hence it lies within the range of validity of the theories under consideration.

The paper is organized as follows. In section 2 we recall the notion of modified GW propagation and GW luminosity distance, that emerges generically in modified theories of gravity. In section 3 we discuss the prediction on modified GW propagation of some of the best studied modified-gravity theories: scalar-tensor theories (with particular emphasis on Horndeski and DHOST theories), infrared non-local modifications of gravity, bigravity, and theories with extra and varying dimensions. We then turn to the study of the capability of LISA to detect modified GW propagation using the coalescence of supermassive black hole binaries as standard sirens. In section 4 we discuss the construction of mock catalogs of events with LISA, and in section 5 we present the results, obtained by running a series of Markov Chain Monte Carlo (MCMC). We present our conclusions in section 6.

2 The gravitational-wave luminosity distance

2.1 Luminosity distance and standard sirens in GR

The standard luminosity distance associated to electromagnetic signals, $d_{L}$, is defined in terms of the energy flux ${\cal F}$ measured in the observer frame, and of the intrinsic luminosity ${\cal L}$ measured at the source frame, as

[TABLE]

From the propagation of electromagnetic signals over an (unperturbed) FRW background, one finds the standard expression for $d_{L}$ as a function of redshift,

[TABLE]

where

[TABLE]

Here $\Omega_{M}=\rho_{M}(t_{0})/\rho_{0}$ is the present matter density fraction (where, as usual, $\rho_{0}=3H_{0}^{2}/(8\pi G)$ is the critical energy density and $t_{0}$ the present value of cosmic time), $\Omega_{R}=\rho_{R}(t_{0})/\rho_{0}$ is the present radiation density fraction, and $\Omega_{\rm DE}(z)=\rho_{\rm DE}(z)/\rho_{0}$, where $\rho_{\rm DE}(z)$ is the dark energy (DE) density in the cosmological model under consideration. In particular, in $\Lambda$CDM, $\Omega_{\rm DE}(z)=\Omega_{\Lambda}$ is constant. We are assuming for simplicity spatial flatness, so $\Omega_{M}+\Omega_{R}+\Omega_{\rm DE}(z=0)=1$. The luminosity distance therefore encodes important information about the cosmological model, and is a prime cosmological observable, that can be measured in particular using type Ia supernovae.

In GR, the amplitude of GWs produced by a binary astrophysical system provides yet another measurement of $d_{L}(z)$. Indeed, introducing a field $\chi_{A}(\eta,{\bf k})$ from

[TABLE]

eq. (1.1) becomes

[TABLE]

For modes well inside the horizon $a^{\prime\prime}/a$ is negligible with respect to $k^{2}$,444In principle, the effect of the $a^{\prime\prime}/a$ term could be included use a WKB approximation, as in [14]. However, the relative size of the $a^{\prime\prime}/a$ and $k^{2}$ terms is of order of the square of the wavelength $\lambda_{\rm GW}$ of the GW over the size of the horizon, $(\lambda_{\rm GW}/H_{0}^{-1})^{2}$, and this correction is therefore not significant for LISA or for ground-based detectors. and we get a solution for $\chi_{A}$ of the form

[TABLE]

with ${\cal A}_{A}$ the amplitude and $\varphi_{A}$ a phase. This shows that, overall, $h_{A}$ scales as $1/a$ in the propagation over cosmological distances. For a coalescing binary, combining this factor with the standard behavior $1/r$ in the near region, as well as with a redshift-dependent factor that arises in transforming frequencies and masses from the rest-frame to the observed one, we obtain the standard result $h_{A}\propto 1/d_{L}(z)$ (see e.g. section 4.1.4 of [20] for a detailed derivation). More precisely, in the so-called restricted post-Newtonian (PN) approximation, where one takes into account the PN corrections to the phase but not to the amplitude, one finds

[TABLE]

where $\Phi(t)$ is the phase, that in general needs to be computed to a high PN order, ${\cal M}_{c}=(1+z)(m_{1}m_{2})^{3/5}(m_{1}+m_{2})^{-1/5}$ is the redshifted chirp mass (i.e. the quantity actually observed in the detector frame), $f(t)$ is the observed GW frequency, that sweeps upward in time, and $\iota$ is the inclination angle of the normal to the orbit with respect to the line of sight.

The chirp mass is accurately determined from the time evolution of $f(t)$ [e.g., to lowest order in the PN expansion, $\dot{f}=(96/5)\pi^{8/3}(G{\cal M}_{c})^{5/3}f^{11/3}$]. Then, as first observed in [21], the amplitude of GWs from a coalescing compact binary provides an absolute measurement of its luminosity distance and in this sense coalescing compact binaries are the GW analogue of standard candles, or standard sirens. As we see from eqs. (2.7) and (2.8), the main uncertainty on the standard siren measurement of $d_{L}(z)$ comes from the partial degeneracy with $\cos\iota$. This can be broken in particular if both polarizations can be measured, or if we have informations on the inclination angle, e.g. from the observation of an electromagnetic jet. Much work as been devoted in the literature to studying how standard sirens can be used for cosmology [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]; see in particular [35, 36, 37] for recent work specifically related to LISA.

2.2 GW luminosity distance in modified gravity theories

2.2.1 GW luminosity distance from modified friction term

Let us now discuss how the situation changes in modified gravity. The expression for the luminosity distance that enters eqs. (2.7) and (2.8) now depends on the equation for propagation of GWs in the theory under consideration. In particular, as recognized in [13, 6, 14, 15, 16] (see also Sect. 19.6.3 of [38]), the function $\delta(\eta)$, that modifies the “friction term” in eq. (1.2) affects the luminosity distance extracted from the observation of GWs from a coalescing binary. Indeed, consider the equation

[TABLE]

where, for the moment, we have only retained the deviations from GR induced by the function $\delta(\eta)$ in eq. (1.2). In this case, to eliminate the friction term, we must introduce $\chi_{A}(\eta,{\bf k})$ from

[TABLE]

where

[TABLE]

Then we get

[TABLE]

Once again, inside the horizon the term $\tilde{a}^{\prime\prime}/\tilde{a}$ is totally negligible, so GWs propagate at the speed of light. However, now the amplitude of $\tilde{h}_{A}$ is proportional to $1/\tilde{a}$ rather than $1/a$. As a result, the GW amplitude measured by coalescing binaries is still given by eqs. (2.7) and (2.8), where now the standard “electromagnetic” luminosity distance $d_{L}(z)$ (that we will henceforth denote by $d_{L}^{\rm em}$) is replaced by a “GW luminosity distance” $d_{L}^{\,\rm gw}$ such that

[TABLE]

By integrating eq. (2.11) (and choosing, without loss of generality, the normalizations $\tilde{a}(t_{0})=a(t_{0})=1$, where $t_{0}$ is the present value of cosmic time), this can be nicely rewritten in terms of $\delta(z)$ as [15, 16]

[TABLE]

2.2.2 Relation to the time variation of the Planck mass

We next discuss the relation between the function $\delta(\eta)$ and the possibility of a time-varying Planck mass. In several explicit models (see section 3), it has been found that the two are related by eq. (1.3). Upon integration, using $d\ln a=-dz/(1+z)$, eq. (1.3) gives

[TABLE]

Comparing with eq. (2.14) we see that

[TABLE]

Equivalently, in terms of the effective Newton constant $G_{\rm eff}(z)=1/M^{2}_{{\rm eff}}(z)$, we have

[TABLE]

The origin of this relation has been discussed in [16], where it has been found that it expresses the conservation of the (comoving) number of gravitons during the propagation. The argument runs as follows. Recall first of all how the usual scaling of the GW energy density $\rho_{\rm gw}\propto a^{-4}$ is deduced from the cosmological evolution in GR (we follow Section 19.5.1 of [38]). From eqs. (2.4) and (2.6) we learn that, for the modes well inside the horizon that we are considering,

[TABLE]

Therefore $h^{\prime}_{A}(\eta,{\bf k})\simeq{\cal A}_{A}k\cos(k\eta+\varphi_{A})/a(\eta)$, since the extra term obtained when the time derivative acts on $1/a(\eta)$ gives a term proportional to $a^{\prime}/a^{2}$, which, for modes well inside the horizon, is negligible compared to $k/a$. Using $\dot{h}_{A}=(1/a)h^{\prime}_{A}$, where the dot is the derivative with respect to cosmic time, we get

[TABLE]

so $\dot{h}_{A}$ scales as $1/a^{2}$. The energy density of the GWs is given by the usual expression

[TABLE]

where the angular bracket denotes an average over several periods. Inserting here eq. (2.19), we find that the term $\cos^{2}(k\eta+\varphi_{A})$, averaged over several periods, gives a factor $1/2$, and it then follows that $\rho_{\rm gw}\propto a^{-4}$. This expresses the fact that the physical number density of gravitons scales as $1/a^{3}$ (i.e. the comoving number density is conserved), and that the energy of each graviton redshifts as $1/a$, leading to the overall $1/a^{4}$ scaling.

In a modified gravity theory in which the comoving number of gravitons remains conserved, the physical number density still scales as $1/a^{3}$. If, furthermore, the graviton dispersion relation remains $E=k$, its energy still scales as $1/a$, so we have again $\rho_{\rm gw}\propto a^{-4}$. However, we have seen that $h_{A}$ now scales as $1/\tilde{a}$, with $\tilde{a}$ in general different from $a$ and, repeating the argument in eqs. (2.18) and (2.19), $\dot{h}_{A}$ scales as $1/(a\tilde{a})$ (since the extra factor of $a$ comes from the transformation between cosmic time and conformal time). On the other hand, $G$ in eq. (2.20) must be replaced by $G_{\rm eff}(t)$. In order to get $\rho_{\rm gw}\propto a^{-4}$, we must then have

[TABLE]

which, upon use of eq. (2.13), is equivalent to eq. (2.17). Thus, eq. (2.17) is valid in any modified gravity theory in which the graviton number is conserved and the graviton dispersion relation is not modified, and therefore in a very broad class of models. In particular, it can be shown [16] that this relation holds in any modified gravity theory described by an action of the form

[TABLE]

which is minimally coupled to matter, where $A$ can be a nontrivial functional of extra fields in the gravitational sector, here denoted collectively as $\phi$, and the dots denote other possible gravitational interaction terms, that can depend on $\phi$ but do not contain terms purely quadratic in the gravitational field $h_{\mu\nu}$ nor interactions with ordinary matter.

It is important to observe that the time variation of the effective Newton constant given in eq. (2.17) [or, equivalently, the time variation of the Planck mass given by eq. (2.16)] can only hold on cosmological scales. At such scales, typical cosmological models predict a time dependence such that, today, $|\dot{G}/G|\,\simeq\,H_{0}$, with a coefficient in general of order one. A scenario in which this result holds down to solar system or Earth-Moon scales would be ruled out by Lunar Laser Ranging experiments, that by now impose a bound |\dot{G}/G|\,\raise 1.72218pt\hbox{<}\kern-8.00003pt\lower 2.6694pt\hbox{\sim}\,10^{-3}H_{0} [39]. If a model predicts a non-trivial value for the ratio $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ then, to be observationally viable, it must either have a screening mechanism at short scales, so that the time-dependent effective Newton constant predicted by eq. (2.17) cannot be extrapolated down to the Earth-Moon scale (this is the case, for instance, for Hordenski and DHOST theories555Very recently it has been pointed out [40] that, for the special case of the cubic Galileon, the screening is not sufficient to guarantee an observable effect in the measurement of (2.17) for the redshifts probed by LISA standard sirens. However, a more general form for the cubic Horndeski term, or the screening provided by DHOST theories, can still provide sizable effects of (2.17) within the reach of LISA as we will see in Section 3.1.2., and bigravity) or else must violate the relation (2.17) and predict $G_{\rm eff}=G$ for all modes well inside the horizon, even if it still predicts a non-trivial result for $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$. As we will see in section 3.2, the latter option is realized in the so-called RT non-local gravity model.

Note also that the relation between the effective Newton constant $G_{\rm eff}$ that governs the gravitational dynamics at cosmological scales and the Newton constant $G$ observed at laboratory or solar-system scales can be model-dependent. In the absence of screening mechanisms, the effective Newton constant at cosmological scales is the same as that at laboratory or solar-system scales, and then $G=G_{\rm eff}(z=0)$. In the presence of screening mechanisms, however, the relation can be more complex.

2.2.3 GW luminosity distance from a time-dependent speed of GWs

We now discuss how the above results for the luminosity distance change if we also include a non-trivial function $c_{T}(\eta)$ in eq. (1.2). In that case we still define $\chi$ from eq. (2.10), with $\tilde{a}$ given by eq. (2.11). Then the “friction term” in eq. (2.11) is again eliminated, and $\chi_{A}$ now satisfies the equation

[TABLE]

For modes inside the horizon (defined now by the condition $|a^{\prime\prime}/a|\ll k^{2}c^{2}_{T}$) we can neglect the term $\tilde{a}^{\prime\prime}/\tilde{a}$. Assuming that the frequency at which $c_{T}(\eta)$ changes is much smaller than the frequency $f=k/(2\pi)$ of the GW [which is an extremely good approximation for GWs detectable with LISA or with ground based interferometers, given that the typical scale of change of $c_{T}(\eta)$ is expected to be fixed by $H(\eta)$], the resulting equation has the WKB-like solution

[TABLE]

where we have normalized the amplitude ${\cal A}_{A}$ so that, for $c_{T}(\eta)$ constant and equal to the present value $c_{T}(\eta_{0})$, we recover eq. (2.6). Thus, a non-trivial speed of GWs, $c_{T}(\eta)$, affects both the phase and the amplitude of the signal from a coalescing binary. Therefore it also affects the GW luminosity distance, so the GR relation $d_{L}(z)=(1+z)a(t_{0})r_{\rm com}(z)$, where $r_{\rm com}$ is the comoving distance to the source, is now modified into

[TABLE]

where $c_{T}(0)\equiv c_{T}(z=0)$ and we used $a(t_{0})=\tilde{a}(t_{0})=1$. Furthermore, the expression of the comoving distance measured with GWs as a function of redshift is also affected by a non-trivial speed of GWs, and is now given by

[TABLE]

Thus, in the end,

[TABLE]

to be compared with

[TABLE]

In App. A we show that the same relation emerges if we define the GW luminosity distance from eq. (2.1), rather than from the $1/d_{L}^{\rm gw}$ behavior of the GW amplitude from a localized source, as we have done here. Also note that the effect can be related to Eq. (2.22) with the conformal factor for scalar-tensor theories corresponding to $M_{\rm eff}^{2}c_{T}$ [6, 41].

The observation of GW170817/GRB 170817A puts a limit $|c_{T}-c|/c<O(10^{-15})$, but this limit only holds for redshifts smaller than the redshift of that source, i.e. for z\,\raise 1.72218pt\hbox{<}\kern-8.00003pt\lower 2.6694pt\hbox{\sim}\,0.01. LISA can measure $c_{T}(z)$ to much higher redshifts. In particular, the simultaneous observation of a GW event and of its electromagnetic counterpart would put extremely strong limits on $|c_{T}(z)-c|/c$ up to the redshift of the source, just as for GW170817/GRB 170817A.

Given the strong observational constraint from GW170817/GRB 170817A, and the lack of explicit models where $c_{T}(z)$ evolves from a value equal to $c$ within 15 digits at $z<0.01$, to a sensibly different value at higher redshift, in the following we will limit our analysis to the case $c_{T}(z)=c$. Note also that, if at higher redshift $c_{T}(z)$ should be sensibly different from $c$, with LISA one would simply not see an electromagnetic counterpart even if it existed, since the time delay of the electromagnetic and gravitational signal, over such distances, would be huge. In that case the analysis of the present paper, that assumes standard sirens with electromagnetic counterpart, would not be applicable, and one would have to resort to statistical methods.666At low-$z$, an alternative way to test an anomalous GW speed at LISA frequencies, $c_{T}(k_{{}_{\rm LISA}})\neq c$, is to measure the phase lag between GW and EM signals of continuous sources such as the LISA verification binaries. This test can constrain the graviton mass [42, 43] as well as the propagation speed [44, 45].

2.2.4 Phenomenological parametrization of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$

In general, in modified gravity, both the cosmological background evolution and the cosmological perturbations are different with respect to GR. It is obviously useful to have phenomenological parametrizations of these effects, that encompass a large class of theories. In modified gravity, the deviation of the background evolution from $\Lambda$CDM is determined by the DE density $\rho_{\rm DE}(z)$ or, equivalently, by the DE equation of state $w_{\rm DE}(z)$. In principle one could try to reconstruct the whole function $w_{\rm DE}(z)$ from cosmological observations, but current results are unavoidably not very accurate (see e.g. fig. 5 of [46]). The standard approach is rather to use a parametrization for this function, that catches the qualitative features of a large class of models. The most common is the Chevallier–Polarski–Linder parametrization [47, 48], which makes use of two parameters $(w_{0},w_{a})$,

[TABLE]

corresponding to the value and the slope of the function at the present time. In terms of redshift,

[TABLE]

One can then analyze the cosmological data adding $(w_{0},w_{a})$ to the standard set of cosmological parameters. Similarly, some standard parametrizations are used for describing the modification from GR in the scalar perturbation sector, in order to compare with structure formation and weak lensing, see e.g. [49, 50]. Here we are interested in tensor perturbations, where the effect is encoded in the non-trivial function $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$. Again, rather than trying to reconstruct this whole function from the data, it is more convenient to look for a simple parametrization that catches the main features of a large class of models in terms of a small number of parameters. We shall adopt the 2-parameter parameterization proposed in Ref. [16],

[TABLE]

which depends on the parameters $\Xi_{0}$ and $n$, both taken to be positive. In terms of the scale factor $a=1/(1+z)$ corresponding to the redshift of the source,

[TABLE]

The value $\Xi_{0}\,=\,1$ corresponds to GR. This parameterization is designed to smoothly interpolate between a unit value

[TABLE]

at small redshifts – where cumulative effects of modified gravity wave propagation have not sufficient time to accumulate differences with respect to GR, see eq. (2.14) – to a constant value $\Xi_{0}$

[TABLE]

at large redshift. Indeed, in the large redshift regime we expect that the effects of modified gravity “turn-off” and $|\delta(z\gg 1)|\ll 1$, since modified gravity should mainly affect late-time evolution (also for ensuring compatibility with CMB observations), in which case the quantity $\Xi(z)$ approaches a constant. This parametrization was originally proposed in [16], inspired by the fact that it fits extremely well the prediction for $\Xi(z)$ obtained from a nonlocal modification of gravity (see section 3.2 and ref. [51] for review), but it was then realized that its features are very general, so that it is expected to fit the predictions from a large class of models. Indeed, in section 3 we will compare this fitting formula to the explicit predictions of several modified gravity models, and confirm that it is appropriate in many situations.777Recently, an analysis of modified GW propagation for ground-based advanced detectors, including ET, has been presented in [52]. That work uses a parametrization for modified GW propagation which corresponds to setting the function $\delta(z)$ in eq. (2.14) to a constant value. As discussed in [16], this is a special case of the parametrization (2.31), with $\Xi_{0}=0$ and $n=\delta$. As we will see in section 3, typical models actually predict variation of $\Xi_{0}$ of only a few percent from the GR value $\Xi_{0}=1$. In any case, the conclusion of [52] that modified GW propagation can be tested at the $1\%$ level at future ground-based detectors such as ET broadly agrees with the finding in [16]. Ref. [53] recently studied the bound on modified GW propagation from the single event GW170817. Again, the parametrization of modified GW propagation is different, with $\delta(z)$ taken to be proportional to $\Omega_{\rm DE}(z)$, but the results are consistent with those presented earlier in [16], where the limit from GW170817 was also computed.

For theories for which eqs. (2.16) and (2.17) hold (see the discussion in section 2.2.2), from eq. (2.31) we can obtain a corresponding parametrization for the time variation of the effective Planck mass or of the effective Newton constant,

[TABLE]

Furthermore, from eq. (2.14), we have

[TABLE]

In this parametrization the quantity $\delta(z)$ indeed goes to zero at large redshifts, as desired: at early times gravity propagates as in GR. At late times, instead, $\delta(z\ll 1)\,=\,n\,(1-\Xi_{0})$. Fig. 1 shows $\Xi(z)$, $\delta(z)$, $M_{\rm eff}(z)/M_{\rm eff}(0)$ and $G_{\rm eff}(z)/G_{\rm eff}(0)$ as function of the redshift, setting for definiteness $n=2.5$ and $\Xi_{0}=0.97$.

It should also be observed that the parametrization (2.31) of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ is more robust than the corresponding parametrization (2.37) of $\delta(z)$. Indeed, even if $\delta(z)$ should have some non-trivial features as a function of redshift, such as a peak, still these features will be smoothed out by the integral in eq. (2.14). Since anyhow $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ must go to one as $z\rightarrow 0$ and we expect that in most models it will go asymptotically to a constant at large $z$, in general the fit (2.31) to $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ will work reasonably well even in cases where the corresponding fit (2.37) is not too good. We will see an explicit example of this behavior in section 3.1.2. Since, for standard sirens, the observable quantity is $d_{L}^{\,\rm gw}(z)$, rather than the function $\delta(z)$ itself, this means that the parametrization (2.31) is quite robust for observational purposes.

Of course, in general, over a large interval of redshifts the simple parametrization (2.31) might no longer be accurate, and in particular the same index $n$ might not fit well the whole range of redshifts of interest. An extreme case are the bigravity theories that we shall discuss in section 3.3, where the effects of graviton oscillations lead to peculiar behaviors for the GW luminosity distance as a function of redshift, see e.g. Fig. 9. Similarly to what is usually done for the parametrization of the DE equation of state, one could further improve the accuracy of the parametrization (2.31) by fixing $n$ at a pivot redshift $z_{p}$ where the observational errors are the lowest,

[TABLE]

Another possibility is to specify the slope of the evolution of $\Xi(z)$ at $z=0$, i.e.,

[TABLE]

However, in considering LISA standard sirens at sufficiently large redshift, where $\Xi(z)$ approaches its asymptotic value, the precise value of $n$ will be of limited importance, and the crucial parameter in the parametrization (2.31) will rather be $\Xi_{0}$.

An alternative parametrization where $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ differs from unit by ${\cal O}(z^{m})$ at small redshifts, rather than by ${\cal O}(z)$, can be obtained by modifying (2.31) into

[TABLE]

In the low-$z$ limit, the expression on the right-hand side approaches $1+(\Xi_{0}-1)z^{m}$, and still goes to $\Xi_{0}$ at large $z$. We will discuss in sections 3.1.2 and 3.4 examples of models where a $1+{\cal O}(z^{m})$ behavior could emerge, as well as other possible parametrizations of this type. Note however that, for a model in which at low redshift $\Xi(z)=1+{\cal O}(z^{m})$, a two-parameter parametrization might not be adequate, since it would connect the index $m$ that determines the behavior at $z\ll 1$ to the index describing the decrease at large $z$, that in general will be given by a different power. Thus, in these cases, it might be necessary to rather use a three-parameter expression such as

[TABLE]

which is obtained by replacing $z\rightarrow z^{m}$ in eq. (2.31). In the low-$z$ limit it approaches $1+n(\Xi_{0}-1)z^{m}$, while at large $z$ it approaches $\Xi_{0}$ with a different power-like behavior.

To conclude this section, let us emphasize the importance of modified GW propagation, as encoded for instance in the parameters $(\Xi_{0},n)$ defined by eq. (2.31), for studies of dark energy and modified gravity at advanced GW detectors. This stems from two important considerations:

We have seen that, in general, a modified gravity model induces deviations from $\Lambda$CDM at the background level (through the DE equation of state), in the scalar perturbation sector, and in the tensor perturbation sector. Concerning the background modifications, as encoded for instance in the $(w_{0},w_{a})$ parameters, several studies [27, 28, 35, 16] have shown that the accuracy that LISA, or third-generation ground-based interferometers such as the Einstein Telescope, could reach on $w_{0}$ is not really better than the measurement, at the level of a few percent, that we already have from Planck in combination with other cosmological probes such as Baryon Acoustic Oscillations (BAO) and Supernovae (SNe). Our analysis in section 5 will indeed confirm these findings. In contrast, $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ is an observable that is only accessible thanks to GW observations. 2. 2.

On top of this, it was realized in [15, 16] that, in a generic modified gravity theory, in which the deviation of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ from 1 is of the same order as the deviation of $w_{\rm DE}(z)$ from $-1$, the effect of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ on standard sirens dominates over the effect of $w_{\rm DE}(z)$. This can be understood from the explicit expression of the standard (“electromagnetic”) luminosity distance given in eqs. (2.2) and (2.3). From this expression one might think, naively, that if one changes the equation of state of DE by, say, $10\%$, this will induce a corresponding relative variation of $d_{L}$ again of order $10\%$. However, this is not true because the cosmological parameters, such as $H_{0}$ and $\Omega_{M}$, that enters in eq. (2.3), are not fixed, but must themselves be determined self-consistently within the model, by comparing with cosmological observations and performing Bayesian parameter estimation, and their best-fit values change if we modify the DE content of the model. The fit to cosmological data basically amounts to requiring that the model reproduces some fixed distance indicators at large redshifts, such as the scale determined by the peaks of the CMB or that from the BAO. Thus, Bayesian parameter estimation has a compensating effect, changing the luminosity distance in a direction opposite to that induced by a change in $\rho_{\rm DE}(z)$, in such a way to keep as small as possible the variation of $d_{L}^{\,\rm em}(z)$ at large redshifts. Thus, after performing Bayesian parameter estimation, a relative change in $w_{\rm DE}$ by, say, $10\%$, would only induce a relative change of order, say, $1\%$ in $d_{L}(z)$. In contrast, modified GW propagation is an extra effect, that is not compensated by degeneracies with other (fitted) cosmological parameters, and therefore dominates over the effect of $w_{\rm DE}$.

This physical argument has been confirmed with an explicit MCMC study for the Einstein Telescope in [16] where it was found that, assuming $10^{3}$ standard sirens with an electromagnetic counterpart, ET, in combination with other cosmological probes, could measure $w_{0}$ with a precision of $3.2\%$ and $\Xi_{0}$ with a precision $0.8\%$. In this paper we will see, from the explicit MCMC study in section 5, that the same pattern holds for LISA, and, using standard sirens in combination with other cosmological probes, $\Xi_{0}$ can be measured with a significantly better accuracy than $w_{0}$.

Therefore $\Xi_{0}$, or more generally modified GW propagation, is a prime observable for dark energy studies with advanced GW detectors.

3 Models of modified gravity

In this section we discuss the prediction for the ratio $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ in various scenarios. We will see that modified GW propagation is a generic phenomenon in modified gravity theories, and we will assess how well the parameterization (2.31) describes the predictions of various models. We examine scalar-tensor theories, non-local models, bigravity and theories with extra dimensions. In our discussion we shall make the simplifying hypothesis that modified gravity affects the propagation, and not the production, of GWs. This is a delicate assumption: on the other hand it is known to hold in specific, interesting examples. For the cubic Galileon model, a scalar-tensor theory of gravity equipped with a Vainshtein screening mechanism (see section 3.1), [54] shown that scalar contributions to GW radiation from merging binaries is very suppressed. Also for the non-local RT model discussed in section 3.2, any effect of modified gravity is expected to be negligible at the (cosmologically) small scales characterizing astrophysical merging events. Moreover, even for systems where our assumption is not strictly valid, we can expect that modified gravity effects occurring when GWs are produced do not influence the redshift dependence of the ratio $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$, and might be readsorbed into appropriate rescalings of (redshift-independent) parameters entering into the expression for the GW luminosity distance.

3.1 Scalar-tensor theories of gravity

The class of scalar-tensor theories of gravity includes some of the best-studied alternatives to GR aimed at explaining late-time cosmic acceleration. Besides a massless spin-2 mode, these theories contain a scalar degree of freedom whose dynamics play an important role for the evolution of the universe. Scalar-tensor theories include the simplest systems that exhibit the following, cosmologically relevant, phenomena:

self-acceleration; it can explain the present day acceleration of the universe without the need to invoke a small positive cosmological constant. See [55] for a comprehensive review. 2. 2.

screening mechanisms; these avoid fifth-force constraints associated with long-range scalar forces by suppressing the effect of the scalar sector on matter at smaller scales. Depending on which screening mechanism is in play, the corresponding suppression relies on specific features of the scalar action and/or its couplings. See [56, 57] for reviews.

In theories where the scalar acquires a non-trivial profile and is non-minimally coupled to gravity, the dynamics of the tensor modes is naturally affected: this may lead to deviations from GR. The physics of GW propagation is a particularly clean probe of modified gravity, the reason being that the tensor sector is typically the one less affected by screening mechanisms. In terms of the evolution equation (1.2), we can expect that scalar-tensor systems lead to a cosmological expansion history different from $\Lambda$CDM, and the gravitational wave equation (1.2) can have $\delta\neq 0$, and $c_{T}\,\neq\,1$. On the other hand, $m_{T}\,=0$ and $\Pi_{ij}\,=\,0$. We examine some aspects of the evolution of tensor modes over cosmological distances in these systems, first by briefly reviewing known results from the perspective of GW physics, then by delivering new predictions for recent degenerate scalar-tensor theories. In light of the GW170817 event, we focus on the case where $c_{T}\,=\,1$.

3.1.1 Horndeski theories of gravity

The most well-known scalar-tensor theories of gravity (Brans–Dicke, $f(R)$, covariant Galileon models) belong to the wide class of Horndeski theories, the most general covariant scalar-tensor theories of gravity leading to second-order equations of motion (see e.g. [58] for a recent review). They are described by the action [59, 60, 61]

[TABLE]

with the Lagrangian densities

[TABLE]

where $X\equiv\partial_{\mu}\phi\partial^{\mu}\phi$; $\phi$ is the scalar; $R$ and $G_{\mu\nu}$ denote the Ricci scalar and Einstein tensor of the Jordan frame metric $g_{\mu\nu}$; and the matter fields $\psi_{\rm m}$ in $\mathcal{L}_{\rm m}$ are minimally coupled with gravity. The $G_{i}$ are arbitrary functions of $\phi$, $X$, and with $G_{iX}$ we denote derivatives along $X$. For convenience, we have set the Planck mass $M_{\rm Pl}$ and the speed of light in vacuum to unity.

In these theories, due to the non-minimal couplings between gravity and the scalar $\phi$, tensor fluctuations typically propagate with a speed different from light, unless the functions $G_{i}$ satisfy some constraints. To ensure $c_{T}=1$, we require that [62, 63, 45]

[TABLE]

where, without loss of generality, $G_{5}$ can be set to zero as should be clear after performing an integration by parts (see also refs. [7, 8, 9, 10, 64] for discussions on fine-tuned solutions, and [65] for further aspects of GW propagation in Horndeski theories of gravity). It follows that Eq. (3.1) can no longer provide an observationally viable acceleration by means of $G_{4},\,G_{5}$ [6, 66]. The function $G_{4}(\phi)$ is associated with the effective Planck mass, $G_{4}=M_{\rm eff}^{2}/2$. As a result, there can be interesting cosmological models where $G_{4}(\phi)$ depends on time, something that can be tested with GW observations. In this set-up, $\delta$ is related to the effective Planck mass according to eq. (1.3) (see ref. [17]):

[TABLE]

as follows from the discussion in section 2.2.2. In this section we use the notation $\alpha_{M}$ [67], usually employed in analyzing scalar-tensor systems. The use of GW standard sirens in order to test the time evolution of $G_{4}$ was already proposed in Ref. [13] and a preliminary forecast at the level $|M_{\rm eff}^{2}(z=0)-1|\lesssim 3.5\times 10^{-3}$ was estimated in Ref. [6] by adapting forecasts on the accuracy that LISA can reach on $H_{0}$. We specify the mapping between a range of Horndeski models and the parametrisation in Eq. (2.31), which will enable us to interpret the constraints on $\Xi_{0}$ and $n$ for given values of the model parameters. The mapping for Horndeski scalar-tensor theories can be generally performed888 Note that one may have to specify $M^{2}(z\rightarrow\infty)$ whenever it does not reduce to $M_{\rm Pl}^{2}$ at early times. The early time matching is usually necessary for the purpose of recovering GR. In the complementary late-time regime, matching may be due to screening effects in the laboratory at $z=0$. by specifying $M^{2}(0)$ and $\alpha_{M0}$ according to

[TABLE]

This overall “dictionary” is summarised in Table 1. Note that we assume the constraint $|\Xi_{0}-1|\ll 1$ (and $n\sim 1$) and that all models recover $M_{\rm eff}(z\rightarrow\infty)=M_{\rm Pl}$, hence, $\Xi_{0}=M_{0}$, and we set $M_{\rm Pl}\equiv 1$ for convenience. In Fig. 2 we illustrate the performance of the fit provided by the parametrisation (2.31), with the values of $\Xi_{0}$ and $n$ given in eqs. (3.8) and (3.9), for two examples embedded in the Horndeski action. We see that the parametrization (2.31) works well.

The first model we shall inspect is $f(R)$ gravity [73], where the Einstein-Hilbert action is generalised by $R\rightarrow R+f(R)$. It can be mapped onto the action (3.1) by defining the scalar field $2G_{4}(\phi)\equiv\phi\equiv 1+f_{R}$ with $f_{R}\equiv{\rm d}f/{\rm d}R$ and $G_{2}\equiv-U(\phi)\equiv\frac{1}{2}[f(R)-f_{R}R]$. Hence, one finds

[TABLE]

[where $f_{R0}\equiv(f_{R})_{0}$ is the present value of $f_{R}$] as well as

[TABLE]

Adopting the particular functional form for $f$ proposed by Hu and Sawicki (HS) [68]

[TABLE]

one obtains

[TABLE]

These relations remain unchanged in a straightforward generalisation of the HS $f(R)$ model to Jordan-Brans–Dicke gravity [70] with appropriate scalar field potential $U(\phi)$ (see, e.g., Refs. [74, 75]). In that case $G_{4}\equiv\phi/2$ and $G_{2}=-U+X\omega(\phi)/\phi$, where $\omega(\phi)$ is the Brans–Dicke function, and $\omega=0$ reduces to the HS $f(R)$ case. We furthermore define $\phi\equiv 1+\delta\phi$. Alternatively to the HS $f(R)$ function, one may, for instance, also consider the designer $f(R)$ model [69] that exactly reproduces the standard cosmological expansion history. It can be parametrised by the Compton wavelength parameter $B_{0}\equiv[(H/H^{\prime})f_{R}^{\prime}/(1+f_{R})]|_{0}$. For simplicity, we make use of the approximation $B_{0}\simeq-2.1\Omega_{m}^{-0.76}f_{R0}$ [76] to derive the expressions in Table 1 (also see Fig. 2).

Next we consider a cubic Galileon model conformally coupled to the Ricci scalar. The model is defined by the following choice of the Horndeski functions

[TABLE]

where $\beta$ regulates the strength of the coupling, with the uncoupled case in the limit $\beta=0$. The coefficients $c_{2},c_{3}$ are dimensionless, and fixed by the scalar field normalization and the value of Galileon energy fraction $\Omega_{\rm gal}$ (see refs. [78, 79] for details). For consistency with expansion probes we consider models where the Galileon is a sub-dominant contribution to the energy density $\Omega_{\rm gal}=0.05$, with $\Omega_{\Lambda}=1-\Omega_{M}-\Omega_{\rm gal}$. The initial condition for the scalar field is $\phi(t_{\rm ini})=0$ to ensure that $M_{\rm eff}^{2}\approx 1$ at early times, and the initial value of the field derivative $\dot{\phi}(t_{\rm ini})$ is fixed by the tracker solution. This leads to a unique solution for $\phi_{0}$. Together with given $\Omega_{\rm gal}$, it fixes $\dot{\phi}_{0}$. It is then straightforward to express $\Xi_{0}$ and $n$ in terms of $\phi$ and its derivative today. We have $\Xi_{0}=\exp\left(\frac{\beta\phi_{0}}{2M_{\rm Pl}}\right)\simeq 1+\frac{\beta\phi_{0}}{2M_{\rm Pl}}$ from which $n=H_{0}^{-1}\frac{\dot{\phi}_{0}}{\phi_{0}}$ using Eq. (3.9). The coupling $\beta\neq 0$ introduces a non-zero value of $\alpha_{M}$ and modifies the ratio of the luminosity distances, as shown in Figure 3. This is a minimal Galileon model that modifies the luminosity distance while remaining compatible with the GW speed.

Due to the wealth of models that can be realised with action (3.1), instead of studying specific models one may adopt a phenomenological parametrisation for the time evolution of $G_{4}$. Common options include a power law for $\alpha_{\rm M}$ in terms of the scale factor $a$ or, motivated by a connection to dark energy, a proportionality to the fractional dark energy density $\Omega_{\Lambda}(a)\equiv 8\pi G\Lambda/(3H^{2})$ as well as a generalisation of that to an evolving effective dark energy contribution [72]. For the power law $\alpha_{M}=\alpha_{M0}a^{\tilde{n}}$, one finds $M(z)=\exp[\alpha_{M0}a^{\tilde{n}}/(2{\tilde{n}})]$ and $\Xi_{0}=e^{\alpha_{M0}/(2{\tilde{n}})}\simeq 1+\alpha_{M0}/(2{\tilde{n}})$ as well as $n\simeq{\tilde{n}}$ (see Fig. 2). Similarly one may also consider a power law for the evolution of the Planck mass instead. For instance, for the EFT function $\Omega=1+\Omega_{+}a^{\tilde{n}}$, corresponding to $\Omega=M^{2}$ for $c_{T}=1$, one immediately finds that $\Xi_{0}=M_{0}\simeq 1+\frac{1}{2}\Omega_{+}$ and $n\simeq\tilde{n}$. Adopting $\tilde{n}=4$, a constraint of $|\Omega_{+}|<1.4\times 10^{-2}$ was inferred from LISA forecasts in Ref. [6], where for cosmic acceleration to be attributed to a running $\Omega$, one requires $\Omega_{+}\lesssim-0.1$. For a scaling of the time evolution as $\alpha_{M}=\alpha_{M0}\Omega_{\Lambda}(a)/\Omega_{\Lambda}=\alpha_{M0}H_{0}^{2}/H^{2}$, one obtains $M^{2}=M_{0}^{2}(a^{3}H^{2}/H_{0}^{2})^{\alpha_{M0}/(3\Omega_{\Lambda})}$ with $\Xi_{0}=M_{0}=\Omega_{m}^{-\alpha_{M0}/(6\Omega_{\Lambda})}\simeq 1-\alpha_{M0}/(6\Omega_{\Lambda})\ln\Omega_{m}$ such that $M^{2}=\Omega_{m}^{-\alpha_{M0}/(3\Omega_{\Lambda})}$ and $n\simeq-3\Omega_{\Lambda}/\ln\Omega_{m}$.

Finally, to more directly address the question whether a time evolution of $G_{4}$ can be made responsible for cosmic acceleration at late times, one can formulate the minimal evolution of the Planck mass required to produce a positive accelerated expansion without the contribution of dark energy or a cosmological constant [66]. The approach defines a threshold $\alpha_{M,min}=C/(aH)-2$, where $C=2H_{0}a_{acc}\sqrt{3(1-\Omega_{m})}$ and $a_{acc}=(\Omega_{m}/(1-\Omega_{m})/2]^{1/3}$. This threshold is, however, in $3\sigma$ tension with cosmological observations, particularly of cross correlations of the CMB temperature with foreground galaxies [66]. A simple rescaling of $\alpha_{M}=\lambda\alpha_{M,min}$ assesses the fraction of cosmic acceleration $\lambda$ that can be attributed to an evolving $G_{4}$. For $\lambda\ll 1$, one finds $\Xi_{0}\simeq(a_{acc})^{\lambda}e^{\lambda C\chi_{acc}/2}\simeq 1+\lambda\left(\ln a_{acc}+C\chi_{acc}/2\right)$ and $n\simeq(C/H_{0}-2)/(\ln a_{acc}^{2}-C\chi_{acc})$, where $\chi_{acc}$ denotes the comoving distance at the acceleration scale $a_{acc}$.

This list is far from being exhaustive of the range of Horndeski models that will be tested with Standard Sirens, and does in fact provide only a small sample. Rather than expanding Table 1, however, we shall address next the constraints that can be inferred on theories that generalise action (3.1) to include higher derivatives in their equations of motion.

3.1.2 DHOST theories of gravity

The Horndeski action is built on the requirement of second order field equations. However, although sufficient, this requirement is not necessary in order to avoid dangerous Ostrogradsky instabilities [80, 81]. This aspect, noticed in [82] using transformations of the metric, led initially to the so called beyond Horndeski theories [83]. Subsequently, it was realised that higher order field equations are actually harmless provided that the Lagrangian is degenerate [84], i.e. there exists an extra primary constraint that removes the Ostrogradsky mode. Exploiting this property, Degenerate Higher Order Scalar-Tensor (DHOST) theories (see [85] for a recent review) have been constructed [86, 87, 88] and they represent (so far) the most general scalar-tensor theories that propagate a single scalar degree of freedom in addition to the helicity-2 mode of a massless graviton.

DHOST theories have been classified up to cubic order in second derivatives of the scalar field, and are divided in 41 classes whose general Lagrangian reads [88]

[TABLE]

where the functions $f_{2}$ and $f_{3}$ depend only on $\phi$ and $X\equiv\nabla_{\mu}\phi\nabla^{\mu}\phi$ and the tensors $C_{(2)}$ and $C_{(3)}$ are the most general tensors constructed with the metric $g_{\mu\nu}$ and the first derivative of the scalar field $\phi_{\mu}\equiv\nabla_{\mu}\phi$. Not all the 41 classes are healthy and many of them feature pathologies such as the absence of tensor modes or gradient instabilities at linear order [89, 90].

Among the healthy classes a further restriction can be imposed if we require that tensor modes propagate with velocity $c_{T}=1$ on cosmological background. The resulting DHOST theory is

[TABLE]

where $A_{3,4,5}$ are functions of $\phi$ and $X$, and the following relations have to be imposed to satisfy the degeneracy conditions that remove the Ostrogradsky ghost

[TABLE]

Notice that we can add a $K$-essence and a cubic Horndeski term to the above action (3.1.2), i.e. $K(\phi,X)+G_{3}(\phi,X)\Box\phi$, without altering the degeneracy conditions, to obtain the most general scalar-tensor theory with $c_{T}=1$. Horndeski theories discussed in the previous subsection are a special subclass of (3.1.2) given by $A_{3,4,5}=f_{2_{X}}=0$.

In order to concretely study the time dependence of the effective tensor Planck mass $M_{\rm eff}$ within DHOST theories, let us focus on a specific model that admits a tracker solution evolving toward a de Sitter fixed point responsible for the present time cosmic acceleration [91] (see also [92]):

[TABLE]

where $M_{\rm Pl}$ is the Planck mass, $c_{2},c_{3},c_{4},\beta$ are free constants and $\Lambda_{3}$ is a strong coupling scale. The tracker solution is characterised by the property that $\dot{\phi}=\xi/H$ ($\xi=$ const.) during both matter domination and the dS phase, each era with its own value of $\xi$. If $\beta=0$ in (3.18), which corresponds to the case of beyond Horndeski theories [83], one recovers a unique value for $\xi$ throughout the evolution of the universe [93]999We do not consider this case here since it is excluded from the constraint derived in [94] to avoid decay of GWs into Dark Energy. This constraint gives $A_{3}=0$ and therefore $\beta=-8c_{4}$. As a consequence, $\beta=0$ would imply also $c_{4}=0$ and the model would simply reduce to GR + K-essence + cubic Galileon, falling in the class of scenarios we discussed in the previous subsection.. Tracker solutions are important because they allow to constrain the model parameters independently of initial conditions (see [95] for other tracker and scaling solutions in DHOST theories).

Working with dimensionless quantities, for this model we have

[TABLE]

and therefore the ratio between the luminosity distances, see eq (2.17), reads

[TABLE]

and the function $\delta(z)$ in equation (1.2) is given by

[TABLE]

Thus, the GW luminosity distance is controlled by the dynamics of the scalar field and of the scale factor over the tracker solution.

In Fig. 4 we show the numerical solution for the two quantities above, eqs (3.20) and (3.21), and the corresponding fit with the phenomenological parametrizations (2.31) and (2.37). In agreement with the discussion in Section 2.2.4 we see that, even if in this case the fit to $\delta(z)$ obtained from eq. (2.37) is quite poor, and misses the presence of a local minimum, still the fit to the directly observed quantity $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ provided by eq. (2.31) is qualitatively acceptable, and catches well the overall properties of the redshift dependent profile for $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$. As a general feature, the ratio of the luminosity distances reaches a plateau in the matter dominated era where the deviation with respect to GR vanishes; instead the function $\delta(z)$ features a local minimum and then goes to zero at large redshifts and in the future time where the universe reaches the de Sitter phase. In order to perform a proper fit with the parametrisation (2.31), a precise determination of the time today $t_{0}$ is required by fitting cosmological parameters. Indeed, it is easy to realise that a change in $t_{0}$ can shift the local minimum at smaller or negative redshifts in order to get a better agreement with (2.31). A similar determination is beyond the scope of this work and the results presented here should only be seen as indicative of possible behaviors of $d_{L}^{\rm gw}(z)/d_{L}^{\rm em}(z)$ and $\delta(z)$ for self-accelerating cosmologies in DHOST theories.

Additional parameterisations

For these cosmologies, at the de Sitter fixed point we have $\delta=0$ and therefore, in order to improve the fit, one can modify the parametrization (2.31) to get $\delta(z\ll 1)\sim 0$ via the addition of some extra terms (with no new parameters) to eq (2.37). A possible choice is

[TABLE]

which has the desired asymptotic values at large and small redshift. By reverse engineering integrating eq. (2.37), we can analytically reconstruct a profile for $\Xi(z)$ as

[TABLE]

which introduces an overall exponential factor to the parametrization (2.31). This profile maintains the same number of free parameters as (2.31) but with the disadvantage of being more complicated. Another alternative parametrization is

[TABLE]

Both the polynomial-exponential parametrization (3.23) and the exponential parametrization (3.24) obey the asymptotic relations (2.33) and (2.34). The three parametrizations are shown in Fig. 5.

To compare the three proposals (2.31), (3.23) and (3.24) and to quantify which fits the theoretical curve $\Xi(z)$ “better”, we perform a simple Bayesian model selection, calculating the Bayes Information Criterion (BIC) [96] or the Akaike Information Criterion (AIC) [97] for the best fits and treating the theoretical curve as data. Our best fits are obtained through a nonlinear fit of the theoretical points with the parametrization curve(s) with the least squares method. Although differences in parametrizations are not especially relevant for LISA standard sirens at high redshifts (e.g., MBH mergers) due to the large error bars, the issue of model reconstruction and the use of Information Criteria may be useful when considering low-$z$ standard sirens (for instance, stellar-mass black holes), or if one uses parametrizations with a different number of free parameters.

The parametrization with smaller BIC or AIC is to be preferred. In the present case, the difference $\Delta:=|\textrm{(IC model 1)}-\textrm{(IC model 2)}|$ is the same for the BIC and the AIC, since the number of data points and the number of free parameters are exactly the same for all parametrizations. According to the classification of [98, 99], one finds weak evidence if $\Delta<2$, positive evidence if $2\leqslant\Delta<6$, strong evidence if $6\leqslant\Delta<10$, and very strong evidence if $\Delta\geqslant 10$. A different number of data points can give rise to different orderings in this triple comparison, a fact to keep in mind when using real data. We make it clear that we are comparing parametrizations with a theoretical curve, which may give different results with respect to when one would use a sparse set of data points, as it will be the case for standard sirens.

Having said that, with very strong evidence ($\Delta\gg 10$) we find that fitting the $\Xi(z)$ data up to $z=5$, the polynomial-exponential parametrization (3.23) is preferred over the vanilla parametrization (2.31), which in turn is preferred over the exponential parametrization (3.24). We represent this situation as POL$>$VAN$>$EXP. Fitting $\Xi(z)$ data up to $z=0.8$, the ordering changes as EXP$>$POL$>$VAN. As one can see from Fig. 6, the exponential parametrization (3.24) is considerably more accurate.

Notice that one could contemplate to use different parametrizations depending on whether one wants to fit what we measure directly (the function $\Xi(z)$) or to reconstruct the function $\delta(z)$ from data. In the second case, if we fit the reconstructed $\delta(z)$ up to $z=5$, then EXP$>$POL$>$VAN and the exponential parametrization (3.24) is the most economic way between the three to achieve that goal. Our conclusion is that:

• •

For $\Xi(z)$, the polynomial-exponential parametrization (3.23) is the second best at small redshift and the best as a global fit, while for $\delta(z)$ it is the second best both at small redshift and as a global fit.

• •

For $\Xi(z)$, the exponential parametrization (3.24) is the best at small redshift but the worst as a global fit, while for $\delta(z)$ it is the best both at small redshift and as a global fit.

While the polynomial-exponential parametrization (3.23) is made ad hoc for DHOST and is the best global fit for $\Xi$, the exponential parametrization (3.24) works as a small-$z$ fit for $\Xi$ not only for DHOST models, but also for models with extra dimensions and models of quantum gravity, as we will see in section 3.4. Therefore, both parametrizations have their own advantages and drawbacks.

3.2 GW propagation in nonlocal infrared modifications of gravity

“Nonlocal infrared modifications of gravity” is a generic denomination for models in which the fundamental theory is local, according to the standard lore of quantum field theory, but non-local terms, relevant in the infrared (IR), emerge at some effective level. A situation where this can happen is when one takes into account quantum fluctuations. In that case, when in the spectrum there are massless particles (such as the graviton in GR) or, more generally, particles that are light with respect to the relevant energy scales, the quantum effective action, which is the quantity whose variation gives the equation of motion of the vacuum expectation value of the fields, develops nonlocal terms. Terms that involve inverse powers of the d’Alembertian are particularly interesting for cosmological applications, since they become important in the IR. The idea that quantum gravity at large distances could induce cosmological effects is quite old, see e.g. [100], and several hints on non-trivial physics at large distances emerge for instance from the study of IR effects in de Sitter space [101, 102, 103, 104, 105, 106, 107, 108]. However, a first-principle understanding of the IR limit of quantum gravity is a highly non-trivial task, and it makes sense to begin with a phenomenological approach, where one studies what sort of nonlocal terms could give a viable and interesting cosmology. In this spirit, a nonlocal gravity model was proposed by Deser and Woodard (DW) [109, 110] (generalizing earlier work in [111]), based on the quantum effective action

[TABLE]

with $f(X)$ a function chosen so to obtain the desired background evolution [112] (see [113] for review). In the DW model the nonlocal term involves a dimensionless function $f(X)$. Another possibility that has been much investigated recently is that the nonlocal term rather involves a mass scale. An explicit example of a model of this class is provided by the Dvali-Gabadadze-Porrati (DGP) model [114] that, when projected onto the four-dimensional brane and linearized around Minkowski space, is equivalent to a nonlocal theory for the fluctuations $h_{\mu\nu}$ whose covariantization, at least to linear order in $h_{\mu\nu}$, reads [115, 116]

[TABLE]

where the mass scale $m=2M_{5}^{3}/M_{4}^{2}$ emerges from a combination of the 5-dimensional and 4-dimensional Planck masses. While DGP is not cosmologically viable on its self-accelerating branch, several other models with nonlocalities associated to a mass term have been studied phenomenologically. In particular, a model of this type was proposed in [117] to introduce the degravitation idea, and is based on the modified Einstein equations

[TABLE]

where $m$ is a mass scale associated to the nonlocal operator $\Box^{-1}$. However, a drawback of eq. (3.28) the energy-momentum tensor is not automatically conserved, since in curved space $[\nabla^{\mu},\Box^{-1}]\neq 0$ and therefore the Bianchi identity $\nabla^{\mu}G_{\mu\nu}=0$ no longer ensures $\nabla^{\mu}T_{\mu\nu}=0$. On the other hand, a symmetric tensor $S_{\mu\nu}$ can always be decomposed as

[TABLE]

where $S_{\mu\nu}^{\rm T}$ is the transverse part of $S_{\mu\nu}$, $\nabla^{\mu}S_{\mu\nu}^{\rm T}=0$. This decomposition was used in [118] to modify eq. (3.28) into

[TABLE]

where the superscript “T” denotes the extraction of the transverse part from of the tensor $\left(\Box^{-1}G_{\mu\nu}\right)$, so that energy–momentum conservation becomes automatic. The cosmological evolution of this model turned out to be unstable, already at the background level, so this model is not phenomenologically viable. Several variants of this idea have then been explored, and it has been found that it is quite difficult to build a model that is cosmologically viable. A rather extensive study of the various possibilities (see [119, 51] for reviews) has finally pinned down two particularly interesting models. The first, proposed in [120], is a variant of eq. (3.30) where the $\Box^{-1}$ operator acts on the Ricci scalar rather than on the Einstein tensor, and is defined by the nonlocal equation of motion

[TABLE]

(where the factor $1/3$ is a useful normalization for the mass scale $m$). We refer to it as the “RT” model, where R stands for the Ricci scalar and T for the extraction of the transverse part. The second, proposed in [121], is instead defined in terms of a quantum effective action

[TABLE]

We refer to it as the “RR” model, after the two occurrences of the Ricci scalar in the nonlocal term. It turns out that the RT and RR models are related by the fact that the equations of motion derived from eq. (3.32), when linearized over Minkowski space, are the same as the linearization of eq. (3.31). However, if one goes beyond linear order, or if one linearizes over a different background such as FRW, the two models are different.

Detailed studies [120, 121, 122, 123, 124, 125, 51] have shown that the RR and the RT models have a viable cosmological background evolution, with the nonlocal term driving accelerated expansion; they have stable cosmological perturbations both in the scalar and tensor sector; tensor perturbations propagate at the speed of light, thus complying with the limit from GW170817/GRB 170817A; both models fit CMB, BAO, SNe and structure formation data at a level statistically equivalent to $\Lambda$CDM (with the same number of parameters, since $m$ replaces $\Lambda$). Very recently, however, it has been found in ref. [126] (elaborating on a discussion in [127]) that the RR model is ruled out by the limit on the time variation of Newton’s constant derived from Lunar Laser Ranging (LLR). This is due to the fact that the RR model predicts a time-varying effective Newton’s constant $G_{\rm eff}(z)$ on cosmological scales, and does not have a screening mechanism on short scales. The analysis in [126] then shows that the same time dependence would appear at the scale of solar system and of the Earth-Moon system, and would then violate the LLR bound. The RT model, in contrast, predicts a constant $G_{\rm eff}=G$ for modes well inside the horizon, and therefore passes the LLR bound. In [126] it is also shown that a screening mechanism that was proposed for the DW model is actually incorrect, and that the DW model is also ruled out by the LLR bound. Thus, it appears that, in this class of nonlocal models, the RT model is the only viable one (see also [128] for an analysis of LLR constraints that rules out other nonlocal gravity models studied in the literature). It is quite remarkable that the whole ensemble of conditions that have been imposed is so strong to single out a single model out of a rather broad class.101010Recently, stimulated by the observation in [126], a new variant of the DW model has been proposed in [129]. This new model solves a part of the problem raised, in the sense that a nonlocal quantity $Y$ used in its construction indeed changes sign between the static and the cosmological solution, in contrast to the original variable $X=\Box^{-1}R$, allowing for a screening mechanism based on setting to zero the function $f(Y)$ that defines the nonlocal model, when $Y<0$. However, the rather general analysis in [127, 126] shows that, for both the RR model and the “old” DW model, once one matches the static solution in the near region to the cosmological solution at large distances, one finds $X_{\rm tot}(t,r)\simeq X_{\rm static}(r)+X_{\rm cosmo}(t)$. Since $X_{\rm cosmo}(t)<0$ and numerically dominates over $X_{\rm static}(r)$, it follows that $X_{\rm tot}(t,r)<0$ independently of the sign of $X_{\rm static}(r)$ (which, furthermore, in this case is also negative), and the screening mechanism proposed for the “old” DW model, based on setting $f(X)=0$ for $X>0$, cannot be realized. In general, the same could happen in the new model for the field $Y=\Box^{-1}(g^{\mu\nu}\partial_{\mu}X\partial_{\nu}X)$, so that $Y_{\rm tot}(t,r)\simeq Y_{\rm static}(r)+Y_{\rm cosmo}(t)$ with $Y_{\rm cosmo}(t)$ positive and numerically dominating over $Y_{\rm static}(r)$. If this is the case, again the model would never be in the regime $Y<0$ where the screening mechanism acts, and would then be ruled out. Further work is then needed to assess whether this model is indeed viable.

These nonlocal gravity models predict modified GW propagation, showing again that this phenomenon is completely general in modified gravity. The RR model prediction is fitted extremely well by the parametrization (2.31), with $\Xi_{0}\simeq 0.970$ and $n\simeq 2.5$; indeed, the prediction of the RR model was the original motivation for introducing this parametrization [16]. The RT model also predicts a modified GW luminosity distance. The equation for the tensor perturbations in this model was computed in [124] and has again the form (2.9), so GWs propagate at the speed of light, with a function $\delta(z)$ determined by the background evolution of some auxiliary field that are introduced to rewrite the model in a local form. Inserting the expression for the background evolution of these auxiliary fields found from the numerical integration of the modified Einstein equation, we get the results shown in Fig. 7 for $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ and $\delta(z)$. In these plots the blue solid line are the predictions of the model, and the magenta dashed lines are the fit to the parametrization (2.31) for $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ [and the corresponding parametrization (2.37) for $\delta(z)$], with the best-fit parameters

[TABLE]

We see that the parametrization works extremely well, and the model predicts a relatively large value of $\Xi_{0}$, which differs from the GR value $\Xi_{0}=1$ by about $6.6\%$.

It should be observed that, in the RT model, for modes well inside the horizon $G_{\rm eff}=G$ [130, 131]. Thus, this is an example of a model where there is a modified GW propagation even if the effective Newton constant, or the effective Planck mass, does not evolve with time, and hence the relation (2.17) does not hold.

3.3 GW propagation in bigravity

The linearised equations controlling the propagation of gravitational waves over cosmological distances can contain couplings with additional fields, for example in scenarios with broken spatial diffeomorphisms or with anisotropic stress. The extra dynamics is typically the result of a dark energy or a dark matter component. In such cases, the evolution of tensor modes can be characterised by new phenomena, such as transitions and oscillations between different states, analogous to neutrino oscillations. One ought to expect that in these scenarios a minimal parametrization for the GW luminosity distance, as the one discussed so far, is not sufficient to fully describe the interplay between different fields. A model-specific analysis is then necessary, and we shall consider one in this subsection. One of the best-developed frameworks that exhibits these extra effects is bigravity, where a massless spin-2 field couples to a massive spin-2 particle in a specifically devised fashion so as to avoid exciting ghostly (i.e. unstable) degrees of freedom. We henceforth concentrate on this theory but stress here that alternative scenarios exist (e.g. dark energy models with vector fields breaking space diffeomorphisms [132, 133]) with similar phenomenological consequences for the propagation of gravitational waves.

A consistent theory of bigravity, free of Ostrogradsky instabilities at the fully non-linear level, has been proposed by Hassan and Rosen (HR) [134], adding an Einstein-Hilbert term for the reference metric of the so-called dRGT theory of massive gravity [135] (see [136, 137, 138] for reviews). Bigravity admits FRW cosmological solutions describing late-time acceleration that differ from $\Lambda$CDM at the level of the background as well as for the dynamics of cosmological fluctuations [139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150]. In this section we study the evolution of tensor modes around homogeneous FRW configurations, developing techniques that will allow us to go beyond existing literature. We will be specifically concerned with GWs propagation (as opposed to generation) and will therefore be allowed to neglect the non-linearities crucial in the strong-gravity GWs generation regime (see e.g. [54] for a recent analysis). Our analysis is focussed on the traceless symmetric part of the two tensor sectors: it will not include the scalar and the (typically decaying) vector degrees of freedom. For other works on GW propagation in bigravity, see for example [151, 152, 153, 154, 155].

We are going to address the question of how the coupling between different modes can affect the graviton propagation in a regime outside a late time de Sitter era, which is the one usually considered when studying oscillation effects in bigravity [156, 157]. This question is particularly relevant for the purposes of this work, since GWs emitted from LISA standard sirens at large redshift can probe phases of cosmological expansion that are not captured by a pure de Sitter space approximation. We provide in Appendix B.1 a brief review of bigravity theory, as well as the equations governing the homogeneous background evolution. Tensor fluctuations $h_{ij}^{(1,\,2)}$ are defined around two FRW line elements, associated to each one of the metrics involved:

[TABLE]

where $c(\tau)$ is the speed of the second tensor fluctuation. In what follows the ratio of scale factors is denoted by

[TABLE]

while $\mathcal{H}=a^{\prime}/a$ is the Hubble parameter corresponding to the first metric. The evolution of the tensor perturbations in bigravity is described by a coupled system of linearised equations for the two tensor modes $h_{1}$ and $h_{2}$ [142]

[TABLE]

where we have dropped the tensorial indices since the propagation is the same for each transverse-traceless polarisation. The constant $\kappa$ controls the relative size of the strength of gravitational interactions in the two sectors, while $m$ sets the scale of the bare graviton mass. The quantity $f_{1}(\tau)$ is a cubic function in $\xi(\tau)$ that depends on the bigravity parameters $a_{i}$:

[TABLE]

The time-dependent coefficients in (3.37) are controlled by the background solutions. We will consider the branch of solutions where scalar and vector modes are not strongly coupled. In this branch, the propagation speed of the second tensor $h_{2}$ is determined by [141]

[TABLE]

and in what follows we assume that matter only couples with the tensor perturbation $h_{1}$. Equation (3.37) implies that the cosmological propagation of GWs in bigravity is characterized by three distinctive effects:

• (i)

the two tensor perturbations $h_{1,2}$ propagate at different speeds ($c\neq 1$ if $\xi^{\prime}\neq 0$),

• (ii)

they have different friction terms,

• (iii)

they mix due to the non-diagonal mass matrix.

In general, given the time dependence of the parameters in (3.37), the propagation of GWs can not be solved analytically, and the system of two tensor modes cannot be diagonalized. On the other hand, in our case we can exploit the fact that the frequency of the GW ($f\sim 10^{-2}$Hz in the LISA band) is much larger than the universe expansion rate, $\mathcal{H}_{0}\sim 10^{-18}$Hz. Thus, the time variation of the parameters is small compared to the frequency of the GW and we can make use of a WKB expansion to obtain approximate analytical solutions for the tensor dynamics. That said, we nevertheless emphasize that for the range of redshifts probed with LISA standard sirens, $z\sim 1-5$, there is always some time dependence left in the parameters introduced by the scale factor $a(z)$. This is in contradistinction to LIGO sources at $z\ll 1$, for which the approximation $a\simeq 1$ can be consistently adopted, as done in previous analysis [156, 157].

The WKB solution for the system of equations (3.37) is described in Appendix B.2, and is summarised by the following expression

[TABLE]

The quantities $\theta_{1,2}$ and $\hat{E}$ correspond respectively to the eigenvalues and to the matrix of eigenvectors that would diagonalise the system in the approximation in which the parameters in eq (3.37) are constant. $\bar{\Phi}_{1,2}$ denote components of a vector controlling the mode amplitudes. The quantities $c_{1,2}$ are constant fixed by initial conditions while $\delta\theta\equiv\theta_{2}-\theta_{1}$.

We relegate technical details on these quantities and how to derive eqs. (3.40)–(3.41) to Appendix B.2. Here we emphasize that the possible mixing between $h_{1}$ and $h_{2}$ is controlled by the non-diagonal elements of the matrix of eigenvectors, i.e. $\hat{E}_{12}$ and $\hat{E}_{21}$. These entries are non-vanishing whenever the mass matrix in (3.37) is non-diagonal. On the other hand, we find that the size of the mixing is controlled by the relative difference of the velocities of each mode. In particular, for the mixing not to be suppressed in the regime of large $k$, one needs to require

[TABLE]

which follows directly from the analytic expression of $\hat{E}$ that can be found in Appendix B.2. If this inequality is satisfied one can have large mixings among modes with interesting phenomenological consequences, such as graviton oscillations. In Appendix B.2 we show that in the small mass regime, $m\sim\mathcal{H}_{0}$, this inequality cannot be satisfied for viable cosmological scenarios. Hence, for the rest of this Section we focus on a large mass limit $m\gg\mathcal{H}_{0}$. In this regime $\xi$ is approximately constant and, as a consequence, $c\simeq 1$ (recall (3.39)). In particular, we find

[TABLE]

where $\Omega_{\Lambda}$ is the density of DE. As a consistency check, we see that in the pure de Sitter limit, the speed is exactly luminal, $c_{dS}=1$. The condition to have mode mixing, eq (3.42), is satisfied for

[TABLE]

For LISA frequencies this leads to a bound $m\gtrsim 10^{8}\mathcal{H}_{0}$, which improves by a couple of orders of magnitude the LIGO detectability range $m\gtrsim 10^{11}\mathcal{H}_{0}$ [157] (recall that $\mathcal{H}_{0}\sim 10^{-33}$eV). We stress here that in the large mass regime the evolution of the universe cannot reproduce the observed accelerated expansion unless we include an effective cosmological constant term compensating the large mass in the Friedmann equation (see Appendix B.1). This means that, by restricting the analysis to the large mass regime, the present work does not fully probe bigravity: it does not capture the region of parameter space where $m\sim\mathcal{H}_{0}$, which is the most propitious for self-acceleration. On the other hand, our analysis can also serve as a proxy for other scenarios supporting oscillations in the luminosity distance at late time. A case in point are models with vector gauge fields [132, 133].

We henceforth restrict the analysis to the regime of (3.44) and investigate its phenomenological consequences. It is convenient to define an effective mass $m_{g}$ and a mixing angle $\theta_{g}$:

[TABLE]

Whenever $\left(c^{2}-1\right)\,k^{2}<m_{g}^{2}a^{2}$, one can expand the phases associated to each tensor mode as

[TABLE]

Analogously, the eigenvectors simplify to

[TABLE]

In order to estimate how large is the correction to the GW amplitude w.r.t. the $c=1$ case, we can parametrize the effective mass with a dimensionless constant $\beta$ via $(m_{g}\,a)\sim\beta\cdot(k\,\mathcal{H}_{0})^{1/2}$, which controls our complying with the large mass regime defined in (3.44). It follows that the largest correction to the matrix of eigenvectors (3.48) scales with $\beta^{-4}$. Therefore, provided that $\beta\gtrsim 5$, we can neglect this correction in the amplitude.

The mixing among modes affect the tensor speed and the luminosity distance. The modification in the propagation speed of the lightest tensor $h_{1}$, defined as $\alpha_{T}\,\equiv\,c^{2}_{T}-1$, is given by

[TABLE]

and scales as $\alpha_{T}\sim\beta^{-2}(k/\mathcal{H}_{0})^{-2}$. For LISA, $\alpha_{T}$ is smaller than $10^{-16}$ in the large mass limit. However, it might still be observable given that this is a cumulative effect over long travel distances. A prompt EM counterpart can give constraints of the order of

[TABLE]

where $\Delta t$ is the difference in the time of arrival and $D$ the distance to the source.

Next, we study how the mixing of the tensor modes $h_{1,2}$ can leave an imprint in the GW luminosity distance. Focussing on the regime $(m_{g}\,a)\gtrsim 5\cdot(k\,\mathcal{H}_{0})^{1/2}$, the only time dependence in the amplitude is introduced by the scale factor $a(\tau)$. Using the definitions (3.45) of the effective mass $m_{g}$ and the mixing angle $\theta_{g}$, the amplitude of the tensor component $h_{1}$ is

[TABLE]

where $h_{0}$ is the amplitude in the $\theta_{g}=0$ case, quantities depend on the time $\tau$, and $\tau_{e}$ is the time of emission. Also, $\theta_{1}^{2}=k^{2}-\mathcal{H}^{2}$, $\theta_{2}^{2}=k^{2}+m_{g}^{2}a^{2}-\mathcal{H}^{2}$, $\delta\theta=m_{g}^{2}a^{2}/(2k)+\mathcal{O}(k^{-3})$. As a result, if $m_{g}^{2}a^{2}\sim k\mathcal{H}_{0}$, the oscillation frequency is of order $\mathcal{H}_{0}$. A similar expression can be obtained for $|a\,h_{2}|$. One should note that while $|a\,h_{1}|$ is never larger than the initial value $h_{0}$, for $\theta_{g}<\pi/4$, $h_{2}^{2}$ can exceed $h_{0}^{2}$ (up to $4h_{0}^{2}$). In the opposite limit, $\theta_{g}>\pi/4$, $h_{2}^{2}$ is always less than $h_{0}^{2}$. This could be observed by comparing the left and right panels of Fig. 8, where the mass is fixed but the mixing angles vary.

We can now compare the luminosity distance of GWs in bigravity with the one of EM radiation, $d_{L}^{\text{gw}}/d_{L}^{\text{em}}$. We focus on the amplitude of the lightest tensor mode, $h_{1}$, which is the one we assume to be coupled to matter. From (3.51), we obtain

[TABLE]

We see that this ratio can become larger than one and display oscillatory patterns. Possible configurations range in principle111111It is actually not possible to push the theory to the asymptotic regions for the following reasons: (i) $\kappa\rightarrow 0$ corresponds to a vanishing kinetic term for the second metric so that it becomes infinitely strongly coupled; (ii) taking the limit $\kappa\rightarrow\infty$, the second metric decouples and gets effectively frozen to a fixed background value. between $\theta_{g}=0$ (corresponding to $\kappa\xi^{2}=0$) and $\theta_{g}=\pi/2$ (corresponding to $\kappa\xi^{2}\rightarrow\infty$), and the maximum mixing occurs at $\theta_{g}=\pi/4$ ($\kappa\xi^{2}=1$). Moreover, in the limit in which $\theta_{2}/\theta_{2}(t_{e})=1$ (which is a good approximation in the high $k$ limit), $d_{L}^{\text{gw}}/d_{L}^{\text{em}}$ is symmetric around $\pi/4$, i.e. $h_{1}^{2}(\pi/4-\varphi)=h_{1}^{2}(\pi/4+\varphi)$.

Given that the ratio (3.52) oscillates as a function of redshift, the two-parameter phenomenological parametrisation of Section 1 is not expected to perform well in this case, and an analysis specific for this model is necessary. We plot the oscillatory behavior of the GW luminosity distance for different masses and mixing angles in Fig. 9. For a qualitative analysis, we include error bars in the measurements of luminosity distances for a representative LISA catalogue (which we will discuss and use in Sections 4 and 5).

Figure 10 shows the projected constraints on GW oscillations in bigravity theories on the parameters $m_{g}$ and $\theta_{g}$. The excluded regions have been obtained comparing the high-mass prediction for the luminosity distance ratio (3.52) with several LISA catalogs, assuming fixed cosmology, and performing a $\chi^{2}$ analysis. Oscillation effects can be observed for masses $m_{g}\gtrsim 2\cdot 10^{-25}\,$eV. The amplitude of the oscillation, and thus its detectability, increases as $\theta_{g}$ approaches to $\pi/4$ from above or below. The mixing angle range where oscillations can be detected is $0.05\pi\lesssim\theta_{g}\lesssim 0.45\pi$. LISA will provide a $\sim 3$ order of magnitude improvement in mass sensitivity over the current LIGO limit, which probes $m_{g}\gtrsim 10^{-22}$eV [157], due to the larger oscillation baseline and the lower detection frequency. We conclude that the use of standard sirens will strengthen existing bounds towards smaller values for the mass.

However, one should note that these bounds are based on standard sirens with a single fiducial frequency, i.e. a monochromatic GW. Since GW oscillation effects depend strongly on the frequency, a coallescing binary would experience a time-dependent modulation of the amplitude as the orbital frequency increases, leading to a distinct signal. Remarkably, no electromagnetic counterparts are necessary to study such effects. As a representative example, we plot in Fig. 11 the strain of a GW signal from a massive BH binary as a function of frequency and its modification in bigravity. The GW oscillations lead to a distinctive frequency profile, which would be interesting to further characterise in order to establish to what extent it can be probed with LISA. Moreover, the initial wave packet emitted might decohere while traveling. In that case the event would be followed by an “echo” signal, which may be detectable if the mixing is sufficiently large [158]. The inclusion of decoherence effects and frequency dependence would provide further means to test GW oscillations. We leave this to future work.

3.4 Models with extra and varying dimensions

Various models of modified gravity predict that gravitational waves propagate in more than four dimensions, while photons and the other Standard Model fields are localized on a four-dimensional brane. Similar considerations apply to models of quantum gravity, where spacetimes with dimensions different from four often enter in theoretical attempts to find quantum theories of gravitational interactions.

Examples of modified gravity scenarios are models in the DGP class [114], where gravitational leakage into an extra fifth dimension has been proposed as a way to explain dark energy, since gravity becomes weaker at very large cosmological scales. In theories involving dimensions different from four, the GW propagation and the luminosity distance relation can differ from GR, and it is interesting to estimate whether GW observations can constrain such systems.

In [14, 159, 160], constraints have been imposed on the parameters entering the following phenomenological expression for the luminosity distance [12], suitable for modified gravity models with extra dimensions:

[TABLE]

with $R_{\rm c}$ is the characteristic length scale beyond which propagating GWs start to leak into the extra dimensions and $n_{\rm c}>0$ is a steepness parameter denoting the transition between GR and the onset of gravitational leakage. For sources emitting GWs accompanied by an electromagnetic counterpart, one can constrain large extra dimensions by comparing the distance inferred from the GW signal with the one inferred from the electromagnetic one [14, 159, 160]. In particular, using the GW170817 detection, the LIGO-Virgo collaboration [160] has recently set $90\%$ upper bounds on the spacetime topological dimension $D$, assuming fixed transition steepness and distance scale. For $R_{\rm c}$ between 1 and 20 Mpc, $4\leqslant D<4.5$. Conversely, fixing $D=5,6,7$, a transition steepness $n_{\rm c}\lesssim 10^{-1}$ is ruled out, while the $10\%$ lower limit on $R_{\rm c}$ is about 20 Mpc for $n_{\rm c}>10$ and between 20 and $10^{4}$ Mpc for $0.5<n_{\rm c}<10$.

Analogous constraints for quantum gravity (QG) can also be investigated. Quantization of the gravitational force intimately affects the dimensionality of spacetime. In fact, in QG approaches the effective dimension of spacetime is not a fixed external parameter, but a dynamical and scale-dependent one, which can change to a value typically smaller (in some theories, larger) than four at very small distances or very high energies [161, 162]. This dimensional flow [163] smoothly spans all scales of the theory, from ultra-short to large scales where a semi-classical, near-GR limit is reached. At an effective level, in a regime where spacetime emerges as a continuum, the dimensional flow is determined by the way volumes scale in position and momentum spaces and by the dispersion relations associated with effective kinetic operators. In turn, these features are encoded into the topological dimension $D$ of spacetime ($D=4$ in physical models) and three effective, scale-dependent “dimensions,” namely the Hausdorff (called also fractal) dimension $d_{\rm H}$ of position space, the Hausdorff dimension $d_{\rm H}^{k}$ of momentum space and the spectral dimension $d_{\rm S}$. In a classical spacetime, $d_{\rm H}=d_{\rm H}^{k}=d_{\rm S}=D$. A general argument based on the scaling of fields and their kinetic terms indicates that correlation functions and distances depend on these dimensions through a geometric parameter $\Gamma:=d_{\rm H}/2-d_{\rm H}^{k}/d_{\rm S}$ that combines them together [164]. The relation between $\Gamma$ and the parameter $\delta$ that expresses the modification of GW propagation (and that, as we have seen in section 2.2.2, in some theories is related to the time variation of the effective Planck mass) is

[TABLE]

while the correction to the GW luminosity distance in a spacetime characterized by dimensional flow is

[TABLE]

where $\gamma\neq 0$ is the value of the parameter $\Gamma$ at small scales, $\ell_{*}$ denotes a length scale at which quantum gravity effects can no longer be neglected (hence it is of order of the Planck scale $\ell_{\rm Pl}$) and, depending on the model, $\varepsilon$ is either a number $O(\gamma-1)$ or a random variable with zero average and unit variance simulating an intrinsic spacetime uncertainty, i.e., a geometry with microscopic fuzziness. Note that (3.55) only relies upon two assumptions. First, that there exists a continuum spacetime limit where the position-space measure, the momentum-space measure and the kinetic term of the tensor mode are deformed in a fairly generic way. Second, that these deformations are such that $d_{\rm S}\neq 0$. Theories where the spectral dimension vanishes in the UV, such as quantum gravity with non-local UV modifications and the low-energy limit of string field theory, must be treated separately (and do not give rise to observable corrections) [165].

The power-law correction in (3.55) of quantum-gravity origin has a strong resemblance with the similar expression (3.53) we discussed above for modified gravity scenarios. Comparing (3.55) with (3.53), we can match the parameters as follows. When the scales $\ell_{*}$ and $R_{\rm c}$ are considered as free parameters, they play the same role. However, in quantum gravity $\ell_{*}$ is expected to be very small or even Planckian, which implies that the power-law correction will be small when $\gamma-1<0$, large when $\gamma-1>0$ and $O(1)$ when $\gamma-1$ is a small positive number. Depending on the theory and the scales at which (3.55) holds, $\gamma$ can take very different values, being either negative, zero, or positive.

In most scenarios, such as string theory in the low-energy limit, group field theory, spin foams, loop quantum gravity, causal dynamical triangulations, asymptotic safety, Stelle gravity, Horava–Lifshitz gravity and various realizations of $\kappa$-Minkowski spacetime, $\gamma<1$ in the deep UV and the correction in (3.55) is Planck-suppressed and negligible, as one would expect from a generic reasoning based on perturbative effective field theory. However, in some QG proposals based on discrete pre-geometries (such as group field theory, spin foams and loop quantum gravity), there is a possibility that $\gamma-1$ is a small positive number at mesoscopic (far from the UV-regime) or near-IR scales for spacetimes described by certain semi-classical states. Hence, in these cases the correction in (3.55) may be non-negligible, even if $\ell_{*}$ is smaller by many orders of magnitude than any cosmological distance.

In analogy with the extra-dimensional case (3.53) studied in [160], one can expect that LIGO-Virgo data, and later LISA data, constrain the free parameters ($\ell_{*},\gamma$) of (3.55), providing a test for QG proposals. This is indeed the case. The model-independent analysis of [164] shows that theories with $\ell_{*}=O(\ell_{\rm Pl})$ and

[TABLE]

can fall into the detectability range of interferometers. Such small non-zero values of $\gamma-1$ can be obtained only in a near-IR regime. Most QG models do not achieve this (because, in this regime, either $\gamma-1<0$ or $\gamma-1$ is positive but too close to zero) and do not affect the luminosity distance. However, QG proposals where spacetime emerges from a discrete structure (group field theory, spin foams and loop quantum gravity), where $0<\gamma-1\ll 1$ ($\gamma$ close to 1), may leave an observable imprint depending on the considered quantum state of geometry (which can be described in terms of geometric indicators, such as the Hausdorff and spectral dimensions). In this latter class of QG models, the spacetime dimension turns out to be highly state-dependent, and hence the effective length scale characterising IR nearly-classical geometries may be several orders of magnitude above the Planck length. Whether such quantum states are realistic and viable within the corresponding theoretical framework remains an open question and a potential subject of future study.

It is worth noting that the quantum correction in (3.55) is not the same as one could find, for instance, in particle-physics or atomic observables, where constraints on new physics are tight to an extreme. In most theories of quantum gravity, these corrections are actually not present, since only the gravitational sector of the dynamics is affected by non-perturbative phenomena such as the dimensional flow. This is indeed the case for two of the three QG proposals that support the parameter space (3.56), namely spin foams and loop quantum gravity. In these two models there is no modification in the Standard Model of particle physics, nor in any fundamental constant of nature. For group field theory though, this is less obvious. However, for theories where particle and atomic physics are also affected, corrections are still of a power-law type $(\ell/\ell_{*})^{\alpha}$ but, typically, with a power $\alpha$ different from the special combination of dimensions (3.54) [164]. Therefore, GW and particle-physics bounds are (at least partially) independent, leading to complementary constraints to the fundamental theories. On the other hand, solar-system tests can yield stronger bounds on $\gamma$ in QG scenarios where also the scalar sector is affected [164]. However, in many QGs models the correction to the Newton potential is not known, while the spin-2 sector is under better control and GWs can be a robust and, for the time being, unique way to constrain such models [164].

To conclude, we comment on the parametrization (2.31) for the systems considered in this subsection. For models that predict relations such as (3.53) and (3.55), the parametrization (2.31) is no longer appropriate. Indeed, (2.31) yields $d_{L}^{\rm gw}(z)/d_{L}^{\rm em}(z)\simeq 1+O(z)$ at small redshift, while (3.53), for $D>4$ and $d_{L}^{\rm em}\ll R_{\rm c}$, gives

[TABLE]

which is of the form $1+{\cal O}(z)$ only for $n_{\rm c}=1$. On the other hand, assuming that the electromagnetic sector is not modified, (3.55) gives

[TABLE]

which is of the form $1+{\cal O}(z)$ only for $\gamma=2$, valid only in two specific models which leave no imprint on the luminosity distance [164]. Moreover, quite importantly, the ratio $d_{L}^{\rm gw}/d_{L}^{\rm em}$ in (3.53) and (3.55) does not tend to a constant at large redshifts – which is the interesting range for the LISA Standard Sirens that we consider in this work – hence it can not be reproduced well with the parametrization (2.31) in the high-$z$ regime.

A possible parametrization inspired by the low-$z$ behaviour (3.57) and (3.58), but such that $d_{L}^{\rm gw}/d_{L}^{\rm em}\simeq{\rm const}$ at large redshift, is given by (3.24). At large $z$, this tends to the constant value $\Xi_{0}$, while in the low-$z$ limit it goes to

[TABLE]

reproducing (3.57) when $\Xi_{0}=1+(D-4)(H_{0}R_{\rm c})^{-n_{\rm c}}/(2n_{\rm c})$ and $n=n_{\rm c}$, and (3.55) when $\Xi_{0}=1+\varepsilon(H_{0}\ell_{*})^{1-\gamma}$ and $n=\gamma-1$. A similar behavior can be obtained with the parametrization (2.40).

4 Construction of catalogs of LISA standard sirens

LISA will detect three classes of binary black holes (BHBs) at cosmological distances [11], namely, massive BHBs (MBHB), extreme mass ratio inspirals (EMRIs) and stellar mass BHBs. These GW sources have different properties and they are expected to be observed in different redshift ranges. In particular, stellar mass BHBs will mainly be detected at $0.01\lesssim z\lesssim 0.1$ [166], EMRIs at $0.1\lesssim z\lesssim 2-3$ [167] and MBHBs at $z\gtrsim 1$ [168], implying that LISA will constitute a cosmological probe able to map the expansion of the universe from local to very large distances [169]. Nevertheless only MBHBs are expected to produce powerful EM counterparts, since they are believed to merge in a gas rich environment that may power EM emission through jets, disk winds or accretion: this fact will allow us to determine precisely the object positions in the sky. MBHBs will thus likely be the main standard siren sources for LISA, though the number of EM counterpart detections might not exceed few tens of events [35]. They will also provide useful data to test modified theories of gravity through a modified cosmological propagation of GWs. In this section we outline our strategy to build simulated catalogs of LISA MBHBs with EM counterparts, while in section 5 we will use these catalogs to constrain modified gravity theories. To produce catalogs of MBHB standard sirens with EM counterpart, we employ a similar approach to the one adopted in [35]. In what follows, we summarize the procedure to construct such catalogs, highlighting the differences with respect to the analysis of [35].

To describe the cosmological evolution of massive black holes, and specifically MBHBs, we use the semi-analytic galaxy formation model of [170] (with later incremental improvements described in [171, 172]). This model was extensively used to predict massive black hole merger rates/stochastic background levels for various LISA/pulsar timing array configurations, respectively in [168, 173] and [174, 175]; to study LISA standard sirens in [35]; to investigate LISA event rates for extreme mass-ratio inspirals in [167]; and to study projected LISA ringdown constraints on the no-hair theorem in [176].

We adopt here the same scenario used in [168, 172]. In more detail, the model follows the evolution of baryonic structures on top of a Dark Matter merger tree produced with an extended Press–Schechter formalism [177], modified to reproduce the results of numerical N-body simulations [178]. The baryonic structures that are evolved include a chemically pristine diffuse inter-galactic medium, which is either shock heated to the virial temperature of Dark Matter halos, or flows into the halos along cold filaments [179]; a cold, chemically enriched inter-stellar medium (in both a disk and bulge component); stellar disks and bulges, forming from the inter-stellar medium (with supernova feedback included in the star formation description); nuclear star clusters, forming from in situ star formation and from the migration of globular clusters to the nuclear region [172, 180]; and massive black holes, accreting from nuclear gas reservoirs (which, in turn, form after major galactic mergers or galactic disk instabilities, which trigger bursts of star formation in the bulge and inflows of gas to the nuclear region). The black holes affect the overall evolution via the feedback of their active galactic nucleus (AGN) phases, which ensures also that the observed local correlations between the properties of massive black holes and those of their galactic hosts are correctly reproduced [170, 171, 181].

Most relevant for our purposes are two ingredients of the model, namely the prescriptions for the initial mass function of the massive black hole seeds, and for the “delays” between galaxy and massive black hole mergers. While we refer to [172] for a complete description of these two aspects of our model, let us here briefly recall their main features.

Regarding the seeding model, we adopt two distinct scenarios: a light seed scenario in which the massive black hole population evolves from the remnants of population III (popIII) stars (with remnant masses $\sim 100M_{\odot}$) at high redshifts $z\gtrsim 15$ [182]; and a heavy seed scenario in which the black hole seeds form, again at high redshifts $z\gtrsim 15$, from bar instabilities of protogalactic disks, with relatively large seed masses $\sim 10^{5}M_{\odot}$ [183]. As for the delay times, our model tracks the evolution of Dark Matter halos driven by dynamical friction, from the moment when they first touch to the final halo/galaxy merger [184], including also environmental effects such as tidal disruption and evaporation [185]. Nevertheless, the evolution of MBHBs that form in the aftermath of a galaxy merger is still poorly understood. Indeed, while MBHBs are generally thought to sink quite efficiently down to separations of a few pc, under the effect of dynamical friction from the background stars and gas, their evolution from pc-separation down to coalescence depends sensitively on the distribution and velocities of stars and gas in the nuclear region (see e.g. [186] for a review). As a result, MBHBs may even stall at pc separations (the so-called final pc problem), at least in the absence of triple massive black hole systems (which, however, naturally form as a result of the hierarchical nature of galaxy evolution [175, 187, 173]). To account for these uncertainties, we adopt a model in which the final pc problem is solved efficiently, i.e. we assume no delays between the merger of galaxies/halos and black hole coalescence; and a model in which we employ a more physical description of the evolution of MBHBs. In more detail, in the latter scenario we account for the interaction of the MBHB with the nuclear gas (if that is present), which may trigger mergers on timescales as short as $10^{7}$ yr, and with the stellar background (“stellar hardening”), which typically induces black hole mergers on timescales of a few Gyr (see e.g. [186] for a review of the physics of these processes). We also account for the formation of triple massive black hole systems, which trigger coalescences on timescales again of a few Gyr, using the simplified model described in [172] (see [188, 189, 190, 175, 187, 173] for more work on triple massive black hole interactions).

The results of these semi-analytical simulations are catalogs of MBHB mergers across the whole cosmic evolution, specifically from $z\simeq 20$ to $z=0$. Using the data of these simulations, in particular the (luminosity) distance and masses of MBHB mergers, we are able to characterize the GW emission of all merger events within the catalogs, and thus to compute the expected GW signal when it reaches the detector. We follow the approach used in [168, 35] to assess the parameter estimation capabilities of LISA for such catalogs.

We use Fourier-domain inspiral-only precessing waveforms using 3.5 post-Newtonian (PN) phase evolution with 2PN spin-spin and 3.5PN spin-orbit couplings [191], and spin precession equations at 2PN spin-spin and 3.5PN spin-orbit orders. Using a set of precessing inpiral-merger-ringdown waveforms, in [168] was shown that the error on the luminosity distance decreased by a factor of the SNR ratio when compared to an inspiral-only analysis, while the errors on the sky location area decreased by a factor of the SNR squared.

In this analysis, we compute Fisher matrices to assess the detector parameter estimation capabilities. We first define the noise-weighted inner product in the space of signals by (see e.g. [192, 193])

[TABLE]

where $\tilde{h}(f)$ is the Fourier transform of the waveform $h(t)$ which folds in the detector response, a star denotes complex conjugation, and $S_{n}(f)$ is the one-sided power spectral density of the detector noise. For the detector noise we use the LISA noise curve 121212In order to obtain our Eq (4.2), we square Eq (3) of [194] and we remove its $20/3$ factor (that comes from sky-averaging, which we do not perform). Eq. (1) of [194] results then equivalent to our Eq. (4.3). Notice also that the requirements reported in the public LISA proposal document [11] are simplified compared to the formulas in [194]. In Section 4.2 of [11] it is proposed to set requirements above $0.1$ mHz, for which low-frequency contributions are not too relevant. There are extra factors that affect the high-frequency acceleration noise and the low-frequency displacement noise, but in general the acceleration noise dominates at low frequencies and the displacement noise at high frequencies. It would be interesting, in future works, to analyse how our results change when changing sensitivity curves, whose features will be modified with a better definition of the instrument properties, and of astrophysical backgrounds. given by [194]

[TABLE]

where the acceleration noise is given by

[TABLE]

The other noise sources combine to

[TABLE]

and the confusion noise from unresolved binaries is approximated by

[TABLE]

with $A=(3/20)3.2665\times 10^{-44}$ Hz4/3, $s_{1}=3014.3$ Hz*-α*, $\alpha=1.183$, $s_{2}=2957.7$ Hz*-1*, and $\kappa=2.0928\times 10^{-3}$ Hz.

We follow [168, 35] and use a low-frequency approximation of the signal in LISA [195], in which the instrument is equivalent to two independent Michelson interferometers. We refer to [168] for details of the implementation.

The SNR of a signal in one interfermometer can be estimated with $\rho=(h|h)^{1/2}$, and the signals from the two interferometers combine as $\rho^{2}=\rho_{I}^{2}+\rho_{II}^{2}$. The Fisher matrix is defined as

[TABLE]

where $\partial h/\partial\theta^{i}$ is the derivative of the wave signal $h$ with respect to parameter $\theta^{i}$. The diagonal elements of its inverse $\Sigma=\Gamma^{-1}$ provide estimates of statistical errors due to the detector noise, as $\Delta\theta^{i}\approx\sqrt{\Sigma^{ii}}$. The Fisher matrices in the two detectors combine as $\Gamma=\Gamma^{I}+\Gamma^{II}$.

For each system in our catalogs, we estimate the error on the luminosity distance using the Fisher matrix method described above using an inspiral-only signal, and rescale it by the ratio of inspiral-merger-ringdown to inspiral-only SNR computed estimated using a phenomC waveform [196], as $\Delta d_{L}=(\rho_{\text{insp}}/\rho_{\text{IMR}})\sqrt{\Sigma^{d_{L}d_{L}}}$. To this experimental error one must add the uncertainties due to weak lensing and peculiar velocities, and we follow the prescriptions outlined in [35] for estimating and incorporating the error budget on $\Delta d_{L}$ of these effects.

Among all the LISA detections with SNR$>8$ we select the events with a sky localization $\Delta\Omega<10\,{\rm deg}^{2}$. As in [35] we characterize the EM emission at merger by assuming the production of an optical accretion-powered luminosity flare and of radio flares and jets, based on results from general-relativistic simulations of merging MBHBs in an external magnetic field [197, 198]. The semi-analytical simulations of MBHB catalogs described above contain all the necessary data (mass in stars, nuclear gas, etc.) needed to compute the magnitude of the EM emission of each MBHB merger in both the optical and radio bands. Following [35] we can thus determine the number of counterparts detected by future EM facilities, specifically LSST [199], SKA [200] and ELT [201]. Note however that over the sky location region of the source there might appear several EM transients that could likely be identified with the EM counterpart of the MBHB merger. As in [35] we will however assume that the true EM counterpart can be efficiently identified (e.g. through timing, EM spectrum characterization and evolution, or information on the host galaxy), and leave a more thorough investigation of this issue to future analyses. Spotting the EM counterpart helps localize the source in the sky, thus allowing one to determine the host galaxy and to improve the errors on the GW parameters by fixing the two sky localization angles.

The identification of the host galaxy in principle allows one to obtain the redshift of the MBHB event using either spectroscopic or photometric techniques. In what follows we distinguish two scenarios depending on the error estimates for these redshift measurements:

• •

Scenario 1 (realistic): Photometric measurements have a relative 1$\sigma$ error on $z$ given by $\Delta z_{\rm photo}=0.05(1+z)$ at best, according to Euclid [202], while spectroscopic measurements have a relative error given by $\Delta z_{\rm spect}=0.01(1+z)^{2}$, assuming we can define a volume-complete spectroscopic sample. However, the flux limit for that sample evolves with redshift. This may explain the redshift dependence, since luminosity distance evolves as approximately the square of $(1+z)$ at $z>1$ [203]. Although extrapolating those expressions beyond $z\sim 3$ might seem unrealistic [204], they are the only guesses we can use at the moment to estimate future redshift measurements at high redshift.

• •

Scenario 2 (optimistic): Photometric measurements have a relative 1$\sigma$ error on $z$ given by $\Delta z_{\rm photo}=0.03(1+z)$, while spectroscopic measurements have a fixed relative error given by $\Delta z_{\rm spect}=0.01$. Moreover here we assume a delensing procedure able to reduce by 50% the luminosity distance uncertainty due to weak lensing. This scenario is equivalent to the “optimistic scenario” considered in [35].

For both scenarios the redshift errors are propagated to the error budget of the luminosity distance using the fiducial cosmology, i.e. $\Lambda$CDM with Planck values for the cosmological parameters.

The total result of the whole procedure described above are thus mock catalogs of MBHB standard sirens, specified by values of the luminosity distance, redshift and total error on the luminosity distance. We simulate 22 4-year catalogs 131313Notice that although the nominal mission lifetime is 4 years (with possible extensions), the duty cycle is only $75\%$ of the total nominal time. It would be interesting to adjust our analysis to properly take into account how frequently data gaps occur: we leave a dedicated study of this point to a future work. for each MBHB astrophysical model (popIII, heavy seeds with and without delays), for a total of 66 catalogs and 264 years of data. For each catalog we then perform a quick Fisher matrix cosmological analysis assuming $\Lambda$CDM. We ranked the 22 catalogs of each MBHB astrophysical model according to the area of the 1$\sigma$ contour ellipses in the $(\Omega_{m},H_{0})$ parameter space of $\Lambda$CDM. This quantity provides us with a rough measure of the cosmological constraining power of each catalog. Finally we select the median catalog among all ranked 22 catalogs as the representative catalog for each MBHB astrophysical model, similarly to the recipe employed in [205]. We perform this procedure for Scenario 1 and then use the same selected catalogs for Scenario 2, but using different assumptions for the redshift errors as explained above. This allows us to directly compare results between the two different scenarios. The final product are thus six catalogs (corresponding to the three different MBHB astrophysical models, each one combined with the two scenarios for the redshift measurement accuracy) to use as mock data for the MCMC analysis presented in the following section. The six selected MBHB catalogs are plotted in Fig. 12.

5 Cosmological parameters and dark energy with LISA

We now use the mock catalogs discussed in the previous section to study the effect of LISA standard sirens on the estimate of the cosmological parameters, with particular emphasis on the DE sector, as described by the parameters $(w_{0},w_{a})$ for the DE equation of state, eq. (2.30), and by the parameters $(\Xi_{0},n)$ that enter eq. (2.31) for modified GW propagation. We adapt to LISA a similar analysis performed in [16] for the Einstein Telescope.

As discussed in section 4, we provide our results for three models for the formation of MBH binaries: heavy seeds and no delays (“hnd” for short), heavy seeds with a delay included and the Toomre parameter (see [35]) set to $Q=3$ (“hQ3”), and light seeds due to pop III stars (“popIII”). For each formation model, we use the two scenarios (1) and (2) of section 4 for estimating the error in the measurement of the redshift: thus, overall, we consider six different possibilities.

From the point of view of the cosmological model, we consider different cases. We will start with $\Lambda$CDM, to investigate how LISA standard sirens help to constrain its parameters. We will then study the extension of the DE sector obtained introducing a non-trivial DE equation of state, parametrized by $w_{0}$ only (without modified GW propagation, i.e. setting $\Xi_{0}=1$); finally, we will introduce modified GW propagation, extending the DE sector through the parameters $(w_{0},\Xi_{0})$ and fixing $w_{a}=0$ and $n$ to a reference value.141414We do not introduce more than two extra parameters at the same time in the DE sector, since the constraining power would decrease and the convergence time of the MCMC chains would become very long. Indeed, already with the two parameters $(w_{0},w_{a})$, we found that, because of the limited number of sources in the LISA catalogs, the MCMC chains do not resolve the degeneracy between these parameters and do not reach a good convergence. For each of these three cosmological models we run MCMCs for the six scenarios describing the formation model and the estimate of the redshift error, as discussed above.

Furthermore, for each of these $3\times 6$ cases we run separate MCMCs to compute the constraints that can be obtained with standard sirens only, and those obtained by combining them with CMB, BAO and SNe data. In that case we use the Planck temperature and polarization power spectra [46], the JLA SNe dataset [206] and a set of isotropic and anisotropic BAO data (see Section 3.3.1 of [51] for details).

To generate the catalogs discussed in section 4 we assume a fiducial cosmological model, that we always take to be $\Lambda$CDM with $H_{0}=67.64$ and $\Omega_{M}=0.3087$, which are the fiducial values obtained from the CMB+BAO+SNe dataset that we use. So, in particular, our fiducial values for the extra parameters in the DE sector are $w_{0}=-1$ and $\Xi_{0}=1$, and we use $\Lambda$CDM and standard perturbation theory within GR as fiducial theoretical framework for treating and analyzing our dataset. Our aim is to evaluate to what accuracy LISA, alone or in combination with other datasets, can reconstruct these fiducial values. We then generate our simulated catalog of events by assuming that, for a source at redshift $z_{i}$, the actual luminosity distance will be $d^{\Lambda\rm CDM}_{L}(z_{i};H_{0},\Omega_{M})$, with the above values of $H_{0}$ and $\Omega_{M}$. The “measured” value of the luminosity distance is then randomly extracted from a Gaussian distribution centered on this fiducial value, and with a width obtained from the estimate of the error on the luminosity distance discussed in section 4.

The relevant cosmological parameters for evaluating the theoretical values of the GW luminosity distance $d_{L}^{\,\rm gw}(z)$ are $\Omega_{M}$ and $H_{0}$ in the $\Lambda$CDM case, with the addition of $w_{0}$ in $w$CDM and the further addition of $\Xi_{0}$ in the $(\Xi_{0},w_{0})$ case. Apart from a constant addend coming from normalization, the logarithm of the likelihood assigned to a given set of values $(H_{0},\Omega_{M},w_{0},\Xi_{0})$ is given by

[TABLE]

where $N_{s}$ is the number of mock sources in the catalog, $d_{L}^{\,\rm gw}(z_{i};H_{0},\Omega_{M},w_{0},\Xi_{0})$ is the theoretical value of the GW luminosity distance for the $i$-th source and $d_{i}$ is the “measured” value of its luminosity distance contained in the catalog. The quantity $\sigma_{i}$ is the error on luminosity distance and it also takes into account the error on redshift determination, which is simply propagated to the luminosity distance using the fiducial $\Lambda$CDM cosmology. The used MCMC code explores the cosmological parameters space, accepting or rejecting the points according to a Metropolis-Hastings algorithm based on the likelihood specified above. The priors assumed on cosmological parameters are Gaussian (or truncated Gaussian) distributions with mean and standard deviation, in this order, given by: $67.8$ and $1.2$ for $H_{0}$ (in $\rm{km}\,\rm{s}^{-1}\rm{Mpc}^{-1}$), $0.02225$ and $0.00028$ for $\omega_{b}$=$\Omega_{b}h^{2}$ (baryon density fraction), $0.1192$ and $0.0027$ for $\omega_{c}=\Omega_{c}h^{2}$ (cold dark matter density fraction), $-1.0$ and $1.0$ for $w_{0}$ (restricted from $-3.0$ to $-0.3$), 1.0 and 0.5 for $\Xi_{0}$ (restricted to positive values).

As we see from Fig. 12, the number of total available standard sirens in our sample is quite limited. For instance, in the specific realization of the catalogs that we use there are 32, 12 and 9 sources in the “hnd”, “hQ3” and “popIII” scenarios, respectively. It is important to realize that, since these numbers are relatively small, the scatter of the mock data around their fiducial values unavoidably induces fluctuations in the reconstruction of the mean value of the cosmological parameters inferred from standard sirens alone, so that the reconstructed mean values in general will not coincide with those obtained from other cosmological observations such as CMB, BAO and SNe. Thus, depending on the specific realization of the catalog, the contours in parameter space of the likelihood obtained from standard sirens could show mild tensions with those obtained from other datasets. This is unavoidable, and will also happen in the actual experimental situation. Thus, in the plots shown below, the overlap (or lack of it) between different contours has little meaning, as it depends on the specific random realization of the catalog. What carries the important information is the relative size of the contours, that will tell us to what extent the addition of standard sirens to other datasets can improve our knowledge of the cosmological parameters. However, when combining contours from different datasets to obtain a combined estimate on the error on a parameter, one must be careful not to use a specific realization where the separate contours are in mild tension among each other.

5.1 $\Lambda$CDM

We begin by studying the effect of standard sirens on cosmological parameter estimation in $\Lambda$CDM, as a benchmark for our subsequent generalizations, and we will then introduce different extensions of the dark energy sector. As baseline $\Lambda$CDM model we use the standard set of six independent cosmological parameters: the Hubble parameter today $H_{0}=100h\,\rm{km}\,\rm{s}^{-1}\rm{Mpc}^{-1}$, the baryon density fraction today $\omega_{b}=\Omega_{b}h^{2}$ , the cold dark matter density fraction today $\omega_{c}=\Omega_{c}h^{2}$, the amplitude $A_{s}$ and tilt $n_{s}$ of the spectrum of primordial scalar perturbations, and the reionization optical depth $\tau_{\rm re}$. We keep the sum of neutrino masses fixed, at the value $\sum_{\nu}m_{\nu}=0.06$ eV, as in the Planck baseline analysis [207]. We then run a series of MCMC, using the CLASS Boltzmann code [208] (or our modification of it, in the case of non-trivial GW propagation).

In $\Lambda$CDM, assuming flatness, $d_{L}(z)$ depends only on $H_{0}$ and $\Omega_{M}$, so these are the parameters for which the inclusion of standard sirens has the most significant impact. Fig. 13 shows the two-dimensional likelihood in the $(\Omega_{M},H_{0})$ plane in $\Lambda$CDM, comparing the contribution from CMB+BAO+SNe (red) to the contribution from LISA standard sirens (gray), and the overall combined contours, for the three formation scenarios and the scenario (1) for the determination of the error on the luminosity distance, while Fig. 14 shows the results for the scenario (2).

In particular in the most favorable scenario (“hnd” seeds and optimistic errors on $d_{L}$), from the corresponding one-dimensional marginalized likelihood we find that, with standard sirens only, the relative error on $H_{0}$ is

[TABLE]

(which raises to $7.7\%$ in the “realistic” scenario with “hnd” seeds) and the one on $\Omega_{M}$ is $\Delta\Omega_{M}/\Omega_{M}=14.7\%$; using our CMB+BAO+SNe dataset we get instead $\Delta H_{0}/H_{0}=0.7\%$ and $\Delta\Omega_{M}/\Omega_{M}=2.1\%$; combining CMB+BAO+SNe+standard sirens we get

[TABLE]

Therefore, from MBH binaries at LISA, we do not find a significant improvement on the accuracy on $H_{0}$, compared to current results from CMB+BAO+SNe. It should however be observed that the measurement from standard sirens is still useful because it has completely different systematic errors from that obtained with CMB, BAO and SNe.

These results should be compared with those found in [35]. To perform the comparison, we should however note that the best results quoted in [35] referred to a LISA configuration with 5 Gm arms, and with all other design specifications to their most optimistic possible choices, while our results use the current LISA configuration with 2.5 Gm arms, given in [11]. The corresponding LISA sensitivity curve used in our study is quite different and generally gives worse results. For example the number of standard sirens used in our analysis is roughly half the number used in the 2016 paper [35].

Furthermore, at the methodological level there are some differences between our analysis and that in [35]. First, as discussed above, to generate our catalogs we have scattered the values $d^{\Lambda\rm CDM}_{L}(z_{i};H_{0},\Omega_{M})$, extracting the “measured” values of the luminosity distance from a Gaussian distribution centered on this fiducial value, and with a width $\Delta d_{L}$ determined by the projected error estimate. Ref. [35] did not consider scattered data, thus the statistical error due to scatter was not included in the analysis. The Fisher matrix analysis performed in [35] was performed on non-scattered data, but the MCMC analysis of the present paper takes also into account the statistical uncertainty due to the scattering and thus gives more realistic results.151515On the other hand, the MCMC analysis is based only on one realization of the dataset, in contrast to the Fisher matrix analysis which can be averaged over many realizations. However, as discussed in section 4, the catalog used is chosen as a typical “average” catalog, between many different realizations. Second, in the present analysis the degeneracies between the cosmological parameters are fully taken into account by our MCMC by freely varying all cosmological parameters of $\Lambda$CDM (and, in the following subsections, of its extensions), while in the Fisher matrix analysis of [35] the most stringent results (that eventually, together with all other assumptions mentioned above, led to the estimate that $H_{0}$ could be measured to $0.5\%$ with standard sirens only) where obtained by varying only $H_{0}$ and assuming a fixed prior on $\Omega_{M}$. The procedure of the present paper therefore gives a larger estimate of the error on the cosmological parameters that can be obtained with LISA standard sirens.

On the other hand, one should also be aware of the fact that, before drawing final conclusions on the sensitivity of LISA to $H_{0}$, as well as to the DE parameters that will be discussed below, much more work is needed. For instance, the counterpart model used in [35] makes a number of assumptions, which will only be validated when (if) LISA electromagnetic counterparts are actually observed. Another caveat is that the error on the localization used in this paper is based on the use of the inspiral waveform only, with corrections based on phenomenological waveforms, and it is still an open issue how this will change when including a full description of the merger and ringdown phases. The estimate of the redshift error also involves other factors, with respect to those that we have discussed, such as the exposure time, that will depend on the availability of telescope time and the pointing accuracy, which are currently difficult to estimate, as well as on the duration of the electromagnetic transient, whose estimation in turn involves the modeling of the counterpart. It should also be observed that the number of sources observed depends strongly on the sensitivity curve at low frequency, where there is potential room for improvement with respect to the sensitivity curve that we have adopted. Finally, in this paper we are only using MBH binaries with counterpart as standard sirens, and we are not considering stellar mass BH binaries and EMRIs, that are not expected to have a counterpart and must be treated with the statistical methods. These issues are currently under active investigations in the LISA Consortium, and are well beyond the scope of this paper. The results of this paper should therefore be understood as the forecast that can be obtained under the specific set of assumptions that we have made.

5.2 $w$CDM

We next add the parameter $w_{0}$, corresponding to the so-called $w$CDM model, where the DE equation of state is taken to be constant, or $w_{a}=0$ in eq. (2.30).161616A further natural extension would be to the $(w_{0},w_{a})$ set of parameters. However we have found that, with the limited number of sources in the LISA catalogs, it is not possible to disentangle the degeneracy between these two parameters, and the MCMC chains do not reach a good convergence. Since the model used to generate the catalog of sources is still $\Lambda$CDM, we are actually asking to what accuracy we can find back the fiducial value $w_{0}=-1$. Figs. 15 and 16 show the two-dimensional likelihood in the $(\Omega_{M},w_{0})$ plane, displaying the combined contribution from CMB + BAO + SNe (red), and the total combined result (blue), for the $3\times 2$ scenarios considered.

Even in the most optimistic scenario that we have considered, we learn from the plots that LISA standard sirens alone do not give any significant constraint on $w_{0}$, and, when combined to CMB+BAO+SNe data, they only induce a very marginal improvement. From the corresponding one-dimensional likelihoods, from CMB+BAO+SNe only, we find that $w_{0}$ can be reconstructed with the accuracy $\Delta w_{0}=0.045$. Combining CMB+BAO+SNe with standard sirens, in the “optimistic hnd” scenario we get

[TABLE]

so the improvement due to SMBH standard sirens is quite negligible. Basically the same results are obtained in all other scenarios considered.

These results significantly degrade the estimates presented in [35]. As discussed in section 5.1, the difference is due to the use of the updated sensitivity curve of LISA, the more realistic assumptions in the construction of the catalogs and estimates of the errors, and the fact that the degeneracies between cosmological parameters are now fully taken into account through a full MCMC. However, as for $H_{0}$, one should be aware of all the assumptions and uncertainties that entered in the construction of the source catalogs, and that could significantly alter the results. In particular these estimates can either be improved by combining the cosmological data collected with SMBHBs with the ones collected from other LISA sources, such as EMRIs and SOBHBs, or by extending the observational period of LISA. In the first case the different redshift ranges where EMRIs and SOBHBs are expected to be observed, will help breaking some degeneracies in the cosmological parameters, such as for example the degeneracy between $H_{0}$ and $\Omega_{M}$ in $\Lambda$CDM. In the second case instead the estimated errors are expected to improve roughly as $\sim\sqrt{N}$ or better, where $N$ is the number of SMBHB merger with EM counterpart observed by LISA, which linearly depends on the observational period. Both these improvements should reduce the error associated with the statistical scatter of the SMBHB data, which was not considered in [35] and in the new analysis performed here, which takes into account more realistic SBHB catalogs and an updated LISA noise sensitivity curve, appears to be the most relevant source of error. Other improvements can be achieved by a more realistic and detailed characterization of the emission and detection of the EM counterparts of SMBHB mergers, which here as in [35] has been modeled only with optical and radio EM emissions, but the consideration of other EM signatures (X-rays, $\gamma$-rays, …) could lead to a higher number of cosmologically useful events.

5.3 $(\Xi_{0},w_{0})$

We next extend the DE sector by introducing the parameter $\Xi_{0}$. In order to keep under control the number of new parameters, which is necessary to ensure the convergence of the MCMC chains, we only take $(\Xi_{0},w_{0})$ as the parameters that describe the DE sector of the theory, fixing $w_{a}=0$ and $n=2.5$; the latter value is of the order of that suggested by the RR and RT nonlocal models, see section 3.2. However, the precise value of $n$ is not very important for the forecasts that we present, since the uncertainty on this quantity is large for $|\Xi_{0}-1|\ll 1$ [16]. Again, we assume $\Lambda$CDM as the fiducial model used to generate the catalog of sources, so the fiducial values for the parameters $\Xi_{0}$ and $w_{0}$ are $\Xi_{0}=1$ and $w_{0}=-1$.

Fig. 17 shows the two-dimensional likelihood in $(\Xi_{0},w_{0})$ plane, for the realistic scenario for the error on the luminosity distance, and the three seed scenarios, while Fig. 18 shows the results obtained using the optimistic scenario for the error on the luminosity distance. We plot the limit from CMB+BAO+SNe and the combined limit from CMB+BAO+SNe+LISA standard sirens while, as before, using only standard sirens, with the addition of these new parameters with respect to $\Lambda$CDM, the MCMC chains fail to converge because of the limited number of sources. Note that CMB, BAO and SNe, as any other electromagnetic probe, are blind to $\Xi_{0}$, and therefore the corresponding contour from CMB+BAO+SNe are flat in the $\Xi_{0}$ direction. Standard sirens, however, lift this flat direction. The errors on $\Xi_{0}$ and $w_{0}$ from the corresponding one-dimensional likelihoods are shown in Tables 2.

Exactly as in the case of $w$CDM discussed above, we see that the accuracy that LISA, combined with CMB+BAO+SNe, can reach on $w_{0}$ is only of about $4.4\%$ (at least using only MBH binaries with counterpart as standard sirens), which is basically entirely determined by the current CMB+BAO+SNe observations.

In contrast, even under the current set of assumptions (see the discussion in sect. 5.1), that have led to relatively large errors on $H_{0}$ and on $w_{0}$ from LISA standard sirens, $\Xi_{0}$ still turns out to be an extremely interesting observable for LISA. First of all it can be observed only with GW experiments, and second, as we already anticipated in the discussion in section 2.2.4, $\Xi_{0}$ can be measured more accurately than $w_{0}$. This is confirmed by the results of our MCMC, which shows that, in the best case, $\Xi_{0}$ can be measured to $1\%$ accuracy, and even in the worst case we still have a $4.4\%$ accuracy, see Table 2. By comparison, for instance, the RT nonlocal model discussed in section 3.2 predicts a deviation from the GR value $\Xi_{0}=1$ at the level of $6.6\%$. Similar values can be obtained for the scalar-tensor theories discussed in Section 3.1: the $f(R)$ and the coupled Galileon models of Section 3.1.1 and the DHOST system of Section 3.1.2 can reach deviations of the $\Xi_{0}$ parameter of order respectively $5\%$ and $2\%$, see Figures, 2, 3 and 4. 171717In sections 2.2.4, 3.1.2 we also discussed more complex parametrizations for the quantity $\Xi(z)$, designed to accommodate the features of specific models: it should be straightforward to apply the methods used so far also to those cases, but we leave this for future works. Thus, the accuracy that LISA can reach on modified GW propagation is extremely interesting for testing modified gravity and dark energy.

It is interesting to compare the results that we have obtained for LISA with the expectations for a third-generation ground based interferometer such as the Einstein Telescope (ET). The Einstein Telescope is expected to detect of order $10^{5}-10^{6}$ binary neutron stars per year, up to redshift $z\sim 2$. For standard sirens with electromagnetic counterpart, the main uncertainty in the analysis comes from the estimate of the fraction of coalescences that will have an observed electromagnetic counterpart. This depends not only on details of the electromagnetic emission, that can in principle be modeled, but, given this very high event rate, also on issues that are currently difficult to predict, such as what will be the network of GRB detectors and of ground-based telescopes that will be available at the time of ET, and what fraction of telescope time will be devoted to the electromagnetic follow up of gravitational events. A rather common working hypothesis is that, out of these $10^{5}-10^{6}$ events per years, in a few years of running one could collect of order $10^{3}$ events with electromagnetic counterpart. Under this hypothesis, it was found in [16] that, by combining ET with the same CMB+BAO+SNe dataset that we are using in this paper, and extending the dark energy sector through the parameters $\Xi_{0}$ and $w_{0}$ (and fixing $n=2.5$, as in this paper), at ET one can get the accuracy $\Delta\Xi_{0}=0.008$ and $\Delta w_{0}=0.032$. Comparing with the results for LISA in Table 2 we see that the accuracy at ET would be slightly better than that in the ‘hnd’ scenario for LISA, although quite comparable. However, this comparison relies heavily on the assumption of having $10^{3}$ sources with counterpart, which at the present stage is only a reasonable working hypothesis. 181818A more detailed study of the accuracy on $(\Xi_{0},w_{0})$ obtained with a realistic modeling of the joint GW-GRB observations between ET and the proposed THESEUS mission will be presented in [209]. One should also observe that ET and LISA observations are complementary because most of the ET events are at redshift z\,\raise 1.72218pt\hbox{<}\kern-8.00003pt\lower 2.6694pt\hbox{\sim}\,1, while, as we see from Fig. 12, LISA has events distributed up to $z\sim 6$. This will allow us to study the ratio $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ in different redshift ranges. In terms of the $(\Xi_{0},n)$ parametrization, this means that LISA will be able to determine the asymptotic value $\Xi_{0}$, without much contamination from the value of $n$ (and, in turn, it will not be able to determine $n$ very accurately), while the lower redshift events of ET will rather be sensitive to a combination of $\Xi_{0}$ and $n$, so the joint LISA and ET data would lift this degeneracy.

6 Conclusions

In this paper we have examined the potential of LISA MBH binaries (with electromagnetic counterparts) for constraining cosmology, and in particular the dark energy sector associated with modified gravity theories. As in [35], we have generated a population of MBHBs using semi-analytic models for the formation and evolution of galaxies, and we have examined three different scenarios for massive black hole formation at high redshifts, namely pop III seeds, and heavy seeds with and without a delay between galaxy merger and the merger of the central MBHs.

Compared to earlier work specifically related to LISA, in particular ref. [35], we have used updated estimates for the LISA configuration and sensitivity curve, we examined a scenario with a more realistic error for the redshift determination of the source and delensing, we produced more realistic mock catalogs by including the scattering of the mock data according to the observational error, and we took fully into account the degeneracies between cosmological parameters by performing a series of MCMCs, both with standard sirens only and with standard sirens combined with current CMB+BAO+SNe data. On the other hand, we have emphasized (see the discussion at the end of sect. 5.1) that many more issues must be addressed before giving a final word of the sensitivity of LISA to cosmological parameters, and this paper should be considered as a contribution toward that goal, under the given assumptions that we have discussed in detail.

One result of this analysis is that, under the assumptions that we have used, the estimate for the accuracy that can be obtained for the Hubble parameter using LISA MBH is revised toward higher values. In particular, in the most favorable formation scenario, which turns out to be the heavy-seed no-delay one, for the relative error on $H_{0}$ from standard sirens we find $3.8\%$ assuming the “optimistic” scenario for the accuracy of redshift measurement and delensing, and $7.7\%$ with possibly more “realistic” assumptions on redshift determination and delensing. We have presented similar results for the DE equation of state, where we found that, under the assumptions we made, the inclusion of MBH LISA standard sirens does not improve significantly the determination of $w_{0}$ compared to current CMB+BAO+SNe observations. Apart from the role of the assumptions that we have used, it is also important to stress that LISA will also see other potential GW standard sirens, in particular stellar mass BH binaries, and extreme mass ratio inspirals (EMRIs). These sources are not expected to have an electromagnetic counterpart, but can still be used as standard sirens by using the statistical method. Work on this is currently in progress within our LISA Cosmology Working Group, and could lead to a significantly better error on $H_{0}$ and $w_{0}$.

The most interesting results concerning the potential of LISA for cosmology studies come however from modified GW propagation, which represents the core of this paper. The fact that a modified friction term in the propagation equation of tensor perturbations over FRW can give rise to a modification of the luminosity distance of standard sirens has been recognized in recent years through the study of explicit models such as scalar-tensor theories and nonlocal gravity [13, 6, 14, 15, 16] (although the possibility of a modification of the luminosity distance due to gravitons leaking in the bulk was already observed in the context of the DGP model [12]). In this paper we analyzed the explicit predictions of some of the best-motivated and most studied modifications of gravity, such as scalar-tensor systems (Horndeski and the more general DHOST family of scalar-tensor theories), nonlocal gravity, bigravity, and theories with extra and varying dimensions. We have found that modified GW propagation is an absolutely generic phenomenon, that takes place in all the theories that we have considered (in spite of the fact that all these theories comply with the limit from GW170817 on the deviation of the speed of GWs from the speed of light).

We have seen that, in most cases, a simple and accurate description of the effect is captured by the $(\Xi_{0},n)$ parametrization of Eq. (2.31) originally introduced in [16]. This parametrization takes into account the fact that the ratio of gravitational to electromagnetic luminosity distances must go to one as the redshift $z\rightarrow 0$, and in most model it saturates to a constant, $\Xi_{0}$, at large redshift (because in typical models dark energy is a recent phenomenon on cosmological scales), and smoothly interpolates between these asymptotic values with a power-law determined by the index $n$. Of course, deviations from this simple behavior can occur in some cases because of specific physical reasons. Most notably, in bigravity we have found a series of oscillations due to the “beatings” between the two metrics. It follows that in some specific cases alternative analytic formulas can be useful, and we have presented some alternative parametrizations.

We have run a series of MCMCs to determine the accuracy that LISA can reach on the parameter $\Xi_{0}$ in the $(\Xi_{0},n)$ parametrization of eq. (2.31).191919The parameter $n$ is less important if the sources are at large redshifts, where the expression for the ratio of gravitational and electromagnetic luminosity distances saturates to its asymptotic value $\Xi_{0}$. In any case, the accuracy that can be reached on it can be estimated analytically as $\Delta n=[\Delta\Xi_{0}/|1-\Xi_{0}|]\,\times[(1+z)^{n}-1]/[\log(1+z)]$ where $z$ is the typical redshift of the sources [16]. In the low-redshift limit this reduces to $\Delta n/n=\Delta\Xi_{0}/|1-\Xi_{0}|$. There are two reasons that make this observable especially important for advanced GW detectors. First of all, modified GW propagation is an observable accessible only to GW observations, and to which electromagnetic observations are blind.202020More precisely, we have seen that in some models (but not in general), modified GW propagation is related to a time dependence of the effective Newton constant, and in these models one could access $\Xi_{0}$ indirectly through the effect of the modified Newton’s constant on structure formation. Second, as discussed in section 2.2.4, in a generic modified gravity theory, in which the deviation of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ from 1 is of the same order as the deviation of $w_{\rm DE}(z)$ from $-1$, the effect of $d_{L}^{\,\rm gw}(z)/d_{L}^{\,\rm em}(z)$ on standard sirens dominates over the effect of the DE equation of state, because the latter is partially compensated by degeneracies with other (fitted) cosmological parameters. As a consequence, the accuracy expected on $\Xi_{0}$ is better than that on $w_{0}$. This argument has been confirmed by our explicit MCMC computations. Combining LISA with CMB+BAO+SNe to reduce the degeneracies with the other parameters, in the best case (heavy-seed no-delay formation scenario and “optimistic” scenario for the accuracy of redshift measurement and delensing) we have found that $\Xi_{0}$ can be measured to an accuracy that reaches $1.1\%$ (to be compared with $4.5\%$ for $w_{0}$) and even in the worst scenario still is $4.4\%$ (see Table 2). Last but not least, in several instances the explicit models that we have considered give predictions for $\Xi_{0}$ larger or equal than these values. For instance the RT model of Section 3.2 predicts for $\Xi_{0}$ a deviation from the GR value of order $6.6\%$. Similar values can be obtained for the scalar-tensor theories: the $f(R)$ and coupled Galileon models of Section 3.1.1 and the DHOST system of Section 3.1.2 can reach deviations of the $\Xi_{0}$ parameter of order $3\%$, $5\%$ and $2\%$ respectively. Due to their complex dynamics, the bigravity set-up studied in section 3.3 and the dimensionally changing systems of section 3.4 cannot be faithfully described in terms of the $(\Xi_{0},n)$ parameterization of Eq. (2.31). We nevertheless quantitatively derived projected constrains on their parameter spaces from standard siren catalogs.

At the theoretical level, we have examined the predictions on modified GW propagation of a large number of the best studied modified gravity theories. In particular, we demonstrated the improved capability of LISA to probe GW oscillations, an analog to neutrino flavour oscillations, present in models with extra tensor interactions such as bigravity, potentially constraining the parameter space of specific theories. We compute predictions for bigravity using a high frequency expansion, focusing on the high mass regime $m_{g}\gg H_{0}$. In this region of parameter space GW oscillations occur in the mHz range, but the theory is not generically able to account for cosmic acceleration. Standard sirens at cosmological distances, as will be observed by LISA, have the potential to constrain the mass range $m_{g}\gtrsim 10^{-25}$eV for most mixing angles, which would improve by three orders of magnitudes the result obtained from the current LIGO-Virgo detection of GW170817. This is a conservative estimate: including frequency-dependent effects on the waveform will improve these bounds. Notice also that the improvement in constraining power in GW oscillations (e.g. $m_{g}$) over LIGO/Virgo will not be matched by standard-siren analyses with third-generation ground detectors, as the extended reach of LISA is driven mainly by the lower frequency range (which increases the sensitivity since the signal is suppressed by $(m_{g}/f)^{2}$), in addition to a longer oscillation baseline (higher source redshift).

The broad conclusion is that modified GW propagation is a prime observable for cosmological studies with advanced GW detectors, and MBH binaries detectable with LISA can be a powerful probe of modified gravity and dark energy. Further significant improvements are expected by the study of LISA standard sirens without electromagnetic counterpart, such as stellar mass BH binaries and EMRIs, that can be used as standard sirens through the statistical method. We leave this for future work.

Acknowledgments

It is a pleasure to thank Chiara Caprini for useful discussions, the LISA Cosmology Working Group Chairs – Robert Caldwell, Chiara Caprini, Germano Nardini – for their input and support during the project, and Vuk Mandic and Alberto Sesana for carefully reading our manuscript and for their useful suggestions. The work of E.B, S.F. and M.M. is supported by the Fonds National Suisse and by the SwissMap National Center for Competence in Research. The MCMCs have been run on the Baobab cluster of the University of Geneva. This project has received funding (to E. Barausse) from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme (grant agreement no. GRAMS-815673; project title “GRavity from Astrophysical to Microscopic Scales”). This work was supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. N.B., D.B., A.G and S. M. acknowledge financial support by ASI Grant 2016-24-H.0. G.C. is supported by the I+D grant FIS2017-86497-C2-2-P of the Spanish Ministry of Science, Innovation and Universities. M.C. is supported by the Labex P2IO and the Enhanced Eurotalents Fellowship. C.D. and L.L. were supported by a Swiss National Science Foundation Professorship grant (No. 170547). J.M.E. is supported by the Spanish FPU Grant No. FPU14/01618, the Research Project FPA2015-68048-03-3P (MINECO-FEDER) and the Centro de Excelencia Severo Ochoa Program SEV- 2016-0597. M.R.F. acknowledges support from STFC grant ST/N000668/1. J.G.B. acknowledges support from the Research Project FPA2015-68048-03-3P [MINECO-FEDER] and the Centro de Excelencia Severo Ochoa Program SEV-2016-0597. A.K. was supported by the Centre National d’Études Spatiales: he also acknowledges support from the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. M.S. is supported in part by STFC grant ST/P000258/1. G.T. is partially supported by STFC grant ST/P00055X/1. M.Z. is supported by the Marie Sklodowska-Curie Global Fellowship Project NLO-CO.

Appendix A GW luminosity distance and the flux-luminosity relation

In this Appendix we show that the GW luminosity distance in modified gravity theories still obeys the standard relation (2.1) between the energy flux ${\cal F}$ measured in the observer frame and the intrinsic luminosity ${\cal L}$ measured at the source frame. Using purely kinematic consideration discussed for example in [20], Section 4.1.4, we derive an expression for the gravity wave luminosity distance in a case where the graviton is massless, $m_{T}=0$, and the non-minimal coupling with other fields are such that $\Pi_{ij}=0$. GW are then produced by events at high redshift associated with standard sirens, and then can propagate freely through space-time. (We examine the effect of a graviton mass and couplings with additional spin-2 fields in Section 3.3.)

In this case, the evolution equation (1.2) can be obtained by varying the effective quadratic action for tensor modes

[TABLE]

Here $\bar{M}$ is a reference mass scale, whose value we will discuss later. Studying the GW dynamics governed by this action is equivalent to studying the propagation of free massless modes on a homogeneous FRW geometry characterised by the line element

[TABLE]

This expression for the homogeneous background metric leads to the following formula for the comoving distance between a source and the observer, computed in terms of the tensor null geodesics $ds\,=\,0$ (we use the relation $d\eta\,=\,dt/a(t)$ connecting physical and conformal time):

[TABLE]

with $t$ an arbitrary time. We denote by $t_{s}$ the time of emission of GW at the source, while $t_{0}$ the observation time, which we assume being today. The comoving distance depends on the speed of the GW. The physical distance results:

[TABLE]

We now choose our reference normalization scale $\bar{M}$ so that in the limit $t\to t_{s}$ the ratio between physical and comoving distance acquires the standard expression

[TABLE]

at the time of the emission from the source. This might be motivated by requiring that nearby the source the effects of modified gravity had not yet have time to develop. Condition (A.5) leads to the definition

[TABLE]

that we will use in what follows.

Suppose now that the observer measures GW signals corresponding to wavecrests emitted at different times from the source. (They travel through the same comoving distance.) Using expression (A.3) for $r_{\rm com}(t_{0})$, at linear order in $\Delta t_{s}$ we can write the relation

[TABLE]

between the time difference of two GW wavecrests as measured at emission and observation times. Formula (A.7), which we more conveniently express in terms of redshift, states that source and observer clocks tick with different rates. We assume that the observer makes its measure today at redshift equal to zero, while the emission occurs at resdshift $z$. Then, from eq. (A.7), we can write

[TABLE]

implying that frequencies measured in the source and observer frames are related by

[TABLE]

We now define the luminosity distance, following [20]. We call ${\cal F}$ the GW energy flux measured by the observer, corresponding to the amount of GW energy per unit time per unit area. ${\cal L}$ is the luminosity of the source, defined as the power it radiates

[TABLE]

Then we define the luminosity distance $d_{L}^{\rm gw}$ as

[TABLE]

Since we measure the energy flux at the observer position, we need to convert $dE_{s}/dt_{s}$ into the observer frame. The energy scales as the frequency, eq. (A.9), while $dt_{s}$ and $dt_{0}$ are related by eq. (A.8). Hence

[TABLE]

The area crossed by the flux of GW which propagate radially from the source is [using eq. (A.6)]

[TABLE]

We now collect the various pieces of information, and obtain the following expression for the luminosity distance [$a(0)=a(z=0)$ is the value of scale factor today]

[TABLE]

Using the fact that

[TABLE]

and eq. (A.3), we can write the following relation for the comoving distance for the GWs:

[TABLE]

The luminosity distance for GW defined from eq. (A.11) can then be written as

[TABLE]

and therefore agrees with eq. (2.27).

Appendix B Technical details on bigravity

B.1 Hassan–Rosen theory of bigravity

The theory known as bigravity [134] is the only known theory of two interactive spin-2 fields that is free of ghosts at the fully non-linear level. See [136, 137, 138] for reviews. It is described by the action

[TABLE]

with $g_{\mu\nu}$ and $\tilde{g}_{\mu\nu}$ the two metric tensors, $M_{\rm Pl}^{2}$ and $\kappa M_{\rm Pl}^{2}$ the corresponding squares of Planck masses and $m$ the graviton mass. Matter is coupled only to the first metric. The interaction potential between the two metrics is indicated by $V$, and it takes the form

[TABLE]

with $a_{n}$ dimensionless parameters, and

[TABLE]

where $\tau_{i}\,=\,{\rm tr}\left[Y^{i}\right]$ with $Y_{\mu}^{\,\,\nu}\,=\,[\sqrt{X}]_{\mu}^{\,\,\nu}$.

This theory admits homogeneous FRW configurations described by two independent metrics

[TABLE]

The ratio of scale factors is denoted with

[TABLE]

and from now on $\mathcal{H}=a^{\prime}/a$ is the Hubble parameter corresponding to the first metric. The Friedmann equation for the first metric reads

[TABLE]

The theory admit two branches of solutions, but only one describes physically interesting cosmological configurations [142]. In this branch the Bianchi identities are realised in the form

[TABLE]

and together with Friedmann equations lead to an algebraic equation for $\xi$:

[TABLE]

where the coefficients $a_{i}$ are the parameters of the bigravity potential (B.2). This information about homogeneous configurations is sufficient for the scope of this work, more details can be found in [139, 140, 141, 142].

B.2 Details on the WKB approximation for bigravity

In this Appendix we spell out details on the calculations using a WKB approximation leading to the solution (3.40–3.41) for the system of equations (3.37) we use in the main text. In order to obtain WKB solutions, it is convenient to work in matrix notation. After absorbing the cosmic friction term and defining a vector containing the two tensor modes,

[TABLE]

we can express the evolution equation (3.37) for $h_{1}$ and $h_{2}$ in bigravity as

[TABLE]

where we have defined

[TABLE]

Next, we introduce a dimensionless expansion parameter $\epsilon$, rescaling time as $d\tau\rightarrow d\tau/\epsilon$, and controlling different orders in a high frequency WKB approximation. We make the following Ansatz for the solution $\Phi$ for tensor modes,

[TABLE]

and aim to solve order by order in the expansion parameter $\epsilon$. In eq. (B.17) we expand over a basis of eigenvectors controlled by the matrix $\hat{E}$, associated with the matrix of eigenfrequencies $\hat{\theta}$, which appears at the exponent of (B.17) and is diagonal by definition. At leading order $\epsilon^{0}$ of our expansion parameter, we obtain the eigenfrequencies

[TABLE]

and the matrix of eigenvectors

[TABLE]

At next to leading order, $\epsilon^{1}$, the amplitude $\Phi_{0}$ can be obtained solving

[TABLE]

where we defined for convenience the matrix $\hat{G}\equiv e^{i\int\hat{\theta}d\tau}$. For general time dependent coefficients, this matrix equation cannot be solved analytically (because the matrices in the parenthesis should commute at any given time). However, within the regime of the WKB, a matrix exponential solution is a very good approximation (that we have checked numerically)

[TABLE]

where $\vec{C}_{0}$ is a vector of constant coefficients to be fixed with the initial conditions. Recall that $\hat{\theta}$ is a diagonal matrix and thus the term $\hat{\theta}^{-1/2}$ in (B.21) is the usual WKB scaling. If there is time dependence, there can be corrections to this scaling, which corresponds to the matrix exponential. The fact that the matrix exponential works as a solution is because the matrix in the exponent is small in this regime and, as a consequence, corrections arising from commutators of this matrix are further suppressed. For higher order corrections in the WKB expansion, one can proceed iteratively and solve $\Phi_{1}$ at order $\epsilon^{2}$ using the solution of $\Phi_{0}$.

Having an approximate analytical solution allows us to understand the role of each parameter. In particular, we note that the speed of the massive mode has a key role in the mixing. This is better seen in the high-$k$ limit where the phases tend to

[TABLE]

Focusing on the non-diagonal terms of the matrix of eigenvectors,

[TABLE]

we notice that in the high-$k$ limit the mode mixing is suppressed, i.e. $\hat{E}_{12}\ll 1$ and $\hat{E}_{21}\ll 1$, if each tensor propagate at a different speed, $c\neq 1$, and there is a large hierarchy between the mass term $m$ and the wavevector $k$.

While in the main text we focussed on a large mass regime, we conclude this Appendix considering a small mass regime, with a mass parameter of the order of the Hubble constant $m\sim\mathcal{H}_{0}$, and show that in this case there is no distinctive observational effect of bigravity in the GW propagation. In this regime, the energy density proportional to $m^{2}$ in the right hand side of the Friedmann equation (B.11) is of the same order of magnitude of the observed vacuum energy. Then, the viable branch of solutions satisfies $\xi\ll 1$, which gives the following value for the speed of the second tensor modes [141]

[TABLE]

This implies that during matter domination $c^{2}\sim 16$, and during radiation domination $c^{2}\sim 25$. Therefore, given that LISA frequency of GWs is much larger than the rate of expansion, $k_{\text{LISA}}\sim 10^{16}\mathcal{H}_{0}$, inequality (3.42) is not satisfied and mixing among different modes is negligible in the small mass regime.

Moreover, this mass range is also still far from being constrained through the modified dispersion relation with LISA. In particular, the bound on an effective mass term $m\,f_{1}^{1/2}$ in this case is [210]

[TABLE]

where we have introduced the expected distance, frequency and signal-to-noise ratio ($\text{SNR}\sim 1/(f\Delta t)$) of a massive BH binary in the LISA band. This is seven orders of magnitude larger than $\mathcal{H}_{0}$, thus far from the small mass regime.

Affiliations

a Département de Physique Théorique and Center for Astroparticle Physics,

Université de Genève, 24 quai Ansermet, CH-1211 Genève 4, Switzerland

b Instituto de Estructura de la Materia, CSIC, Serrano 121, 28006 Madrid, Spain

c Institut de physique théorique, Université Paris Saclay CEA, CNRS, 91191 Gif-sur-Yvette, France

d AIM, CEA, CNRS, Univ. Paris-Saclay, Univ. Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

e Laboratoire de Physique Théorique, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91405 Orsay, France

f Center for Theoretical Astrophysics and Cosmology,

Institute for Computational Science, University of Zürich, CH-8057 Zürich, Switzerland

g Instituto de Física Teórica UAM-CSIC, Universidad Autonóma de Madrid, Cantoblanco, 28049 Madrid, Spain

h Institute of Cosmology and Gravitation, University of Portsmouth,

Dennis Sciama Building, Portsmouth, PO1 3FX, United Kingdom

i Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova,

and INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy

j Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am Mühlenberg 1, 14476 Potsdam-Golm, Germany

k Department of Physics, Swansea University, Swansea, SA2 8PP, UK

l Berkeley Center for Cosmological Physics, LBNL and University of California at Berkeley, Berkeley, California 94720, USA

m Institut d’Astrophysique de Paris, CNRS and Sorbonne Universités, UMR 7095, 98 bis bd Arago, 75014 Paris, France

n Scuola Internazionale di Studi Superiori Avanzati (SISSA), via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste

o INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122 Padova, Italy

p Gran Sasso Science Institute, Viale F. Crispi 7, I-67100 L’Aquila, Italy

q Theoretical Particle Physics and Cosmology Group, Physics Department, King’s College London, University of London, Strand, London WC2R 2LS, United Kingdom

r Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom