All these methods have been first proposed to measure \(H_0\)  in the context of GR. However, they can be extended to the context of modified gravity. In this case, the GW observations measure \(d_L^{\,\rm gw}(z)\)  and, in the parametrization (REF ), the parameter space is therefore enlarged to include \((\Xi _0,n)\) . The most general approach is to perform a joint inference over all cosmological and astrophysical parameters; a simpler approach, for what concerns the cosmological parameters, is to assume a value (or a narrow prior) for \(H_0\) , such as that given by Planck  [1]}, and focus on the accuracy that can be obtained on \(\Xi _0\) .The parameter \(n\)  is less important since it only determines the precise shape of the function that interpolates between \(d_L^{\,\rm gw}(z)/d_L^{\,\rm em}(z)=1\)  at \(z=0\)  and \(d_L^{\,\rm gw}(z)/d_L^{\,\rm em}(z)\simeq \Xi _0\)  at large \(z\) . One could then fix it to typical values suggested by explicit model, such as the value \(n\simeq 1.91\)  obtained in the RT nonlocal model. Alternatively, one can marginalize also over \(n\) , which results in limits on \(\Xi _0\)  larger by about a factor of 2 [2]}. In this paper we will include also \(n\)  among the parameters determined by the Bayesian inference. Several results have already been obtained with current data.
Using GW170817 as a standard siren with counterpart,
a first limit was obtained in [3]}, using the electromagnetic luminosity distance of the galaxy hosting the counterpart, obtained from surface brightness fluctuations (see also [4]}, and the discussion in App. , for other analysis using GW170817). Because of the very small redshift
of GW170817, this is really a limit on \(\delta (z=0)\) , independently of the parametrization used, and the result from [3]} is
\(\delta (0)=-7.8^{+9.7}_{-18.4}\, .\) 
