  We present a detailed study of the methodology for correlating `dark sirens'
(compact binaries coalescences without electromagnetic counterpart) with galaxy
catalogs. We propose several improvements on the current state of the art, and
we apply them to the published LIGO/Virgo gravitational wave (GW) detections,
performing a detailed study of several sources of systematic errors that, with
the expected increase in statistics, will eventually become the dominant
limitation. We provide a measurement of $H_0$ from dark sirens alone, finding
as the best result $H_0=75^{+25}_{-22}\,\,{\rm km}\, {\rm s}^{-1}\, {\rm
Mpc}^{-1}$ ($68\%$ c.l., for a flat prior in the range $[20,140] \,\,{\rm km}\,
{\rm s}^{-1}\, {\rm Mpc}^{-1}$) which is, currently, the most stringent
constraint obtained using only dark sirens. Combining dark sirens with the
counterpart for GW170817 we find $H_0=70^{+11}_{-7} \,{\rm km}\, {\rm s}^{-1}\,
{\rm Mpc}^{-1}$. We also study modified GW propagation, which is a smoking gun
of dark energy and modifications of gravity at cosmological scales, and we show
that current observations of dark sirens already start to provide interesting
limits. From dark sirens alone, our best result for the parameter $\Xi_0$ that
measures deviations from GR (with $\Xi_0=1$ in GR) is
$\Xi_0=1.88^{+3.83}_{-1.10}$. We finally discuss limits on modified GW
propagation under the tentative identification of the flare ZTF19abanrhr as the
electromagnetic counterpart of the binary black hole coalescence GW190521, in
which case our most stringent result is $\Xi_0=1.6^{+1.0}_{-0.6}$.
  We release the publicly available code $\tt{DarkSirensStat}$, which is
available under open source license at
https://github.com/CosmoStatGW/DarkSirensStat.
