Probing observational bounds on scalar-tensor theories from standard sirens

TL;DR

This paper demonstrates how future gravitational wave observations from the Einstein telescope can tightly constrain scalar-tensor theories and $f(R)$ gravity models, improving our understanding of cosmic acceleration.

Contribution

It introduces a simulated catalog of standard sirens for the Einstein telescope and uses it to set new bounds on scalar-tensor theories and $f(R)$ gravity, highlighting the potential of future GW data.

Findings

New bounds on the running of the Planck mass.

Potential to test $f(R)$ gravity deviations up to 95% confidence.

Future GW observations can significantly constrain modified gravity theories.

Abstract

Standard sirens are the gravitational wave (GW) analog of the astronomical standard candles, and can provide powerful information about the dynamics of the Universe. In this work, we simulate a catalog with 1000 standard siren events from binary neutron star mergers, within the sensitivity predicted for the third generation of the ground GW detector called Einstein telescope. After correctly modifying the propagation of GWs as input to generate the catalog, we apply our mock data set on scalar-tensor theories where the speed of GW propagation is equal to the speed of light. As a first application, we find new observational bounds on the running of the Planck mass, when considering appropriate values within the stability condition of the theory, and we discuss some consequences on the amplitude of the running of the Planck mass. In the second part, we combine our simulated standard sirens catalog with other geometric cosmological tests (Supernovae Ia and cosmic chronometers measurements) to constrain the Hu-Sawicki $f(R)$ gravity model. We thus find new and non-null deviations from the standard $\Lambda$CDM model, showing that in the future the $f(R)$ gravity can be tested up to 95\% confidence level. The results obtained here show that the statistical accuracy achievable by future ground based GW observations, mainly with the ET detector (and planed detectors with a similar sensitivity), can provide strong observational bounds on modified gravity theories.

Tables

Table 1. Table 1: Summary of the MCMC results for the cases α M ​ 0 > 0 subscript 𝛼 𝑀 0 0 \alpha_{M0}>0 and α M ​ 0 < 0 subscript 𝛼 𝑀 0 0 \alpha_{M0}<0 . The upper and lower values next to the mean value of each parameter denote the 68% and 95% CL errors.

Stability Conditions	$H_0$	${\Omega }_{m0}$	${\alpha }_{M0}$
	${67.466}_{-0.067\left(0.179\right)}^{+0.036\left(0.143\right)}$	${0.328}_{-0.014\left(0.028\right)}^{+0.015\left(0.028\right)}$	$-{0.100}_{-0.043\left(0.085\right)}^{+0.051\left(0.092\right)}$
${\alpha }_{M0}>0$	${67.390}_{-0.050\left(0.095\right)}^{+0.047\left(0.098\right)}$	${0.297}_{-0.044\left(0.072\right)}^{+0.029\left(0.083\right)}$	${0.199}_{-0.097\left(0.167\right)}^{+0.069\left(0.178\right)}$

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\sum _{i=2}^5{M}_P^2{\mathcal{L}}_i+{\mathcal{L}}_m\right],$$

$${\mathcal{L}}_2=$$

$${\mathcal{L}}_3=$$

$${\mathcal{L}}_4=$$

$${\mathcal{L}}_5=$$

$$+2{\varphi }_{;\mu \nu }{\varphi }^{;\mu \sigma }{\varphi }_{;\sigma }^{;\nu }-3{\varphi }_{;\mu \nu }{\varphi }^{;\mu \nu }\square \varphi ].$$

$$h_{ij}^{\prime \prime }+(2+{\alpha }_M)\mathcal{H}{h}_{ij}^{\prime }+k^2{h}_{ij}=0,$$

$${\alpha }_M=\frac{1}{H{M}_{\ast }^2}\frac{d{M}_{\ast }^2}{dt},$$

$$M_{\ast }^2=2(G_4-2X{G}_{4X}+X{G}_{5\varphi }-\dot{\varphi }HX{G}_{5X}),$$

$$h=e^{-\mathcal{D}}{h}_{GR},$$

$$\mathcal{D}=\frac{1}{2}{\int }^{\tau }{\alpha }_M\mathcal{H}d{\tau }^{\prime }.$$

$$d_L^{GW}(z)=d_L^{EM}(z)\exp \left[\frac{1}{2}{\int }_0^z\frac{d{z}^{\prime }}{1+z^{\prime }}{\alpha }_M(z^{\prime })\right],$$

$$d_L^{EM}(z)=(1+z){\int }_0^z\frac{d{z}^{\prime }}{H(z^{\prime })}.$$

$$\tilde{h}(f)=Q\mathcal{A}{f}^{-7/6}{e}^{i\Phi (f)},$$

$$\mathcal{A}=\sqrt{\frac{5}{96}}\frac{{\mathcal{M}}_c^{5/6}}{{\pi }^{2/3}{d}_L^{GW}}\left(\sum _{i=0}^6{A}_i(\pi f{)}^{i/3}\right),$$

$$Q^2=F_{+}^2(1+co{s}^2(\iota ){)}^2+2F_{\times }^{2}co{s}^2(\iota ),$$

$$\Phi (f)=2\pi f{t}_c-{\varphi }_c-\frac{\pi }{4}+\frac{3}{128\eta v^5}\left[1+\sum _{i=2}^7{\alpha }_i{v}^i\right],$$

