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.


INTRODUCTION
The ΛCDM model is the currently adopted standard model of cosmology.Great effort has been put into the estimation of one of its key parameters, the Hubble constant  0 , the current expansion rate of the Universe.The ΛCDM model calibrated with data from the Planck mission, that is, from early-Universe physics, predicts the Hubble constant to 1% precision: 67.4 ± 0.5 km s −1 Mpc −1 (Planck Collaboration et al. 2020) (we quote medians and 68% credible intervals).However,  0 can also be empirically measured locally ( < 1), in the late-time Universe.The latter kind of measurements, such as from SH0ES (Supernovae  0 for the Equation of State, Riess et al. 2019) and H0liCOW ( 0 lenses in COSMOGRAIL's Wellspring, Wong et al. 2019), favour larger values of  0 : 74.0 ± 1.4 km s −1 Mpc −1 and 73.3 ± 1.8 km s −1 Mpc −1 , respectively.Thus, the early-Universe data seem to be consistently predicting a low value of  0 , while the late-time Universe data a higher one, leading to a more than 3 discrepancy (see for an extensive discussion Verde et al. 2019).
A way to solve this discrepancy is to measure the Hubble constant through an independent method, using, for example, gravitational ★ E-mail: giulia.gianfagna@inaf.itwaves (GWs), where the distance is directly estimated from fitting the waveform, relying only on the general theory of relativity.This estimation does not depend on cosmic distance ladders.GWs can determine the Hubble constant if the redshift is provided by an electromagnetic (EM) counterpart, for example by a kilonova (Taylor et al. 2012;Feeney et al. 2021).This is the so-called standard sirens method (Holz & Hughes 2005;Nissanke et al. 2010).Even when a unique counterpart cannot be identified, the redshifts of all the potential host candidates can be incorporated in the analysis, when the localization volume is sufficiently small.This is not as constraining as the first scenario, but it is still informative, once many detections are available.In this case, more than 50 binary neutron stars are needed to reach a 6%  0 measurement (Chen et al. 2018).The same holds also for GWs emitted by binary black holes, even if the localization volumes are usually much larger than for binary neutron stars.In this case ∼ 500 events are needed to reach a precision < 7% on  0 (Chen et al. 2022;Bom & Palmese 2023).
The main problem of the standard sirens method is the degeneracy between the luminosity distance and inclination (the angle between the total angular momentum and the line of sight) estimated from GWs.They are measured from the amplitude of the two GW polarizations.At small inclinations, the cross and plus polarizations have nearly the same amplitude, but the larger the inclination, the more they decrease and start to differ (Usman et al. 2019).This means that the GW signal is strongest at small inclinations (face-on or face-off), but, in these cases, we cannot measure distance and inclination separately.Therefore, associated EM observations can lead to a tighter measurement of  0 by providing additional constraints on the inclination.
The EM information on the inclination derived from the afterglow and the relativistic jet motion of GW170817 allow us to improve the Hubble constant measurement for the reason stated above (see also Bulla et al. 2022, for a review).The common practice is to use the GW analysis results (posterior) for inclination and luminosity distance, and apply these as a priori information on the inclination obtained by fitting the EM data sets (or the other way around).This can be done using Bayesian analysis.The results retrieved in this way run from low values such as  0 = 66.2 +4.4  −4.2 km s −1 Mpc −1 , from Dietrich et al. (2020), who fit the kilonova emission and the jet centroid motion, and  0 = 69.5 ± 4 km s −1 Mpc −1 , from Wang & Giannios (2021), who use the afterglow emission, to high values such as  0 = 75.5 +11.6  −9.6 km s −1 Mpc −1 , by Guidorzi et al. (2017), who fit the afterglow up to 40 days from the merger.Wang et al. (2023) estimate  0 = 71.80+4.15  −4.07 km s −1 Mpc −1 , modelling the jet with hydrodynamic simulations, including also a sub-relativistic kilonova outflow.Palmese et al. (2023) use the same model as Wang et al. (2023), fitting the afterglow, but including a priori information on the Lorentz factor from the jet centroid motion.They find  0 = 75.46+5.34  −5.39 km s −1 Mpc −1 .Hotokezaka et al. (2018) fit the afterglow and the jet centroid motion, finding  0 = 68.9+4.7 −4.6 km s −1 Mpc −1 .In general, the smaller is the viewing angle, the higher is the luminosity distance (because of their degeneracy), the lower is  0 .
In general, the best approach to estimate  0 is to include all the available information about the analyzed event, so, in case of GW170817, the GW, the afterglow light curve, the jet centroid motion and the kilonova.However, even analyses including only some of these are useful, especially to quantify possible systematics, that will count as reference for future events.In the case of GW170817, the afterglow light curve and the centroid motion tend to give contrasting results, as will be shown in this work, mainly because of the late time data points in the light curve, which seem to be showing a flux excess.Therefore, the use of the complete available information is even more important.
At present, the aforementioned EM-informed  0 measurements are a factor 2 more precise than the first standard-siren measurement for GW170817 that fitted GW data only (Abbott et al. 2017b).For this reason, this method is very compelling.However, there are potential systematics that should be addressed.An open issue, as outlined in Nakar & Piran (2021), is the sensitivity of EM-derived parameters, as the inclination, on the assumed jet structure.A related problem is the presence of deviations from the assumed model due to a possible flux excess at late times.In this work, we assess these issues with a comprehensive approach.We estimate the Hubble constant exploiting the GW, the broad-band afterglow and the centroid motion of the relativistic jet of the GW170817 event.We test the sensitivity of the results on the jet structure, and check for potential biases, both due to the jet model assumption and to the possible presence of an excess at late times in the afterglow.We fit the GW, the afterglow and the centroid motion data sets simultaneously using a Bayesian approach (Gianfagna et al. 2023) and compare the results obtained fitting only afterglow and GW data, and then including also the centroid motion of the relativistic jet.We focus on the degeneracy between the viewing angle and the luminosity distance.As already stated above, because of this degeneracy, the standard sirens method at present cannot give a Hubble constant estimation at the Planck or SH0ES level of precision.Here we study how we can break this degeneracy including different types of EM messengers.We also test the robustness of the derived  0 by implementing the different jet models on real data and simulations.
In Section 2 we present the data sets used in this work, analysed following the method presented in Section 3. In Section 4 we show the results that we obtained, both for the energetics, microphysics and geometry of the event, and for the Hubble constant.Finally in Section 5 we summarize our conclusions.

DATA
This work used three data sets pertaining to the GW170817 event and analyzes them simultaneously.These are the broad-band afterglow emission, the centroid position of the jet as a function of time, and the GW strain timeseries.
We also include the centroid motion of the relativistic jet, visible in optical and radio images (Mooley et al. 2018;Ghirlanda et al. 2019;Mooley et al. 2022).For this analysis, we use the positions and uncertainties of the data points from VLBI (Very Long Baseline Interferometry) at 75, 206, 230 days reported in Mooley et al. (2018); Ghirlanda et al. (2019), and from HST at 8 days (Mooley et al. 2022).For the latter we use the positions (right ascension, RA, and declination, Dec) and their statistical uncertainties, to which we add From bottom to top, red points refer to the X-ray observations by Chandra and XMM at 5 keV, orange ones to observations by HST, F606W filter, in the optical band, and blue ones to observations in the radio band from VLA (Very Large Array) at 3 GHz.The continuous and dotted lines represent the fit of the GW, broad-band afterglow, and centroid motion (GW+AG+C) data and of the GW and afterglow (GW+AG) data, respectively.For sake of simplicity, the fit for the radio band is plotted only for the observations at 3 GHz, but it is not limited to this single frequency.Bottom panel.Centroid motion of the relativistic jet from HST and VLBI images at 8 (negative RAs), 75, 206 and 230 days (Mooley et al. 2018;Ghirlanda et al. 2019).The blue dots represents the positions predicted by the model, the blue contours represent the 68% probability region.in quadrature the two systematic uncertainty contributions to take into account the different reference frame of the optical and radio images (as in Mooley et al. 2022).
The GW data of GW170817 are publicly available at the GW Open Science Center1 (Abbott et al. 2021).We use the cleaned version of the strain data, where the glitch discussed in Abbott et al. (2017a) has been removed.

JOINT ANALYSIS OF ELECTROMAGNETIC AND GRAVITATIONAL-WAVE DATA
We use the Bayesian inference to process the GW and EM data.The two domains can be joined in one analysis as the models describing the two emissions (EM and GW) have parameters in common, namely the viewing angle and the luminosity distance.We perform three fits, one including only GWs, one folding GWs and the afterglow emission, and one including also the jet centroid motion data set.In this Section we describe the GW and the EM models, along with the joint fit method (see also Gianfagna et al. 2023).

Afterglow light curve and centroid motion
We model the broad-band afterglow light curve using the python package afterglowpy (Ryan et al. 2020).The observer frame flux of synchrotron radiation is estimated for various jet geometries.In this work we we use a Gaussian structured jet model, where the energy drops according to  () =  0 exp(− 2 /2 2  ), up to a truncating angle   . 0 ,   and   are free parameters in the fit, representing the on-axis isotropic equivalent kinetic energy of the blast wave, the jet opening angle and the jet total angular width.We also use a power law jet model, where the energy is given by , where  is the power law index.The electrons are shock-accelerated and emit synchrotron radiation, with an energy distribution given by a power law with slope −, the fraction of their post-shock internal energy is   , while the fraction of post-shock internal energy in the magnetic field is denoted by   .Furthermore, the circumburst medium number density  0 , the viewing angle   , between the jet axis and the line of sight, and the luminosity distance   are also free parameters.The participation fraction   is fixed to 1.0.
In order to model also the jet centroid motion, we use an extended version of afterglowpy (Ryan et al. 2023;van Eerten & Ryan 2023), where the afterglow image centroid position and sizes can be estimated.The imaging plane is perpendicular to the line of sight of the observer, and the centroid position and sizes are computed as an intensity-weighted quantity.The outputs of the model that we use in this work are the centroid position in the sky (RA and Dec) and the flux expected at each particular time.At 8 days, in the optical, only the kilonova emission is visible, and not the afterglow.For this reason, we place a 5 upper limit for the optical flux of 4×10 −5 mJy.The parameters are the same as above, with an extra three parameters: RA 0 , Dec 0 , which represent the jet origin in the sky image, and the position angle PA, which is the orientation of the jet direction in the image.
The prior probability distributions are reported in Table 1.The prior for the viewing angle   is isotropic, meaning a sinusoidal distribution from 0°to 90°(uniform in cosine).For the luminosity distance, we use a uniform-in-volume prior (∝  2  ) from 1 to 75 Mpc, which distributes mergers uniformly throughout a Euclidean universe.

Gravitational waves
We use the IMRPhenomPv2_NRTidal waveform approximant (Hannam et al. 2014;Dietrich et al. 2017Dietrich et al. , 2019a,b) ,b) to model GWs from binary neutron star mergers.The intrinsic parameters (the source physical parameters that shape the emitted signal) used by this model refer to the masses, the spins and the tidal deformabilities of the two Table 1.Prior probability distributions for the shared, the EM and GW fitted parameters.

