1 Introduction

It is well known that measuring the relative luminosity distance of different astrophysical sources, or the relative time-of-flight of different ultra-relativistic particles emitted by the same source (see e.g. [1]), we can obtain important information and constraints on the fundamental cosmological parameters. Recently, the simultaneous detection of gravitational wave (GW) and electromagnetic (e.m.) signals with the same origin—namely, the binary neutron star merger GW170817 [2] and the associate \(\gamma \)-ray burst GRB170817A [3,4,5,6]—has opened the possibility of including also GW observations among the useful tools for testing the standard \(\varLambda \hbox {CDM}\) model and its possible modifications. Particularly useful (as noted long ago [7]) to that purpose are GW sources like coalescing compact binaries located at different cosmic distances, which can play the role of “standard sirens” analogous to the role of the “standard candles” typical of e.m. signals.

In such a context it is important to note that, by comparing the propagation of GW signals received from a given source and the propagation of other ultra-relativistic signals emitted by the same source, we have also the possibility of testing alternative theories of gravity on a cosmic scale of distance, as recently stressed and discussed in [8,9,10,11] (see also [12,13,14,15]).

We should recall, in this respect, that in the standard Einstein’s theory of gravity both GW and e.m. signals propagate with the same speed along the light cone of the given background geometry, thus defining exactly the same luminosity distance \(d_L\) for the emitting source. In models of gravity different from General Relativity, on the contrary, the propagation of GW and e.m. signals may be different, and not only because the metric satisfies modified Einstein equations, whose solution may define a background geometry different from the standard one. In addition, in fact, the propagation of massless signals may deviate from the usual light-cone propagation, and the deviation of GW signals may be different from the deviation of e.m. signals.

Such a modified propagation may be associated, in particular, to a different speed and/or to a different variation of the amplitude of the signal with the distance from the source (for what concerns a possible modified speed of GW signals, however, there are strong constraints imposed by present observations [5], see also [16,17,18,19]). The modified variation of the amplitude, in its turn, produces an effective variation of the associated luminosity distance. Hence, when computing the relevant luminosity distance in the context of a generalized theory of gravity, we may have to face the (in principle interesting) situation in which \(d_L^{GW} \not = d_L^{em} \not =d_L^{GR}\), where the last value of \(d_L\) refers to the one computed for the cosmological scenario based on the General Relativity equations. Comparing the above (possibly different) three expressions of \(d_L\) may thus represent an efficient tool for testing alternative theories of gravity on cosmological scales of distance [8,9,10,11].

In this paper we will concentrate on a simple generalization of the Einstein theory which contains only one additional (constant) parameter, controlling a quadratic curvature correction to the Einstein action. The model is also known as “Starobinsky gravity” [20], and can be regarded as a truncated higher-curvature expansion of the gravitational action. Also, it corresponds to a simple type of f(R)-gravity and can be related to the largest class of Horndeski theories, widely used to test gravity through modified GW propagation (see in particular [21] and [11, 13, 14]). It was shown [22] that such a model can fit the Hubble diagram of Supernovae, consistently with other recent cosmic data, even in the absence of an explicit dark-energy contribution. Here we show that such a simple modification of Einstein’s gravity is already enough to produce an interesting difference between the propagation of e.m. and GW signals, and to predict possibly testable effects when estimating the luminosity distance of standard sirens.

The paper is organized as follows. In Sect. 2 we recall some general result concerning the propagation of tensor metric perturbations in modified theories of gravity and the corresponding definition of the modified luminosity distance. In Sect. 3 we apply such result to a gravitational action with quadratic curvature corrections. In Sect. 4, by assuming that the deviations from the Einstein theory are small, we obtain an approximate expression for the GW luminosity distance by expanding to first order the solution of the modified cosmological equations around the standard-model solution. In Sect. 5 we numerically evaluate the parameters of such a modified cosmological solution by fitting the UNION 2 data set [23] of Supernovae. In Sect. 5.1 we compare our results with the constraints imposed by the standard siren GW170817 and its e.m. partner GRB170817A. In Sect. 5.2 we briefly discuss the differences between the GW and e.m. luminosity distance and the standard luminosity distance \(d_L^{GR}\), based on the Einstein equations and on exact light-cone propagation. Finally, Sect. 6 is devoted to our conclusive remarks. In the “Appendix” we present an explicit computation of the propagation equation for tensor metric perturbations for our model with higher-curvature corrections.

2 GW luminosity distance in generalized models of gravity

Let us start by recalling, for the reader’s convenience, a few well know results concerning the propagation of tensor metric perturbations in a spatially flat, four-dimensional geometry of the Friedmann–Lemaitre–Robertson–Walker (FLRW) type, described in conformal time \(\eta \) by the background metric \(g_{00}= a^2(\eta )\), \(g_{ij}= -a^2(\eta ) \delta _{ij}\). In the so-called transverse-traceless gauge, each polarization component h of the first order metric perturbation \(h_{\mu \nu }\) satisfies the linearized propagation equation

$$\begin{aligned} {d^2h\over d\eta ^2}+ 2 \mathcal{H}{d h\over d\eta }-\nabla ^{2}h=0, \end{aligned}$$
(1)

where we have defined \(\mathcal{H}=a^{-1}da/d\eta \) (we have assumed no anisotropic stress in the matter sources). The corresponding canonical variable [24] \(v=ah\) satisfies

$$\begin{aligned} {d^2v\over d\eta ^2}-\left( \nabla ^{2}+{1\over a}{d^2 a\over d\eta ^2}\right) v=0. \end{aligned}$$
(2)

By expanding the solution in Fourier modes, let us consider those modes well “inside” our present horizon, satisfying the condition \(\left| a^{-1} d^2a /d\eta ^2 v \right| \ll \left| \nabla ^2 v\right| \), like (in particular) the modes relevant to the sensitivity frequency band of present interferometric detectors. For such modes the equation for v is unaffected by the background geometry, so that the amplitude of a gravitational signal \(h=v/a\) emitted at the time \(\eta \) and received at the time \(\eta _0>\eta \) decreases, during the propagation from the source to the observer, by an amount determined by the standard redshift factor, i.e. by the factor \(a(\eta )/a_0 \equiv (1+z)^{-1}\), where \(a_0\equiv a(\eta _0)\), and z is the redshift of the emitting source.