$$\Delta {\theta }^i=\sqrt{{\Sigma }^{ii}},$$

$${\Gamma }_{ij}=\left(\frac{\partial \tilde{h}}{\partial {\theta }^i}\mid \frac{\partial \tilde{h}}{\partial {\theta }^j}\right).$$

$$({\tilde{h}}_1\mid {\tilde{h}}_2)\equiv 2{\int }_{f_{low}}^{f_{upper}}\frac{{\tilde{h}}_1{\tilde{h}}_2^{\ast }+{\tilde{h}}_1^{\ast }{\tilde{h}}_2}{S_n(f)}df,$$

$${\mathrm{SNR}}^2\equiv 4Re{\int }_{f_{low}}^{f_{upper}}\frac{|h(f){|}^2}{S_n}df.$$

$${\sigma }_{d_L}^2={\left(\frac{\partial \tilde{h}(f)}{\partial d_L},\frac{\partial \tilde{h}(f)}{\partial d_L}\right)}^{-1}.$$

$${\sigma }_{d_L}$$

$$=\sqrt{{\left(\frac{2d_L(z)}{\mathrm{SNR}}\right)}^2+(0.05z{d}_L(z){)}^2}.$$

$$P(z)\propto \frac{4\pi d_C^2(z)r(z)}{H(z)(1+z)},$$