Parameter
Prior functional form Bounds neutron stars.The two component masses  1 and  2 , for which we follow the common convention  1 ≥  2 , will be quoted as the chirp mass (Finn & Chernoff 1993;Blanchet et al. 1995;Cutler & Flanagan 1994;Poisson & Will 1995), and the mass ratio  =  2 / 1 ≤ 1.The components of the dimensionless spin angular momenta of each neutron star, a 1 and a 2 , constitute six additional parameters, which are:  1 and  2 , the dimensionless spin magnitudes;  1 and  2 , the tilt angles between the spins and the orbital angular momentum;  1,2 , the azimuthal angle separating the spin vectors; and    , the opening angle of the cone of precession of the orbital angular momentum about the system's total angular momentum.The tidal deformability of each star is described by the dimensionless parameters Λ 1 and Λ 2 .
The extrinsic parameters (that further shape the observed GW signal) in this model are the RA and Dec of the source (i.e., its sky position), the luminosity distance   of the source, the inclination angle  JN between the total angular momentum of the binary and the line of sight from the source to the observer, the polarization angle , and the phase and time of coalescence.In this work, we fix the sky-position (RA and Dec) of the source to the one of AT 2017gfo (Abbott et al. 2017c).The GW RA and Dec correspond to the RA 0 , Dec 0 parameters of the centroid motion model (Section 3.1.1);however, the precision on RA and Dec from the GW data does not reach the mas level, as instead do RA 0 and Dec 0 , so the analysis does not benefit from promoting the RA and Dec of the GW and the EM models to common, free parameters.Moreover, we do not report the time of coalescence in the results, as this is of little interest in the context of our study, and we marginalize over the phase of coalescence.The latter marginalization is justified by the small spin magnitudes (see Table 1), and hence the negligible precession effects Romero-Shaw et al. (2020).
Given that the GRB jet develops around the total angular momentum, the inclination angle  JN and the viewing angle   introduced in Sec.3.1.1are essentially the same quantity, and thus a common parameter of the GW and EM domains.More precisely, in the case of GW170817, the two angles are supplementary (see Eq.( 1) in Gianfagna et al. 2023).The other parameter shared by the GW and EM domains is the luminosity distance.This implies that there are 23 parameters when addressing the GW and EM domains in our approach.
The priors for the intrinsic and extrinsic GW parameters are set as in Romero-Shaw et al. (2020), in the case of the "Low Spin" analysis, see Table 1.

