  Gravitational-wave (GW) radiation from a coalescing compact binary is a
standard siren as the luminosity distance of each event can be directly
measured from the amplitude of the signal. One possibility to constrain
cosmology using the GW siren is to perform statistical inference on a
population of binary black hole (BBH) events. In essence, this statistical
method can be viewed as follows. We can modify the shape of the distribution of
observed BBH events by changing cosmological parameters until it eventually
matches the distribution constructed from an astrophysical population model,
thereby allowing us to determine the cosmological parameters. In this work, we
derive the Cram\'er-Rao bound for both cosmological parameters and those
governing the astrophysical population model from this statistical dark siren
method by examining the Fisher information contained in the event distribution.
Our study provides analytical insights and enables fast yet accurate
estimations of the statistical accuracy of dark siren cosmology. Furthermore,
we consider the bias in cosmology due to unmodeled substructures in the merger
rate and the mass distribution. We find a $1\%$ deviation in the astrophysical
model can lead to a more than $1\%$ error in the Hubble constant. This could
limit the accuracy of dark siren cosmology when there are more than $10^4$ BBH
events detected.