$$r(z)=\{\begin{array}{rl}& 1+2z,z\le 1\\ & (15-3z)/4,1<z<5\\ & 0,z\ge 5\end{array}$$

$${\mathcal{L}}_{\text{GW}}\propto \exp \left[-\frac{1}{2}\sum _{i=1}^{1000}\left(\frac{d_L^{obs}(z_i)-d_L^{th}(z_i)}{{\sigma }_{d_L,i}}\right)\right].$$

$${\alpha }_M=\frac{\dot{G_4}}{H{G}_4}.$$

$$S=\int d^{4}x\sqrt{-g}\frac{M_P^2}{2}f(R)+S_m,$$

$$3F{H}^2=8\pi G\left({\rho }_m+{\rho }_r\right)+\frac{1}{2}\left(FR-f\right)-3H\dot{F},$$

$$\displaystyle-2F\dot{H}=8\pi G\left(\rho_{m}+P_{m}+\rho_{r}+P_{r}\right)+\ddot{F}-H\dot{F}\$$

$$R=6\left(2H^2+\dot{H}\right).$$

$$f(R)=R-2\Lambda y(R,b),$$

Full text

Modified Gravity, Gravitational Waves Standard Sirens

I Introduction

After 20 years of research, the nature of the physical mechanism responsible for accelerating the Universe at late times is still an open question, and a large variety of cosmological models have been and are continually proposed in the literature to explain such observations (see DE_review ; MG_review01 ; MG_review03 for review). This is essentially due to the difficulty of discriminating among different scenarios that respond to the observations in the same way, leading to a theoretical degeneracy.

The observation of new astrophysical sources, through a direct manifestation of gravitational effects, can provide rich physical information about the nature of gravity, which should play a key role to probe new (or rule out) additional gravitational degree(s) of freedom, or exotic forms of energy such as dark energy. The gravitational waves (GWs) issued by binary systems, such as binary black hole (BBH) and/or binary neutron star (BNS) detected by LIGO/VIRGO, certainly open a new window to investigate fundamental physics in this direction. At present, catalogues of GWs from 10 BBH mergers and 1 BNS merger are available LIGO01 . The latter, the GW170817 event GW170817 , observed at $z\simeq 0.009$, has imposed strong constraints on modified gravity/dark energy models GW_MG01 ; GW_MG02 ; GW_MG03 ; GW_MG04 ; GW_MG041 ; GW_MG05 . Also, GW170817 was the first standard siren (the gravitational-wave analog of an astronomical standard candle) event to be catalogued, once its electromagnetic counterpart (GRB170817) was measured. These observations were also used to measure $H_{0}$ at the 12% accuracy, assuming a fiducial $\Lambda$CDM cosmology SS01 . We refer the reader to H0_GW1 ; H0_GW2 ; H0_GW3 for proposals to use standard siren to measure $H_{0}$ with more accuracy.

Given the central importance of GW astronomy, beyond the present performance of the LIGO and Virgo interferometers, plans for construction of several GW observatory interferometers (on earth and in space) are currently in preparation, such like LIGO Voyager CE , Cosmic Explore CE , Einstein Telescope (ET) ET_design01 ; ET_design02 , LISA LISA , DECIGO DECIGO , TianQin TianQin , among others, to observe GWs in the most diverse frequencies bands and different types of GW sources. In this paper, we are particularly interested to use the sensitivity predicted for the ET ET_design01 ; ET_design02 , which is a third-generation ground-based detector and it is envisaged to be several times more sensitive in amplitude than the advanced ground-based detectors in operation, covering the frequency band range from 1 Hz to $10^{4}$ Hz. Also, the ET is expected to have signal-to-noise ratio (SNR) for BBH and BNS mergers several times larger than the current measures, as well as to observe hundreds, thousands of these events throughout the whole operational time. Several works have been done using the ET sensitivity to simulate GWs standard siren in order to investigate diverse aspects in cosmology ET01 ; ET02 ; ET03 ; ET04 ; ET05 ; ET06 ; ET07 ; ET08 ; ET09 ; ET10 ; ET11 ; ET12 ; ET13 ; ET14 ; ET15 ; ET16 ; ET17 ; ET18 ; ET19 ; ET20 ; ET21 ; ET22 ; ET23 .

In this work, we generate a simulated catalog with 1000 standard siren events from BNS mergers, from the ET power spectral density noise, in order to evaluate forecasting observational constraints on scalar-tensor theories where the speed of GW propagation is equal to speed of light. First, assuming a well-known parametric model for the running of the Planck mass, and assuming appropriate stability conditions on the theory, we find new observational bounds on the amplitude of the running of the Planck mass and we discuss its possible implications. In the second part, we apply our simulated standard siren data on $f(R)$ gravity given by the Hu-Sawicki model in order to find new observational limits on such model. In both analysis, we find that the parameters that characterize deviations from General Relativity (GR) may be non-null, within some statistical borders.

The manuscript is organized as follows. In Section II, we set our theoretical framework to show how the GW propagation is modified from scalar-tensor theories. In Section III, we describe our methodology to generate standard siren mock catalogs. In Section IV and V, we present our main results. Finally, in Section VI, we outline our final considerations and future perspectives.

Throughout the text, we use units such that $c=\hbar=1$, and $M_{P}=1/\sqrt{8\pi G}$ is the Planck mass. Moreover, we adopt the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, $ds^{2}=-dt^{2}+a(t)^{2}\delta_{ij}dx^{i}dx^{j}$, where $a$ is the scale factor, normalized to unity today. As usual notation, we denote by a subscript ‘0’ physical quantities evaluated at the present time, and by the prime and dot symbols the derivatives with respect to the conformal time ($\tau$) and cosmic time ($t$), respectively, related by $dt=ad\tau$. Lastly, we express the Hubble constant $(H_{0})$ results in units of km/s/Mpc.

II Modified Gravitational Wave Propagation in Scalar-Tensor Gravity

The Horndeski theories of gravity Horndeski ; Deffayet (see Horndeski_revisao01 ; Horndeski_revisao02 for a review) are the most general Lorentz invariant scalar-tensor theories with second-order equations of motion and where all matter is universally coupled to gravity. They include, as a sub set, the archetypal modifications of gravity such as metric and Palatini $f(R)$ gravity, Brans-Dicke theories, galileons, among others. The Horndeski action reads

[TABLE]

where $g$ is the determinant of the metric tensor, and

[TABLE]

Here, $G_{i}$ ($i$ runs over 2, 3, 4, 5) are functions of a scalar field $\phi$ and the kinetic term $X\equiv-1/2\nabla^{\nu}\phi\nabla_{\nu}\phi$, and $G_{i,X}\equiv\partial G_{i}/\partial X$. For $G_{2}=\Lambda$, $G_{4}=M^{2}_{P}/2$ and $G_{3}=G_{5}=0$, we recover GR with a cosmological constant. For a general discussion on the model varieties for different $G_{i}$ choices after GW170817, see Horndeski_revisao02 .

Recently, the GW170817 event together with the electromagnetic counterpart showed that the speed of GW, $c_{T}$, is very close to the speed of light for $z<0.01$, i.e. $|c_{T}/c-1|<10^{-15}$ GW170817 . In the context of Horndeski gravity, in order to explain these constraints, the only option is to consider $G_{4,X}\approx 0$ and $G_{5}\approx constant$ in the action above. It is important to note that this restriction applies only to the local Universe ($\lesssim$ 40 Mpc). Thus, in principle, nothing prevents from considering the presence of these terms at redshifts larger than $z=0.01$. In fact, only future measurements at high $z$ can confirm whether $c_{T}=c$. Here, we assume $c_{T}=c$, without loss of generality in the analysis we are going to develop. Under this condition, the GW propagation obeys the equation of motion Saltas04

[TABLE]

where $h_{ij}$ is the metric tensor perturbation and $\mathcal{H}\equiv a^{\prime}/a$ is the Hubble rate in conformal time. Moreover, $\alpha_{M}$ is the running of the Planck mass, which enters as a friction term responsible for modifying the amplitude of GWs acting as a damping term:

[TABLE]

where $M_{*}$ is the effective Planck mass:

[TABLE]

and $H\equiv\dot{a}/a$ is the Hubble parameter. Following the methodology presented in Atsushi01 (see also Atsushi02 ; Atsushi03 ), we can write a generalized GW amplitude propagation for scalar-tensor theories as

[TABLE]

where

[TABLE]

Note that due the condition $c_{T}=c$, that is, $G_{4,X}\approx 0$ and $G_{5}\approx constant$, we do not have phase corrections in Eq. (10). As the GW amplitude is inversely proportional to the distance, one can interpret the amplitude modification in Eq. (10) as a correction to the luminosity distance, defining an effective luminosity distance, or equivalently, an effective amplitude correction as Atsushi02 ; MG_amplitude01 ; MG_amplitude02 :

[TABLE]

where $d_{L}^{EM}$ is the standard electromagnetic luminosity distance as a function of the redshift111The redshift is defined as $z=a^{-1}-1$.:

[TABLE]

This generalization has been recently investigated in some contexts of modified gravity (see, e.g., Atsushi01 ; Atsushi02 ; Atsushi03 ; MG_amplitude01 ; MG_amplitude02 ; MG_dL01 ; MG_dL02 ; MG_dL03 ; MG_dL04 ; MG_dL05 ; MG_dL06 ; MG_dL07 ; MG_dL08 ).

It is usual to choose phenomenologically motivated functional forms for the functions $\alpha_{i}$ (see, e.g., Bellini ; alphai_01 ; alphai_02 ; alphai_03 ). Typically, their evolution are tied to the scale factor $a(t)$ or to the dark energy density $\Omega_{de}(a)$ raised to some power $n$. On the other hand, an important point within Horndeski gravity are the stability conditions of the theory. Appropriate values of the free parameters functions must be considered in order to have a stable theory throughout the evolution of the Universe (see alphai_02 and reference therein). Following alphai_02 , we adopt the parametrization $\alpha_{M}=\alpha_{M0}a^{n}$, so that the stability conditions can be summarized as follows:

$n>\dfrac{5}{2}$: stable for $\alpha_{M0}<0$; 2. 2.

$0<n<1+\dfrac{3\Omega_{m0}}{2}$: stable for $\alpha_{M0}>0$.

Here, $\Omega_{m0}$ is the present dimensionless matter density. Under these considerations, we can note from Eq. (12) that the changes in the GW amplitude propagation will be sensitive to the sign of $\alpha_{M0}$. Possible corrections with $\alpha_{M0}>0$ or $\alpha_{M0}<0$ will induce $d_{L}^{GW}>d_{L}^{EM}$ and $d_{L}^{GW}<d_{L}^{EM}$, respectively. We quantify these effects in Figure 1. We note that variations on $\alpha_{M0}>0\,(<0)$ can produce changes up to 10% (5%), respectively, on the effective GW amplitude, for reasonable values of the running of the Planck mass today.

III Methodology and GW standard sirens data set

In order to move on, we need to define the GW signal $h_{GR}$. In modeling the gravitational waves form, given a GW strain signal $h(t)=A(t)\cos[\Phi(t)]$, we can obtain its Fourier transform $\tilde{h}(f)$ using the stationary phase approximation for the orbital phase of inspiraling binary system. For a coalescing binary system with component masses $m_{1}$ and $m_{2}$, we have

[TABLE]

where $\mathcal{A}$ is the GW inspiral amplitude computed perturbatively within the so-called post-Newtonian (PN) formalism up until 3 PN corrections,

[TABLE]

where $d_{L}^{GW}$ is the modified luminosity distance as in Eq. (12), and the coefficients $A_{i}$ are given in Appendix A. The function $Q$ is expressed by

[TABLE]

where $\iota$ is the inclination angle of the binary orbital angular momentum with respect to the line of sight, and $F^{2}_{+}$, $F^{2}_{\times}$ are the two antenna pattern functions. In Eq. (14), the function $\Phi(f)$ is the inspiral phase of the binary system:

[TABLE]

where the coefficients $\alpha_{i}$ are the corrections up to the 3.5 PN corrections. In Appendix A, we also provide the expressions for these coefficients. In the above equation, we have defined $v\equiv(\pi Mf)^{1/3}$, $M\equiv m_{1}+m_{2}$, $\eta\equiv m_{1}m_{2}/(m_{1}+m_{2})^{2}$, and $\mathcal{M}_{c}\equiv(1+z)M\eta^{3/5}$ to be the inspiral reduced frequency, total mass, symmetric mass ratio, and the redshifted chirp mass, respectively. The quantities $t_{c}$ and $\phi_{c}$ are the time and phase of coalescence, respectively.

After having defined the modified GW signal for compact binaries, in what follows we summarize the already know methodology used to estimate $d_{L}(z)$ measures from GW standard sirens. We refer to Schutz ; Holz for pioneer works in this regard.

For a high enough signal-to-noise ratio (SNR) and a given waveform model, $h(f,\theta_{i})$, with free parameters $\theta_{i}$, we can use the Fisher matrix analysis to provide upper bounds for the free parameters of the models by means of the Cramer-Rao bound Fisher01 ; Fisher02 . We refer the reader to Fisher03 ; Fisher04 ; Fisher05 ; Fisher06 ; Fisher07 ; Fisher08 for a discussion on the Fisher analysis to estimate parameters in binary systems for a given GW signal. Once the waveform model is defined, the root-mean-squared error on any parameter is determined by

[TABLE]

where $\Sigma^{ij}$ is the covariance matrix, i.e, the inverse of the Fisher matrix, $\Sigma^{ij}=\Gamma_{ij}^{-1}$. The Fisher matrix is given by

[TABLE]

The inner product between two waveform models is defined as

[TABLE]

where the ‘star’ stands for complex conjugation, and $S_{n}(f)$ is the detector spectral noise density. With this definition of the inner product, the SNR is defined as

[TABLE]

In what follows, we consider the Einstein telescope (ET) detector power spectral density noise. The ET is a third-generation ground-based detector of GWs and it is envisaged to be ten times more sensitive in amplitude than the advanced ground-based detectors in operation nowadays, covering the frequency range $1-10^{4}$ Hz. Unlike the current detectors, from the ET conceptual design study, the expected rates of BNS detections per year are of the order of $10^{3}-10^{7}$ ET02 . However, we can expect only a small fraction ($\sim 10^{-3}$) of them accompanied with the observation of a short $\gamma$-ray burst. If we assume that the detection rate is in the middle range around $\mathcal{O}(10^{5})$, we can expect to see $\mathcal{O}(10^{2})$ events with short $\gamma$-ray bursts per year.

Thus, let us consider in our simulations a mock GW standard sirens data set composed by 1000 BNS merger events. Assuming that the errors on $d_{L}$ are uncorrelated with errors on the remaining GW parameters, we have

[TABLE]

Since $\tilde{h}(f)\propto{(d_{L}^{GW}})^{-1}$, then $\sigma_{d_{L}}=d_{L}/{\rm SNR}$. However, when we estimate the practical uncertainty of the measurements of $d_{L}$, we should take the orbital inclination into account. The maximal effect of the inclination on the SNR is a factor of 2 (between $\iota=0^{\circ}$ and $\iota=90^{\circ}$). Therefore, we add this factor to the instrumental error for a conservative estimation. Thus, the estimate of the instrumental error is given by $\sigma_{d_{L}}=2d_{L}/{\rm SNR}$. On the other hand, GWs are lensed in the same way as the electromagnetic waves, resulting into a weak lensing effect error, which we model as $\sigma^{\rm lens}_{d_{L}}=0.05\ z\ d_{L}(z)$ ET02 ; ET_lens . In our study, we do not consider possible errors induced from the peculiar velocity due to the clustering of galaxies. Since we are interested in simulating events at high $z$ mainly, we can neglect such contributions, which are significant only for $z\ll 1$. In fact, at high $z$, the dominant source of uncertainty is the one due to weak lensing. Therefore, the total uncertainty $\sigma_{d_{L}}$ on the luminosity distance measurements associated to each event is obtained by combining the instrumental and weak lensing uncertainties as

[TABLE]

The redshift distribution of the BNS sources is taken to be of the form

[TABLE]

where $d_{C}(z)$ is the comoving distance, and $r(z)$ describes the time-evolution of the burst rate:

[TABLE]

The distribution of the neutron star masses is chosen to be randomly sampled from uniform distributions within $[1-2]~{}{\rm M}_{\odot}$, also under the condition $m_{1}\gtrsim m_{2}$ and $\eta<0.25$. In this case, in the mock data generation we take $\chi_{1}=\chi_{2}=0$, where $\chi_{1},\chi_{2}$ are the associated spin magnitudes on each mass component. Then, we simulate BNS mergers up to $z=2$, which represents the maximum distance at which these events can be observed from the power spectral density noise from the ET ET02 . Also, we checked that, beyond $z=2$, the SNR presents low values, which also can limit the use of the Fisher information for mock data. Now, in order to realistically generate a mock catalog using modified gravity, we shall consider non-zero values for the parameters $\alpha_{M0}$ and $n$, which are compatible with the current cosmological observation as well as with the stability criteria of the theory. Lastly, when generating our mock GW standard sirens data set, we only consider BNS mergers with ${\rm SNR}>8$.

In order to estimate the observational constraints on the free parameters of the models, we apply Markov Chain Monte Carlo (MCMC) method through the Metropolis-Hastings algorithm Metropolis-Hastings , where the likelihood function for the GW standard sirens mock data set is built in the form

[TABLE]

Here, $d_{L}^{obs}(z_{i})$ are the 1000 simulated BNS merger events with their associated uncertainties $\sigma_{d_{L},i}$, while $d_{L}^{th}(z_{i})$ is the theoretical prediction on each $i$th event.

IV Constraints on the running of the Planck mass

In this section, we present and discuss our results regarding the future observational constraints that GW standard sirens can impose on a possible time variation of the Planck mass within of the ET sensitivity. The running of the Planck mass is an important physical quantity, which essentially is present in any and all modified gravity models. To generate a simulated $d_{L}(z_{i})$ catalog using modified gravity, we assume realistic values for the pair ($\alpha_{M0}$, $n$), on each triplet ($z_{i},d_{L}(z_{i}),\sigma_{d_{L}(z_{i})}$) evaluated at each point $i$, as follows.

We first note that the parameter $n$ is statistically degenerate. This fact is already well known and expected to happen. In the literature, it is usual to assume $n=1$, but here we follow the stability conditions discussed in Section II and weakly generate random values for $n$, within the range of stability of the theory: $i$) for the case $\alpha_{M0}>0$, we randomly sampled from uniform distributions: $\alpha_{M0}\in[0,0.5]$ and $n\in[0,1.40]$; $ii$) for the case $\alpha_{M0}<0$, we randomly sampled $\alpha_{M0}\in[-0.5,0]$ and $n\in[2.5,3.5]$. We found that different prior ranges on $n$ change the simulated catalogs very weakly. Only very different prior ranges on $\alpha_{M0}$ can significantly change the pair ($d_{L}(z_{i}),\sigma_{d_{L}(z_{i})}$). The range assumed on $\alpha_{M0}$ is fully compatible with current constraints Horndeski_constraints_01 ; Horndeski_constraints_02 ; Horndeski_constraints_03 ; Horndeski_constraints_04 ; Horndeski_constraints_05 ; Horndeski_constraints_06 ; Horndeski_constraints_07 ; Horndeski_constraints_08 ; Horndeski_constraints_09 ; Horndeski_constraints_10 ; Horndeski_constraints_11 . We used as input values $H_{0}=67.4$ km/s/Mpc and $\Omega_{m0}=0.31$ for the Hubble constant and matter density parameter, respectively, in agreement with the most recent Planck CMB data Planck2018 . Hence, these values are reasonable for our purpose to generate GW standard sirens mock data.