Joint fit
We use Bayesian inference to analyse jointly the data from the GW and EM domains, which we denote as   and    , respectively (the same methodology is presented in Gianfagna et al. 2023).The three main components are: a prior distribution, which models the available knowledge about a given parameter before data collection in a statistical distribution; the likelihood function, which encloses the information about the parameter from observed data; the posterior distribution, which combines the prior distribution and the likelihood function using the Bayes theorem.Thus, the multi-dimensional posterior probability distribution for our set of parameters ì  is: where L EM+GW ( EM ,  GW | ì ) is the likelihood function that folds the EM and GW domains, ( ì ) is the multi-dimensional prior probability distribution for our parameters, and Z ì  is the Bayesian evidence.This is obtained by marginalizing the joint likelihood over the GRB and GW parameters: When the two data sets are independent, as is the case here, the likelihood L EM+GW (see also Fan et al. 2014;Biscoveanu et al. 2020) is simply given by the product of the EM and GW likelihoods The EM and GW likelihoods are both Normal distributions.The GW likelihood function is defined in, e.g., Finn (1992); Romano & Cornish (2017); Romero-Shaw et al. (2020); in this likelihood, both the data and the model are expressed in the frequency domain.
In the EM case, when only the afterglow is folded with the GW data (GW+AG), the likelihood function is proportional to exp(− 2 /2), where  2 is given by the comparison between the expected flux and the entire broadband set of afterglow data.
In the case of the fit of the afterglow, centroid motion and GW strain (GW+AG+C), we assume the afterglow and the centroid motion data sets to be independent, being taken with different telescopes and at different times.The centroid data set includes the positions (RA and Dec) at each time and their respective fluxes, we take the data point at 8 days as the position of the merger.We assume the likelihood function to be a multivariate Normal distribution, where the expected centroid positions and fluxes from the model are compared with the three offset positions (RA and Dec) and the corresponding flux measurements.Moreover, we assume the covariance matrix to be diagonal (see Ryan et al. 2023 for more details).We place the centre of the centroid motion reference system at the positions corresponding to the observations at 75 days, as in Mooley et al. (2018); Ghirlanda et al. (2019); Mooley et al. (2022).
We use the Bayesian inference library bilby (Ashton et al. 2019;Smith et al. 2020) and the dynamic nested sampling package dynesty (Speagle 2020) to simultaneously fit the EM and GW data sets.We use 2000 live points and multiple bounding ellipsoids as bounding strategy.The corner plots are created with the corner package (Foreman-Mackey 2016).

Hubble constant estimation
At small redshifts, as in the GW170817 case, the luminosity distance does not depend on the cosmological model, so the Hubble constant can estimated from where   is the local "Hubble flow" velocity, in this case at the position of GW170817, and   is the luminosity distance to the source.We follow the same procedure as Abbott et al. (2017b), assuming a Normal distribution for   = 3017 ± 166 km s −1 .

