Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues

In the absence of numerous gravitational-wave detections with confirmed electromagnetic counterparts, the "dark siren" method has emerged as a leading technique of gravitational-wave cosmology. The method allows redshift information of such events to be inferred statistically from a catalogue of potential host galaxies. Due to selection effects, dark siren analyses necessarily depend on the mass distribution of compact objects and the evolution of their merger rate with redshift. Informative priors on these quantities will impact the inferred posterior constraints on the Hubble constant ($H_0$). It is thus crucial to vary these unknown distributions during an $H_0$ inference. This was not possible in earlier analyses due to the high computational cost, restricting them to either excluding galaxy catalogue information, or fixing the gravitational-wave population mass distribution and risking introducing bias to the $H_0$ measurement. This paper introduces a significantly enhanced version of the Python package GWCOSMO, which allows joint estimation of cosmological and compact binary population parameters. This thereby ensures the analysis is now robust to a major source of potential bias. The gravitational-wave events from the Third Gravitational-Wave Transient Catalogue are reanalysed with the GLADE+ galaxy catalogue, and an updated, more reliable measurement of $H_0=69^{+12}_{-7}$ km s$^{-1}$ Mpc$^{-1}$ is found (maximum a posteriori probability and 68% highest density interval). This improved method will enable cosmological analyses with future gravitational-wave detections to make full use of the information available (both from galaxy catalogues and the compact binary population itself), leading to promising new independent bounds on the Hubble constant.