Potential biases and prospects for the Hubble constant estimation via electromagnetic and gravitational-wave joint analyses

GW170817 is a binary neutron star merger that exhibited a gravitational wave (GW) and a gamma-ray burst, followed by an afterglow. In this work, we estimate the Hubble constant (H₀) using broad-band afterglow emission and relativistic jet motion from the Very Long Baseline Interferometry and HST images of GW170817. Compared to previous attempts, we combine these messengers with GW in a simultaneous Bayesian fit. We probe the H₀ measurement robustness depending on the data set used, the assumed jet model, the possible presence of a late time flux excess. Using the sole GW leads to a 20 per cent error ($77^{+21}_{-10}$ $\rm km\, s^{-1}\, Mpc^{-1}$, medians, 16th-84th percentiles), because of the degeneracy between viewing angle (θ_v) and luminosity distance (d_L). The latter is reduced by the inclusion in the fit of the afterglow light curve, leading to $H_0=96^{+13}_{-10}$ $\rm km\, s^{-1}\, Mpc^{-1}$, a large value, caused by the fit preference for high viewing angles due to the possible presence of a late-time excess in the afterglow flux. Accounting for the latter by including a constant flux component at late times brings $H_0=78.5^{+7.9}_{-6.4}$$\rm km\, s^{-1}\, Mpc^{-1}$. Adding the centroid motion in the analysis efficiently breaks, the d_L - θ_v degeneracy and overcome the late-time deviations, giving $H_0 = 69.0^{+4.4}_{-4.3}$ $\rm km\, s^{-1}\, Mpc^{-1}$ (in agreement with Planck and SH0ES measurements) and $\theta _{\rm v} = 18.2^{+1.2}_{-1.5}$°. This is valid regardless of the jet structure assumption. Our simulations show that for next GW runs radio observations are expected to provide at most few other similar events.