RESULTS AND DISCUSSION
We assume a Gaussian jet profile throughout the work, with the exception of Sec.4.2, where we assume a power law profile, in order to test the sensitivity of the results on the jet model.
The parameter medians and 16th-84th percentiles are collected in Table 2.The second column reports the results of the GW-only fit, while the third and fourth column refer to the fit including the broadband afterglow and GW (GW+AG), and to the complete fit that also includes the centroid (GW+AG+C), respectively.
The results from the GW fit are in agreement with previous works Abbott et al. (2017bAbbott et al. ( , 2019)); Romero-Shaw et al. (2020), the  0 value that we retrieve from the GW-only fit is  0 = 77 +21 −10 km s −1 Mpc −1 (median, 16th-84th percentiles), see Fig. 2, bottom panel, and Fig. 3.As we already pointed out above, one of the main sources of uncertainty in the GW measurement of the inclination and of the distance (and  0 ) is due to their degeneracy, see the light blue contours in Fig. 2, top panel.This means that it is hard to distinguish whether a source is further away with the binary orbit facing Earth (face-on or face-off), or closer but highly inclined (edge-on, Usman et al. 2019).If we assume to have inclinations from 0 to 90 deg (like in our case),   is a decreasing function of the inclination (viewing angle,   ).Another independent messenger is needed to break this degeneracy, which, in this case, comes from the afterglow.
The afterglow light curve alone, however, is not enough to efficiently break this degeneracy.Including it in the fit, only helps in shrinking the degeneracy region, see green filled contours in Fig. 2, top panel, the uncertainty on the viewing angle is reduced by a factor ∼ 3 (from    = 146 +16 −18 deg to 129.9 +5.1 −5.4 deg), the one on the distance by a factor ∼ 2 (from   = 39.2 +5.4  −8.6 Mpc to 31.3 +3.0 −3.6 Mpc).However, this is not an accurate measurement, in fact the medians are on the high-  -low-  end of the GW 1 region, leading to quite a low distance (and large viewing angle,   = 50.1 +5.1 −5.4 deg), which is however within 3 from the generally accepted value of ∼ 40 Mpc.Our  0 value from the GW+AG fit is quite high: we retrieve  0 = 96 +13 −10 km s −1 Mpc −1 (median, 16th-84th percentiles), see the green filled contours in Fig. 2 −9.6 km s −1 Mpc −1 , assuming a Top Hat jet and fitting the afterglow data up to 40 days from the merger.The latter is the reason why the  0 uncertainties are larger with respect to more recent works, their   posterior distribution peaks at ∼30 deg.However, the Top Hat jet is not the best choice for GW170817 light curve, as it cannot reproduce the slope before the peak.
In Section 4.3 we account for the possible late time excess in the GW+AG fit with a constant flux component.In this case, our results are in agreement with the aforementioned works, see the complete analysis in Section 4.3.It is also interesting to note that the similar results of these works are obtained using different jet structures.We will come back to this point and explore how  0 changes depending on the jet structure in Section 4.2.
From these results, we find that limiting the analysis to GW+AG domains could be subject to possible systematics in the  0 determination due to the detections at late times in the afterglow light curve.This adds up to the degeneracy between   and   , proper of a Gaussian modelling of the jet, that is evident in the left panel of Fig. 4.Here we show the marginalised, 2D posterior probability distributions for the jet opening angle   and the viewing angle   in the cases of the the joint GW+AG fit (in red contours).
For the reasons stated above, in order to break both the   −   GW degeneracy and the   −   EM one, we have to include not only the afterglow light curve, but also the centroid motion in the analysis.We note that the sole centroid motion is not enough to break the   −   degeneracy, being itself subjected to some level of degeneracy between these two parameters, see Appendix A for a more detailed discussion.
The results for the GW+AG+C fit, using both the afterglow and the centroid, are written in Table 2, fourth column.This fit not only shifts the viewing angle to lower values, but also shrinks further the degeneracy between   and   , see left panel of Fig. 4, blue contours.This happens because the relativistic jet motion strongly constrain the viewing angle.In Fig. 1, bottom panel, the jet positions are well reconstructed by the model, within 1, while in top panel, we see that the GW+AG+C model does not fit well the late time light curve, especially in the X-rays and in the radio bands, unlike the GW+AG fit, recognising it as a possible excess not due to the afterglow emission.We try to account for this adding a constant flux component at late times.The latter helps in fitting that part of the light curve, but results in very similar posteriors to the GW+AG+C fit without it (see the full results in Section 4.3).This shows that, adding the afterglow centroid motion in the analysis, provides robustness to the fit.The discrepancy in the fit of the light curve at late time when including the jet centroid motion is due to the fits preference for two different viewing angles and types of jet: the GW+AG fit prefers a Table 2. Fit results for GW170817 for a Gaussian (GJ) and a power law jet (PLJ).We report the medians and the 16th-84th percentiles.In the second column we report the results for the GW-only fit, in the third and fifth columns the results of the fit of the broad-band afterglow and the GW, in the fourth and sixth columns the results of the joint fit of broad-band afterglow, centroid motion and GW.The two last columns provide the results of the GW+AG+C fit with a constant component of the type  , = 10  modelling the late time emission.simulations), and use a prior on the jet break Lorentz factor from the centroid measurements, which acts also on the jet opening angle and on the viewing angle.They find  0 = 75.46+5.34 −5.39 km s −1 Mpc −1 .Leaving the Lorentz factor free leads to an opening angle of around 7 deg, which is instead consistent with our GW+AG results for   .

Parameter
Our values of   and   from the GW+AG+C fit are in agreement with other works that included the centroid motion in their analysis.Ghirlanda et al. (2019) predicts   = 3.1 ± 1 deg, with a viewing angle of about 15 deg, Mooley et al. (2018Mooley et al. ( , 2022) )

About the difference between GW+AG and GW+AG+C fits
The GW+AG and GW+AG+C produce quite different results, not only regarding the Hubble constant, the luminosity distance and the viewing angle, but also the energetics and microphysics of the jet.Indeed, while the GW+AG fit results in a large   , a broader jet profile and less energy on the jet axis, the GW+AG+C fit results in a small   , a highly collimated jet with a large energy on the jet axis and a less dense circumburst medium.The viewing angle values are about 5 away, which is quite singular, considering that the event is the same.This, as we stated above, is due to the light curve data points at late times, which are well captured by the GW+AG fit, but  not by the GW+AG+C fit, see Fig. 1.In particular, when including the afterglow centroid motion in the fit, the latter prevails over the light curve data points at late times, resulting in a low viewing angle and a fit of the light curve that, at late times, shows deviations from the flux observations.In the GW+AG+C, the centroid motion is able to constrain very well   to 18.2 +1.2 −1.5 deg, which then translates into a constraint also on   = 2.85 +0.24  −0.20 deg.This happens because of the degeneracy between the two angles, proper of the Gaussian jet light curve (see Fig. 4).In fact, its rising slope depends on their ratio, which, in this fit, is about 6.4.In the GW+AG, instead, there are no constraints on   or   individually, but just on their ratio, from the rising slope of the light curve (see also Ryan et al. 2020;Nakar & Piran 2021).This still leads to the same ratio of about 6.5, but   = 50.1 +5.1 −5.4 deg and   = 7.73 +0.86  −0.80 deg, about 5 sigma away from the GW+AG+C case.The GW+AG fit, not being constrained by the centroid data set, is free to account for the mild decay of the light curve at late times by anticipating the non-relativistic phase.In particular, we estimate the non-relativistic time (Ryan et al. 2020) to be    = 880 +290 −210 days (GW+AG), with respect to    = 13000 +2700 −2400 days (GW+AG+C).Therefore, at late times, according to the parameters of the GW+AG fit, the jet is non-relativistic.The anticipation of the relativistic phase is obtained mainly by acting on the  0 ,  0 parameters.However, the one order of magnitude lower energy and two orders of magnitude higher circumburst density would shift the flux at low values and the break at earlier times, since   ∝ ( 0 / 0 ) 1/3 (  + 1.24  ) 8/3 (Ryan et al. 2020).This is balanced by the fit with higher values of   and   , in order to bring back the jet break (the peak) at about 130 days, and to adjust the rising slope of the light curve.This influences also the early decreasing slope (before the non-relativistic phase), as a higher   (wider jet) provides a larger surface area, so the jet is brighter and the flux is higher.The parameters   ,   ,   have mainly the role of shifting the flux.The  parameter stays the same in the two fits, as it is constrained by the spectrum.  is unconstrained in both fits, however it is better constrained in the GW+AG fit, mainly because   is larger, and, being a Gaussian jet,   has to be lower than   .
In other words, the good fit in the GW+AG case is provided by a combined effect of the high   (in the decreasing slope right after the peak) and the anticipation of the non-relativistic phase (in the slope at late times).