On the other hand, if we are interested in GW signal emitted by coalescing binaries (like the so-called standard sirens), it is well known [25, 26] that such a “weakening” of the amplitude (due to the propagation across the expanding cosmic geometry) leads to introduce the so-called luminosity distance of the emitting source defined in terms of the received flux of radiation, and thus inversely proportional to the amplitude h of the received signal, namely \(h \sim d_L^{-1}(z)\). Here \(d_L\) is the standard luminosity distance referred to signals propagating on the light-cone of our FLRW geometry, and given, for a metric with scale factor a and Hubble factor \(H= a^{-2} da/d\eta \), by

$$\begin{aligned} d_{L}(z) = (1+z)\int ^z_0 \frac{dz'}{H(z')} \end{aligned}$$
(3)

(see also [8,9,10,11, 13, 14] for a detailed discussion of how the GW amplitude h is related to the GW luminosity distance).

Let us now consider, following [8,9,10,11], a generalized model of gravity different from Einstein’s theory and leading, in particular, to a different set of cosmological equations. Their perturbation gives us a modified equation for the tensor modes h which, even in the same FLRW metric as before, takes a form different from Eq. (1):

$$\begin{aligned} {d^2h\over d\eta ^2}+ 2 \mathcal{H}{d h\over d\eta } \left[ 1- \delta (\eta )\right] -c_s^2(\eta ) \nabla ^{2}h=0. \end{aligned}$$
(4)

We may have, in general, a modified speed parameter, \(c_s\not =1\), and a modified geometric friction term with \(\delta \not =0\). Both parameters \(c_s\) and \(\delta \) may be in principle time-dependent. Even in that case, however, we can still introduce the analogous of a canonical variable v defined by \(v\equiv \widetilde{a} h\), where \(\widetilde{a}\) satisfies

$$\begin{aligned} {1\over \widetilde{a}} \frac{d\,\widetilde{a}}{d\eta }=\mathcal{H}\left[ 1-\delta (\eta )\right] , \end{aligned}$$
(5)

and obtain the corresponding propagation equation in canonical form (but with a modified “pump field”, \(a \rightarrow \widetilde{a}\)) as follows:

$$\begin{aligned} {d^2v\over d\eta ^2}-\left( c_s^2\,\nabla ^{2}+{1\over \widetilde{a}}{d^2 \widetilde{a}\over d\eta ^2}\right) v=0. \end{aligned}$$
(6)

Let us concentrate, again, on modes well inside the horizon, and on a theory of gravity which implies \(c_s=1\). In that case, the only difference from the previous case described by the perturbation of the Einstein equations is that the amplitude \(h= v/\widetilde{a}\) of the GW signal, while traveling from the source to the observer, is now decreased (or increased) by the modified factor \(\widetilde{a}(\eta )/\widetilde{a}_0 \equiv (1+\widetilde{z})^{-1}\) (instead of the standard redshift factor \(a(\eta )/a_0\)). Hence, the received amplitude turns out to be inversely proportional to a modified luminosity distance, which we may call \(d_L^{GW}\), and which we may relate to the standard luminosity distance \(d_L\) as follows [8,9,10,11]:

$$\begin{aligned} {1\over d_L^{GW}(z)}= & {} \left( 1+z\over 1+ \widetilde{z}\right) {1\over d_L(z)}\equiv {\widetilde{a} (z)\over \widetilde{a}_0} {a_0\over a(z)}{1\over d_L(z)}= {\widetilde{a}(z)\over a(z)} {1\over d_L(z)} \end{aligned}$$
(7)

Here \(d_L\) is the distance given by Eq. (3) and associated with the standard light-cone propagation in the geometry with scale factor a. Also, we have normalized the solution of Eq. (5) for \(\widetilde{a}\) in such a way that \(\widetilde{a}(\eta _0) = a(\eta _0)\), so as to recover the obvious result \(d_L^{GW}(0)=d_L(0)\) for a source at zero cosmic distance. With such a normalization, a simple integration of Eq. (5), rewritten as \(d(\log a/\widetilde{a})/d\eta = a^{-1}(da/d\eta ) \delta \), leads to express \(\widetilde{a}\) in terms of a as follows [8,9,10]:

$$\begin{aligned} {a\over \widetilde{a}}(z) = {d_L^{GW}\over d_L}(z)= \exp \left[ -\int _0^z {dz'\over 1+z'}\delta (z')\right] . \end{aligned}$$
(8)

In the following sections we will apply the above results to a generalized model of gravity with quadratic curvature corrections.

3 A simple model of quadratic gravity

The generalized model of gravity we shall consider here is a particularly simple case of the so-called f(R)-theories (see e.g. [27,28,29] for comprehensive reviews), with a function \(f(R)=R+c_2R^2\), which we may also regard as an expansion in power series of the scalar curvature R, truncated to second order. Our generalized action is then

$$\begin{aligned} S= -{M_{\mathrm{P}}^2\over 2 } \int d^4 x \sqrt{-g} \left( R+ c_2 R^2 \right) + S_m, \end{aligned}$$
(9)

where \(S_m\) is the matter action (assumed to describe perfect fluid sources), \(M_{\mathrm{P}}=1/\sqrt{8\pi G}\) is the Planck mass, and \(c_2\) an unknown constant parameter with dimensions of length squared (we are using units \(\hbar =c=1\)).

By varying with respect to the metric \(g_{\mu \nu }\), and by adding to the action (9) an appropriate generalized version of the York–Gibbons–Hawking action (see e.g. [22]), so as to get rid of the unwanted boundary contributions (as in General Relativity), we obtain the following generalized Einstein equations

$$\begin{aligned}&8 \pi G \,T_\mu ^\nu =R_\mu ^\nu -{1\over 2} \delta _\mu ^\nu R +c_2 \left[ 2RR_\mu ^\nu + \delta _\mu ^\nu \left( 2 \Box R-{1\over 2} R^2 \right) -2 \nabla _\mu \nabla ^\nu R\right] ,\nonumber \\ \end{aligned}$$
(10)

which are fourth-order differential equations for the components of the metric tensor. We have defined \(\Box \equiv \nabla _\mu \nabla ^\mu \), and \(T_{\mu \nu }\) is the source stress tensor, given by the variation with respect to \(g_{\mu \nu }\) of the matter part of the action. It can be easily checked that by taking the covariant divergence of both sides of the above equation, and using the contracted Bianchi identity, one consistently recovers the standard result

$$\begin{aligned} \nabla _\nu T_\mu \,^\nu =0, \end{aligned}$$
(11)

which guarantees the covariant conservation of the stress tensor of the gravitational sources.

Let us specify now the above equations to the case of a spatially flat FLRW metric with scale factor a, denoting with a dot the derivative with respect to the cosmic time coordinate t, related to conformal time by \(dt=ad\eta \). We shall adopt the following conventions: \(g_{\mu \nu }= \mathrm{diag} (+---)\), \(R_{\nu \alpha }= R_{\mu \nu \alpha }\,^\mu \), and \(R_{\mu \nu \alpha }\,^\beta = \partial _\mu \varGamma _{\nu \alpha }\,^\beta - \cdots \). We then find, in particular,

$$\begin{aligned} R_0^0=-3\left( \dot{H} +H^2\right) , ~~~ R_i^j= -\left( \dot{H} +3H^2\right) \delta _i^j, \end{aligned}$$
(12)

and

$$\begin{aligned} R=-6\left( \dot{H} +2H^2\right) , \end{aligned}$$
(13)

where \(H=\dot{a}/a\). For a perfect fluid source, with \(T_\mu ^\nu = \mathrm{diag}(\rho , -p\delta _i^j)\), the time component of Eq. (10) thus gives

$$\begin{aligned} 3H^{2} +18 c_{2}\left[ {\dot{H}}^{2}-6H^{2}{\dot{H}}-2H{\ddot{H}}\right] = 8 \pi G \rho . \end{aligned}$$
(14)

The space-like (diagonal) components of Eq. (10) lead to the condition

$$\begin{aligned} 2{\dot{H}}+3H^{2} -6c_{2}\left[ 9{\dot{H}}^2+18{\dot{H}}H^{2}+12H{\ddot{H}}+ 2\dddot{H}\right] =- 8 \pi Gp. \end{aligned}$$
(15)

The mixed components of Eq. (10) are identically satisfied. Finally, the conservation Eq. (11) takes the usual form

$$\begin{aligned} {\dot{\rho }} + 3H(\rho +p)=0. \end{aligned}$$
(16)

To obtain the evolution equation for the first-order tensor perturbations of the metric, defined by \(g_{\mu \nu } \rightarrow g_{\mu \nu }+\delta g_{\mu \nu }\), \(\delta g_{\mu \nu }=h_{\mu \nu }\), we shall work in the transverse-traceless (TT) gauge \(\nabla _\nu h_\mu ^\nu =0=h_\mu \,^\mu \), which in our context takes the simple form \(h_{\mu 0}=0\), \(\nabla _jh_i\,^j\equiv \partial _j h_i\,^j=0=g^{ij}h_{ij}\). By directly perturbing the background equations (10) (see “Appendix” for an explicit computation) we then obtain for each polarization mode of \(h_i\,^j\) a linearized equation which, expressed in conformal time, can be written as

$$\begin{aligned} {d^2 h\over d\eta ^2}+2\mathcal{H}\left[ 1+ {c_2 (dR/d\eta )\over \mathcal{H}(1+2c_2R)}\right] {dh\over d\eta }-\nabla ^2 h=0, \end{aligned}$$
(17)

where R is the scalar curvature (13). Comparing with the modified Eq. (4) we then find, for our model, a standard propagation velocity \(c_s=1\), but a non-standard friction term which can be written in the form (4) with the time-dependent parameter

$$\begin{aligned} \delta (\eta )=- {c_2 (dR/d\eta )\over \mathcal{H}(1+2c_2R)}. \end{aligned}$$
(18)

It is important to recall, at this point, that the absence of modifications of the standard speed parameter is a peculiar property of the model of gravity we are considering. More general quadratic-curvature corrections to the Einstein–Hilbert action can lead indeed to a propagation Eq. (4) with an effective speed parameter \(c_s \not =1\). An important example of this type is provided by the Gauss–Bonnet two-form typically predicted in a string-theory context, and generated by the high-curvature \(\alpha '\) corrections to the String-frame gravi-dilaton action (see [30] for an explicit computation of the modified propagation of tensor perturbations in the presence of the quadratic Gauss–Bonnet term in the gravi-dilaton action).

It may be useful to check, also, that the previous results (17), (18) can be exactly reproduced starting with the modified action (9), and expanding the action up to terms quadratic in the metric fluctuations \(h_i\,^j\) (so as to obtain an effective Lagrangian for the first-order perturbations, from which to deduce their equations of motion through the standard variational formalism).

In fact, working again in the TT gauge, we have the following quadratic contributions to the perturbed action (see also [30,31,32]):

$$\begin{aligned} \delta ^{(2)}\sqrt{-g}= & {} -{1\over 4} \sqrt{-g}\, {\mathrm{Tr}} \,(h^2), \end{aligned}$$
(19)
$$\begin{aligned} \delta ^{(2)}(\sqrt{-g}R)= & {} a^{3}{\mathrm{Tr}} \left[ \left( \frac{3}{2}{\dot{H}}+3H^{2}\right) h^{2}+h\ddot{h} \, \right. \nonumber \\&+ \left. 4H h {\dot{h}}+\frac{3}{4}{\dot{h}}^{2}-\frac{1}{4}h\frac{\nabla ^{2}}{a^{2}}h\right] , \end{aligned}$$
(20)
$$\begin{aligned} \delta ^{(2)}(\sqrt{-g}R^{2})= & {} R \, a^{3}{\mathrm{Tr}} \left[ \left( \frac{3}{2}{\dot{H}}+3H^{2}\right) h^{2}+2h\ddot{h} \, \right. \nonumber \\&+ \left. 8H h {\dot{h}}+\frac{3}{2}{\dot{h}}^{2}-\frac{1}{2}h\frac{\nabla ^{2}}{a^{2}}h\right] , \end{aligned}$$
(21)

where \({\mathrm{Tr}}(h^{2})={h_{i}}^{j}{h_{j}}^{i}\), \({\mathrm{Tr}}(h{\dot{h}})={h_{i}}^{j}\dot{ h_{j}}^{i}\), and so on. Summing up all contributions, integrating by parts the terms in \(h{\dot{h}}\) and \(h\ddot{h}\), and neglecting a total derivative which does not contribute to the equations of motion, it turns out that the coefficient of the \(h^2\) term appearing in the quadratic perturbed action is identically vanishing, thanks to the background equation (15). By decomposing the fluctuations into the two physical polarization modes \(h_+\) and \(h_\times \), i.e. \({\mathrm{Tr}}(h^{2})={h_{i}}^{j}{h_{j}}^{i}= 2 (h_+^2 + h_\times ^2)\), and switching again to the confromal time coordinate where \(\dot{h}= a^{-1}dh/d\eta \), we are eventually led (for each polarization mode) to the following quadratic effective action:

$$\begin{aligned} \delta ^{(2)}S=\frac{M_{\mathrm{P}}^2}{4}\int d^{3}xd\eta \,a^{2}\left( 1+2c_{2}R\right) \left[ \left( dh\over d\eta \right) ^{2}+h\nabla ^{2}h\right] . \end{aligned}$$
(22)

We can now easily recover the standard form of the Einstein perturbed action by defining a modified scale factor \(\widetilde{a}\), related to a (and satisfying the correct normalization \(\widetilde{a}_0=a_0\), see Sect. 2) in such a way that

$$\begin{aligned} \widetilde{a}(\eta )= a(\eta ) \left[ 1+2c_2 R(\eta )\over 1+2 c_2 R(0)\right] ^{1/2}. \end{aligned}$$
(23)

By replacing a with \(\widetilde{a}\) in the action (22), and applying to h the standard variational formalism, we are then led to the equation of motion for tensor perturbations in the form

$$\begin{aligned} {d^2 h\over d\eta ^2}+ {2\over \widetilde{a}}{d \widetilde{a} \over d\eta }{dh\over d\eta }-\nabla ^{2}h=0. \end{aligned}$$
(24)

Finally, by computing \(\widetilde{a}'/\widetilde{a}\) from the definition (23), and comparing the result with Eq. (17), we can immediately recover the expression of the parameter \(\delta (\eta )\) anticipated in Eq. (18).

To conclude this section let us stress that by following this second procedure, based on the computation of the quadratic perturbed action, we can directly obtain from Eq. (23) the important ratio \(d_L/d_L^{GW}= \widetilde{a}/a\) [see Eq. (8)], with no need of performing the explicit integration of the parameter \(\delta (\eta )\) along the travel of the signal from the source to the observer.

4 First-order corrections to the standard cosmological dynamics

In the previous section we have introduced a modified theory of gravity where the propagation of GW signals, and the associated luminosity distance \(d_L^{GW}\), are different from those predicted by the standard Einstein theory. In the following sections we will discuss the possibility that such a modified GW dynamics may be compatible with (and possibly useful to interpret) present cosmological observations, by concentrating our attention, in particular, on the distance scales typical of the Supernovae and of the GW sources known as standard sirens.

In order to compare the model predictions with the observational data in that range of distances, what we need, first of all, is to express both the friction correction \(\delta \) and the ratio \(\widetilde{a}/a=d_L/d_L^{GW}\) as a function of the redshift parameter z. To this purpose, by recalling the definition \((1+z)=a_0/a\), we first note that \(\dot{z}= - (1+z)H\), so that

$$\begin{aligned} {d R\over d\eta }= a \dot{R}= a \dot{z}{dR\over dz} = -(1+z) \mathcal{H} {dR\over dz}. \end{aligned}$$
(25)

We shall adopt, from now on, the convenient notation in which a prime denotes differentiation with respect to z. We can then rewrite Eqs. (18) and (23), respectively, as

$$\begin{aligned}&\delta (z)=c_{2}\frac{(1+z)R'}{1+2c_{2}R}, \end{aligned}$$
(26)
$$\begin{aligned}&\frac{d_{L}^{GW}}{d_{L}}(z)=\left[ 1+2c_{2}R(0)\over 1+2c_{2}R(z)\right] ^{1/2}, \end{aligned}$$
(27)

where, from Eq. (13),

$$\begin{aligned}&R(z)=-6H\left[ 2H-(1+z)H'\right] , \end{aligned}$$
(28)
$$\begin{aligned}&R'(z)=-6\left[ 3HH'-(1+z)\left( H^{\prime 2} + H H''\right) \right] . \end{aligned}$$
(29)

It may be noted that, in spite of the fact that the exact definition of the modified luminosity distance (27) is formally equivalent to the GW luminosity distance defined in the presence of a time-varying effective Planck mass (considered e.g. in [11]), the particular z-dependence predicted by the above equations is quite different, in principle, from the z-dependence of \(d_L^{GW}\) assumed in [11] for their phenomenological discussion.

For an explicit evaluation of the modified luminosity distance (27) we need now H(z), which has to be determined by solving the modified cosmological equations (14), (15). In addition, and for a numerical estimate of the differences with the standard results, we should determine a possible allowed range of values for the parameter \(c_2\). To this second purpose we will follow the procedure already suggested in a previous paper [22], namely fitting a set of Supernovae data with an expression of the luminosity distance \(d_L\) computed according to the standard definition (3), but with an Hubble parameter H(z) which satisfies the modified cosmological equations of this paper.

It is important to stress, in this respect, that in spite of the modified gravitational action, the matter part of our model (9) keeps unchanged. This implies, in particular, that the e.m. field is still minimally coupled to the background geometry (like in General Relativity), and that the propagation equation of e.m. signals keeps the same form as in the standard Einstein theory. Hence, the e.m. luminosity distance, \(d_L^{em}\), is still given by Eq. (3), with the only difference that H(z) satisfies the generalized cosmological equations (14), (15) with \(c_2 \not =0\). A fit of the Supernovae data with such a generalized form of \(d_L(z)\) can thus provide a reliable estimate of the possibly allowed value of \(c_2\).

In order to obtain the required solution for H(z) we shall work under the assumption that the higher-curvature corrections to the standard equations are small, and that the solution of the modified equations can be expanded in power series of \(c_2\) around the zeroth-order solution of the standard cosmological equations.

Let us consider, in particular, Eqs. (14), (16) as our independent equations, and use the standard (\(\varLambda \hbox {CDM}\)) model of perfect fluid source with two non-interacting components, namely \(\rho = \rho _m+\rho _\varLambda \), with \(p_m=0\) and \(p_\varLambda = - \rho _\varLambda \). The conservation equation (16), expressed in terms of z, thus reduces to

$$\begin{aligned} \rho _\varLambda '=0~, ~~~~~~~~~~ (1+z) \rho _m' =3 \rho _m ~, \end{aligned}$$
(30)

and can be solved exactly to give \(\rho _\varLambda = \rho _{\varLambda 0}\) and \(\rho _m= \rho _{m 0} (1+z)^3\), where \(\rho _{\varLambda 0}\) and \( \rho _{m 0}\) are integration constants. By inserting such solution into the generalized Friedmann equation (14) we then find that the exact equation for H(z) takes the form:

$$\begin{aligned}&3H^2+ 18 c_2 H^2 \left[ 4(1+z) H H' - (1+z)^2 (H^{\prime 2} +2H H'') \right] \nonumber \\&\quad = 8 \pi G \left[ \rho _{m 0} (1+z)^3 + \rho _\varLambda \right] . \end{aligned}$$
(31)

Let us now expand the solution of the above equation, to first order in \(c_2\), as follows:

$$\begin{aligned} H(z)= H^{(0)}(z)+ c_2 H^{(1)}(z)+ \cdots , \end{aligned}$$
(32)

where \(H^{(0)}\) is the solution of the standard equation with \(c_2=0\), namely

$$\begin{aligned} H^{(0)}(z) = H_0 \left[ \varOmega _{m0}(1+z)^3+\varOmega _\varLambda \right] ^{1/2}, \end{aligned}$$
(33)

and where we have defined \(\varOmega _{m0}= 8 \pi G \rho _{m0}/3H_0^2\) and \(\varOmega _\varLambda = 8 \pi G \rho _{\varLambda 0}/3H_0^2\). To first order in \(c_2\) we then obtain a simple algebraic equation for \(H^{(1)}\), which can be easily solved and leads to:

$$\begin{aligned} H^{(1)}(z) = - {27\over 4} H_0^3 \varOmega _{m0}^2 (1+z)^6 \left[ \varOmega _{m0}(1+z)^3+\varOmega _\varLambda \right] ^{-1/2}. \end{aligned}$$
(34)

In order to give to \(\varOmega _{m0}\) the same physical meaning as in the standard scenario (namely, to interpret \(\varOmega _{m0}\) as the present fraction of pressureless matter in units of the today value of the critical energy density, \(\rho _c(t_0)= 3H^2(t_0)/ 8 \pi G\)), we have now to impose that \(H_0\) represents the today value of our generalized Hubble parameter, namely that

$$\begin{aligned} H_0= \lim _{z \rightarrow 0} H(z) = \lim _{z \rightarrow 0} \left[ H^{(0)}(z)+ c_2 H^{(1)}(z) \right] . \end{aligned}$$
(35)

By applying this condition to Eqs. (33) and (34) we obtain that the four parameters \(H_0\), \(c_2\), \(\varOmega _{m0}\) and \(\varOmega _\varLambda \) cannot be independent, but must be related (to first order in \(c_2\)) by the non-standard condition

$$\begin{aligned} \varOmega _\varLambda =1 - \varOmega _{m0}+ {27\over 2} c_2 H_0^2 \varOmega _{m0}^2 ~. \end{aligned}$$
(36)

By inserting this condition into Eqs. (33), (34), and combining all zeroth-order and first-order contributions to H, we can eventually write our approximate solution as follows:

$$\begin{aligned} H(z)= & {} H^{(0)}+ c_2 H^{(1)}= H_0 \left[ \varOmega _{m0}(1+z)^3+1- \varOmega _{m0} \right] ^{1/2} \nonumber \\&\quad \times \left\{ 1+ {27\over 4} c_2 H_0^2 \varOmega _{m0}^2{1-(1+z)^6\over [\varOmega _{m0}(1+z)^3+1- \varOmega _{m0}}\right\} . \end{aligned}$$
(37)

Let us stress again that the parameters \(H_0\) and \( \varOmega _{m0}\) represent the same physical quantities as in the standard \(\varLambda \hbox {CDM}\) scenario (namely, the today value of the Hubble parameter and of the critical fraction of pressureless matter). However, they are related to the critical fraction of dark energy \(\varOmega _\varLambda \) by the non-standard relation (36), imposed (to first order in \(c_2\)) by our generalized cosmological equation (14).

Given the above (approximated) expression for H(z) it is finally useful, for the application of this paper, to explicitly evaluate also the modified friction parameter and luminosity distance of Eqs. (26) and (27). By using Eqs. (28) and (29) we easily obtain, to first order in \(c_2\):

$$\begin{aligned} \delta (z) \simeq c_2 (1+z) R^{\prime (0)} = - 9 c_2 (1+z)^3 H_0^2 \varOmega _{m0}, \end{aligned}$$
(38)

and

$$\begin{aligned} \frac{d_{L}^{GW}}{d_{L}}(z)\simeq & {} 1+ c_2 \left[ R^{(0)}(0)- R^{(0)}(z)\right] = 1+ 3 c_2 H_0^2 \varOmega _{m0}\left[ (1+z)^3 -1 \right] . \end{aligned}$$
(39)

5 Numerical results using Supernovae data

By inserting the explicit form of our modified Hubble parameter, Eq. (37), into the luminosity distance (3), we can now fit the Supernovae data to obtain a numerical estimate of the standard cosmological parameters, as well as of the new parameter \(c_2\).

The method is the same as the one already applied in [22], with the important difference that in this paper we have included the standard dark-energy contribution into the lowest-order background solution. In addition, we recall that the fit presented in [22] was performed at fixed given value of \(H_0\), using for \(H_0\) the value suggested by Planck data. Here we prefer to keep \(H_0\) as a free parameter (after all, the local value of \(H_0\) in the redshift range typical of Supernovae, , might be different from the corresponding value at much larger scales of distance). However, we will perform the fit at constant given value of \(\varOmega _{m0}\), using to this purpose the recent results of Planck + BAO data (see e.g. [33]).

We have fitted, in particular, the UNION 2 data set [23] consisting of 557 measurements of redshift/distance-modulus of Type Ia Supernovae, with redshift between 0 and 1.5. We have used for the theoretical curve the standard expression of the so-called distance modulus \(\mu (z)\), with the normalization usually adopted for SNIa and defined by

$$\begin{aligned} \mu (z,c_2,H_0)= 5 \log _{10}\left[ 10^5 d_L(z)\over 1 \mathrm{Mpc}\right] , \end{aligned}$$
(40)

where \(d_L\) is given by Eq. (3) and H(z) is given by the generalized solution (37). Given the experimental points \(\mu ^{\mathrm{obs}}(z_i) \pm \varDelta \mu (z_i)\), where \(\varDelta \mu \) is the relative error of the distance modulus for the i-th Supernova, we have performed a standard \(\chi ^2\)-analysis, with

$$\begin{aligned} \chi ^2= \sum _{i=1}^{557} \left[ \mu ^{\mathrm{obs}}(z_i)- \mu (z_i, c_2, H_0)\over \varDelta \mu (z_i)\right] ^2. \end{aligned}$$
(41)

By inserting into Eq. (37) the standard result [33] \(\varOmega _{m0} = 0.311 \pm 0.006\), and minimizing the above \(\chi ^2\) expression, we have then obtained the following best-fit values:

$$\begin{aligned} H_0= & {} (69.84 \pm 0.30) \,{\mathrm{km}} \,\mathrm{s}^{-1}\, \mathrm{Mpc}^{-1}, \nonumber \\ c_2= & {} (1.128\pm 0.670)\times 10^{-6} \,{\mathrm{km}}^{-2} \,\mathrm{s}^{2} \,{\mathrm{Mpc}}^{2}, \end{aligned}$$
(42)

with a goodness of fit \(\chi ^2/{\mathrm{d.o.f}} = 0.98\). Correspondingly, the critical fraction of dark energy turns out to be fixed, according to Eq. (36), to the value \(\varOmega _\varLambda =0.696 \pm 0.100\). Note that the above value of \(c_2\) is compatible with possible constraints arising from other cosmic data [34].

Fig. 1
figure 1

Plot of the friction parameter \(\delta (z)\) of Eq. (38) (black dashed curve), with the associated error band (grey region) controlled by the best-fit parameters of Eq. (42)

With the numerical values of Eq. (42) we are now in the position of plotting the friction term (38) and the modified GW luminosity distance (39), consistently with present observational results on the chosen scale of distances. The results are shown in Figs. 1 and 2, with the associated error bands obtained by propagating into Eqs. (38) and (39) the one sigma best-fit errors of Eq. (42). It should be mentioned that we shall henceforth evaluate all errors under the conservative assumption that the uncertainties on the given parameters are to be summed with their absolute value. For the particular case in which there are no cross-correlations among the errors of the single parameters we are certainly overestimating the resulting error bars. However, the difference turns out to be very small and, in any case, our choice of using upper bound errors does not forbid a significant comparison with observations (and with standard predictions).

As shown by the plot of \(\delta (z)\) in Fig. 1, the deviations from the standard Einstein predictions tend to increase with the distance, but they are very small in the limit \(z \rightarrow 0\), where we find

$$\begin{aligned} \delta (0)\simeq - 0.015 \pm 0.009. \end{aligned}$$
(43)

From the plot of Fig. 2 we can check the trivial identity \(d_L^{GW}(0)= d_L(0)\) at \(z=0\), consistently with the chosen normalization of the modified scale factor [see Eqs. (7) and (8)]. In the limit of large z we can see that the modification of the luminosity distance tends to grow monotonically, as expected (at least in the considered range of distances), and it is important to note that \(d_L^{GW} > d_L\) at all z. This means that, in the model of gravity we are considering, the received GW signals turn out to be weaker with respect to the value predicted by the standard Einstein theory (and thus more difficult to be detected, given the same distance and detection sensitivity).

Fig. 2
figure 2

Plot of the ratio (39) between GW and e.m. luminosity distance (black dashed curve), with the associated error band (grey region) controlled by the best-fit parameters of Eq. (42)

Finally, comparing our results with those based on the gravitational model discussed in [8, 9] (the so-called RR model [35, 36]), it may be important to note that we find a different qualitative behavior for \(\delta \) and \(d_L^{GW}\) as a function of z (actually, the opposite behavior).

In the RR model (in particular, in its minimal version where one of its parameters is set to zero) one finds \(\delta (0) \simeq 0.062\), which is definitely larger than our result (43). In the large z limit, however, the ratio \(d_L^{GW} /d_L\) tends to stabilize on the value \(\simeq 0.970\) (see e.g. Figs. 2 and 3 of [9]), which corresponds to a deviation from the standard Einstein result of about \(3 \%\). In our case, instaed, we expect a deviation of about \(8 \%\) already at \(z=1.5\) (see Fig. 2). It should be taken into account, however, that our plots and the plots of [9] are performed not only with a different functional dependence on z, but also with a different \(H_0\) parameter.

5.1 Comparison with GW170817 data

It should be checked, at this point, whether the results we have presented for \(\delta \) and \(d_L^{GW}\) may be compatible or not with other numerical estimates of such parameters, obtained from presently available data. Let us consider, in particular, the constraints on modified GW propagation following from the observation of the gravitational source GW170817 [2] and of its e.m. partner GRB170817A [3,4,5,6].

In that case the redshift parameter is fixed at \(z=0.00980\) (with an error which is negligible with respect to other uncertainties), and since this value is very small, we can use Eq. (8) expanded to first order as

$$\begin{aligned} \frac{d_{L}^{GW}}{d_{L}}(z)\simeq 1-z\delta (0) + \cdots ~. \end{aligned}$$
(44)

Following the discussion of [9], we can then compare our previous result for \(\delta (0)\) with the ones obtained from the standard siren GW170817 in two ways: by extracting \(\delta (0)\) from the measured value of the Hubble parameter or from the measured luminosity distance.

Concerning the first approach, what is needed is the local Hubble constant \(H_0^{GW}\), obtained from GW170817/GRB170817A measurements (updated to take into account the correct inclination of the source and the peculiar velocity of the host galaxy [37,38,39]), together with the local value of \(H_0\) obtained with the e.m. measurements of a set of standard candles at very small redshift [40, 41]. By using such values (reported in [9]), by expanding to first order the luminosity-redshift relation as \(d_L \simeq z/H_0+ \cdots \), by identifying, to this order, the measured e.m. distance with our \(d_L\), one then obtains from Eq. (44) the same result already presented in [9]:

$$\begin{aligned} \delta (0)\simeq \frac{H_0^{GW}-H_0}{zH_0^{GW}}\simeq 2.7^{+15.4}_{-12.8}. \end{aligned}$$
(45)

As expected, this is not a very stringent constraint, but it is well compatible with our previous estimate (43).

Conversely, we may note that using our result (43) for \(\delta (0)\), and the previously quoted [40, 41] local value \(H_0= (73.48 \pm 1.66)\) km \(\hbox {s}^{-1}\) \(\hbox {Mpc}^{-1}\) for the e.m. determined Hubble constant, we obtain from Eq. (45):

$$\begin{aligned} H_0^{GW} \simeq H_0 \left[ 1+z \delta (0)\right] \simeq (73.47 \pm 1.67) \,{\mathrm{km}} \,\,\mathrm{s}^{-1} {\mathrm{Mpc}}^{-1}. \end{aligned}$$
(46)

This agrees, within the experimental errors, with the estimate of \(H_0^{GW}\) obtained from the updated GW170817/GRB170817A measurements [37,38,39], and given by \(H_0^{GW}= 75.50^{+11.60}_{-9.60}\) km \(\hbox {s}^{-1}\) \(\hbox {Mpc}^{-1}\).

Finally, let us consider the second approach based on the direct comparison of the gravitational luminosity distance \(d_L^{GW}\) of the standard siren GW170817 [2] and the e.m. luminosity distance \(d_L^{em}\) of the host galaxy NGC4993 [42, 43]. By working as before to lowest order in z, by identifying to this order \(d_L^{em}\) with our \(d_L\), and using the expression (44), we are led to the numerical result also reported in [9],

$$\begin{aligned} \delta (0)\simeq \frac{d_L-d_L^{GW}}{zd_L} \simeq -7.8^{+9.7}_{-18.4}. \end{aligned}$$
(47)

Again, this is not a stringent limit, and is clearly compatible with our previous estimate (43).

Conversely, using the numerical value (43) for \(\delta (0)\), using for \(d_L\) the e.m. luminosity distance of the host galaxy [42, 43], \(d_L= (40.7\pm 2.4)\) Mpc, and inverting Eq. (47), we find no predicted difference for \(d_L^{GW}\), namely \(d_L^{GW}\simeq d_L\) . This is because of the very small deviations from the standard theory induced by our modified model of gravity in the small redshift limit. This results also agrees with the estimated gravitational distance of the standard siren [2], \(d_L^{GW} = 43.8^{+2.9}_{-6.9}\) Mpc.

5.2 Comparison with the Einstein cosmological equations

In the standard cosmological scenario, based on the equations of the General Relativity theory, the luminosity distance of a source at redshift z, embedded in a spatially flat FLRW geometry, is given by the well-known expression:

$$\begin{aligned} d_{L}^{GR}(z) = \frac{1+z}{H_0} \int ^z_0 \frac{dz'}{\left[ \varOmega _{m0} (1+z')^3+ (1-\varOmega _{m 0})\right] ^{1/2}} \end{aligned}$$
(48)

(where GR denotes General Relativity). It depends on two observational parameters, \(H_0\) and \(\varOmega _{m0}\), and is the same for both e.m. and GW signals (in general, it is the same for all signals propagating on the light-cone of the given FLRW metric background).

In a cosmological scenario based on the modified gravitational equations (10) we have instead two different types of luminosity distance.

We have still the luminosity distance \(d_L^{em}\) associated with light-cone propagation (and typical, for instance of e.m. signals): it is defined by Eq. (3), where, however, the function H(z) is now a solution of the generalized gravitational equations and is thus different, in general, from the standard expression of H(z) appearing in Eq. (48). In addition, we have a modified luminosity distance \(d_L^{GW}\) associated with the propagation of GW signals (distorted by the presence of a non standard friction coefficient), and related to \(d_L^{em}\) by Eq. (27). The aim of this subsection is to compare the generalized expressions of both \(d_L^{em}\) and \(d_L^{GW}\) with the standard luminosity distance \(d_L^{GR}\) of Eq. (48).

To this purpose, and for a more consistent comparison, let us first evaluate the parameters controlling the standard expression of \(d_L^{GR}\) in the same way used for our model of modified gravity: namely, in the same redshift range, by fitting the same experimental data as before, the UNION 2 data set [23], and at the same fixed value of the parameter \(\varOmega _{m0}\), namely [33] \(\varOmega _{m0}=0.311 \pm 0.006\). The standard best-fit procedure, starting with the expression (40) of the distance modulus written for the luminosity distance (48), then gives:

$$\begin{aligned} H_0=(69.53 \pm 0.24) \,{\mathrm{km}} \,\,\mathrm{s}^{-1}{\mathrm{Mpc}}^{-1}, \end{aligned}$$
(49)

with \(\chi ^2/{\mathrm{d.o.f}}=0.98\). This completely fixes (with the relative uncertainty) the standard model luminosity (48) as a function of z.

We are now in the position of comparing the above expression of \(d_L^{GR}\) with the e.m. distance \(d_L^{em}\) of our model—given by Eq. (3) in terms of the generalized H(z) of Eq. (37)—and with the GW distance \(d_L^{GW}\), related to \(d_L^{em}\) by Eq. (39). Of course, the modified luminosity distances are to be computed (as emphasized also in [9]) with the parameters determined by the associated fitting procedures [i.e. the ones reported in Eq. (42)], while the standard distance (48) refers to the parameter (49). The plots of the corresponding fractional corrections,

$$\begin{aligned} \varDelta _{em}(z)= { d_L^{em}- d_L^{GR}\over d_L^{GR}}, ~~~~~~~~~ \varDelta _{GW}(z)= { d_L^{GW}- d_L^{GR}\over d_L^{GR}}, \end{aligned}$$
(50)

together with the associated error bars, are presented in Figs. 3 and 4.

Fig. 3
figure 3

Plot of \(\varDelta _{em}(z)\) (black dashed curve) with the relative error band (grey region), controlled by the best fit parameters of Eq. (49) for \(d_L^{GR}\) and of Eq. (42) for \(d_L^{em}\). In both cases we have set \(\varOmega _{m0}=0.311 \pm 0.006\). The use of different parameters for the two luminosity distances leads to \(\varDelta _{em}(0)\not =0\) (see the text)

Fig. 4
figure 4

Plot of \(\varDelta _{GW}(z)\) (black dashed curve) with the relative error band (grey region), controlled by the best fit parameters of Eq. (49) for \(d_L^{GR}\) and of Eq. (42) for \(d_L^{GW}\). In both cases we have set \(\varOmega _{m0}=0.311 \pm 0.006\). The use of different parameters for the two luminosity distances leads to \(\varDelta _{GW}(0)\not =0\) (see the text)

We can easily check that the behavior of the plotted functions is in agreement (at all z) with the property \(d_L^{GW}>d_L^{em}\), already emphasized by Fig. 2 and due, as previously stressed, to a faster decrease with the distance of the GW signals with respect to e.m. signals. We can also see, from Fig. 3, that \(\varDelta _{em}\) tends to be positive at large enough z, and this implies that also the e.m. signals, in our modified-gravity scenario, tend to be weaker than in the standard scenario, at large enough distances. In addition, by comparing the plots of Figs. 3 and 4, we find that \(|\varDelta _{GW}|>|\varDelta _{em}|\), at any z (at least in the considered redshift range). This seems to suggest that detecting GW signals might represent a more efficient tool (than detecting e.m. signals) to test generalized model of gravity on the given scale of cosmic distances.

Finally, it is important to recall that in the limit \(z \rightarrow 0\) the two modified distances tend (by definition) to coincide, i.e. \(d_L^{em}(0)=d_L^{GW}(0)\) (as previously stressed, and as illustrated by Fig. 2). This explains why, for \(z=0\), the two fractional corrections of Figs. 3 and 4 go to the same (non-vanishing) value

$$\begin{aligned} \varDelta _{em}(0)=\varDelta _{GW}(0)\simeq -0.004 \pm 0.008. \end{aligned}$$
(51)

Note that the difference \(\varDelta (0)\) (still compatible with zero, to one sigma, in the context of our numerical analysis), in principle can always be expected to be non-vanishing, because we have plotted the standard and the modified luminosity distances using, for each curve, their peculiar best-fit parameters: that of Eq. (49) for the Einstein gravitational equations, and those of Eq. (42) for the modified equations. Let us note, in this respect, that the importance of using for each model its own set of estimated parameters has been repeatedly stressed also in [9]. Only in this way we obtain for \(\varDelta _{em}(z)\) and \(\varDelta _{GW}(z)\) a relevant quantity to be compared with observations [9]. This implies, however, that the results we are presenting cannot be reliably extrapolated outside the range of z that we have considered to fit the data, and to obtain the corresponding estimates.

6 Conclusion

In this paper we have studied the possible differences arising between the propagation of GW signals and e.m. signals—as well as the differences between the related luminosity distances—in the context of a generalized theory of gravity with quadratic curvature corrections.

Considering in particular the range of distances typical of Supernovae (), we have found that the deviations from the standard expression of the luminosity distance is larger for GW signals than for e.m. signals (at least, in the range of redshift we have considered). This seems to confirm previous suggestions [9, 10] that GW detectors might play a crucial (future) role in discriminating between different models of gravity on cosmic scales of distances. For the case discussed in this paper, however, the amplitude of GW signals received by a distant source, after their modified propagation, tends to be weaker than in the context of the standard Einstein theory.

It should be stressed, finally, that the main purpose of this paper is not only to propose possible observational constraints on the particular model of gravity we have considered, but also to present and discuss a general approach, possibly useful for future applications to other (more general and motivated) non-standard theories of gravity.

The particular choice of the model discussed in this paper has been motivated by previous applications of the same model in an astrophysical context to study deviations from the standard cosmological dynamics [44], as well as the detailed mechanisms of neutron star formations [45, 46]. The validity of the analysis presented in this paper is limited to the redshift range . We are planning to discuss the compatibility of this model with other available cosmological data, and in a larger redshift range, in a future paper.