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_0$) using broad-band afterglow emission and relativistic jet motion from the Very Long Baseline Interferometry and Hubble Space Telescope images of GW170817. Compared to previous attempts, we combine these messengers with GW in a simultaneous Bayesian fit. We probe the $H_0$ 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\%$ error ($77^{+21}_{-10}$ km/s/Mpc, medians, 16th-84th percentiles), because of the degeneracy between viewing angle ($\theta_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}$ km/s/Mpc, 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}$ km/s/Mpc. Adding the centroid motion in the analysis efficiently breaks the $d_L-\theta_v$ degeneracy and overcome the late-time deviations, giving $H_0 = 69.0^{+4.4}_{-4.3}$ km/s/Mpc (in agreement with Planck and SH0ES measurements) and $\theta_v = 18.2^{+1.2}_{-1.5}$ deg. 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.