In the realization of the MCMC analysis, the sampling has been done assuming the following uniform priors for the cosmological parameters: $H_{0}\in[55,90]$, $\Omega_{m0}\in[0,1]$, and $\alpha_{M0}\in[-1,0]$, $\alpha_{M0}\in[0,1]$ for each case. Due to the large statistical degeneracy on $n$, as commented above, we fixed $n=3$ and $n=1$ for the cases of $\alpha_{M0}>0$ and $\alpha_{M0}<0$, respectively. Table 1 summarizes the constraints at the 68% and 95% confidence levels (CL). In Figure 2, we show the parametric space and the one-dimensional marginalized distribution for the parameters $\Omega_{m0}$, $H_{0}$ and $\alpha_{M0}$, in both $\alpha_{M0}>0$ and $\alpha_{M0}<0$ cases. We note from both analyses that the parameter $\alpha_{M0}$ is non-null at the 68% CL. Assuming the stability condition where $\alpha_{M0}$ is negative, we find the new lower limit $\alpha_{M0}>-0.2$ at the 95% CL. On the other hand, under the condition that the running of the Planck mass is positive defined, we find that $\alpha_{M0}$ can be non-null up to 95% CL, more specifically $0.03\lesssim\alpha_{M0}\lesssim 0.38$. For a qualitative comparison, we show also the constraints from the $\Lambda$CDM model in both cases. It is important to note that the $\Lambda$CDM model is a particular case of our simulated data, where both mock catalogs are mainly controlled by possible corrections on the amplitude of the running of the Planck mass as explained above. Therefore, extended estimates on the pair ($d_{L}(z_{i}),\sigma_{d_{L}(z_{i})}$) beyond the $\Lambda$CDM prediction will be dependent on the $\alpha_{M0}$ correction, where we have considered two different scenarios $\alpha_{M0}>0$ and $\alpha_{M0}<0$. This should induce a minimal bias when analyzing the parameters within the $\Lambda$CDM scenario due to the $\alpha_{M0}$ prior range in the mock data generation. In the case $\alpha_{M0}<0$, our minimal baseline, i.e. $\Omega_{m0}$ and $H_{0}$, is completely compatible for the two scenarios. On the other hand, when analyzing $\alpha_{M0}>0$, we can note a minimal bias manifested on $\Omega_{m0}$, resulting in a minimum tension at 1$\sigma$ CL on this parameter, but still making the model compatible beyond that statistical significance.