Changing the structure of the jet
In the case of a power law jet, the degeneracy between   and   is not as strong as for the Gaussian geometry, the rising phase slope is a function of ,   and   (see Eq.( 33) of Ryan et al. 2020).This can be seen in the right panel of Fig. 4, where the GW+AG fit is represented in red contours.The GW+AG and GW+AG+C are written in the fifth and sixth column of Table 2, while the fits of the afterglow light curve and centroid motion are in Fig. 5. Also for this jet structure, the GW+AG and GW+AG+C produce quite different results, and the reasoning in Section 4.1 is still valid.The majority of the parameters from the GW+AG+C and GW+AG fits, assuming a power law model, are in agreement within 1 with the Gaussian jet model, this is probably due to the fact that, at early times, the afterglow light curve rises, so  has to be large.At the same time, the larger is , the more the Gaussian and power law structures are similar.For example, in the case of the GW+AG+C fit, for  = 10.8, the  () of a Gaussian and power law structures are very similar within ∼ 3  , after which the decay is shallower for the power law structure.
The GW+AG fit produces a larger   = 62.7 +5.0 −4.3 deg, a smaller   = 5.57 +0.69  −0.62 deg and smaller   = 23.7 +3.8  −3.4 Mpc, than the Gaussian jet, these parameters are, however, in agreement within 2 with the latter.The 2D posterior for   and   are represented in Fig. 2, top panel, in green dashed contours.The microphysics and the energetics are in agreement within 1 with the Gaussian jet results.
In the GW+AG+C fit, the parameters are in agreement within 1 with the Gaussian jet model, except for   = 2.18 +0.20  −0.16 deg, which is within 2.The   and   2D posteriors for the power law jet are represented in the right panel of Fig. 4, we can see that at 3 there are samples with large viewing angles, usually preferred instead from the GW+AG fits.The results for   is 19.8 +1.3  −1.8 deg and for   = 43.0+1.4  −1.4 Mpc.The 2D posterior distributions of   and   are in Fig. 2, top panel, in black dashed contours, almost superimposed to the Gaussian jet results (purple and yellow coloured contours).Also in this case there is a small region of the parameter space at 3 at large   and small   , which gets cancelled when estimating  0 , bottom panel.
The  0 values retrieved in these fits are 70.2 +4.6 −4.4 km s −1 Mpc −1 for the GW+AG+C fit, and 127 +22 −19 km s −1 Mpc −1 (medians, 16th-84th percentiles) for the GW+AG fit, the  0 posteriors are represented in the central panel of Fig. 3, in blue and red respectively.For this event, if we use the complete data set (the most robust fit), the change in jet model does not significantly influence  0 , the power law jet predicts an  0 larger than 1.2 km s −1 Mpc −1 , which is a 2% difference, with respect to a Gaussian jet, but still in agreement within the uncertainties.In the GW+AG there is a 30% difference, but the two  0 s are compatible within 1.
In order to assess if the unknown jet structure leads to systematics in the estimation of  0 , we simulate an afterglow light curve and centroid movement using a Gaussian jet, then we fit them twice, assuming a Gaussian and a power law structure.To keep this simulation as similar to GW170817 as possible, we keep the GW170817 detection times, errors and frequencies for the afterglow light curve and centroid motion, but we adopt fluxes and positions predicted by the model, with a Gaussian variation.We simulate the EM data sets  The GW-only fit is represented in black (same distribution as Fig. 3), the GW+AG+C assuming a Gaussian jet in the fit is represented in violet, while the GW+AG+C assuming a power law jet is represented in green.The vertical dashed lines represent the 16th and 84th percentiles of each distribution.The magenta and yellow shaded regions represent the 1 interval of the Planck and SH0ES measurements respectively.
assuming a Gaussian jet and the parameters in Table 2, medians in the fourth column.In this way, we do not include the excess in the flux at late times, which we are not interested in, as we are focusing on the influence of the jet structure.These EM data sets are then fitted with GW two times, one assuming a Gaussian jet, and the other assuming a power law jet.
For both jet structures, we retrieve the parameters of the GW, the energetics and microphysics in agreeement within 1 with the median values in Table 2, fourth column.Focusing on the distance and the geometry of the system, assuming a Gaussian jet, we retrieve   = 19.3+1.5  −1.7 deg,   = 3.01 +0.28 −0.25 deg and   = 43.8+1.5 −1.7 Mpc, while for a power law jet   = 20.2+1.6  −1.8 deg,   = 2.40 +0.24 −0.21 deg and   = 43.6 +1.5 −1.5 Mpc.As for the case of GW170817, the power law jet tends to give a slightly higher (lower) viewing angle (jet opening angle), which is in agreement within 2 with the simulated values.This, however, does not influence much the luminosity distance.The  0 posteriors that we retrieve from these fits are represented in Fig. 6, with purple (Gaussian jet fit) and green (power law jet fit) colors.It seems that the Hubble constant, as   , is not influenced by the different structure, resulting in  0 = 68.8+4.5 −4.4 km s −1 Mpc −1 for a Gaussian jet and  0 = 69.4+4.5 −4.4 km s −1 Mpc −1 for a power law jet (medians, 16th-84th percentiles).This is a less than 1% difference, which is well inside the 1 range, nonetheless is also at the same level of the Planck uncertainty on  0 .For this reason, in the future, with a larger number of events, it could be important to assess if this, at the moment negligible difference, is just a statistical fluctuation or a real fluctuation due to the changing jet structure.

