  Gravitational waves (GWs) from binary neutron stars (BNSs) are expected to be
accompanied by electromagnetic (EM) emissions, which help to identify the host
galaxy. Since GW events directly measure their luminosity distances, joint
GW-EM observations from BNSs help to study cosmology, particularly the Hubble
constant, unaffected by cosmic distance ladder systematics. However, detecting
EM counterparts from BNS mergers is not always possible. Additionally, the
tidal deformations of BNS components offer insights into the neutron star (NS)
equation of state (EoS). In such cases, the tidal parameters of NSs, combined
with the knowledge of the NS EoS, can break the degeneracy between mass
parameters and redshift, allowing for the inference of the Hubble constant.
Several efforts have aimed to infer the Hubble constant using the tidal
parameters of BNSs, without EM counterparts, termed dark BNSs. Moreover, some
studies have focused on the joint estimation of population and NS EoS for
unbiased NS EoS estimation. However, none of the works consistently combined
the uncertainties of population, cosmology, and NS EoS within a Bayesian
framework. In this study, we propose a novel Bayesian analysis to jointly
constrain the NS EoS, population, and cosmological parameters using a
population of dark BNSs detected through GW observations. This method can well
constrain the Hubble constant with as few as $5$ BNS observations using
current-generation detectors. This level of precision is unattainable without
incorporating the NS EoS, especially when observing BNS mergers without EM
counterpart information.