It is interesting to compare our results with others already obtained in the literature. For instance, in Horndeski_constraints_11 , using measurements of the growth rate of structures from DESI survey, it is observed that the amplitude of the running of the Planck mass (quantity physically analogous to our $\alpha_{M0}>0$) can be detected up to 99% CL. In Horndeski_constraints_08 , a 95% CL upper limit of 0.015 is found from CMB data. In MG_dL01 , analyzing the standard siren GW17081 event, the authors found the amplitude of the running of the Planck mass to be $\in$ [-80, 28] at the 95% CL. In Atsushi02 , the amplitude damping $\alpha_{M0}<0$ appear to be preferentially at low $z$ from GWs observations. Other analyzes can be found in Horndeski_constraints_01 ; Horndeski_constraints_02 ; Horndeski_constraints_03 ; Horndeski_constraints_04 ; Horndeski_constraints_05 ; Horndeski_constraints_06 ; Horndeski_constraints_07 ; Horndeski_constraints_08 ; Horndeski_constraints_09 ; Horndeski_constraints_10 ; Horndeski_constraints_11 . We note that the new borders on the amplitude of the running of the Planck mass derived in this work may have also an impact on the modified propagation primordial gravitational waves spectrum Horndeski_constraints_12 .

Now, we shall briefly discuss the consequences of our results. Based on the arguments developed in Section II, we can write the running of the Planck mass as