Adding a constant component in the flux at late times
In the case of GW170817, the high viewing angle preference mainly arise at late times, where there is a flux excess.This is either due to some missing emission at late times in the jet model itself, or due to a new component becoming visible, like a kilonova afterglow or the emission from a long-lived pulsar, in the former case we expect to see a rising flux in future observations, in the latter a constant flux (Troja   et al. 2020;Hajela et al. 2019;Piro et al. 2019).If a flux additive component is included in the fit, indeed the jet viewing angle slightly decreases, see for example Troja et al. (2021); Balasubramanian et al. (2021); Hajela et al. (2022); Wang et al. (2023); Ryan et al. (2023).
We fit the same data set in Fig. 1, adding a constant flux component of the type   =  ,agpy + 10  , where  ,agpy is the flux predicted by afterglowpy and  is a parameter in the fit.This is done only at late times and at all frequencies.The  parameter has three possible values, depending on the frequency:  radio , with a uniform prior in [-3.5,-2],  optical with a uniform prior in [-5.5, -4.5] and  X−rays , with uniform prior in [-8, -7].
The results of the GW+AG+C and GW+AG fit are written in Table 2, last two columns, while the fit of the broad-band afterglow light curve and centroid motion are in Fig. 7.
This model can well fit the afterglow light curve and centroid motion, in both cases.The parameters values from the GW+AG+C fit are in agreement within 1 with the ones from the simple Gaussian jet model (Table 2, fourth column).In the GW+AG+C the viewing angle   = 17.2 +1.1 −1.2 deg is lower with respect to the simple Gaussian jet model, so the distance   = 44.3+1.4  −1.3 Mpc is slightly larger (see, for example, Fig. 8, top panel).The viewing angle and the jet opening angle are better constrained, but the error on the distance is unvaried with respect to the previous analysis.This fit leads to an  0 of 68.0 +4.4  −4.2 km s −1 Mpc −1 (see bottom panel of Fig. 3 and bottom panel of Fig. 8), which is in agreement with the value from the GW+AG+C fit with a simple Gaussian jet.The inclusion of a constant component that accounts for the latetime behaviour does not significantly influence the parameter posteriors with respect to the model without it, so we can say that the GW+AG+C fit and model are robust.
In the case of the GW+AG fit, the addition of the constant component brings some improvements in the results.The jet parameters are compatible at most within 2 with the GW+AG+C (with constant) fit, except for   = 5.37 +0.97  −0.87 deg and   = 35.2+5.7 −6.2 deg, which are within 3.Thanks to the inclusion of the constant component at late times, the viewing angle decreases with respect to the GW+AG fit with the simple Gaussian jet model, and the error on the distance is about 2 times better that the GW fit, cutting part of the tails of the   −   degeneracy, see Fig. 8, top panel.Indeed, the Hubble constant value that we retrieve is 78.5 +7.9 −6.4 km s −1 Mpc −1 , see the bottom panels of Fig. 3 and Fig. 8, that is compatible within 1 with the GW+AG+C fit, including a constant component.Moreover, this result is compatible within 1 with Guidorzi et al. (2017); Wang & Giannios (2021); Wang et al. (2023).

Prospects for jet centroid observations
Using GW, afterglow light curve and centroid motion, leads to  0 = 69.0+4.4  −4.3 km s −1 Mpc −1 .However, the precision of this measure is not at the level of SH0ES or Planck, in order to reach these, we would need at least ∼ 10 events and ∼ 60 events respectively.In this Section, we estimate the likelihood that a new GW event, followed by the detection of the afterglow light curve and the measurement of the afterglow centroid motion, is seen in the next GW Observing runs O4 and O5.
From the GW simulations of Petrov et al. (2022), we generate the EM counterparts of more than a thousand binary neutron star events detectable in O4 (Singer 2021a) and in O5 (Singer 2021b), assuming a Gaussian jet.Each GW event is characterized by an inclination and a luminosity distance, which we use to generate their afterglow light curve and centroid motion.For inclinations larger than 90 deg, we convert them in EM viewing angles as explained at the end of Section 3.1.2and in Eq. (1) of Gianfagna et al. 2023.We assume all the other parameters to be the same as GW170817 (see Table 2, fourth column).Moreover, we assume that all the events are well localized and easy to be followed up by the radio telescopes.This will lead to very optimistic rates.We adopt VLBI as the reference radio facility, both for O4 and O5, so we assume a sensitivity in the radio band of 24Jy (the observations of GW170817 afterglow centroid motion reached an RMS of about 8Jy), and a resolution of 1.5 mas (Ghirlanda et al. 2019).These performances can be achieved also, for example, with the European VLBI Network (EVN) 2 .
The centroid data set is composed of the same detection times of GW170817, but we adopt fluxes and positions predicted by the model.We assume that the afterglow centroid motion is visible if the offset between two data points is above the assumed resolution.Regarding the afterglow light curve, we define an event as detectable if its afterglow peak is above the sensitivity.
In the case of O4 (operating from 2023 to 2025), the GW rate of events is 34 +78 −25 yr −1 (Petrov et al. 2022).We find that 7% of the total have detectable flux, resulting in a rate of 2.4 +5.5  −1.8 yr −1 (for the whole sky).Regarding jet centroid observations, we find that only 0.13% of events has a detectable afterglow flux and centroid, see red dots in Fig. 9, in agreement with Mastrogiovanni et al. (2021).This translate into a rate of 0.05 +0.11 −0.03 yr −1 , therefore it is very unlikely that the jet centroid will be measured again during O4.
In the case of the O5 run, which is due after 2027, the predicted GW rate is 190 +410 −130 yr −1 .We find that 6% have a peak flux above the sensitivity, resulting in a rate of 11 +25 −8 yr −1 .The jet centroid motion is visible in 0.09% of the cases leading to a rate of 0.17 +0.36  −0.12 yr −1 .The latter is still a very low rate, with a slightly lower event fraction than O4, due to the fact that O5 will probe larger distances, which very unlikely will have a detectable afterglow centroid motion.The event rate for the GW, afterglow light curve and centroid motion is slightly larger than O4, despite the same number of events at distances lower than 100 Mpc (∼ 1yr −1 both for O4 and O5).For this reason, we can say that the rate fluctuation is just due to the small number of events.
As is shown in Fig. 9, at large distances we mainly see on-axis or almost on-axis events (with a small   ), these events will not have a visible jet motion, as the observer is within (or just outside) the jet's opening angle.This results in a small or null offset, which is hardly detected with sensitivities of the order of the mas.However, Figure 9.The dots represent GW events simulated by (Petrov et al. 2022), in the case of the O4 run (Singer 2021a).Depending on their   and   , we highlight in blue the ones that have a detectable afterglow counterpart in the radio band and in red the ones that have also a detectable afterglow centroid motion.
if the jet has a large viewing angle, the peak of the afterglow will be at low fluxes, not reaching the VLBI sensitivity.Indeed, for the O4 run, the events that have a coincident detectable GW, afterglow light curve and centroid motion are very similar to GW170817 (at small distances and with   ∼ 20 deg, red dots in Fig. 9).

