  Gravitational wave (GW) standard sirens are well-established probes with
which one can measure cosmological parameters, and are complementary to other
probes like the cosmic microwave background or supernovae standard candles.
Here we focus on dark GW sirens, specifically binary black holes (BBHs) for
which there is only GW data. Our approach relies on the assumption of a source
frame mass model for the BBH distribution, and we consider four models that are
representative of the BBH population observed so far. In addition to inferring
cosmological and mass model parameters, we use dark sirens to test modified
gravity theories. These theories often predict different GW propagation
equations on cosmological scales, leading to a different GW luminosity distance
which in some cases can be parametrized by variables $\Xi_0$ and $n$. General
relativity (GR) corresponds to $\Xi_0= 1$. We perform a joint estimate of the
population parameters governing mass, redshift, the variables characterizing
the cosmology, and the modified GW luminosity distance. We use data from the
third LVK observation run and find - for the four mass models and for three SNR
cuts - that GR is consistently the preferred model to describe all observed BBH
GW signals. Furthermore, all modified gravity parameters have posteriors that
are compatible with the values predicted by GR at the 90% confidence interval.
We show that there are strong correlations between cosmological, astrophysical
and modified gravity parameters. If GR is the correct theory of gravity, and
assuming narrow priors on the cosmological parameters, we forecast an
uncertainty of the modified gravity parameter $\Xi_0$ of 51% with $\sim 90$
detections at O4-like sensitivities, and $\Xi_0$ of 20% with an additional
$\sim 400$ detections at O5-like sensitivity. We also consider how these
forecasts depend on the current uncertainties of the BBH population.