[TABLE]

One of the surviving classes of models under the condition $c_{T}=c$ are the non-minimal theories in which the scalar field $\phi$ is coupled with the curvature scalar $R$ in the form $G_{4}(\phi)R$. This class includes the metric $f(R)$ gravity and the Brans-Dicke theory Brans61 . The original Brans-Dicke theory, for instance, is obtained by setting $G_{4}=\phi$. By substituting this in Eq. (27), it is possible to obtain $\phi$ as a function of the cosmic time. In Figure 3, we show a reconstruction for the evolution of the field $\phi/\phi_{0}$ in Planck mass units, where $\phi_{0}$ is $\phi(z=0)$. Also, it is important to remember that in such a theory the gravitational constant is not presumed to be constant, but $G_{N}(\phi)\varpropto 1/\phi$. This fact is physically encoded in $\alpha_{M}(a)$, which measures the gravity strength. We can note that at late times, the gravity strength, $G_{N}$, is greater (smaller) than predicted by GR when $\alpha_{M0}>0\,(<0)$, up to 20% (30%) at $z=0$. On the other hand, when $z\gg 0$, GR is recovered and we do not expect to have significant variations at early times.

V Constraints on parametric $f(R)$ gravity

In this section, we briefly review $f(R)$ cosmology and show new observational constraints on the Hu-Sawicki (HS) model from our standard siren mock data set. We refer to fR1 ; fR2 ; Capozziello11 ; Nojiri17 for reviews on $f(R)$ gravity.