CONCLUSIONS
The estimation of the Hubble constant  0 exploiting GW, also known as standard sirens method, is a very powerful tool to try to solve the Hubble tension.However, its main issue is the degeneracy between the viewing angle and the luminosity distance of the event, which precludes reaching the level of precision of Planck and SH0ES.In this work, we use this method to estimate  0 , with additional constraints that help in breaking this degeneracy.Using Bayesian analysis, we fit simultaneously the EM and GW domains for the event GW170817.The electromagnetic data set includes the broadband afterglow and the centroid motion of the relativistic jet from HST and VLBI observations.From here, we estimate the Hubble constant and we test its robustness depending on the data set used, on the assumed structure of the jet and on the presence of a possible late time flux excess in the afterglow light curve.
A GW-only fit leads to an Hubble constant value of  0 = 77 +21 −10 km s −1 Mpc −1 (median, 16th-84th percentiles).The almost 20% error is due to the degeneracy stated above.The latter can be broken exploiting independent EM messengers, like the afterglow light curve and centroid motion, at least in the case of GW170817.
In GW+AG analysis, we join the GW and the afterglow light curve.This fit reduces the   −   degeneracy, but gives  0 = 96 +13 −10 km s −1 Mpc −1 .This high value follows from the low value of distance (  = 31.3+3.0 −3.6 Mpc) and high value of viewing angle (  = 50.1 +5.1 −5.4 deg).This behaviour is caused by a possible late time excess in the afterglow flux, these data points are well modelled and are driving the result of the fit.Therefore, for the specific case of GW170817, using only the afterglow as EM counterpart is not enough to get a reliable measurement of  0 .
The GW+AG+C fit, instead, joining the GW, the afterglow light curve and the centroid motion, breaks the   −   degeneracy and results in  0 = 69.0+4.4  −4.3 km s −1 Mpc −1 , which is in agreement with other estimations of this parameter using GW170817 and is about 3 times more precise than the GW-only  0 measurement.This is because of the very strong constraint on the viewing angle given by the afterglow centroid data set.As a consequence, the latter model does not fit well (even if the residuals are ≤ 3.5) the late time flux data points.The viewing angle is   = 18.2 +1.2 −1.5 deg and the distance is   = 43.7 +1.4  −1.4 Mpc.Thus, in the GW+AG+C fit a small value of   , and consequently a highly collimated jet and a large energy on the jet axis, is preferred.In the GW+AG fit a large   , with a broader profile and less energy on the jet axis, is preferred instead.
The possible excess in the afterglow light curve at late times can be explained as either something missing in the jet model, or as a new emission becoming visible (Troja et al. 2020;Hajela et al. 2019;Piro et al. 2019).In either cases, adding a constant flux component to the GW+AG+C model at late times leads to posterior probabilities that are in agreement within 1 with the fit without this constant component ( 0 = 68.0+4.4  −4.2 km s −1 Mpc −1 ), but helps in better fit the late times data.This shows that the model and the GW+AG+C results are robust.Instead, adding this constant flux component to the fit of GW+AG leads to more acceptable values of viewing angle, luminosity distance and Hubble constant:   = 35.2+5.7 −6.2 deg,   = 38.6 +2.5 −3.0 Mpc and  0 = 78.5 +7.9 −6.4 km s −1 Mpc −1 .The latter is compatible within 1 with the GW+AG+C fit.
Finally, it seems that the Hubble constant is not influenced by the assumption on the structure of the jet (either Gaussian or power law), at the present level of precision.
We also note that other systematic uncertainties in the Hubble constant estimation can arise from the estimation of the peculiar velocity of the host galaxy (Hjorth et al. 2017;Nicolaou et al. 2020).The latter, in this work (see Eq. ( 5)), is included in the Hubble flow velocity (Abbott et al. 2017b).In our analysis, a shift in the NGC4993 peculiar velocity of ∼ 140km/s leads to a shift in  0 of ∼ 4km s −1 Mpc −1 , in agreement with Nicolaou et al. (2020).This is still inside the  0 precision reached in this work, but, for a larger number of events, will become one of the main sources of uncertainty.
The best  0 precision reached with this method is 4km s −1 Mpc −1 , in the case of GW+AG+C fit.This is not good enough to prefer either the Planck or the SHoES  0 , yet.More events are needed to reach their level of precision.However, in the future, we do not expect many events that have coincident detections of GW, afterglow light curve and centroid motion.Using the GW simulations from Petrov et al. (2022), we generate the EM counterparts of more than a thousand binary neutron star events detected in O4 (Singer 2021a) and O5 (Singer 2021b).With the VLBI image resolution and sensitivity, we estimate that, both for O4 and O5, the rate of GW, afterglow light curve and centroid motion joint detections is a fraction of event per year.
To conclude, by introducing additional constraints based on astronomical observations, there is the potential to introduce systematic biases, that could affect the standard sirens measurements (Chen 2020; Nicolaou et al. 2020;Govreen-Segal & Nakar 2023).As we show in this work, the viewing angle in the EM modelling is affected by the type of data set used.For this reason, it is fundamental to include in the analysis all the messengers available, in order to have robust results.At the moment, the uncertainty of the standard sirens methods is still too large with respect to the early or late-time Universe  0 s, but in the future attention should be taken, to avoid biases.For example, in a GW170817-like case, measurements at very late times could confirm or qualify as a systematic the milder decrease of the flux at late times.Regarding jet centroid studies, measures at both early and late times will be important to constrain its motion and its viewing angle.For these reasons, highly sensitive instruments are needed, like Athena (Piro et al. 2022) in the X-rays or SKA (Square Kilometre Array, Braun et al. 2019) in the radio band.In the distant future (mid 2030s), facilities like Next Generation VLA (ngVLA) will reach the mas resolution (or lower, Beasley et al. 2019), increasing the chances of detecting the motion of the relativistic jet.