The $f(R)$ gravitational theories consist of extending the Einstein-Hilbert action in the form

[TABLE]

where $f(R)$ is a function of the Ricci scalar, and $S_{m}$ is the action for matter fields. For $f(R)=R$, the GR case is recovered.

Let us consider a spatially flat FLRW Universe dominated by pressureless matter (baryonic plus dark matter) and radiation with energy densities $\rho_{m}$, $\rho_{r}$ and pressures $P_{m}$, $P_{r}$, respectively. The modified Friedmann equations in the metric formalism are given by fR1 ; fR2

[TABLE]

Moreover, one obtains the following useful relation:

[TABLE]

It can be shown that, through the transformation of the scalar degree of freedom, $\phi=M_{P}df(R)/dR$, the metric $f(R)$ gravity is equivalent to the Brans-Dicke theory (with $w_{BD}=0$). Thus, when considering a parametric function $f(R)$, given our constraints on $\phi$, some bounds can also be found on $f(R)$ gravity. For example, adopting the formalism presented in Basilakos ; Rafael (and reference therein), one can write

[TABLE]

where the function $y(R,b)$ quantifies the deviation from Einstein gravity, i.e. the effect of the $f(R)$ modification, through the distortion parameter $b$. Then, the scalar field can be expressed as

[TABLE]

Interpreting $\phi_{0}=M_{P}$ and given a function $y(R,b)$, we can use our limits on $\phi/\phi_{0}$ and place observational bounds on $b$. We note that the stability conditions assumed in Section II and used in the development of this section are completely in agreement with $f(R)$ gravity alphai_02 . Clearly, more direct observational boundaries on $f(R)$ gravity can be obtained by modifying appropriately Eq. (12) to include a function $f(R)$. We are motivated to present a more detailed study in this sense in a future communication.

Without loss of generality, we can use our simulated standard siren catalog from BNS mergers, within a parametric limit and certain stability conditions, to model the $f(R)$ gravity dynamic. Thus, in what follows, we consider a parametric $f(R)$ gravity scenario and investigate the observational bounds that our standard siren mock data set can impose on the free parameters of the theory.

Let us consider viable models that have up to two parameters, where $f(R)$ function is given by Eq. (32). This methodology has been used earlier to investigate the observational constraints on $f(R)$ gravity in Basilakos ; Rafael . In this respect, one of the most well-known models in the modified gravity theory literature is the Hu-Sawicki (HS) model Hu:2007nk , which satisfies all the dynamics conditions required for a given $f(R)$ function. The function $y(R,b)$ for the HS model is given by

[TABLE]

where $n$ is an intrinsic parameter of the model. In what follows, we assume $n=1$ and refer to Basilakos ; Rafael for details.

As a direct application of the standard siren events from BNS mergers, within the sensitivity predicted for the ET, we proceed to constrain the HS model using our mock data set generated by the condition $\alpha_{M0}<0$. Together with the standard siren data, in the present analysis we also employed the type Ia Supernova (SN) Pantheon data Scolnic18 and the cosmic chronometers (CC) measurements Moresco12 in order to obtain tighter constraints on the free parameters of the HS scenario (see Appendix B for the details on the SNe and CC data.) In the MCMC analysis, we assumed uniform priors for the cosmological parameters: $H_{0}\in[55,90]$, $\Omega_{m0}\in[0,1]$, and $b\in[0,1]$. In Figure 4, we show the parametric space and the one-dimensional posterior distributions for the HS baseline model from SN + CC + GW and SN + CC analysis. Thus, we can quantify how much the addition of the GWs standard siren data can improve the constraints and break the degeneracy on the parametric space of the model, in particular on $b$. For comparison, we show also the constraints for the $\Lambda$CDM model. In particular, we found the following mean values with the relative 68% CL and 95% CL errors from the joint analysis (SNe + CC + GWs):

[TABLE]

We note that $b$ tends to be non-null up to the 95% CL with the addition of GWs data, therefore breaking the degeneracy enough to have a non-null value on $b$ up to 2$\sigma$. This demonstrates the potential of the future standard sirens catalogs in joint analysis with another geometric probes. As a comparison, in Rafael it was found $b<0.50$ and $b<0.13$ at the 95% CL from CC + $H_{0}$ and JLA + BAO + CC + $H_{0}$, respectively, while in Basilakos it was found $b<0.25$ at the 68% CL. The data used in these works are from different physical nature and accuracy and, thus, a direct comparison seems not to be appropriate. However, all these constraints are compatible with each other within the 95% CL. Also, we can note how much the addition of GWs data can improve the constraints on $\Omega_{m0}$ and $H_{0}$ parameters. A detailed analysis from GW standard sirens on other viable $f(R)$ models as well as full discussion will be presented in a future communication.

Finally, it is interesting to analyze the consequences of our results on the effective dark energy equation of state parameter. This can be expressed as Basilakos

[TABLE]

where $E(a)=H(a)/H_{0}$ and

[TABLE]

In the case of the standard $\Lambda$CDM model, $w_{DE}=-1$ throughout the entire cosmological evolution. In Figure 5, we show the $1\sigma$ and $2\sigma$ reconstructions of $w_{DE}(z)$ for the Hu-Sawicki model compared to the prediction of the standard $\Lambda$CDM model. The non-vanishing value of $b$ makes the effective dark energy term behave as quintessence at late times for $z<0.4\ (0.1)$ at the 68% (95%) CL. Beyond this range, the effective equation of state is compatible $w_{DE}=-1$ at the 95% CL. The green line represents the evolution of $w_{DE}$ from the mean values of the MCMC analysis on SNe + CC + GWs.

VI Final remarks

Assuming an effective luminosity distance for GWs, which can be physically interpreted as possible corrections on the GW amplitude propagation between source and detector, we obtained new observational constraints on the running of the Planck mass under appropriate stability conditions (see summary in Table 1), from 1000 standard siren events from binary neutron star mergers within the ET power spectral density noise. We found that the GW amplitude damping correction can be preferentially non-zero at the 68% CL and 95% CL, for $\alpha_{M0}<0$, $\alpha_{M0}>0$, respectively.

Furthermore, we combined our simulated GW data with the latest available SNe and CC measurements to constrain the parametric Hu-Sawicki gravity model. In doing so, we found that the deviation from GR may be evidenced as non-null at the 95% CL. These results demonstrate the statistical accuracy that can be achieved by future ground based GW observatory as ET detector.

The underlying assumption of our analysis is $c_{T}=c$ at all redshifts. This hypothesis is strongly motivated from the GW170817 constraints to hold locally. On the other hand, the evidence for $c_{T}=c$ might not be falsified due to the lack the gravitational wave observations at high $z$, even in the near future, since within the LIGO/VIRGO sensitivity, we expect to measure standard siren events from BNS only at very low $z$. Thus, as significant deviations from GR are expected only at moderate to high $z$ and/or large scale, it may be interesting to relax the condition $c_{T}/c=1$ and perform forecast analyses on the ratio $c_{T}/c$ from binary systems at high $z$ within the sensitivity of the future detectors such as ET and DECIGO, and check for possible consequences on modified gravity phenomenology. It would be also interesting to apply a similar methodology as developed here to investigate other well-motivated classes of modified gravity models, as well as to perform joint analysis between standard siren mock events and CMB data (present and future missions). In this way, we will be able to determine the level of deviations from GR as suggested by the combination of future astronomical observations in the next decade.

Acknowledgements.

The authors are grateful to Salvatore Capozziello and Jose C. N. de Araujo for useful comments and discussions. R.C.N. would like to thank FAPESP for financial support under the project # 2018/18036-5. The authors also wish to thank the referee for his/her constructive comments.

Appendix A: PN coefficients

For the convenience of the reader, we list below the PN coefficients TF2 we used for the waveform model Eq. (15) and Eq. (17). The individual masses and spin parameters, $m_{i}$ and $\chi_{i}$ $(i=1,2)$, are encoded in the following parameter combinations:

[TABLE]

The PN amplitude expansion coefficients are:

[TABLE]

The phase $\Phi(f)$ expansion coefficients are:

[TABLE]

Appendix B: SNe and CC data set

We provide here some details of the low-redshift cosmological observables that we used to complement the GW mock data in the statistical analysis on the HS model.

The first data set is the Pantheon sample Scolnic18 , composed of 1048 Supernovae (SNe) Ia in the redshift range $0.01<z<2.3$. In this compilation, all the SNe are standardized through the SALT2 light-curve fitter, in which the distance modulus is modelled as follows SALT2 :

[TABLE]

where $m_{b}$ is the $B$-band apparent magnitude of each SN and $M$ is its absolute magnitude, while $\Delta_{M}$ and $\Delta_{B}$ account for the host-mass galaxy and the distance bias corrections, respectively. Moreover, $x_{1}$ and $C$ are the stretch and color parameters of each SN light-curve, respectively, with their relative coefficients $\alpha$ and $\beta$. On the other hand, the distance modulus predicted by a cosmological model is given as

[TABLE]

As shown in Riess18 , under the assumption of a flat universe one can compress the full SN sample into a set of 6 cosmological model-independent measurements of $E(z)^{-1}$, where $E(z)\equiv H(z)/H_{0}$ is the dimensionless Hubble parameter. This approach allows us to properly marginalize over the SN nuisance parameters in the fit. Thus, taking into account the correlations among the $E^{-1}(z)$ measurements, we can write the likelihood function associated to the SN data as DAgostino19

[TABLE]

where $\mathbf{A}\equiv E^{-1}_{obs,i}-E^{-1}_{th}(z_{i})$ quantifies the difference between the measured values and the values predicted by a cosmological model, and $\mathbf{C}_{\text{SN}}$ is the covariance matrix resulting from the correlation matrix given in Riess18 .

The second data set is built upon the differential age approach Jimenez02 , which represents a model-independent method to characterize the expansion of the universe up to $z<2$. In this technique, passively evolving red galaxies are used as cosmic chronometers (CC) to measure the age difference $(dt)$ of the universe at two close redshifts $(dz)$. Thus, one can estimate the Hubble parameter as

[TABLE]

In our analysis, we used the compilation of 31 $H(z)$ uncorrelated measurements collected in Capozziello18 (see references therein). We can then write the likelihood function relative to the CC data as

[TABLE]

where $H_{obs,i}$ are the observed measurements with their relative uncertainties $\sigma_{H,i}$, while $H_{th}(z_{i})$ are the theoretical values of the Hubble parameter estimated from using a specific cosmological model.