Figure 1 .
Figure 1.Top panel.Broad-band afterglow of GW170817: data and fits.From bottom to top, red points refer to the X-ray observations by Chandra and XMM at 5 keV, orange ones to observations by HST, F606W filter, in the optical band, and blue ones to observations in the radio band from VLA (Very Large Array) at 3 GHz.The continuous and dotted lines represent the fit of the GW, broad-band afterglow, and centroid motion (GW+AG+C) data and of the GW and afterglow (GW+AG) data, respectively.For sake of simplicity, the fit for the radio band is plotted only for the observations at 3 GHz, but it is not limited to this single frequency.Bottom panel.Centroid motion of the relativistic jet from HST and VLBI images at 8 (negative RAs), 75, 206 and 230 days(Mooley et al. 2018;Ghirlanda et al. 2019).The blue dots represents the positions predicted by the model, the blue contours represent the 68% probability region.
, bottom panel, and Fig.3, top panel.As explained in more details in Sections 4.1 and 4.3, this result is mostly driven by the possible presence of a late time additional component, which can be seen in the top panel of Fig.1.The GW+AG model (dotted line) fits very well the light curve, especially the data points at late time.The latter force the model to prefer a high   with respect to the fit including also the jet centroid motion, GW+AG+C, represented with a solid line.Indeed,Wang & Giannios (2021), using the same messengers but limiting the light curve data up to ∼ 300 days (when no flux deviation is present yet), retrieve  0 = 69.5 ± 4 km s −1 Mpc −1 , with a   = 43.4± 1 Mpc and   = 22 ± 1 deg.The jet structure model they use is from 3-dimensional general-relativistic megnetohydrodynamical simulations.AlsoWang et al. (2023) fit the afterglow light curve and include an additional component at late times, a sub-relativistic kilonova outflow.They estimate  0 = 71.80+4.15−4.07 km s −1 Mpc −1 .The kilonova component helps in the fit of the light curve, keeping the viewing angle around 30 deg.Also in this case, they model the jet using hydrodynamic simulations.Guidorzi et al. (2017) get  0 = 75.5+11.6 an opening angle of < 5 deg, and a viewing angle < 24 deg, while Ren et al. (2020) find   = 3.1 deg and   = 17.4 deg.

Figure 3 .
Figure 3. Histograms of the Hubble constant  0 posterior from the GW-only fit, in black, the GW+AG in red and the GW+AG+C in blue.The vertical dashed lines represent the 16th and 84th percentiles of each distribution.The magenta and yellow shaded regions represent the 1 interval of the Planck and SH0ES measurements respectively.Top panel: Gaussian jet.Central panel: power law jet.Bottom panel: Gaussian jet with the addition of a constant component at late times.

Figure 4 .
Figure 4. Contour plots of the viewing angle and jet opening angle.The contours represent the 68%, 95%, 99.7% probabilities.The blue contour lines represent the result from the joint GW+AG+C fit, while the red ones represent the result from the GW+AG fit.Left panel: Gaussian jet.Right panel: power law jet.

Figure 5 .
Figure 5. Same as Fig. 1, but assuming a power law structure for the jet.

Figure 6 .
Figure 6.Histogram of the  0 posterior distribution for a simulated event.The GW-only fit is represented in black (same distribution as Fig.3), the GW+AG+C assuming a Gaussian jet in the fit is represented in violet, while the GW+AG+C assuming a power law jet is represented in green.The vertical dashed lines represent the 16th and 84th percentiles of each distribution.The magenta and yellow shaded regions represent the 1 interval of the Planck and SH0ES measurements respectively.

Figure 7 .
Figure 7. Same as Fig. 1, but including an additional constant flux component in the model at late times.

Figure 8 .
Figure 8. Same as Fig. 2, but using a Gaussian jet with the addition of a constant flux component at late times.
Top panel.Contour plot of the viewing angle and luminosity distance for the GW, GW+AG and GW+AG+C fits.The contours represent the 68%, 95%, 99.7% probability regions.Bottom panel.The same contour plot as above, but switching to  0 , instead of   .The magenta and yellow regions represent the 1  of the Planck and SHoES measurements respectively.The Gaussian jet results are represented with filled contours, while the power law jet with empty contours and dashed lines.